text:
module ReactiveMPModelsGMMTest

using Test, InteractiveUtils
using Rocket, ReactiveMP, GraphPPL, Distributions
using BenchmarkTools, Random, Plots, Dates, LinearAlgebra
using StableRNGs

## Model definition
## -------------------------------------------- ##
@model [ default_factorisation = MeanField() ] function univariate_gaussian_mixture_model(n)

    s ~ Beta(1.0, 1.0)

    m1 ~ NormalMeanVariance(-2.0, 1e3)
    w1 ~ GammaShapeRate(0.01, 0.01)

    m2 ~ NormalMeanVariance(2.0, 1e3)
    w2 ~ GammaShapeRate(0.01, 0.01)

    z = randomvar(n)
    y = datavar(Float64, n)

    for i in 1:n
        z[i] ~ Bernoulli(s)
        y[i] ~ NormalMixture(z[i], (m1, m2), (w1, w2))
    end

    return s, m1, w1, m2, w2, z, y
end

@model function multivariate_gaussian_mixture_model(rng, L, nmixtures, n)

    z = randomvar(n)
    m = randomvar(nmixtures)
    w = randomvar(nmixtures)

    basis_v = [ 1.0, 0.0 ]

    for i in 1:nmixtures
        # Assume we now only approximate location of cluters's mean
        approximate_angle_prior = ((2π + + rand(rng)) / nmixtures) * (i - 1)
        approximate_basis_v  = L / 2 * (basis_v .+ rand(rng, 2))
        approximate_rotation = [ cos(approximate_angle_prior) -sin(approximate_angle_prior); sin(approximate_angle_prior) cos(approximate_angle_prior) ]
        mean_mean_prior     = approximate_rotation * approximate_basis_v
        mean_mean_cov       = [ 1e6 0.0; 0.0 1e6 ]

        m[i] ~ MvNormalMeanCovariance(mean_mean_prior, mean_mean_cov)
        w[i] ~ Wishart(3, [ 1e2 0.0; 0.0 1e2 ])
    end

    s ~ Dirichlet(ones(nmixtures))

    y = datavar(Vector{Float64}, n)

    means = tuple(m...)
    precs = tuple(w...)

    for i in 1:n
        z[i] ~ Categorical(s) where { q = MeanField() }
        y[i] ~ NormalMixture(z[i], means, precs) where { q = MeanField() }
    end

    return s, z, m, w, y
end
## -------------------------------------------- ##
## Inference definition
## -------------------------------------------- ##
function inference_univariate(data, n_its)
    n = length(data)
    model, (s, m1, w1, m2, w2, z, y) = univariate_gaussian_mixture_model(n);

    mswitch = keep(Marginal)
    mm1 = keep(Marginal)
    mm2 = keep(Marginal)
    mw1 = keep(Marginal)
    mw2 = keep(Marginal)

    fe = keep(Float64)

    m1sub = subscribe!(getmarginal(m1), mm1)
    m2sub = subscribe!(getmarginal(m2), mm2)
    w1sub = subscribe!(getmarginal(w1), mw1)
    w2sub = subscribe!(getmarginal(w2), mw2)
    switchsub = subscribe!(getmarginal(s), mswitch)
    fesub = subscribe!(score(Float64, BetheFreeEnergy(), model), fe)

    setmarginal!(s, vague(Beta))
    setmarginal!(m1, NormalMeanVariance(-2.0, 1e3))
    setmarginal!(m2, NormalMeanVariance(2.0, 1e3))
    setmarginal!(w1, vague(GammaShapeRate))
    setmarginal!(w2, vague(GammaShapeRate))

    for i in 1:n_its
        update!(y, data)
    end

    unsubscribe!((fesub, switchsub, m1sub, m2sub, w1sub, w2sub))

    return mswitch, mm1, mm2, mw1, mw2, fe
end

function inference_multivariate(rng, L, nmixtures, data, viters)
    n = length(data)

    model, (s, z, m, w, y) = multivariate_gaussian_mixture_model(rng, L, nmixtures, n)

    means_estimates  = keep(Vector{Marginal})
    precs_estimates  = keep(Vector{Marginal})
    switch_estimates = keep(Marginal)
    fe_values        = keep(Float64)

    s_sub = subscribe!(getmarginal(s), switch_estimates)
    m_sub = subscribe!(getmarginals(m), means_estimates)
    p_sub = subscribe!(getmarginals(w), precs_estimates)
    f_sub = subscribe!(score(Float64, BetheFreeEnergy(), model), fe_values)

    setmarginal!(s, vague(Dirichlet, nmixtures))

    basis_v = [ 1.0, 0.0 ]

    for i in 1:nmixtures
        # Assume we now only approximate location of cluters's mean
        approximate_angle_prior = ((2π + + rand(rng)) / nmixtures) * (i - 1)
        approximate_basis_v  = L / 2 * (basis_v .+ rand(rng, 2))
        approximate_rotation = [ cos(approximate_angle_prior) -sin(approximate_angle_prior); sin(approximate_angle_prior) cos(approximate_angle_prior) ]
        mean_mean_prior     = approximate_rotation * approximate_basis_v
        mean_mean_cov       = [ 1e6 0.0; 0.0 1e6 ]

        setmarginal!(m[i], MvNormalMeanCovariance(mean_mean_prior, mean_mean_cov))
        setmarginal!(w[i], Wishart(3, [ 1e2 0.0; 0.0 1e2 ]))
    end

    for i in 1:viters
        update!(y, data)
    end

    unsubscribe!((s_sub, m_sub, p_sub, f_sub))

    return switch_estimates, means_estimates, precs_estimates, fe_values
end

@testset "Gaussian Mixture Model" begin

    @testset "Univariate" begin
        ## -------------------------------------------- ##
        ## Data creation
        ## -------------------------------------------- ##
        n = 150

        rng = StableRNG(12345)

        switch = [ 1/3, 2/3 ]
        z      = rand(rng, Categorical(switch), n)
        y      = Vector{Float64}(undef, n)

        μ1 = -10.0
        μ2 = 10.0
        w  = 1.777

        dists = [
            Normal(μ1, sqrt(inv(w))),
            Normal(μ2, sqrt(inv(w))),
        ]

        for i in 1:n
            y[i] = rand(rng, dists[z[i]])
        end
        ## -------------------------------------------- ##
        ## Inference execution
        mswitch, mm1, mm2, mw1, mw2, fe = inference_univariate(y, 10);
        ## -------------------------------------------- ##
        ## Test inference results
        @test length(mswitch) === 10
        @test length(mm1) === 10
        @test length(mm2) === 10
        @test length(mw1) === 10
        @test length(mw2) === 10
        @test length(fe) === 10 && all(filter(e -> abs(e) > 1e-3, diff(getvalues(fe))) .< 0)
        @test abs(last(fe) - 280.3276104) < 0.01

        ms = mean(last(mswitch))

        @test ((abs(ms - switch[1]) < 0.1) || (abs(ms - switch[2]) < 0.1))

        ems = sort([ last(mm1), last(mm2) ], by = mean)
        rms = sort([ μ1, μ2 ])

        foreach(zip(rms, ems)) do (r, e)
            @test abs(r - mean(e)) < 0.2
        end

        ews = sort([ last(mw1), last(mw2) ], by = mean)
        rws = sort([ w, w ])

        foreach(zip(rws, ews)) do (r, e)
            @test abs(r - mean(e)) < 0.2
        end
        ## -------------------------------------------- ##
        ## Form debug output
        base_output = joinpath(pwd(), "_output", "models")
        mkpath(base_output)
        timestamp        = Dates.format(now(), "dd-mm-yyyy-HH-MM")
        plot_output      = joinpath(base_output, "univariate_gmm_model_plot_$(timestamp)_v$(VERSION).png")
        benchmark_output = joinpath(base_output, "univariate_gmm_model_benchmark_$(timestamp)_v$(VERSION).txt")
        ## -------------------------------------------- ##
        ## Create output plots
        dim(d) = (a) -> map(r -> r[d], a)
        mp = plot(mean.(mm1), ribbon = var.(mm1) .|> sqrt, label = "m1 prediction")
        mp = plot!(mean.(mm2), ribbon = var.(mm2) .|> sqrt, label = "m2 prediction")
        mp = plot!(mp, [ μ1 ], seriestype = :hline, label = "real m1")
        mp = plot!(mp, [ μ2 ], seriestype = :hline, label = "real m2")

        wp = plot(mean.(mw1), ribbon = var.(mw1) .|> sqrt, label = "w1 prediction", legend = :bottomleft, ylim = (-1, 3))
        wp = plot!(wp, [ w ], seriestype = :hline, label = "real w1")
        wp = plot!(wp, mean.(mw2), ribbon = var.(mw2) .|> sqrt, label = "w2 prediction")
        wp = plot!(wp, [ w ], seriestype = :hline, label = "real w2")

        swp = plot(mean.(mswitch), ribbon = var.(mswitch) .|> sqrt, label = "Switch prediction")

        swp = plot!(swp, [ switch[1] ], seriestype = :hline, label = "switch[1]")
        swp = plot!(swp, [ switch[2] ], seriestype = :hline, label = "switch[2]")

        fep = plot(fe[2:end], label = "Free Energy", legend = :bottomleft)

        p = plot(mp, wp, swp, fep, layout = @layout([ a b; c d ]), size = (1000, 700))
        savefig(p, plot_output)
        ## -------------------------------------------- ##
        ## Create output benchmarks
        benchmark = @benchmark inference_univariate($y, 10) seconds=15
        open(benchmark_output, "w") do io
            show(io, MIME("text/plain"), benchmark)
            versioninfo(io)
        end
        ## -------------------------------------------- ##
    end

    @testset "Multivariate" begin
        ## -------------------------------------------- ##
        ## Data creation
        ## -------------------------------------------- ##
        rng = StableRNG(43)

        L         = 50.0
        nmixtures = 3
        n_samples = 500

        probvec = ones(nmixtures)
        probvec = probvec ./ sum(probvec)

        switch = Categorical(probvec)

        gaussians = map(1:nmixtures) do index
            angle      = 2π / nmixtures * (index - 1)
            basis_v    = L * [ 1.0, 0.0 ]
            rotationm  = [ cos(angle) -sin(angle); sin(angle) cos(angle) ]
            mean       = rotationm * basis_v
            covariance = Matrix(Hermitian(rotationm * [ 10.0 0.0; 0.0 20.0 ] * transpose(rotationm)))
            return MvNormal(mean, covariance)
        end

        z = rand(rng, switch, n_samples)
        y = Vector{Vector{Float64}}(undef, n_samples)

        for i in 1:n_samples
            y[i] = rand(rng, gaussians[z[i]])
        end
        ## -------------------------------------------- ##
        ## Inference execution
        s, m, w, fe = inference_multivariate(rng, L, nmixtures, y, 25);
        ## -------------------------------------------- ##
        ## Test inference results
        @test length(s) === 25
        @test length(m) === 25
        @test length(w) === 25
        @test length(fe) === 25
        @test all(filter(e -> abs(e) > 1e-3, diff(getvalues(fe))) .< 0)
        @test abs(last(fe) - 3442.4015524445967) < 0.01

        ems = sort(mean.(last(m)), by = x -> atan(x[2] / x[1]))
        rms = sort(mean.(gaussians), by = x -> atan(x[2] / x[1]))

        foreach(zip(ems, rms)) do (estimated, real)
            @test norm(normalize(estimated) .- normalize(real)) < 0.1
        end
        ## -------------------------------------------- ##
        ## Form debug output
        base_output = joinpath(pwd(), "_output", "models")
        mkpath(base_output)
        timestamp        = Dates.format(now(), "dd-mm-yyyy-HH-MM")
        plot_output      = joinpath(base_output, "multivariate_gmm_model_plot_$(timestamp)_v$(VERSION).png")
        benchmark_output = joinpath(base_output, "multivariate_gmm_model_benchmark_$(timestamp)_v$(VERSION).txt")
        ## -------------------------------------------- ##
        ## Create output plots
        sdim(n) = (a) -> map(d -> d[n], a)

        pe = plot(xlim = (-1.5L, 1.5L), ylim = (-1.5L, 1.5L))

        rp = scatter(y |> sdim(1), y |> sdim(2))

        pe = scatter!(pe, y |> sdim(1), y |> sdim(2))

        e_means = mean.(m[end])
        e_precs = mean.(w[end])

        for (e_m, e_w) in zip(e_means, e_precs)
            gaussian = MvNormal(e_m, Matrix(Hermitian(inv(e_w))))
            pe = contour!(pe, range(-2L, 2L, step = 0.25), range(-2L, 2L, step = 0.25), (x, y) -> pdf(gaussian, [ x, y ]), levels = 7, colorbar = false)
        end

        pfe = plot(fe[2:end], label = "Free Energy")

        p = plot(rp, pe, pfe, size = (600, 600), layout = @layout([ a b; c ]))

        savefig(p, plot_output)
        ## -------------------------------------------- ##
        ## Create output benchmarks
        benchmark = @benchmark inference_multivariate(StableRNG(123), $L, $nmixtures, $y, 25) seconds=15
        open(benchmark_output, "w") do io
            show(io, MIME("text/plain"), benchmark)
            versioninfo(io)
        end
        ## -------------------------------------------- ##
    end

end

end


text:
# script for Fig.1, Fig.2, Fig.4, Fig.6

using MultiscaleGraphSignalTransforms, LightGraphs, Plots, LaTeXStrings, MultivariateStats
pyplot(dpi = 200)

Nx, Ny = 7, 3
G = LightGraphs.grid([Nx, Ny]); N = nv(G);
L = Matrix(laplacian_matrix(G))
Q = incidence_matrix(G; oriented = true)
𝛌, 𝚽 = eigen(L); 𝚽 = 𝚽 .* sign.(𝚽[1, :])';

#################### Fig. 1(a) eigenvectors by nondecreasing eigenvalue ordering
plot(layout = Plots.grid(3, 7))
for i in 1:N
    heatmap!(reshape(𝚽[:, i],(Nx, Ny))', c = :viridis, cbar = false,
                clims = (-0.4,0.4), frame = :none, ratio = 1,
                title = latexstring("\\phi_{", i-1, "}"), titlefont = 12,
                subplot = i)
end
plt = current()
# savefig(plt, joinpath(@__DIR__, "../paperfigs/grid7x3_evsp_title.png"))

#################### Fig. 1(b) eigenvectors by natural frequency ordering
# find correct 2D index
grid2eig_ind = [1,2,3,6,8,12,15,4,5,7,9,13,16,18,10,11,14,17,19,20,21];
eig2grid_ind = sortperm(grid2eig_ind);
eig2dct = Array{Int64,3}(undef, Nx, Ny, 2);
for i = 1:Nx; for j = 1:Ny; eig2dct[i,j,1] = i-1; eig2dct[i,j,2] = j-1; end; end
eig2dct = reshape(eig2dct, (N, 2)); eig2dct = eig2dct[eig2grid_ind, :];

plot(layout = Plots.grid(3, 7))
for i in 1:N
    k = grid2eig_ind[i]
    heatmap!(reshape(𝚽[:,k],(Nx,Ny))', c = :viridis, cbar = false,
                clims = (-0.4,0.4), frame = :none, ratio = 1,
                title = latexstring("\\varphi_{", string(eig2dct[k,1]),
                ",", string(eig2dct[k,2]), "}"), titlefont = 12, subplot = i)
end
plt = current()
# savefig(plt, joinpath(@__DIR__, "../paperfigs/grid7x3_dct_title2.png"))

# DAG pseudo-metric
distDAG = eigDAG_Distance(𝚽, Q, N)

# MDS embedding into R²
D = distDAG
E = transform(fit(MDS, D, maxoutdim=2, distances=true))

# set up all heatmap plots' positions
dx = 0.01; dy = dx;
xej = zeros(Nx, N); yej=zeros(Ny, N);
a = 5.0; b = 9.0;
for k = 1:N
    xej[:,k] = LinRange(E[1,k] - Ny * a * dx, E[1, k] + Ny * a * dx, Nx)
    yej[:,k] = LinRange(E[2,k] - a * dy, E[2, k] + a * dy, Ny)
end

#################### Fig. 2
plot()
for k=1:N
    heatmap!(xej[:, k], yej[:, k], reshape(𝚽[:, k], (Nx, Ny))', c = :viridis,
                colorbar = false, ratio = 1, annotations = (xej[4, k],
                yej[3, k] + b*dy, text(latexstring("\\varphi_{",
                string(eig2dct[k, 1]), ",", string(eig2dct[k, 2]), "}"), 10)))
end
plt = plot!(xlim = [-1.4, 1.3], ylim = [-1.4, 1.3], grid = false, clims = (-0.4, 0.4))
# savefig(plt, joinpath(@__DIR__, "../paperfigs/Grid7x3_DAG_MDS.png"))

#################### Fig. 4
# first level partition
p1x = [-0.2, 1.0, NaN]; p1y = [1.3, -1.0, NaN];
plot!(p1x, p1y, c = :red, legend = false, width = 3)
# second level partition
p2x = [-1.0, 0.2, NaN, 0.4, 1.2, NaN]; p2y = [-0.8, 0.45, NaN, 0.25, 0.2, NaN];
plot!(p2x, p2y, c=:orange, legend = false, width = 2)
plt = current()
# savefig(plt, joinpath(@__DIR__, "../paperfigs/Grid7x3_DAG_2levels_partition.png"))

## Build Dual Graph
Gstar_Sig = dualgraph(distDAG)
GP_dual = partition_tree_fiedler(Gstar_Sig; swapRegion = false)
GP_primal = pairclustering(𝚽, GP_dual)

@time VM_NGWP = vm_ngwp(𝚽, GP_dual)

## level 2 VM-NGWP vectors
j = 3; W_VM = VM_NGWP[:, j, :]'

wav_kl = [[0 0];[0 1];[0 2];[1 0];[1 1];[1 2];[1 3];[2 0];[2 1];[2 2];[3 0];
            [3 1];[3 2];[3 3];[2 3];[2 4];[2 5];[2 6];[3 4];[3 5];[3 6]];
wav_kl = wav_kl[eig2grid_ind,:];

# reorder_ind = [2,3,1,5,7,4,6, 16,17,15,  9,11,8,10, 18,20,21,19, 13,14,12]
reorder_ind = [1,3,2,5,7,4,6, 9,10,8,16,18,15,17, 11,13,14,12,20,21,19]
W_VM = W_VM[:,reorder_ind[eig2grid_ind]];
sgn = ones(N); sgn[grid2eig_ind[[4,6,8,10,14]]] .= -1; W_VM = W_VM * Diagonal(sgn);

#################### Fig. 6
plot()
for k=1:N
    heatmap!(xej[:,k],yej[:,k],reshape(W_VM[:,k],(Nx,Ny))',c=:viridis,colorbar=false,ratio=1,annotations=(xej[4,k], yej[3,k]+b*dy, text(latexstring("\\psi_{", string(wav_kl[k,1]), ",", string(wav_kl[k,2]), "}"),10)))
end
plot!(aspect_ratio = 1, xlim = [-1.4, 1.3], ylim = [-1.4, 1.3], grid = false, clims=(-0.34,0.34))
# first level partition
p1x = [-0.2, 1.0, NaN]; p1y = [1.3, -1.0, NaN]; plot!(p1x, p1y, c = :red, legend = false, width = 3)
# second level partition
p2x = [-1.0, 0.2, NaN, 0.4, 1.2, NaN]; p2y = [-0.8, 0.45, NaN, 0.25, 0.2, NaN]; plot!(p2x, p2y, c=:orange, legend = false, width = 2)
plt = current()
# savefig(plt, joinpath(@__DIR__, "../paperfigs/Grid7x3_DAG_VM_NGWP_lvl2_wavelets.png"))

ext: jl


text:
@compat const DefaultRefType = UInt32

## Pools

# V is always set to CategoricalValue{T}
# This workaround is needed since this type not defined yet
# See JuliaLang/julia#269
if VERSION >= v"0.6.0-dev.2643"
    include_string("""
        type CategoricalPool{T, R <: Integer, V}
            index::Vector{T}
            invindex::Dict{T, R}
            order::Vector{R}
            levels::Vector{T}
            valindex::Vector{V}
            ordered::Bool

            function CategoricalPool{T, R, V}(index::Vector{T},
                                              invindex::Dict{T, R},
                                              order::Vector{R},
                                              ordered::Bool) where {T, R, V}
                pool = new(index, invindex, order, index[order], V[], ordered)
                buildvalues!(pool)
                return pool
            end
        end
    """)
else
    @eval begin
        type CategoricalPool{T, R <: Integer, V}
            index::Vector{T}
            invindex::Dict{T, R}
            order::Vector{R}
            levels::Vector{T}
            valindex::Vector{V}
            ordered::Bool

            function CategoricalPool{T, R}(index::Vector{T},
                                           invindex::Dict{T, R},
                                           order::Vector{R},
                                           ordered::Bool)
                pool = new(index, invindex, order, index[order], V[], ordered)
                buildvalues!(pool)
                return pool
            end
        end
    end
end

immutable LevelsException{T, R} <: Exception
    levels::Vector{T}
end

## Values

immutable CategoricalValue{T, R <: Integer}
    level::R
    pool::CategoricalPool{T, R, CategoricalValue{T, R}}
end

## Arrays

@compat abstract type AbstractCategoricalArray{T, N, R} <: AbstractArray{CategoricalValue{T, R}, N} end
@compat AbstractCategoricalVector{T, R} = AbstractCategoricalArray{T, 1, R}
@compat AbstractCategoricalMatrix{T, R} = AbstractCategoricalArray{T, 2, R}

immutable CategoricalArray{T, N, R <: Integer} <: AbstractCategoricalArray{T, N, R}
    refs::Array{R, N}
    pool::CategoricalPool{T, R, CategoricalValue{T, R}}
end
@compat CategoricalVector{T, R} = CategoricalArray{T, 1, R}
@compat CategoricalMatrix{T, R} = CategoricalArray{T, 2, R}

## Nullable Arrays

@compat abstract type AbstractNullableCategoricalArray{T, N, R} <: AbstractArray{Nullable{CategoricalValue{T, R}}, N} end
@compat AbstractNullableCategoricalVector{T, R} = AbstractNullableCategoricalArray{T, 1, R}
@compat AbstractNullableCategoricalMatrix{T, R} = AbstractNullableCategoricalArray{T, 2, R}

immutable NullableCategoricalArray{T, N, R <: Integer} <: AbstractNullableCategoricalArray{T, N, R}
    refs::Array{R, N}
    pool::CategoricalPool{T, R, CategoricalValue{T, R}}
end
@compat NullableCategoricalVector{T, R} = NullableCategoricalArray{T, 1, R}
@compat NullableCategoricalMatrix{T, R} = NullableCategoricalArray{T, 2, R}

## Type Aliases

@compat CatArray{T, N, R} = Union{CategoricalArray{T, N, R}, NullableCategoricalArray{T, N, R}}
@compat CatVector{T, R} = Union{CategoricalVector{T, R}, NullableCategoricalVector{T, R}}


text:
# This file is a part of RRRMC.jl. License is MIT: http://github.com/carlobaldassi/RRRMC.jl/LICENCE.md

module PSpin3

using Random
using ExtractMacro
using ..Interface
using ..Common

export GraphPSpin3

import ..Interface: energy, delta_energy, neighbors, allΔE, update_cache!

struct GraphPSpin3{K} <: DiscrGraph{Int}
    N::Int
    #allA::Vector{Int} # aux contiguous memory
    A::Vector{Vector{Int}} # neighbors
    U::Vector{Vector{Int}} # unique neighbors
    #J::Vector{Vector{Int}} # currently unused (all 1)
    cache::LocalFields{Int}
    function GraphPSpin3(N::Integer, K::Integer)
        @assert K ≥ 1

        N % 3 == 0 || throw(ArgumentError("N must be divisible by 3, given: $N"))
        A = Vector{Int}[zeros(Int, 2*K) for i = 1:N]
        #allA = zeros(Int, 2 * K * N)
        #A = Array{Vector{Int}}(undef, N)
        #for i = 1:N
            #A[i] = pointer_to_array(pointer(allA, (i-1) * 2K + 1), 2K)
        #end

        for k = 1:K
            l = randperm(N)
            for i = 1:3:(N-2)
                v1, v2, v3 = l[i], l[i+1], l[i+2]
                A[v1][2*(k-1) + 1] = v2
                A[v1][2*(k-1) + 2] = v3
                A[v2][2*(k-1) + 1] = v1
                A[v2][2*(k-1) + 2] = v3
                A[v3][2*(k-1) + 1] = v1
                A[v3][2*(k-1) + 2] = v2
            end
        end

        @assert all(a->!any(x->x==0, a), A)

        U = Vector{Int}[unique(a) for a in A]

        cache = LocalFields{Int}(N)

        return new{K}(N, A, U, cache)
    end
end

@doc """
    GraphPSpin3(N::Integer, K::Integer) <: DiscrGraph{Int}

A `DiscrGraph` implementing a \$p\$-spin regular graph with \$p=3\$. `N` is the number of spins, and must
be divisible by \$3\$; `K` is the connectivity. All interactions are set to \$J=1\$.
""" -> GraphPSpin3(N, K)

function energy(X::GraphPSpin3{K}, C::Config) where {K}
    @assert X.N == C.N
    @extract C : s
    @extract X : A cache #J
    @extract cache : lfields lfields_last
    n = 0
    for x = 1:length(A)
        #Jx = J[x]
        σx = 2 * s[x] - 1
        Ax = A[x]
        p = 1
        lf = 0
        for k = 1:K
            y = Ax[p]
            z = Ax[p+1]
            σy = 2 * s[y] - 1
            σz = 2 * s[z] - 1
            #Jxyz = Jx[k]

            lf -= σy * σz
            p += 2
        end
        lf *= σx
        n += lf
        lfields[x] = 2lf
    end
    @assert n % 3 == 0
    n ÷= 3
    cache.move_last = 0
    fill!(lfields_last, 0)
    return n
end

function update_cache!(X::GraphPSpin3{K}, C::Config, move::Int) where {K}
    @assert X.N == C.N
    @assert 1 ≤ move ≤ C.N
    @extract C : N s
    @extract X : A U cache

    @extract cache : lfields lfields_last move_last
    if move_last == move
        @inbounds begin
            Ux = U[move]
            for y in Ux
                lfields[y], lfields_last[y] = lfields_last[y], lfields[y]
            end
            lfields[move] = -lfields[move]
            lfields_last[move] = -lfields_last[move]
        end
        return
    end

    @inbounds begin
        Ux = U[move]
        for y in Ux
            lfields_last[y] = lfields[y]
        end
        #Jx = J[move]
        σx = 2 * s[move] - 1
        Ax = A[move]
        p = 1
        for k = 1:K
            y = Ax[p]
            z = Ax[p+1]
            σy = 2 * s[y] - 1
            σz = 2 * s[z] - 1
            #Jxyz = Jx[k]
            Δ = 4 * σx * σy * σz
            lfields[y] -= Δ
            lfields[z] -= Δ
            p += 2
        end
        lfm = lfields[move]
        lfields_last[move] = lfm
        lfields[move] = -lfm
    end
    cache.move_last = move

    #lfields_bk = copy(lfields)
    #energy(X, C)
    #@assert lfields_bk == lfields

    return
end

function delta_energy(X::GraphPSpin3, C::Config, move::Int) where {K}
    @assert X.N == C.N
    @assert 1 ≤ move ≤ C.N
    #@extract C : s
    @extract X : cache # A #J
    @extract cache : lfields

    @inbounds Δ = -lfields[move]
    return Δ

    # @inbounds begin
    #     Δ = 0
    #     #Jx = J[move]
    #     σx = 2 * s[move] - 1
    #     Ax = A[move]
    #     p = 1
    #     for k = 1:K
    #         y = Ax[p]
    #         z = Ax[p+1]
    #         σy = 2 * s[y] - 1
    #         σz = 2 * s[z] - 1
    #         #Jxyz = Jx[k]
    #         Δ += σy * σz
    #         p += 2
    #     end
    #     Δ *= 2 * σx
    # end
    # return Δ
end

neighbors(X::GraphPSpin3, i::Int) = return X.U[i]
@generated allΔE(::Type{GraphPSpin3{K}}) where {K} =
    iseven(K) ? Expr(:tuple, ntuple(d->4*(d-1), K÷2+1)...) :
                Expr(:tuple, ntuple(d->2*(2d-1), (K+1)÷2)...)

end


text:
#=
  Original Author: Anastasia Yendiki

  Copyright © 2022 The General Hospital Corporation (Boston, MA) "MGH"

  Terms and conditions for use, reproduction, distribution and contribution
  are found in the 'FreeSurfer Software License Agreement' contained
  in the file 'LICENSE' found in the FreeSurfer distribution, and here:

  https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense

  Reporting: freesurfer@nmr.mgh.harvard.edu
=#

export cart2pol, pol2cart, cart2sph, sph2cart

"""
    cart2pol(x, y)

    Transform Cartesian coordinates (`x`, `y`) to polar coordinates (`φ`, `ρ`),
    where `φ` is in radians.
"""
function cart2pol(x, y)

  φ = atan.(y, x)
  ρ = hypot.(x, y)

  return φ, ρ
end

"""
    pol2cart(φ, ρ)

    Transform polar coordinates (`φ`, `ρ`) to Cartesian coordinates (`x`, `y`),
    where `φ` is in radians.
"""
function pol2cart(φ, ρ)

  x = ρ .* cos.(φ)
  y = ρ .* sin.(φ)

  return x, y
end

"""
    cart2sph(x, y, z)

    Transform Cartesian coordinates (`x`, `y`, `z`) to spherical coordinates
    (`φ`, `θ`, `ρ`), where `φ` and `θ` are in radians.
"""
function cart2sph(x, y, z)

  hypotxy = hypot.(x, y)
  ρ = hypot.(hypotxy, z)
  θ = atan.(z, hypotxy)
  φ = atan.(y, x)

  return φ, θ, ρ
end

"""
    sph2cart(φ, θ, ρ)

    Transform spherical coordinates (`φ`, `θ`, `ρ`) to Cartesian coordinates
    (`x`, `y`, `z`).

    `φ` and `θ` must be in radians.
"""
function sph2cart(φ, θ, ρ)

  z = ρ .* sin.(θ)
  ρcosθ = ρ .* cos.(θ)
  x = ρcosθ .* cos.(φ)
  y = ρcosθ .* sin.(φ)

  return x, y, z
end

"""
    ang2rot(φ, θ)

    Convert polar `φ` and azimuthal `θ` angles to a rotation matrix.

    `φ` and `θ` must be in radians.
"""
function ang2rot(φ, θ)

  # Rotate by φ around the z-axis
  Rz = [cos(φ)  -sin(φ)  0
        sin(φ)   cos(φ)  0
        0         0      1]

  # Rotate by θ around the y-axis
  Ry = [cos(θ)   0   sin(θ)
        0        1   0
       -sin(θ)   0   cos(θ)]

  R =  Rz * Ry

  return R
end

max_stars_repo_name: freesurfer/jupysurf

max_issues_repo_name: freesurfer/jupysurf

max_forks_repo_name: freesurfer/jupysurf


text:
include("utils.jl")

struct TouchstoneSnP
    frequency_unit::String
    parameter_type::String
    parameter_format::String
    reference_resistance::Int # [Ohm]
    number_of_ports::Int
    number_of_frequencies::Int
    number_of_noise_frequencies::Int
    reference::Array{Float64}
    two_port_data_order::String
    matrix_format::String
    network_data::Dict#{Array{Float64}, 1}
    noise_data::Dict#{Array{Float64}, 1}

    function TouchstoneSnP(snp::AbstractString)

        sentences = map(uppercase, readlines(snp))
        # Default parameters
        frequency_unit = "GHZ"
        parameter_type = "S"
        parameter_format = "MA"
        reference_resistance = 50 # [Ohm]
        number_of_ports = 2 # Required
        number_of_frequencies = 2 # Required
        number_of_noise_frequencies = 0 # Required if noise data is present
        reference = []
        two_port_data_order = ""
        matrix_format = "FULL"
        network_data = Dict()
        noise_data = Dict()
        t_network = []
        t_noise = []

        i = 1
        i = skip_comments(sentences, i)
        m = match(r"^\s{0,}\[VERSION\]\s(2\.0)", sentences[i])[1]
        if m != "2.0"
            throw("Only version 2.0 is currently supported")
        end
        i = skip_comments(sentences, i+1)

        # Parse option line
        freq_unit = match(r"((?:GHZ)|(?:KHZ)|(?:MHZ)|(?:HZ))", sentences[i])
        if !(freq_unit === nothing)
            frequency_unit = freq_unit[1][1:end-1]*"z"
        end

        param_type = match(r"\s(S|Y|Z|H|G)\s", sentences[i])
        if !(param_type === nothing)
            parameter_type = param_type[1]
        end

        param_format = match(r"((?:DB)|(?:MA)|(?:RI))", sentences[i])
        if !(param_format === nothing)
            parameter_format = param_format[1]
        end

        ref_resistance = match(r"R (\d+)", sentences[i])

        if !(ref_resistance === nothing)
            try
                reference_resistance = parse(Int, ref_resistance[1])
            catch
                throw("Incorrect number of reference resistance")
            end
        end

        i = skip_comments(sentences, i+1)

        # Parse keywords that don't appear in a specific order
        while true
            keyword = match(r"^\s*\[(.+)\]", sentences[i])[1]
            if keyword == "NUMBER OF PORTS"
                number_of_ports = parse(Int, match(r"^\[NUMBER OF PORTS\]\s+(\d+)", sentences[i])[1])
            elseif keyword == "REFERENCE"
                reference_resistance = parse(Int, match(r"(\d+)", sentences[i])[1])
            elseif keyword == "TWO-PORT ORDER"
                two_port_data_order = match(r"^\[TWO-PORT DATA ORDER\]\s+((?:12_21)|(?:21_12))", sentences[i])[1]
            elseif keyword == "NUMBER OF FREQUENCIES"
                number_of_frequencies = parse(Int, match(r"^\[NUMBER OF FREQUENCIES\]\s+(\d+)", sentences[i])[1])
            elseif keyword == "NUMBER OF NOISE FREQUENCIES"
                number_of_noise_frequencies = parse(Int, match(r"^\[NUMBER OF NOISE FREQUENCIES\]\s+(\d+)", sentences[i])[1])
            elseif keyword == "REFERENCE"
                # TODO: Think about a more elegant way to do this
                reference = []
                while true
                   for k in split(split(sentences[i], "!")[1])
                        if k != "[REFERENCE]"
                            append!(reference, parse(Float64, k))
                        end
                    end
                    if length(reference) >= number_of_ports
                        break
                    end
                    i = skip_comments(sentences, i+1)
                end
            elseif keyword == "MATRIX FORMAT"
                matrix_format = match(r"((?:FULL)|(?:LOWER)|(?:UPPER))", sentences[i])[1]
            elseif keyword == "NETWORK DATA"
                # This should always come at the end along with network noise data
                network_data_len = number_of_frequencies * ((matrix_format=="FULL") ? 2*number_of_ports^2 + 1 : number_of_ports^2 + number_of_ports + 1)
                i = skip_comments(sentences, i+1)
                while true
                    t_network = vcat(t_network, parse.(Float64, split(split(sentences[i], "!")[1])))
                    (length(t_network) < network_data_len) || break
                    i = skip_comments(sentences, i+1)
                end
            elseif keyword == "NOISE DATA"
                # This should always come at the end along with network data
                noise_data_len = number_of_noise_frequencies * ((matrix_format=="FULL") ? 2*number_of_ports^2 + 1 : number_of_ports^2 + number_of_ports + 1)

                i = skip_comments(sentences, i+1)
                while true
                    t_noise = vcat(t_noise, parse.(Float64, split(split(sentences[i], "!")[1])))
                    (length(t_noise) < noise_data_len) || break
                    i = skip_comments(sentences, i+1)
                end
            end
        try
            i = skip_comments(sentences, i+1)
        catch
            break
        end

        end

        network_data = transform_vector_to_matrix_mode(t_network, String(matrix_format), number_of_frequencies)
        if number_of_noise_frequencies > 0
            noise_data = transform_vector_to_matrix_mode(t_noise, String(matrix_format), number_of_noise_frequencies)
        end

        new(frequency_unit, parameter_type, parameter_format, reference_resistance, number_of_ports,
        number_of_frequencies, number_of_noise_frequencies, reference, two_port_data_order, matrix_format,
        network_data, noise_data)
    end
end

"""
Converts Network/Noise data to complex matrix
...
# Arguments
- `snp::TouchstonerSnP`: provided Touchstone file
- `frequency::Int`: target frequency (must be specified in snp)
- `noise_data::Bool` if set to `true` returns results for Noise Data. Else,
    Network Data are assumed
...
"""
function get_ri_form(snp::TouchstoneSnP, frequency::Float64, noise_data::Bool=false)
data = noise_data ? snp.noise_data : snp.network_data
if !haskey(data, frequency)
    throw("Frequency not specified in data")
end

parameters = data[frequency]
n, m = size(parameters)
if snp.parameter_format == "RI"
    return [to_complex([i, 2j - 1:2j]) for i = 1:n, j = 1:n]
elseif snp.parameter_format == "MA"
   return [to_complex(ma_to_ri(parameters[i, 2j - 1:2j])) for i = 1:n, j = 1:n]
elseif snp.parameter_format == "DB"
    return [to_complex(ma_to_ri(db_to_ma(parameters[i, 2j - 1:2j]))) for i = 1:n, j = 1:n]
end
end

"""
Converts Network/Noise parameters to specified parameters
...
# Arguments
- `snp::TouchstonerSnP`: provided Touchstone file
- `frequency::Int`: target frequency (must be specified in snp)
- `parameter_type`::Str : parameter typ to convert to (source type being the  one specified in SnP)
- `noise_data::Bool` if set to `true` returns results for Noise Data. Else,
    Network Data are assumed
...
"""
function convert_param_type(snp::TouchstoneSnP, frequency::Float64, parameters_type::String, noise_data::Bool=false)
    data = noise_data ? snp.noise_data : snp.network_data
    if !haskey(data, frequency)
        throw("Frequency not specified in data")
    end

    # Convert to R-I format for easier conversion
    parameters = get_ri_form(snp, frequency, noise_data)
    # todo: Mixed-port parameters
    # Also I don't like this at all :P
    if parameters_type == "S"
        if snp.parameter_type == "S"
            return parameters
        elseif snp.parameter_type == "Y"
            return admittance_params_to_scatter(parameters, snp.reference_resistance)
        elseif snp.parameter_type == "Z"
            return impedance_params_to_scatter(parameters, snp.reference_resistance)
        end
    elseif parameters_type == "Y"
        if snp.parameter_type == "S"
            return scatter_params_to_admittance(parameters, snp.reference_resistance)
        elseif snp.parameter_type == "Y"
            return parameters
        elseif snp.parameter_type == "Z"
            return inv(parameters)
        end
    elseif parameters_type == "Z"
        if snp.parameter_type == "S"
            return scatter_params_to_impedance(parameters, snp.reference_resistance)
        elseif snp.parameter_type == "Y"
            return inv(parameters)
        elseif snp.parameter_type == "Z"
            return parameters
    end
end

end


text:
module Math2D

import Base: size, eltype, isempty, min, max, union, intersect, zero, one, typemin, typemax, +, -, *, /

export Vec2, Box, Rect
export size, eltype, len, len2, min, max, isempty, union, intersect, zero, one, normalize, orthogonal, dot, rect, typemin, typemax

struct Vec2{T}
	x::T
	y::T
end

size(::Type{Vec2{T}}) where T = (2,)
size(v::Vec2) = size(typeof(v))
eltype(::Type{Vec2{T}}) where T = T
eltype(v::Vec2) = eltype(typeof(v))

dot(v1::Vec2, v2::Vec2) = v1.x*v2.x + v1.y*v2.y

len2(v::Vec2) = v.x*v.x + v.y*v.y
len(v::Vec2) = sqrt(len2(v))

function normalize(v::Vec2{T}) where T
	l = len(v)
	l < eps(T) ? v : v / l
end

orthogonal(v::Vec2) = Vec2(v.y, -v.x)

zero(::Type{Vec2{T}}) where T = Vec2(zero(T), zero(T))
zero(v::Vec2) = zero(typeof(v))
one(::Type{Vec2{T}}) where T = Vec2(one(T), one(T))
one(v::Vec2) = one(typeof(v))
typemin(::Type{Vec2{T}}) where T = Vec2(typemin(T), typemin(T))
typemin(v::Vec2) = typemin(typeof(v))
typemax(::Type{Vec2{T}}) where T = Vec2(typemax(T), typemax(T))
typemax(v::Vec2) = typemax(typeof(v))

(+)(v::Vec2) = v
(-)(v::Vec2) = Vec2(-v.x, -v.y)

(+)(v1::Vec2, v2::Vec2) = Vec2(v1.x + v2.x, v1.y + v2.y)
(-)(v1::Vec2, v2::Vec2) = Vec2(v1.x - v2.x, v1.y - v2.y)
(*)(v1::Vec2, v2::Vec2) = Vec2(v1.x * v2.x, v1.y * v2.y)
(/)(v1::Vec2, v2::Vec2) = Vec2(v1.x / v2.x, v1.y / v2.y)

(*)(v::Vec2, n::Number) = Vec2(v.x * n, v.y * n)
(*)(n::Number, v::Vec2) = Vec2(n * v.x, n * v.y)
(/)(v::Vec2, n::Number) = Vec2(v.x / n, v.y / n)
(/)(n::Number, v::Vec2) = Vec2(n / v.x, n / v.y)

broadcast(::typeof(<), v1::Vec2, v2::Vec2) = v1.x < v2.x && v1.y < v2.y

min(v::Vec2) = v
min(v1::Vec2, v2::Vec2) = Vec2(min(v1.x, v2.x), min(v1.y, v2.y))
min(v::Vec2, vs::Vec2...) = min(v, min(vs...))

max(v::Vec2) = v
max(v1::Vec2, v2::Vec2) = Vec2(max(v1.x, v2.x), max(v1.y, v2.y))
max(v::Vec2, vs::Vec2...) = max(v, max(vs...))

mutable struct Box{V}
	min::V
	max::V
end

size(b::Box) = b.max - b.min
isempty(b::Box) = b.max .< b.min

union(bs::Box...) = Box(min([b.min for b in bs]), max([b.max for b in bs]))
intersect(bs::Box...) = Box(max([b.min for b in bs]), min([b.max for b in bs]))

const Rect{T} = Box{Vec2{T}}

rect(T, minX, minY, maxX, maxY) = Box(Vec2{T}(minX, minY), Vec2{T}(maxX, maxY))
rect(T) = Box(typemax(Vec2{T}), typemin(Vec2{T}))

end


text:
"""
    Third experimental scenario.

    In this experimental scenario, we varied each node threshold proportionally
    - from 0.2 to 0.8 - to the degree of the node,
    and fixed each hyperedge threshold to a random value between 1 and its size.
"""

using Pkg, Distributed
#addprocs(4)
Pkg.activate(".")
@everywhere using Distributed, Pkg
@everywhere Pkg.activate(".")
using LTMSim, SimpleHypergraphs, Random, Serialization
@everywhere using LTMSim, SimpleHypergraphs, Random, Serialization

#
#
#
project_path = dirname(pathof(LTMSim))
data_path = joinpath(project_path, "..", "data", "hgs")
res_path = joinpath(project_path, "..", "res", "journal", "propV_randE_w_wo_opt.data")

hg_files = readdir(data_path)
hgs = [prune_hypergraph!(hg_load(joinpath(data_path, hg_file))) for hg_file in hg_files]
hg_names = [split(fname, ".")[1] for fname in hg_files]

algorithms = [
    "BinarySearch(H)", "Greedy([H]₂)", "Greedy(H)", "SubTSS(H)",
    "BinarySearch(H)-noOpt", "Greedy([H]₂)-noOpt", "Greedy(H)-noOpt", "SubTSS(H)-noOpt"
]

data = Dict{String, Dict{String, Vector{Vector{Int}}}}()

for algo in algorithms
    for hg in hg_names
        push!(
            get!(data, algo, Dict{String, Vector{Vector{Int}}}()),
            hg => Vector{Vector{Int}}()
        )
    end
end

nvalues = range(0.2, stop=0.8, step=0.1)
runs = 50

for (index, h) in enumerate(hgs)
    println("Index=$index, $(hg_files[index])")

    for n=nvalues
        println("n=$n")

        results = @distributed (append!) for run=1:runs
            println("run=$run at proc $(myid())")
            Random.seed!(run)

            metaV = proportionalMetaV(h, n)
            metaE = randMetaE(h)

            # with optimization
            r1 = greedy_tss_2section(h,metaV,metaE; opt=true)
            r2 = bisect(h,metaV,metaE; opt=true)
            r3 = greedy_tss(h,metaV,metaE; opt=true)
            r4 = sub_tss(h,metaV,metaE; opt=true)

            # without optimization
            r5 = greedy_tss_2section(h,metaV,metaE; opt=false)
            r6 = bisect(h,metaV,metaE; opt=false)
            r7 = greedy_tss(h,metaV,metaE; opt=false)
            r8 = sub_tss(h,metaV,metaE; opt=false)

            [(r1,r2,r3,r4,r5,r6,r7,r8)]
        end

        push!(data["Greedy([H]₂)"][hg_names[index]], [r[1] for r in results])
        push!(data["BinarySearch(H)"][hg_names[index]], [r[2] for r in results])
        push!(data["Greedy(H)"][hg_names[index]], [r[3] for r in results])
        push!(data["SubTSS(H)"][hg_names[index]], [r[4][1] for r in results])

        push!(data["Greedy([H]₂)-noOpt"][hg_names[index]], [r[5] for r in results])
        push!(data["BinarySearch(H)-noOpt"][hg_names[index]], [r[6] for r in results])
        push!(data["Greedy(H)-noOpt"][hg_names[index]], [r[7] for r in results])
        push!(data["SubTSS(H)-noOpt"][hg_names[index]], [r[8][1] for r in results])

        println()
    end
    println()
end

serialize(res_path, data)


text:
## Select keyword arguments from lists
const KEYS_GEOM_ATTRIBUTES = [:accelerate, :algorithm, :alpha, :baseline, :clabels, :color, :headsize,
    :horizontal, :label, :linecolor, :linewidth, :markercolor, :markersize, :shadelines, :spec, :stair_position, :xform]
const KEYS_PLOT_ATTRIBUTES = [:backgroundalpha, :backgroundcolor, :colorbar, :colormap, :location, :hold, :overlay_axes, :radians, :ratio, :scheme, :subplot, :title,
    :xflip, :xlabel, :xlim, :xlog, :xticklabels, :yflip, :ylabel, :ylim, :ylog, :yticklabels, :zflip, :zlabel, :zlim, :zlog]

geom_attributes(; kwargs...) = filter(p -> p.first ∈ KEYS_GEOM_ATTRIBUTES, kwargs)
plot_attributes(; kwargs...) = filter(p -> p.first ∈ KEYS_PLOT_ATTRIBUTES, kwargs)

_setargs_default(f, args...; kwargs...) = (args, kwargs)

"""
    @plotfunction(fname, options...)

Macro to create plotting functions. E.g. `@plotfunction plot` creates two
functions:

    * `plot!(f::Figure, args...; kwargs...)`
    * `plot(args...; kwargs...)`

The first of those functions (the one whose name ends with an exclamation)
edits the figure given as first argument, replacing its last plot by a new
one. The second function (the one without exclamation) creates the plot in the
current figure. How those functions work depends on the options that are passed
after the function name to the macro. Those options are expressed in the fashion
of keyword argments, i.e. as `key = value`, and they can be the following:

* **`geom`**: a `Symbol` with the name of the kind of the `Geometry` that is created.
* **`axes`**: a `Symbol` with the name of the kind of the `Axes` that are created.
* **`plotkind`**: a `Symbol` with the name of the plot kind (only needed as meta-data).
    If this option is not given, the name of the function is used by default.
* **`setargs`**: a function that takes the positional and keyword arguments that are
    passed to the functions, and returns: (a) a tuple of positional arguments
    to be passed to the function [`geometries`](@ref)), and
    (b) the set of keyword arguments that are passed to the constructor of
    geometries, axes, and the plot object. If `setargs` is not defined, the
    positional and keyword arguments are returned untransformed.
* **`kwargs`**: a named tuple with extra keyword arguments that are passed to
    the constructors of geometries, axes and the plot object.
* **`docstring`**: the documentation string that will be assigned to the plotting function.
"""
macro plotfunction(fname, options...)
    # Parse options - minimum geom and axes
    dict_op = Dict{Symbol, Any}()
    for op in options
        if typeof(op) <: Expr && op.head ∈ (:(=), :kw)
            dict_op[op.args[1]] = op.args[2]
        end
    end
    if !haskey(dict_op, :geom)
        throw(ArgumentError("`geom` not specified"))
    end
    if !haskey(dict_op, :axes)
        throw(ArgumentError("`axes` not specified"))
    end
    # Define functions
    geom_k = dict_op[:geom]
    axes_k = dict_op[:axes]
    setargs_fun = get(dict_op, :setargs, _setargs_default)
    plotkind = get(dict_op, :kind, Symbol(fname))
    def_kwargs = get(dict_op, :kwargs, NamedTuple())
    fname! = Symbol(fname, :!)
    expr = quote
        function $(fname!)(f::GRUtils.Figure, args...; kwargs...)
            kwargs = (; $(def_kwargs)..., kwargs...)
            p = GRUtils.currentplot(f)
            if haskey(kwargs, :hold)
                holdstate = kwargs[:hold]
            else
                holdstate = get(p.attributes, :hold, false)
            end
            if holdstate
                # Keep all attributes
                kwargs = (; p.attributes..., kwargs...)
                args, kwargs = $setargs_fun(f, args...; kwargs...)
                geoms = [p.geoms; GRUtils.geometries(Symbol($geom_k), args...; GRUtils.geom_attributes(;kwargs...)...)]
            else
                # Only keep previous subplot
                kwargs = (subplot = p.attributes[:subplot], kwargs...)
                args, kwargs = $setargs_fun(f, args...; kwargs...)
                geoms = GRUtils.geometries(Symbol($geom_k), args...; GRUtils.geom_attributes(;kwargs...)...)
            end
            axes = GRUtils.Axes(Symbol($axes_k), geoms; kwargs...)
            GRUtils.makeplot!(f.plots[end], axes, geoms; kind=$plotkind, GRUtils.plot_attributes(; kwargs...)...)
            return f
        end
        $fname(args...; kwargs...) = $fname!(GRUtils.gcf(), args...; kwargs...)
    end
    # Add docstrings if available
    if haskey(dict_op, :docstring)
        push!(expr.args, quote @doc $(dict_op[:docstring]) $fname end)
    end
    esc(expr)
end

function _setargs_line(f, args...; kwargs...)
    if typeof(args[end]) <: AbstractString
        kwargs = (spec=args[end], kwargs...)
        args = args[1:end-1]
    end
    return (args, kwargs)
end

@plotfunction(plot, geom = :line, axes = :axes2d, setargs=_setargs_line, kind = :line, docstring="""
    plot(x[, y, spec; kwargs...])
    plot(x1, y1, x2, y2...; kwargs...)
    plot(x1, y1, spec1...; kwargs...)

Draw one or more line plots.

Lines are defined by the `x` and `y` coordinates of the connected points, given as
numeric vectors, and optionally the format string `spec` that defines the line
and marker style and color as in
[matplotlib](https://matplotlib.org/3.1.1/api/_as_gen/matplotlib.pyplot.plot.html).

The `y` vector can be replaced by a callable that defines the Y coordinates as a
function of the X coordinates.

Multiple lines can be defined by pairs of `x` and `y` coordinates (and optionally
their format strings), passed sequentially as arguments of `plot`.
Alternatively, if various lines have the same X coordinates, their Y values can
be grouped as columns in a matrix.

If no `spec` is given, the series will be plotted as solid lines with a
predefined sequence of colors.

Additionally, specifications of lines and markers can be defined by keyword arguments:

* `linewidth`: line width scale factor.
* `markersize`: marker size scale factor.
* `linecolor`: hexadecimal RGB color code for the line.
* `markercolor`: hexadecimal RGB color code for the markers.

This function can receive a single numeric vector or matrix, which will be
interpreted as the Y coordinates; in such case the X coordinates will be a
sequence of integers starting at 1.

# Examples

```julia
$(_example("plot"))
```
""")

function _setargs_stair(f, args...; kwargs...)
    stair_position_str = get(kwargs, :where, "mid")
    if stair_position_str == "mid"
        stair_position = 0.0
    elseif stair_position_str == "post"
        stair_position = 1.0
    elseif stair_position_str == "pre"
        stair_position = -1.0
    else
        throw(ArgumentError("""`where` must be one of `"mid"`, `"pre"` or `"post"`"""))
    end
    return _setargs_line(f, args...; stair_position=stair_position, kwargs...)
end

@plotfunction(stair, geom = :stair, axes = :axes2d, setargs=_setargs_stair, docstring="""
    stair(x[, y, spec; kwargs...])
    stair(x1, y1, x2, y2...; kwargs...)
    stair(x1, y1, spec1...; kwargs...)

Draw one or more staircase or step plots.

The coordinates and format of the stair outlines are defined as for line plots
(cf. [`plot`](@ref)).

Additionally, the keyword argument `where` can be used to define where the "stairs"
(vertical discontinuities between Y values) shoud be placed:

* `where = "pre"` to make the steps stop at each point (`x[i]`, `y[i]`),
    starting at the previous `x` coordinate except for the first point.
* `where = "post"` to make the steps start at each point (`x[i]`, `y[i]`),
    stopping at the next `x` coordinate except for the last point.
* `where = "mid"` (default) to make the steps go through each point (`x[i]`, `y[i]`)
    starting and ending in the middle of the surrounding x-intervals,
    except for the first and last points.

# Examples

```julia
$(_example("stair"))
```
""")

function _setargs_stem(f, args...; baseline=0.0, kwargs...)
    args, kwargs = _setargs_line(f, args...; kwargs...)
    x = collect(args[1])
    y = collect(args[2])
    args = args[3:end]
    prepend!(x, [minimum(x), maximum(x)])
    prepend!(y, float.([baseline, baseline]))
    return ((x, y, args...), kwargs)
end

@plotfunction(stem, geom = :stem, axes = :axes2d, setargs=_setargs_stem, docstring="""
    stem(x[, y, spec; kwargs...])
    stem(x1, y1, x2, y2...; kwargs...)
    stem(x1, y1, spec1...; kwargs...)

Draw a stem plot

The coordinates and format of the stems and markers are defined as for line plots
(cf. [`plot`](@ref)).

Additionally, the keyword argument `baseline` can be used to define the
Y coordinate where stems should start from.

# Examples

```julia
$(_example("stem"))
```
""")

function _setargs_errorbar(f, x, y, args...; kwargs...)
    if typeof(args[end]) <: AbstractString
        kwargs = (; spec=args[end], kwargs...)
        return _setargs_errorbar(f, x, y, args[1:end-1]...; kwargs...)
    end
    # Define bar size
    if length(args) == 0
        throw(ArgumentError("errorbar sizes not defined"))
    else
        low = args[1]
        high = (length(args)==1) ? low : args[2]
    end
    if length(low) == 1
        low = repeat(low, length(x))
    end
    if length(high) == 1
        high = repeat(high, length(x))
    end
    # Cap width
    horizontal = get(kwargs, :horizontal, false)
    if haskey(kwargs, :capwidth)
        w = kwargs[:capwidth]/2
    elseif horizontal
        w = 0.015 * (maximum(y) - minimum(y))
    else
        w = 0.015 * (maximum(x) - minimum(x))
    end
    # Coordinates
    if horizontal
        x3 = (x' .+ [-vec(low)'; zeros(1, length(x)); vec(high)'])[:]
        y3 = (y' .+ [-w; 0.0; w])[:]
    else
        x3 = (x' .+ [-w; 0.0; w])[:]
        y3 = (y' .+ [-vec(low)'; zeros(1, length(y)); vec(high)'])[:]
    end
    return ((x3, y3), kwargs)
end

@plotfunction(errorbar, geom = :errorbar, axes = :axes2d, setargs=_setargs_errorbar, docstring="""
    errorbar(x, y, err[, spec; kwargs...])
    errorbar(x, y, errlow, errhigh[, spec; kwargs...])

Draw a series of error bars.

Error bars are defined by their `x` and `y` coordinates, and the size of the
error bars at either of each `(x, y)` point. For symmetric error bars,
only a vector `err` is required, such that their total size will be `2 .* err`.
For asymmetric error bars, two vectors `errlow` and `errhigh` are required,
such that the size of the error bars will be `errlow .+ errhigh`.

The optional format string `spec` defines the style and color of the lines
of error bars and the markers at `(x, y)`, as in
[matplotlib](https://matplotlib.org/3.1.1/api/_as_gen/matplotlib.pyplot.plot.html).
If no `specs` are given, the error bars will be plotted as solid lines with a
predefined sequence of colors, without markers.

Additionally, the following keyword arguments can be used to modify the aspect
of the error bars:

* `linewidth::Float64`: line width scale factor.
* `markersize::Float64`: marker size scale factor.
* `linecolor`: hexadecimal RGB color code for the line.
* `horizontal::Bool`: set it to `true` to draw horizontal error bars).
* `capwidth`: fixed value of the width of the bar "caps", in units of
    the X axis (or Y axis if `horizontal` is `true`). If it is not given,
    the cap width will be automatically adjusted to 0.3 times the mean
    separation between data points.

# Examples

```julia
$(_example("errorbar"))
```
""")

function _setargs_polar(f, x, y, args...; kwargs...)
    if !get(kwargs, :radians, true)
        x = collect(x) .* (π/180)
    end
    _setargs_line(f, x, y, args...; kwargs...)
end

@plotfunction(polar, geom = :polarline, axes = :polar, setargs=_setargs_polar,
kwargs = (ratio=1.0,), docstring="""
    polar(angle, radius[, spec; kwargs...])

Draw one or more polar plots.

The first coordinate the represents the angle, and the second the
radius of the line points. The rest is defined as for line plots,
except that the first variable (`angle`) is always required.
(cf. [`plot`](@ref)).

The first variable is by default considered to be radians, and the angular
labels of the grid are shown as factors of π. Use the keyword argument
`radians = false` to pass and show angles in degrees.

!!! note

    Logarithmic and reversed scales ar disabled in polar plots

# Examples

```julia
$(_example("polar"))
```
""")

# Recursive call in case of multiple x-y pairs
for fun = [:plot!, :stair!, :stem!, :polar!]
    @eval function $fun(f::Figure, x, y, u, v, args...; kwargs...)
        holdstate = get(currentplot(f).attributes, :hold, false)
        if typeof(u) <: AbstractString
            $fun(f, x, y, u; kwargs...)
            hold!(currentplot(f), true)
            $fun(f, v, args...; kwargs...)
        else
            $fun(f, x, y; kwargs...)
            hold!(currentplot(f), true)
            $fun(f, u, v, args...; kwargs...)
        end
        hold!(currentplot(f), holdstate)
        return f
    end
end

@plotfunction(scatter, geom = :scatter, axes = :axes2d, kwargs=(colorbar=true,),
docstring="""
    scatter(x[, y, size, color; kwargs...])

Draw a scatter plot.

Points are defined by their `x` and `y` coordinates, given as numeric vectors.
Additionally, values for markers' `size` and `color` can be provided.
Size values will determine the marker size in percent of the regular size,
and color values will be used in combination with the current colormap.

The last variable can be replaced by a callable that defines it as a
function of the previous variables.

This function can receive a single numeric vector or matrix, which will be
interpreted as the Y coordinates; in such case the X coordinates will be a
sequence of integers starting at 1.

# Examples

```julia
$(_example("scatter"))
```
""")

# horizontal and vertical coordinates of bar edges
function barcoordinates(heights; barwidth=0.8, baseline=0.0, kwargs...)
    n = length(heights)
    halfw = barwidth/2
    wc = zeros(2n)
    hc  = zeros(2n)
    for (i, value) in enumerate(heights)
        wc[2i-1] = i - halfw
        wc[2i]   = i + halfw
        hc[2i-1] = baseline
        hc[2i]   = value
    end
    (wc, hc)
end

function barcoordinates(heights::AbstractMatrix; barposition="grouped", kwargs...)
    if barposition == "grouped"
        groupedbars(heights; kwargs...)
    elseif barposition == "stacked"
        stackedbars(heights; kwargs...)
    else
        throw(ArgumentError("""`barposition` must be `"grouped"` or `"stacked"`"""))
    end
end

# stacked bar coordinates
function stackedbars(heights; barwidth=0.8, baseline=0.0, kwargs...)
    n, m = size(heights)
    halfw = barwidth/2
    wc = zeros(2n, m)
    hc  = zeros(2n, m)
    bases = repeat([baseline], n)
    for c=1:m, r=1:n
        value = heights[r, c]
        wc[2r-1, c] = r - halfw
        wc[2r, c]   = r + halfw
        hc[2r-1, c] = bases[r]
        hc[2r, c]   = bases[r] + value
        bases[r] = hc[2r, c]
    end
    (wc, hc)
end

# grouped bar coordinates
function groupedbars(heights; barwidth=0.8, baseline=0.0, kwargs...)
    n, m = size(heights)
    halfw = barwidth/2m
    wc = zeros(2n, m)
    hc  = zeros(2n, m)
    for c=1:m, r=1:n
        value = heights[r, c]
        offset = (2c - 1 - m) * halfw
        wc[2r-1, c] = r - halfw + offset
        wc[2r, c]   = r + halfw + offset
        hc[2r-1, c] = baseline
        hc[2r, c]   = value
    end
    (wc, hc)
end

function _setargs_bar(f, labels, heights; fillcolor=nothing, horizontal=false, kwargs...)
    if fillcolor ≠ nothing # deprecate?
        kwargs = (; color=fillcolor, kwargs...)
    end
    wc, hc = barcoordinates(heights; kwargs...)
    if horizontal
        args = (hc, wc)
        tickoptions = (yticks = (1,1), yticklabels = string.(labels))
    else
        args = (wc, hc)
        tickoptions = (xticks = (1,1), xticklabels = string.(labels))
    end
    return (args, (; tickoptions..., kwargs...))
end

function _setargs_bar(f, heights; kwargs...)
    n = length(heights)
    _setargs_bar(f, string.(1:n), heights; kwargs...)
end

@plotfunction(barplot, geom = :bar, axes = :axes2d, setargs=_setargs_bar, docstring="""
    bar(labels, heights; kwargs...)
    bar(heights; kwargs...)

Draw a bar plot.

If no specific labels are given, the bars are labelled with integer
numbers starting from 1.

If `heights` is a matrix, each column is taken as a different set of data,
which are represented as bars of different colors.
Use the keyword argument `barposition` with the values `"grouped"` (default)
or `"stacked"` to control if the bars are positioned side by side or
stacked on top of the previous series.

The keyword arguments `barwidth`, `baseline` and `horizontal`
can also be used to modify the aspect of the bars, which by default is:

* `barwidth = 0.8` (80% of the separation between bars).
* `baseline = 0.0` (bars starting at zero).
* `horizontal = false` (vertical bars)

The color of the bars is selected automatically, unless
a specific hexadecimal RGB color code is given through
the keyword argument `color`.

# Examples

```julia
$(_example("barplot"))
```
""")

# Coordinates of the bars of a histogram of the values in `x`
function hist(x, nbins=0, baseline=0.0)
    if nbins <= 1
        nbins = round(Int, 3.3 * log10(length(x))) + 1
    end

    xmin, xmax = extrema(x)
    edges = range(xmin, stop = xmax, length = nbins + 1)
    counts = zeros(nbins)
    buckets = Int[max(2, min(searchsortedfirst(edges, xᵢ), length(edges)))-1 for xᵢ in x]
    for b in buckets
        counts[b] += 1
    end
    wc = zeros(2nbins)
    hc  = zeros(2nbins)
    for (i, value) in enumerate(counts)
        wc[2i-1] = edges[i]
        wc[2i]   = edges[i+1]
        hc[2i-1] = baseline
        hc[2i]   = value
    end
    (wc, hc)
end

function _setargs_hist(f, x; nbins = 0, fillcolor=nothing, horizontal = false, kwargs...)
    if fillcolor ≠ nothing # deprecate?
        kwargs = (; color=fillcolor, kwargs...)
    end
    # Define baseline - 0.0 by default, unless using log scale
    if get(kwargs, :ylog, false) || horizontal && get(kwargs, :xlog, false)
        baseline = 1.0
    else
        baseline = 0.0
    end
    wc, hc = hist(x, nbins, baseline)
    args = horizontal ? (hc, wc) : (wc, hc)
    return (args, kwargs)
end

@plotfunction(histogram, geom = :bar, axes = :axes2d, kind = :hist, setargs = _setargs_hist,
docstring="""
    histogram(data; kwargs...)

Draw a histogram of `data`.

The following keyword arguments can be supplied:

* `nbins`: Number of bins; by default, or if a number smaller than 1 is given,
    the number of bins is computed as `3.3 * log10(n) + 1`,  with `n` being the
    number of elements in `data`.
* `horizontal`: whether the histogram should be horizontal (`false` by default).
* `color`: hexadecimal RGB color code for the bars.

!!! note

    If the vertical axis (or the horizontal axis if `horizontal == true`) is set
    in logarithmic scale, the bars of the histogram will start at 1.0.

# Examples

```julia
$(_example("histogram"))
```
""")

function _setargs_polarhist(f, x; kwargs...)
    if get(kwargs, :fullcircle, false)
        minval, maxval = extrema(x)
        x = (collect(x) .- minval) .* (2π / (maxval - minval))
    elseif !get(kwargs, :radians, true)
        x = collect(x) .* π / 180
    end
    _setargs_hist(f, x; kwargs..., horizontal=false)
end

@plotfunction(polarhistogram, geom = :polarhist, axes = :polar, kind = :polarhist,
setargs = _setargs_polarhist, kwargs = (ratio=1.0, overlay_axes=true), docstring="""
    polarhistogram(data; kwargs...)

Draw a polar histogram of angle values contained in `data`.

The following keyword arguments can be supplied:

* `nbins`: Number of bins; by default, or if a number smaller than 1 is given,
    the number of bins is computed as `3.3 * log10(n) + 1`,  with `n` being the
    number of elements in `data`.
* `radians`: Set this argument to `false` to pass and show the angles as degrees.
    By default, `data` is assumed to be radians and the angular labels of the
    grid are presented as factors of π.
* `fullcircle`: Set this argument to `true` to scale the angular coordinates of
    the histogram and make the bars span over the whole circle.
* `color`: hexadecimal RGB color code for the bars.

!!! note

    Logarithmic and reversed scales ar disabled in polar plots

# Examples

```julia
$(_example("polarhistogram"))
```
""")

@plotfunction(plot3, geom = :line3d, axes = :axes3d, kwargs = (ratio=1.0,), setargs=_setargs_line, docstring="""
    plot3(x, y, z[, spec; kwargs...])
    plot3(x1, y1, z1, x2, y2, z2...; kwargs...)
    plot3(x1, y1, z1, spec1...; kwargs...)

Draw one or more three-dimensional line plots.

Lines are defined by the `x`, `y` and `z` coordinates of the connected points,
given as numeric vectors, and optionally the format string `spec` that defines
the line and marker style and color as in
[matplotlib](https://matplotlib.org/3.1.1/api/_as_gen/matplotlib.pyplot.plot.html).

The `z` vector can be replaced by a callable that defines the Z coordinates as a
function of the X and Y coordinates.

Multiple lines can be defined by triplets of `x`, `y` and `z` coordinates (and
optionally their format strings), passed sequentially as arguments of `plot3`.

If no `specs` are given, the series will be plotted as solid lines with a
predefined sequence of colors.

Additionally, specifications of lines and markers can be defined by keyword arguments:

* `linewidth`: line width scale factor.
* `markersize`: marker size scale factor.
* `linecolor`: hexadecimal RGB color code for the line.
* `markercolor`: hexadecimal RGB color code for the markers.

# Examples

```julia
$(_example("plot3"))
```
""")

function plot3!(f::Figure, x, y, z, u, v, args...; kwargs...)
    holdstate = get(currentplot(f).attributes, :hold, false)
    if typeof(u) <: AbstractString
        plot3!(f, x, y, z, u; kwargs...)
        hold!(currentplot(f), true)
        plot3!(f, v, args...; kwargs...)
    else
        plot3!(f, x, y, z; kwargs...)
        hold!(currentplot(f), true)
        plot3!(f, u, v, args...; kwargs...)
    end
    hold!(currentplot(f), holdstate)
    return f
end

# Quiver:

# Arrow head corners for a given relative size
# The arrow flanks are 22.5º away from (u,v)
# The point of the arrow is at (0, 0)
# Arrows are scaled by xf, yf in the X, Y axes
function arrowhead(u, v, hsize, xf, yf)
    x1 = -0.1848 * hsize * u - 0.0765 * hsize * v *xf
    x2 = x1 + 0.1531 * hsize * v *xf
    y1 = -0.1848 * hsize * v + 0.0765 * hsize * u * yf
    y2 = y1 - 0.1531 * hsize * u * yf
    return (x1, x2, y1, y2)
end

# Calculate the x and y factors taking into account
# the ranges of the area covered by the arrows
function quiver_xyfactors(x, y, u, v)
    xrange = max(maximum(x), maximum(x.+u)) - min(minimum(x), minimum(x.+u))
    yrange = max(maximum(y), maximum(y.+v)) - min(minimum(y), minimum(y.+v))
    xrange = max(xrange, 0.01*yrange)
    yrange = max(yrange, 0.01*xrange)
    xf = xrange/yrange
    yf = yrange/xrange
    return xf, yf
end

function _setargs_quiver(f, x, y, u, v, spec=""; arrowscale=1.0, kwargs...)
    n = length(x)
    xa = zeros(2n)
    ya = zeros(2n)
    for i = 1:n
        j = 2i - 1
        xa[j] = x[i]
        xa[j+1] = x[i] + arrowscale*u[i]
        ya[j] = y[i]
        ya[j+1] = y[i] + arrowscale*v[i]
    end
    if !isempty(spec)
        kwargs = (; spec=spec, kwargs...)
    end
    return ((xa, ya), kwargs)
end

@plotfunction(quiver, geom = :quiver, axes = :axes2d, setargs = _setargs_quiver, docstring="""
    quiver(x, y, u, v[, spec; kwargs...])

Draw a vector field at points `(x,y)` with arrows
of size `(u,v)`.

The style and color of the arrow lines and -- optionally --
the markers drawn at the `(x,y)` points can be configured through the
format string `spec` and keyword arguments, as in [`plot`](@ref) and other functions.

Use the keyword arguments `arrowscale` and `headsize` to give a scale factor
for the size of the arrows and their heads.

# Examples

```julia
$(_example("quiver"))
```
""")

function _setargs_quiver3(f, x, y, z, u, v, w, spec=""; arrowscale=1.0, kwargs...)
    n = length(x)
    xa = zeros(2n)
    ya = zeros(2n)
    za = zeros(2n)
    for i = 1:n
        j = 2i - 1
        xa[j] = x[i]
        xa[j+1] = x[i] + arrowscale*u[i]
        ya[j] = y[i]
        ya[j+1] = y[i] + arrowscale*v[i]
        za[j] = z[i]
        za[j+1] = z[i] + arrowscale*w[i]
    end
    if !isempty(spec)
        kwargs = (; spec=spec, kwargs...)
    end
    return ((xa, ya, za), kwargs)
end

@plotfunction(quiver3, geom = :quiver3, axes = :axes3d, setargs = _setargs_quiver3,
kwargs = (ratio=1.0,), docstring="""
    quiver3(x, y, z, u, v, w[, spec; kwargs...])

Draw a three-dimensional vector field at points `(x,y,z)` with arrows
of size `(u,v,w)`.

The style and color of the arrow lines and -- optionally --
the markers drawn at the `(x,y,z)` points can be configured through the
format string `spec` and keyword arguments, as in [`plot`](@ref) and other functions.

Use the keyword argument `arrowscale` to give a scale factor for the size of the arrows.

!!! note

    Unlike in the case of 2d [`quiver`](@ref) plots, the 3-D arrows of
    `quiver3` do not have heads.

# Examples

```julia
$(_example("quiver3"))
```
""")

_setargs_scatter3(f, x, y, z; kwargs...) = ((x,y,z), kwargs)
_setargs_scatter3(f, x, y, z, c; kwargs...) = ((x,y,z,c), (;colorbar=true, kwargs...))

@plotfunction(scatter3, geom = :scatter3, axes = :axes3d, setargs = _setargs_scatter3,
kwargs = (ratio=1.0,), docstring="""
    scatter3(x, y, z[, color; kwargs...])

Draw a three-dimensional scatter plot.

Points are defined by their `x`, `y` and `z` coordinates, given as numeric vectors.
Additionally, values for markers' `color` can be provided, which will be used in
combination with the current colormap.

The last variable can be replaced by a callable that defines it as a
function of the previous variables.

# Examples

```julia
$(_example("scatter3"))
```
""")

# Contour arguments for different inputs:
# Coordinates (x, y, z) and countour line levels
function _setargs_contour(f, x, y, z, h; kwargs...)
    if length(x) == length(y) == length(z)
        x, y, z = GR.gridit(vec(x), vec(y), vec(z), 200, 200)
    else
        z = z'
    end
    if get(kwargs, :colorbar, true)
        majorlevels = get(kwargs, :majorlevels, 0)
        clabels = float(1000 + majorlevels)
        kwargs = (; colorbar = true, clabels = clabels, kwargs...)
    else
        majorlevels = get(kwargs, :majorlevels, 1)
        clabels = float(majorlevels)
        kwargs = (; clabels = clabels, kwargs...)
    end
    if majorlevels ≠ 0
        kwargs = (; kwargs..., ratio = 1.0)
    end
    return ((vec(x), vec(y), vec(z), vec(h)), kwargs)
end

# Coordinates (x, y, z) with contour lines automatically calculated
function _setargs_contour(f, x, y, z; levels = 20, kwargs...)
    (x, y, z, _), kwargs = _setargs_contour(f, x, y, z, []; kwargs...)
    levels = Int(levels)
    zmin, zmax = get(kwargs, :zlim, (_min(z), _max(z)))
    hmin, hmax = GR.adjustrange(zmin, zmax)
    h = range(hmin, stop = hmax, length = levels + 1)
    return ((x, y, z, h), kwargs)
end

# z values are calculated from a function
function _setargs_contour(f, x, y, fz::Function, args...; kwargs...)
    z = fz.(transpose(vec(x)), vec(y))
    _setargs_contour(f, x, y, z, args...; kwargs...)
end

# z values alone, without x,y coordinates
function _setargs_contour(f, z; kwargs...)
    sz = size(z)
    _setargs_contour(f, 1:sz[2], 1:sz[1], z; kwargs...)
end

function _setargs_contour(f, z, h; kwargs...)
    sz = size(z)
    _setargs_contour(f, 1:sz[2], 1:sz[1], z, h; kwargs...)
end

@plotfunction(contour, geom = :contour, axes = :axes3d, setargs = _setargs_contour,
kwargs = (rotation=0, tilt=90), docstring="""
    contour([x, y,] z; kwargs...)

Draw a contour plot.

The current colormap is used to display a either a series of
points or a two-dimensional array as a contour plot. It can receive one
of the following:

- `x` values, `y` values and `z` values.
- *M* sorted values of the `x` axis, *N* sorted values of the `y` axis,
    and a set of `z` values on a *N*×*M* grid.
- *M* sorted values of the `x` axis, *N* sorted values of the `y` axis,
    and a callable to determine `z` values.
- `z` values on a *N*×*M* grid, with the *x* and *y* axes
    defined as the indices of the columns and rows of the grid, respectively.

If a series of points is passed to this function, their values will be
interpolated on a grid. For grid points outside the convex hull of the
provided points, a value of 0 will be used.

Contour lines are colored by default, using the current colormap as color scale.
Colored lines can be disabled by removing the color bar (with keyword argument
`colorbar == false`).

The following keyword arguments can be provided to set the number and aspect of
the contour lines:

* `levels::Int`, the number of contour lines that will be fitted to the data
    (20 by default).
* `majorlevels::Int`, the number of levels between labelled contour lines
    (no labels by default for colored lines, all lines labelled by default
    if color is removed).

# Examples

```julia
$(_example("contour"))
```
""")

@plotfunction(contourf, geom = :contourf, axes = :axes3d, setargs = _setargs_contour,
kwargs = (rotation=0, tilt=90, tickdir=-1), docstring="""
    contourf([x, y,] z; kwargs...)

Draw a filled contour plot.

The current colormap is used to display a either a series of
points or a two-dimensional array as a filled contour plot. It can receive one
of the following:

- `x` values, `y` values and `z` values.
- *M* sorted values of the `x` axis, *N* sorted values of the `y` axis,
    and a set of `z` values on a *N*×*M* grid.
- *M* sorted values of the `x` axis, *N* sorted values of the `y` axis,
    and a callable to determine `z` values.
- `z` values on a *N*×*M* grid, with the *x* and *y* axes
    defined as the indices of the columns and rows of the grid, respectively.

If a series of points is passed to this function, their values will be
interpolated on a grid. For grid points outside the convex hull of the
provided points, a value of 0 will be used.

The following keyword arguments can be provided to set the number and aspect of
the contour lines:

* `levels::Int`, the number of contour lines that will be fitted to the data
    (20 by default).
* `majorlevels::Int`, the number of levels between labelled contour lines
    (no labels by default).

# Examples

```julia
$(_example("contourf"))
```
""")

_setargs_tricont(f, x, y, z, h; kwargs...) = ((x, y, z, h), kwargs...)

function _setargs_tricont(f, x, y, z; levels = 20, kwargs...)
    levels = Int(levels)
    zmin, zmax = get(kwargs, :zlim, (_min(z), _max(z)))
    hmin, hmax = GR.adjustrange(zmin, zmax)
    h = range(hmin, stop = hmax, length = levels)
    return ((x, y, z, h), kwargs)
end

function _setargs_tricont(f, x, y, fz::Function, args...; kwargs...)
    z = fz.(transpose(vec(x)), vec(y))
    _setargs_tricont(f, x, y, z, args...; kwargs...)
end

@plotfunction(tricont, geom = :tricont, axes = :axes3d, setargs = _setargs_tricont,
kwargs = (colorbar=true, rotation=0, tilt=90), docstring="""
    tricont(x, y, z; kwargs...)

Draw a triangular contour plot.

The current colormap is used to display a series of points
as a triangular contour plot. `z` values are interpolated between `x` and `y`
values through [Delaunay triangulation](http://mathworld.wolfram.com/DelaunayTriangulation.html).

The number of contour lines can be set by the keyword argument `levels`
(by default `levels = 20`).

!!! note

    If the series of points is concave, there may be interpolation artifacts on
    the edges of the plot, as the interpolation may occur in very acute triangles.

# Examples

```julia
$(_example("tricont"))
```
""")

function _setargs_surface(f, x, y, z; accelerate = true, kwargs...)
    if length(x) == length(y) == length(z)
        x, y, z = GR.gridit(vec(x), vec(y), vec(z), 200, 200)
    else
        z = z'
    end
    accelerate = Bool(accelerate) ? 1.0 : 0.0
    ((vec(x), vec(y), vec(z), vec(z)), (; accelerate = accelerate, kwargs...))
end

# z values are calculated from a function
function _setargs_surface(f, x, y, fz::Function, args...; kwargs...)
    z = fz.(transpose(vec(x)), vec(y))
    _setargs_surface(f, x, y, z, args...; kwargs...)
end

# z values alone, without x,y coordinates
function _setargs_surface(f, z; kwargs...)
    sz = size(z)
    _setargs_surface(f, 1:sz[2], 1:sz[1], z; kwargs...)
end

function _setargs_surface(f, z, h; kwargs...)
    sz = size(z)
    _setargs_surface(f, 1:sz[2], 1:sz[1], z, h; kwargs...)
end

@plotfunction(surface, geom = :surface, axes = :axes3d, setargs = _setargs_surface,
kwargs = (ratio=1.0, colorbar=true, accelerate=true), docstring="""
    surface([x, y,] z; kwargs...)

Draw a three-dimensional surface plot.

Either a series of points or a two-dimensional array is drawn as
a surface plot, colored according to the Z coordinates and the
current colormap. It can receive one of the following:

- `x` values, `y` values and `z` values.
- *M* sorted values of the `x` axis, *N* sorted values of the `y` axis,
    and a set of `z` values on a *N*×*M* grid.
- *M* sorted values of the `x` axis, *N* sorted values of the `y` axis,
    and a callable to determine `z` values.
- `z` values on a *N*×*M* grid, with the *x* and *y* axes
    defined as the indices of the columns and rows of the grid, respectively.

If a series of points is passed to this function, their values will be
interpolated on a grid. For grid points outside the convex hull of the
provided points, a value of 0 will be used.

# Examples

```julia
$(_example("surface"))
```
""")

@plotfunction(wireframe, geom = :wireframe, axes = :axes3d, setargs = _setargs_surface,
kwargs = (ratio=1.0,), docstring="""
    wireframe([x, y,] z; kwargs...)

Draw a three-dimensional wireframe plot.

Either a series of points or a two-dimensional array is drawn as
a wireframe plot. It can receive one of the following:

- `x` values, `y` values and `z` values.
- *M* sorted values of the `x` axis, *N* sorted values of the `y` axis,
    and a set of `z` values on a *N*×*M* grid.
- *M* sorted values of the `x` axis, *N* sorted values of the `y` axis,
    and a callable to determine `z` values.
- `z` values on a *N*×*M* grid, with the *x* and *y* axes
    defined as the indices of the columns and rows of the grid, respectively.

Also use the attributes `color` and `linecolor` to set the color of the
surface and lines of the mesh, as RGB hexadecimal color values.

If a series of points is passed to this function, their values will be
interpolated on a grid. For grid points outside the convex hull of the
provided points, a value of 0 will be used.

# Examples

```julia
$(_example("wireframe"))
```
""")

@plotfunction(trisurf, geom = :trisurf, axes = :axes3d, setargs = _setargs_tricont,
kwargs = (colorbar=true, ratio=1.0), docstring="""
    tricont(x, y, z; kwargs...)

Draw a triangular surface plot.

Either a series of points or a two-dimensional array is drawn as a
triangular surface plot. `z` values are interpolated between `x` and `y`
values through [Delaunay triangulation](http://mathworld.wolfram.com/DelaunayTriangulation.html).

!!! note

    If the series of points is concave, there may be interpolation artifacts on
    the edges of the plot, as the interpolation may occur in very acute triangles.

# Examples

```julia
$(_example("tricont"))
```
""")

function _setargs_heatmap(f, args...; kwargs...)
    if length(args) == 1
        data = args[1]
        h, w = float.(size(data))
        x = [w]
        y = [h]
        xl = (0.0, w)
        yl = (0.0, h)
    else
        x, y, data = args
        @assert size(data) == (length(y)-1, length(x)-1) "inconsistent dimensions of input data"
        xl = extrema(x)
        yl = extrema(y)
    end
    if get(kwargs, :xflip, false)
        data = reverse(data, dims=1)
    end
    if get(kwargs, :yflip, false)
        data = reverse(data, dims=2)
    end
    kwargs = (; xlim = xl, ylim = yl, kwargs...)
    ((x, y, emptyvector(Float64), vec(data')), kwargs)
end

@plotfunction(heatmap, geom = :heatmap, axes = :axes2d, setargs = _setargs_heatmap,
kwargs = (colorbar=true, tickdir=-1), docstring="""
    heatmap([x, y,] data; kwargs...)

Draw a heatmap.

The current colormap is used to display a two-dimensional array `data` as a heatmap.

If `data` is an *N*×*M* array, the cells of the heatmap will be plotted in
a grid of rectangular cells, whose edges are defined by the coordinates
`x` (a sorted vector with `M+1` values) and `y` (a sorted vector with `N+1` values).

If `x` and `y` are not given, a uniform grid is drawn spanning the interval
`[1, M+1]` in the X-axis, and `[1, N+1]` in the Y-axis.

The array is drawn with its first value in the bottom left corner, so in some
cases it may be neccessary to flip the axes.

By default column and row indices are used for the x- and
y-axes, respectively, so setting the axis limits is recommended. Also note that the
values in the array must lie within the current z-axis limits so it may
be neccessary to adjust these limits or clip the range of array values.

# Examples

```julia
$(_example("heatmap"))
```
""")

@plotfunction(polarheatmap, geom = :polarheatmap, axes = :polar, setargs = _setargs_heatmap,
kwargs = (colorbar=true, overlay_axes=true, ratio=1.0), docstring="""
    polarheatmap(data; kwargs...)

Draw a polar heatmap

The current colormap is used to display a two-dimensional array `data` as a polar heatmap.

If `data` is an *N*×*M* array, the cells of the matrix will be plotted in
a circle divided in *M* angular sectors and *N* rings, such that
the values of the first row will be concentrated in the center of the circle,
and the following rows will be drawn in increasing concentric rings.
The columns of the array are drawn as radii of the circle in a
counterclockwise order, starting and ending in horizontal axis pointing to the right.

# Examples

```julia
$(_example("polarheatmap"))
```
""")

_setargs_hexbin(f, x, y; kwargs...) = ((x, y, emptyvector(Float64), [0.0, 1.0]), kwargs)

@plotfunction(hexbin, geom = :hexbin, axes = :axes2d, setargs = _setargs_hexbin,
kwargs = (colorbar=true, nbins=40), docstring="""
    hexbin(x, y; kwargs...)

Draw a hexagon binning plot.

Hexagonal binning and the the current colormap are used to display a bi-dimensional
series of points given by `x` and `y`.
The number of bins is 40 by default; use the keyword argument `nbins` to set
it as a different number.

# Examples

```julia
$(_example("hexbin"))
```
""")

function _setargs_imshow(f, data, limits=(0,1); kwargs...)
    if isa(data, AbstractString)
        w, h, rgbdata = GR.readimage(data)
    else
        sz = size(data)
        if length(sz) == 2
            h, w = sz
            data_t = data'
            minval, maxval = set_limits(limits, extrema(data))
            valrange = maxval - minval
            GR.setcolormap(get(kwargs, :colormap, COLOR_INDICES[:colormap]))
            rgbdata = zeros(UInt32, w, h)
            for j=1:h, i=1:w
                value = min(max(data_t[i,j], minval), maxval)
                rgbdata[i,j] = to_rgba((value - minval)/valrange)
            end
        elseif length(sz) == 3
            h, w, c = sz
            rgbdata = switchbytes.(color.(data[:,:,1], data[:,:,2], data[:,:,3]))
            if c > 3
                rgbdata .+= round.(UInt32, data[:,:,4] * 255) .<< 24
            else
                rgbdata .+= 0xff000000
            end
            rgbdata = transpose(rgbdata)
        end
    end
    if get(kwargs, :xflip, false)
        rgbdata = reverse(rgbdata, dims=2)
    end
    if get(kwargs, :yflip, false)
        rgbdata = reverse(rgbdata, dims=1)
    end
    kwargs = (; xlim = (0.0, float(w)), ylim = (0.0, float(h)), ratio = w/h, kwargs...)
    ((1.0:w, 1.0:h, emptyvector(Float64), float.(rgbdata[:])), kwargs)
end

@plotfunction(imshow, geom = :image, axes = :axes2d, setargs = _setargs_imshow,
kwargs = (xticks=NULLPAIR, yticks=NULLPAIR, noframe=true), docstring="""
    imshow(img; kwargs...)
    imshow(img, (minval, maxval); kwargs...)

Draw an image, defined by any of the following:
* A string with a valid file name of an image.
* A 3-dimensional *H*×*W*×*C* array, which will be drawn as an
  image *W* pixels wide and *H* pixels high, with *C* channels
  (3 or 4), corresponding to RGB or RGBA values between 0 and 1.
* A matrix of values, which will be drawn with a hue corresponding to
  their relative position in the current colormap. The range of values
  mapped in the colormap is by default `[0, 1]`, but can be set to
  an arbitrary range between `minval` and `maxval`. Setting either limit
  to `nothing` will cause it to be automatically determined based on the data.

# Examples

```julia
$(_example("imshow"))
```
""")

function _setargs_isosurf(f, v, isovalue; skincolor=nothing, kwargs...)
    if skincolor ≠ nothing # deprecate?
        kwargs = (; color=skincolor, kwargs...)
    end
    values = round.((v .- _min(v)) ./ (_max(v) .- _min(v)) .* (2^16-1))
    dimensions = float.(collect(size(v)))
    isoval_norm = (isovalue - _min(v)) / (_max(v) - _min(v))
    # x = dimensions, y = isovalue, z = values, c = color
    ((dimensions, [isoval_norm], values[:]), kwargs)
end

@plotfunction(isosurface, geom = :isosurf, axes = :axes3d, setargs = _setargs_isosurf,
kwargs = (xticks=NULLPAIR, yticks=NULLPAIR, zticks=NULLPAIR, ratio=1.0, gr3=true), docstring="""
    isosurface(data, isovalue; kwargs...)

Draw an isosurface determined by the region of the three-dimensional array `data`
around a given `isovalue`.

The isosurface is calculated so that values in `data` greater than `isovalue` are
considered to be outside the surface, and the values lower than `isovalue` are
inside the surface.

The color of the isosurface can be chosen with the keyword argument
`color`, with the hexadecimal RGB color code.

# Examples

```julia
$(_example("isosurface"))
```
""")

const XFORMS = Dict(
    "boolean"=>0, "linear"=>1, "log"=>2,
    "loglog"=>3, "cubic"=>4, "equalized"=>5
)

function _setargs_shade(f, x, y; kwargs...)
    # Determine type of footprint
    default_footprint = hasnan(x) || hasnan(y) ? "lines" : "points"
    footprint = get(kwargs, :footprint, default_footprint)
    if footprint == "lines"
        kwargs = (; shadelines = 1.0, kwargs...)
    elseif footprint == "points"
        kwargs = (; shadelines = 0.0, kwargs...)
    else
        throw(ArgumentError("""`footprint` must be either `"lines"` or `"points"`"""))
    end
    # Transformation
    if haskey(kwargs, :xform)
        xf = lookup(kwargs[:xform], XFORMS)
        kwargs = (; kwargs..., xform=float(xf))
    end
    return ((x, y), kwargs)
end

@plotfunction(shade, geom = :shade, axes = :axes2d, kwargs = (tickdir=-1,),
setargs = _setargs_shade, docstring="""
    shade(x, y; kwargs...)

Draw a point- or line-based heatmap.

The current colormap is used to display the footprint left by the pairs of
`x`, `y` values. If the data contain `NaN` or `missing`, the footprints
will be based on lines separated by those values. Otherwise the footprints
will be based on points. The type of footprint can be enforced regardless
of the input by the keyword argument `footprint = "lines"` or
`footprint = "points"`.

The value of that footprint is determined by a transformation that can be
adjusted by the keyword argument `xform` --- a number or a string from the
following table:

| # |String       |description                  |
|:-:|:------------|:----------------------------|
| 0 |`"boolean"`  |boolean                      |
| 1 |`"linear"`   |linear                       |
| 2 |`"log"`      |logarithmic                  |
| 3 |`"loglog"`   |double logarithmic           |
| 4 |`"cubic"`    |cubic                        |
| 5 |`"equalized"`|histogram equalized (default)|

# Examples

```julia
$(_example("shade"))
```
""")

const ALGORITHMS = Dict("emission"=>0, "absorption"=>1, "mip"=>2)

function _setargs_volume(f, v::Array{T, 3}; kwargs...) where {T}
    (nx, ny, nz) = size(v)
    # Algorithm
    if haskey(kwargs, :algorithm)
        alg = lookup(kwargs[:algorithm], ALGORITHMS)
        kwargs = (; kwargs..., algorithm=float(alg))
    end
    (([nx], [ny], [nz], vec(v)), kwargs)
end

@plotfunction(volume, geom = :volume, axes = :axes3d, setargs = _setargs_volume,
kwargs = (colorbar=true, ratio=1.0), docstring="""
    volume(v; kwargs...)

Draw a the three-dimensional array `v`, using volume rendering.

The volume data is reduced to a two-dimensional image using
an emission or absorption model, or by a maximum intensity projection.
After the projection the current colormap is applied to the result.

The method to reduce volume data can be defined by the keyword argument
`algorithm` --- a number or a string from the
following table:

| # |String        |description                 |
|:-:|:-------------|:---------------------------|
| 0 |`"emission"`  |emission model (default)    |
| 1 |`"absorption"`|absorption model            |
| 2 |`"mip"`       |maximum intensity projection|

# Examples

```julia
$(_example("volume"))
```
""")

function oplot!(f::Figure, args...; kwargs...)
    p = currentplot(f)
    kwargs = (; p.attributes..., kwargs...)
    if typeof(args[end]) <: AbstractString
        kwargs = (spec=args[end], kwargs...)
        args = args[1:end-1]
    end
    append!(p.geoms, geometries(:line, args...; geom_attributes(; kwargs...)...))
    return f
end

"""
    oplot(args...; kwargs...)

Draw one or more line plots over another plot.

Equivalent to calling [`plot`](@ref) after holding the current plot,
except that the axes limits are not re-adjusted to the new data.

# Examples

```julia
$(_example("oplot"))
```
"""
oplot(args...; kwargs...) = oplot!(gcf(), args...; kwargs...)

function _setargs_annotation(f, x, y, s::AbstractString; kwargs...)
    # only works in axes2d
    p = currentplot(f)
    p.axes.kind ≠ :axes2d && throw(ErrorException("annotations are only available for 2D plots"))
    # Get coordinates of text box in world coordinates
    GR.savestate()
    GR.setwindow(p.axes.ranges[:x]..., p.axes.ranges[:y]...)
    GR.setviewport(p.viewport.inner...)
    width, height = stringsize(s, charheight(p.viewport.inner), true)
    GR.restorestate()
    halign = get(kwargs, :halign, "left")
    valign = get(kwargs, :valign, "bottom")
    if halign == "center"
        x -= 0.5 * width
    elseif halign == "right"
        x -= width
    end
    if valign == "center"
        y -= 0.5 * height
    elseif valign == "top"
        y -= height
    end
    (([x, x+width], [y, y+height]), (; kwargs..., label=s))
end

@plotfunction(_annotation, kind=:text, geom=:text, axes=:axes2d, setargs=_setargs_annotation)

function annotations!(f::Figure, x::Real, y::Real, s::AbstractString; kwargs...)
    holdstate = get(currentplot(f).attributes, :hold, false)
    hold!(currentplot(f), true)
    _annotation!(f, x, y, s; kwargs...)
    hold!(currentplot(f), holdstate)
    draw(f)
end

function annotations!(f::Figure, x::AbstractArray, y::AbstractArray, s::AbstractArray{<:AbstractString}; kwargs...)
    p = currentplot(f)
    kwargs_ext = (; p.attributes..., kwargs...)
    for i = 1:length(s)-1
        args, kwargs_ext = _setargs_annotation(f, x[i], y[i], s[i]; kwargs_ext...)
        append!(p.geoms, GRUtils.geometries(:text, args...; GRUtils.geom_attributes(;kwargs_ext...)...))
    end
    annotations!(f, x[end], y[end], s[end]; kwargs...)
end

"""
    annotations(x, y, s)

Add one ore more text annotations to the current plot.

`x` and `y` can be scalars or vectors of coordinates, and `s` a string
or a vector of strings, respectively. Annotations are only available
for 2-D plots.

By default the coordinates indicate the lower left corner of the text box.
This can be changed by the following keyword arguments:

* `halign` for the horizontal alignment (`"left"`, `"center"` or `"right"`).
* `valign` for the vertical alignment (`"bottom"`, `"center"` or `"top"`).

# Examples
```julia
$(_example("annotations"))
```
"""
annotations(x, y, s; kwargs...) = annotations!(gcf(), x, y, s; kwargs...)

"""
    savefig(filename[, fig])

Save a figure to a file.

Draw the current figure in a file of the given name.
Which file types are supported depends on the installed workstation types,
but GR usually is built with support for
.png, .jpg, .pdf, .ps, .gif and various other file formats.

If no figure is given (optional argument `fig`), the current figure
is used.

# Examples

```julia
# Create a simple plot
x = 1:100
plot(x, 1 ./ (x .+ 1))
# Save the figure to a file
savefig("example.png")
```
"""
function savefig(filename, fig=gcf())
    GR.beginprint(filename)
    draw(fig)
    GR.endprint()
end


text:
# This file contains functions for the Revenue Efficiency Additive DEA model
"""
    RevenueDEAModel
An data structure representing a revenue additive DEA model.
"""
struct RevenueAdditiveDEAModel <: AbstractRevenueDEAModel
    n::Int64
    m::Int64
    s::Int64
    weights::Symbol
    monetary::Bool
    rts::Symbol
    disposX::Symbol
    dmunames::Union{Vector{String},Nothing}
    eff::Vector
    lambda::SparseMatrixCSC{Float64, Int64}
    techeff::Vector
    alloceff::Vector
    normalization::Vector
    Xtarget::Matrix
    Ytarget::Matrix
end

"""
    dearevenueadd(X, Y, P, model)
Compute revenue efficiency using additive data envelopment analysis for
inputs `X`, outputs `Y` and price of outputs `P`.

Model specification:
- `:Ones`: standard additive DEA model.
- `:MIP`: Measure of Inefficiency Proportions. (Charnes et al., 1987; Cooper et al., 1999)
- `:Normalized`: Normalized weighted additive DEA model. (Lovell and Pastor, 1995)
- `:RAM`: Range Adjusted Measure. (Cooper et al., 1999)
- `:BAM`: Bounded Adjusted Measure. (Cooper et al, 2011)
- `:Custom`: User supplied weights.

# Optional Arguments
- `rts=:VRS`: chooses variable returns to scale. For constant returns to scale choose `:CRS`.
- `rhoY`: matrix of weights of outputs. Only if `model=:Custom`.
- `disposal=:Strong`: chooses strong disposal of inputs. For weak disposal choose `:Weak`.
- `names`: a vector of strings with the names of the decision making units.

# Examples
```jldoctest
julia> X = [1; 1; 1; 1; 1; 1; 1; 1];

julia> Y = [7 7; 4 8; 8 4; 3 5; 3 3; 8 2; 6 4; 1.5 5];

julia> P = [1 1; 1 1; 1 1; 1 1; 1 1; 1 1; 1 1; 1 1];

julia> dearevenueadd(X, Y, P, :Ones)
Revenue Additive DEA Model
DMUs = 8; Inputs = 1; Outputs = 2
Orientation = Output; Returns to Scale = VRS
Weights = Ones
───────────────────────────────────
   Revenue  Technical    Allocative
───────────────────────────────────
1      0.0        0.0   0.0
2      2.0        0.0   2.0
3      2.0        0.0   2.0
4      6.0        6.0   0.0
5      8.0        8.0   0.0
6      4.0        2.0   2.0
7      4.0        4.0  -8.88178e-16
8      7.5        7.5   0.0
───────────────────────────────────

```
"""
function dearevenueadd(X::Union{Matrix,Vector}, Y::Union{Matrix,Vector},
    P::Union{Matrix,Vector}, model::Symbol = :Default;
    rts::Symbol = :VRS,
    rhoY::Union{Matrix,Vector,Nothing} = nothing,
    dispos::Symbol = :Strong, monetary::Bool = false,
    names::Union{Vector{String},Nothing} = nothing,
    optimizer::Union{DEAOptimizer,Nothing} = nothing)::RevenueAdditiveDEAModel

    # Check parameters
    nx, m = size(X, 1), size(X, 2)
    ny, s = size(Y, 1), size(Y, 2)

    np, sp = size(P, 1), size(P, 2)

    if nx != ny
        throw(DimensionMismatch("number of rows in X and Y ($nx, $ny) are not equal"));
    end
    if np != ny
        throw(DimensionMismatch("number of rows in P and Y ($np, $ny) are not equal"));
    end
    if sp != s
        throw(DimensionMismatch("number of columns in P and Y ($sp, $s) are not equal"));
    end

    if dispos != :Strong && dispos != :Weak
        throw(ArgumentError("`dispos` must be :Strong or :Weak"));
    end

    # Default behaviour
    if model == :Default
        # If no weights are specified use :Ones
        if rhoY === nothing
            model = :Ones
        else
            model = :Custom
        end
    end

    # Get weights based on the selected model
    if model != :Custom
        # Display error if both model and weights are specified
        if rhoY !== nothing
            throw(ArgumentError("`rhoY` not allowed if model != :Custom"));
        end

        # Get weights for selected model
        rhoX, rhoY = deaaddweights(X, Y, model, orient = :Output)
    end

    if size(rhoY) != size(Y)
        throw(DimensionMismatch("size of rhoY and Y ($(size(rhoY)), $(size(Y))) are not equal"));
    end

    # Default optimizer
    if optimizer === nothing
        optimizer = DEAOptimizer(:LP)
    end

    # Get maximum revenue targets and lambdas
    n = nx

    Xtarget = X[:,:]
    Ytarget, rlambdaeff = deamaxrevenue(X, Y, P, rts = rts, dispos = dispos, optimizer = optimizer)

    # Revenue, technical and allocative efficiency
    refficiency  = vec(sum(P .* Ytarget, dims = 2)  .- sum(P .* Y, dims = 2))
    normalization = vec(minimum(P ./ rhoY, dims = 2) )

    if model == :Custom
        rhoX = ones(size(X))
        techefficiency = efficiency(deaadd(X, Y, model, orient = :Output, rts = rts, rhoX = rhoX, rhoY = rhoY, disposX = dispos, optimizer = optimizer))
    else
        techefficiency = efficiency(deaadd(X, Y, model, orient = :Output, rts = rts, disposX = dispos, optimizer = optimizer))
    end

    if monetary
        techefficiency = techefficiency .* normalization
    else
        refficiency = refficiency ./ normalization
    end

    allocefficiency = refficiency - techefficiency

    return RevenueAdditiveDEAModel(n, m, s, model, monetary, rts, dispos, names, refficiency, rlambdaeff, techefficiency, allocefficiency, normalization, Xtarget, Ytarget)

end

function Base.show(io::IO, x::RevenueAdditiveDEAModel)
    compact = get(io, :compact, false)

    n = nobs(x)
    m = ninputs(x)
    s = noutputs(x)
    disposX = x.disposX
    dmunames = names(x)

    eff = efficiency(x)
    techeff = efficiency(x, :Technical)
    alloceff = efficiency(x, :Allocative)

    if !compact
        print(io, "Revenue Additive DEA Model \n")
        print(io, "DMUs = ", n)
        print(io, "; Inputs = ", m)
        print(io, "; Outputs = ", s)
        print(io, "\n")
        print(io, "Orientation = Output")
        print(io, "; Returns to Scale = ", string(x.rts))
        print(io, "\n")
        print(io, "Weights = ", string(x.weights))
        print(io, "\n")

        if disposX == :Weak println(io, "Weak disposability of inputs") end

        show(io, CoefTable(hcat(eff, techeff, alloceff), ["Revenue", "Technical", "Allocative"], dmunames))
    end

end