text:
"Units in plots, for MechanicalUnits.jl>Unitful.jl"
module MechGluePlots
import Unitfu
using Unitfu: Quantity, unit, ustrip, ∙, numtype
import Plots
using Plots: default
using RecipesBase
using RecipesBase: _debug_recipes
struct SeriesUnitInfo
    serno
    unit_foo
    label
    relevant_axis_guide
    letter
    sertyp
end

include("internal_functions.jl")

#########
# Arguments that are spit into series to handle units by series
# In a few combinations of arguments, we need to split before evaluating
# in order to capture the unit information. E.g. in a case like
# plot( xlims =(-5s, 5s), [f_q_q, x-> 0.5f_q_q(1.5x)], ribbon = [x-> f_q_q(0.25x), x-> f_q_q(0.5x)])
#########

@recipe function f(fs::AbstractArray{F}, xmin::T, xmax::T) where {F <: Function, T<:Quantity}
    _debug_recipes[1] && printstyled(color=:172, "\n    MechGluePlots fs::AbstractArray{F}, xmin::T, xmax::T) where {F <: Function, T<:Quantity}\n")
    ribbon = get(plotattributes, :ribbon, nothing)
    for (serno, f) in enumerate(fs)
        @series begin
            if !isnothing(ribbon)
                r = extract_series(ribbon, serno)
                ribbon := r
            end
            series_plotindex := serno
            f, xmin, xmax
        end
    end
    nothing
end

@recipe function f(fs::VecOrMat{F}, y::T) where {F<:Function, T <: AbstractArray{<:Union{Missing,<:Quantity}}}
    _debug_recipes[1] && printstyled(color=:172, "\n    MechGluePlots fs::AbstractArray{F}, y::T) where {F<:Function, T <: AbstractArray{<:Union{Missing,<:Quantity}}}\n")
    ribbon = get(plotattributes, :ribbon, nothing)
    for (serno, f) in enumerate(fs)
        @series begin
            if !isnothing(ribbon)
                r = extract_series(ribbon, serno)
                ribbon := r
            end
            ys = extract_series(y, serno)
            series_plotindex := serno
            f, ys
        end
    end
end

@recipe function f( x::T, fs::VecOrMat{F}) where {F<:Function, T <: AbstractArray{<:Union{Missing,<:Quantity}}}
    _debug_recipes[1] && printstyled(color=:172, "\n    MechGluePlots  where {F<:Function, T <: AbstractArray{<:Union{Missing,<:Quantity}}}\n")
    ribbon = get(plotattributes, :ribbon, nothing)
    for (serno, f) in enumerate(fs)
        @series begin
            if !isnothing(ribbon)
                r = extract_series(ribbon, serno)
                ribbon := r
            end
            xs = extract_series(x, serno)
            series_plotindex := serno
            xs, f
        end
    end
end

#########
# Arguments that are evalated before quantities are dropped
# from axes. For example,
# plot(xlims = (0, 10)s, sin(t /s))
# needs to be evaluated before 'xlims' is converted to a unitless tuple.
#########

@recipe function f(foo::F, x::T)  where {F<:Function, T<:AbstractArray{<:Union{Missing,<:Quantity}}}
    _debug_recipes[1] && printstyled(color=:yellow, "\n    MechGluePlots foo::F, x::T)  where {F<:Function, T<:AbstractArray{<:Union{Missing,<:Quantity}}}\n")
    x, map(foo, x) # Note the swap. For series with values, the ordinat is placed second.
end

@recipe function f(x::T, foo::F) where {F<:Function, T<:AbstractArray{<:Union{Missing,<:Quantity}}}
    _debug_recipes[1] && printstyled(color=:yellow, "\n    MechGluePlots x::T, foo::F) where {F<:Function, T<:AbstractArray{<:Union{Missing,<:Quantity}}}\n")
    # Without recipes, plot(0:10, sin) and plot(sin, 0:10) produce the same plot, sin on y. Same here:
    x, map(foo, x)
end

#########
# Arguments that include quantities.
# In the pipeline, functions are evaluated using units.
# If the evaluated functions include units, the result wil be stripped of
# units here. The information is not lost, but stored in 'plotattributes[:unitinfo].
# Data that pass through here are temporarily assigned type :postfixseries.
# For such series, the stored unit info is added to axis guides after all series
# have been parsed.
#########

@recipe function f(::Type{T}, x::T) where T <: AbstractArray{<:Union{Missing,<:Quantity}}
    # while plotting, change a vector of quantities to a unitless vector,
    # but pass the lost info on to the later in the pipeline.
    _debug_recipes[1] && printstyled(color=:red, "\n    MechGluePlots ::Type{T}, x::T) where T <: AbstractArray{<:Union{Missing,<:Quantity}} \n")
    _debug_recipes[1] &&  println("    size(x) = $(size(x))")
    _debug_recipes[1] &&  print_prettyln(plotattributes)
    _debug_recipes[1] &&  println(" ")
    serno = get(plotattributes, :series_plotindex, 1)
    letter = plotattributes[:letter]
    # Prior, if any, unit and series type info
    vui = get(plotattributes, :unitinfo, SeriesUnitInfo[])
    # How can we tell if there are more series or not?
    # this is how: Multiple column x, but one row vector: The folloing columns are ignored.
    # But with one column vector, everything is supposed to belong to this series.
    # So passing one column in, everything belong to the first series.

    ribbon = get(plotattributes, :ribbon, nothing)
    # Prepare unitless output values
    vals = similar(x, Float64)
    local modrib, u
    # When in a type recipe, we can't actually split in two series,
    # but we can take note of units, strip them, and passively let the split be done by other
    # parts of the pipeline.
    # For the 'ribbon', this is a bit more complicated but essentially we do the same.
    for colno in 1:size(x, 2)
        thisserno = serno + colno - 1

        sx = x[:, colno]
        u = unit(first(sx))
        _debug_recipes[1] &&  println("        serno = $serno, u = $u")
        vui  = accumulate_unit_info(vui, plotattributes, u, thisserno)
        # Numeric type
        nut = numtype(first(sx))
        sux = map(nut, ustrip(x[:, colno] ./ u))
        vals[:, colno] = sux
        if !isnothing(ribbon)
            modrib = modified_ribbon(colno == 1 ? ribbon : modrib, colno, letter, u)
            if _debug_recipes[1]
                    println("    Original ribbon : $ribbon")
                    println("    Modified ribbon : $modrib")
                    if letter == :y && modrib isa Function
                        println("    Modified ribbon with unitless input: $modrib(0.5) = $(modrib(0.5))")
                    end
            end
        end
    end
    # Pass the info stripped from 'x' to dictionary 'plotattributes' in the pipeline
    unitinfo := vui
    if !isnothing(ribbon)
        ribbon := modrib
    end
    # Strip units off limits.
    li = get(plotattributes, Symbol(letter, :lims), nothing)
    if !isnothing(li)
        ulli = ustrip.(li)
        if letter == :x
            xlims := ustrip.(li)
        elseif letter == :y
            ylims := ustrip.(li)
        elseif letter == :z
            zlims := ustrip.(li)
        end
    end
    # For dispatch to series recipe
    seriestype := :postfixseries
    vals
end

#########
# Functions with arguments that include quantities.
# We do the evaluation here, since the order in which
# limits are determined by Plots.jl varies depending on
# the order of arguments.
#########

#f(y), at user_recipe.jl:36
@recipe function f(foo::F, arg1::T, arg2::T) where {F <: Function, T <: AbstractArray{<:Union{Missing, <:Quantity}}}
    _debug_recipes[1] && printstyled(color=:172, "\n    foo::F, arg1::T, arg2::T) where {F <: Function, T <: AbstractArray{<:Union{Missing,<:Quantity}}}\n")
    map(foo, arg1), arg1, arg2
end

#f(x)
@recipe function f(arg1::T, foo::F, arg2::T) where {F <: Function, T <: AbstractArray{<:Union{Missing,<:Real}}}
    _debug_recipes[1] && printstyled(color=:172, "\n    arg1::T, foo::F, arg2::T) where {F <: Function, T <: AbstractArray{<:Union{Missing,<:Real}}}\n")
    # Prior, if any, unit and series type info
    vui = get(plotattributes, :unitinfo, SeriesUnitInfo[])
    isnothing(vui) && error("f(x): Wrong  MechGluePlots recipe triggered. Missing unit info.")
    # Evaluate foo with correct units, although they have been stripped from arg1 and arg2.
    u1, serno = axis_unit_series_no(:x, vui)
    # We'll pass on inlcuding units. If necessary, the pipeline will remove units and store the information.
    arg1, map(foo, arg1∙u1), arg2
end

#########
# Lastly, find the info about stripped units again.
# This is triggered once per series, but we can now access
# a storage that is common to all series (we couldn't before!).
# When the current series is the last series, the accumulated
# info is placed on the axis guides.
#########

@shorthands postfixseries
@recipe function f(::Type{Val{:postfixseries}}, plt::AbstractPlot)
    _debug_recipes[1] && printstyled(color=:blue, "\n    MechGluePlots ::Type{Val{:postfixseries}}, plt)\n")
    # Info stored earlier in the pipeline for this series. There may be more. In that case,
    # this will be called several times for the same plot.
    unitinfo = plotattributes[:unitinfo]
    serno = plotattributes[:series_plotindex]

    _debug_recipes[1] &&  print_prettyln(plotattributes)

    # Pick the interesting attributes dictionary from which to read existing guides
    read_from_plot_object = length(plt.attr.explicit) > 0
    read_attr = read_from_plot_object ? plt.attr.explicit : plotattributes
    if read_from_plot_object
        _debug_recipes[1] &&  println("    Previously stored attributes in plot object:")
        _debug_recipes[1] &&  print_prettyln(read_attr)
    end

    # Preserve existing guides and units, add bracketed units from this series
    gux = axis_guide_with_units(:x, read_attr, unitinfo)
    guy = axis_guide_with_units(:y, read_attr, unitinfo)
    guz = axis_guide_with_units(:z, read_attr, unitinfo)

    if serno == plt.n
        # Last series in plot, write to (local) plotattributes
        # and delete the temporary storage we used in 'plt'
        a = plt.attr.explicit
        if gux != ""
            xguide := gux
            haskey(a, :xguide) && pop!(a, :xguide)
        end
        if guy != ""
            yguide := guy
            haskey(a, :yguide) && pop!(a, :yguide)
        end
        if guz != ""
            zguide := guz
            haskey(a, :zguide) && pop!(a, :zguide)
        end
        _debug_recipes[1] && print_prettyln(plt.attr.explicit)
    else
        # Preserve
        gux != "" && push!(plt.attr.explicit, (:xguide => gux))
        guy != "" && push!(plt.attr.explicit, (:yguide => guy))
        guz != "" && push!(plt.attr.explicit, (:zguide => guz))
    end
    # Signal to the pipeline that we're through with this series
    seriestyp = sertype(serno, unitinfo)
    seriestype := seriestyp
    plt
end

end


text:
# This file contains functions for the Profit Efficiency Enhanced Russell Graph Slack Based Measure DEA model
"""
    ProfitERGDEAModel
An data structure representing a profit Enhanced Russell Graph Slack Based Measure DEA model.
"""
struct ProfitERGDEAModel <: AbstractProfitDEAModel
    n::Int64
    m::Int64
    s::Int64
    monetary::Bool
    dmunames::Union{Vector{String},Nothing}
    eff::Vector
    lambda::SparseMatrixCSC{Float64, Int64}
    techeff::Vector
    alloceff::Vector
    normalization::Vector
    Xtarget::Matrix
    Ytarget::Matrix
end

"""
    deaprofiterg(X, Y, W, P)
Compute profit efficiency using data envelopment analysis Enhanced Russell Graph Slack Based Measure
model for inputs `X`, outputs `Y`, price of inputs `W`, and price of outputs `P`.

# Optional Arguments
- `names`: a vector of strings with the names of the decision making units.

# Examples
```jldoctest
julia> X = [2; 4; 8; 12; 6; 14; 14; 9.412];

julia> Y = [1; 5; 8; 9; 3; 7; 9; 2.353];

julia> W = [1; 1; 1; 1; 1; 1; 1; 1];

julia> P = [2; 2; 2; 2; 2; 2; 2; 2];

julia> deaprofiterg(X, Y, W, P)
Enhanced Russell Graph Slack Based Measure Profit DEA Model
DMUs = 8; Inputs = 1; Outputs = 1
Returns to Scale = VRS
──────────────────────────────────
     Profit  Technical  Allocative
──────────────────────────────────
1  4.0        0.0         4.0
2  0.5        0.0         0.5
3  0.0        0.0         0.0
4  0.166667   0.0         0.166667
5  0.8        0.6         0.2
6  0.571429   0.52381     0.047619
7  0.285714   0.142857    0.142857
8  1.2706     0.8         0.4706
──────────────────────────────────
```
"""
function deaprofiterg(X::Union{Matrix,Vector}, Y::Union{Matrix,Vector},
    W::Union{Matrix,Vector}, P::Union{Matrix,Vector};
    monetary::Bool = false,
    names::Union{Vector{String},Nothing} = nothing,
    optimizer::Union{DEAOptimizer,Nothing} = nothing)::ProfitERGDEAModel

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

    nw, mw = size(W, 1), size(W, 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 nw != nx
        throw(DimensionMismatch("number of rows in W and X ($nw, $nx) are not equal"));
    end
    if np != ny
        throw(DimensionMismatch("number of rows in P and Y ($np, $ny) are not equal"));
    end
    if mw != m
        throw(DimensionMismatch("number of columns in W and X ($mw, $m) are not equal"));
    end
    if sp != s
        throw(DimensionMismatch("number of columns in P and Y ($sp, $s) are not equal"));
    end

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

    # Get maximum profit targets and lambdas
    n = nx

    Xtarget, Ytarget, plambdaeff = deamaxprofit(X, Y, W, P, optimizer = optimizer)

    # Profit, technical and allocative efficiency
    maxprofit = sum(P .* Ytarget, dims = 2) .- sum(W .* Xtarget, dims = 2)
    pefficiency  = vec(maxprofit .- ( sum(P .* Y, dims = 2) .- sum(W .* X, dims = 2)))
    techergmodel = deaerg(X, Y, rts = :VRS, optimizer = optimizer)
    normalization = vec( minimum(hcat(m .* W .* X, s .* P .* Y), dims = 2) .* (1 .+ 1 ./ s .* sum(slacks(techergmodel, :Y) ./ Y, dims = 2)) )
    techefficiency = 1 .- efficiency(techergmodel)

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

    allocefficiency = pefficiency - techefficiency

    return ProfitERGDEAModel(n, m, s, monetary, names, pefficiency, plambdaeff, techefficiency, allocefficiency, normalization, Xtarget, Ytarget)

end

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

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

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

    if !compact
        print(io, "Enhanced Russell Graph Slack Based Measure Profit DEA Model \n")
        print(io, "DMUs = ", n)
        print(io, "; Inputs = ", m)
        print(io, "; Outputs = ", s)
        print(io, "\n")
        print(io, "Returns to Scale = VRS")
        print(io, "\n")
        show(io, CoefTable(hcat(eff, techeff, alloceff), ["Profit", "Technical", "Allocative"], dmunames))
    end
end

ext: jl


text:
export get_all_rutes_from_home;
export solve_with_very_brute_force;

"""
Get all the possibles ways to visit all the cities from home.
"""
function get_all_rutes_from_home(ncities; home::Int = 1)::Vector{Vector{Int}}
    return pushfirst!.(all_possible_rutes(copy_and_remove(collect(1:ncities),home)), home);
end

"""
Solve teh TSP using brute force. It finds all the rutes and then pick the cheapest.
Returns a CheapestRute struct with the rute, the total rute cost and the
cost per interval for the computed rute.
"""
function solve_with_very_brute_force(costmatrix::Matrix{T}; home::Int = 1)::CheapestRute where T
    checkmatrix(costmatrix);
    cheapest_rute = pick_cheaper(costmatrix, get_all_rutes_from_home(size(costmatrix,1); home = home));
    return CheapestRute(cheapest_rute, rute_total_cost(costmatrix, cheapest_rute), rute_stretches_costs(costmatrix, cheapest_rute));
end


text:
@doc raw"""
The Gauss nodes are given by the roots of the shifted Legendre polynomial
$P_s (2x-1)$ with $s$ the number of stages.
"""
function get_gauss_nodes(::Type{T}, s) where {T}
    sort(T.(Polynomials.roots(_shifted_legendre(s,T))))
end

get_gauss_nodes(s) = get_gauss_nodes(BigFloat, s)

@doc raw"""
The Gauss weights are given by the following integrals
```math
b_i = \bigg( \frac{dP}{dx} (c_i) \bigg)^{-2} \int \limits_0^1 \bigg( \frac{P(x)}{x - c_i} \bigg)^2 dx ,
```
where $P(x)$ denotes the shifted Legendre polynomial
$P(x) = P_s (2x-1)$ with $s$ the number of stages.
"""
function get_gauss_weights(::Type{T}, s) where {T}
    c = get_gauss_nodes(T,s)
    P = _shifted_legendre(s,T)
    D = Polynomials.derivative(P)

    inti(i) = begin
        I = Polynomials.integrate( ( P ÷ Polynomial(T[-c[i], 1]) )^2 )
        I(1) - I(0)
    end

    b = [ inti(i) / D(c[i])^2  for i in 1:s ]
end

get_gauss_weights(s) = get_gauss_weights(BigFloat, s)

@doc raw"""
The Gauss coefficients are implicitly given by the so-called simplifying assumption $C(s)$:
```math
\sum \limits_{j=1}^{s} a_{ij} c_{j}^{k-1} = \frac{c_i^k}{k}  \qquad i = 1 , \, ... , \, s , \; k = 1 , \, ... , \, s .
```
"""
function get_gauss_coefficients(::Type{T}, s) where {T}
    solve_simplifying_assumption_c(get_gauss_nodes(T,s))
end

get_gauss_coefficients(s) = get_gauss_coefficients(BigFloat, s)

"""
Gauss tableau with s stages

```julia
TableauGauss(::Type{T}, s)
TableauGauss(s) = TableauGauss(Float64, s)
```

The constructor takes the number of stages `s` and optionally the element type `T` of the tableau.

References:

    John C. Butcher.
    Implicit Runge-Kutta processes.
    Mathematics of Computation, Volume 18, Pages 50-64, 1964.
    doi: 10.1090/S0025-5718-1964-0159424-9.

    John C. Butcher.
    Gauss Methods.
    In: Engquist B. (eds). Encyclopedia of Applied and Computational Mathematics. Springer, Berlin, Heidelberg. 2015.
    doi: 10.1007/978-3-540-70529-1_115.

"""
function TableauGauss(::Type{T}, s) where {T}
    Tableau{T}(Symbol("Gauss($s)"), 2s, get_gauss_coefficients(s), get_gauss_weights(s), get_gauss_nodes(s))
end

TableauGauss(s) = TableauGauss(Float64, s)


text:
#################### Posterior Plots ####################

#################### Generic Methods ####################

function draw(p::Array{Plot}; fmt::Symbol=:svg, filename::AbstractString="",
              width::MeasureOrNumber=8inch, height::MeasureOrNumber=8inch,
              nrow::Integer=3, ncol::Integer=2, byrow::Bool=true,
              ask::Bool=true)

  fmt in [:pdf, :pgf, :png, :ps, :svg] ||
    throw(ArgumentError("unsupported draw format $fmt"))

  f(args...) = fmt == :pdf ? PDF(args...) :
               fmt == :pgf ? PGF(args...) :
               fmt == :png ? PNG(args...) :
               fmt == :ps  ? PS(args...)  : SVG(args...)

  isexternalfile = length(filename) > 0
  addextension = isexternalfile && something(findfirst(isequal('.'), filename), 0) == 0
  args = (width, height)

  pp = nrow * ncol               ## plots per page
  ps = length(p)                 ## number of plots
  np = ceil(Int, ps / pp)        ## number of pages

  mat = Array{Context}(undef, pp)
  for page in 1:np
    if ask && page > 1 && !addextension
      println("Press ENTER to draw next plot")
      readline(stdin)
    end

    if isexternalfile
      fname = filename
      if addextension
        fname = string(fname, '-', page, '.', fmt)
      end
      args = (fname, width, height)
    end
    img = f(args...)

    nrem = ps - (page - 1) * pp
    for j in 1:pp
      if j <= nrem
        mat[j] = render(p[(page - 1) * pp + j])
      else
        mat[j] = context()
      end
    end
    result = byrow ? permutedims(reshape(mat, ncol, nrow), [2, 1]) :
                     reshape(mat, nrow, ncol)

    draw(img, gridstack(result))
  end

end

function plot(c::AbstractChains, ptype::Vector{Symbol}=[:trace, :density];
              legend::Bool=false, args...)
  n = length(ptype)
  p = Array{Plot}(undef, n, size(c, 2))
  for i in 1:n
    showlegend = legend && i == n
    p[i, :] = plot(c, ptype[i]; legend=showlegend, args...)
  end
  p
end

function plot(c::AbstractChains, ptype::Symbol; legend::Bool=false, args...)
  ptype == :autocor      ? autocorplot(c; legend=legend, args...) :
  ptype == :bar          ? barplot(c; legend=legend, args...) :
  ptype == :contour      ? contourplot(c; args...) :
  ptype == :density      ? densityplot(c; legend=legend, args...) :
  ptype == :mean         ? meanplot(c; legend=legend, args...) :
  ptype == :mixeddensity ? mixeddensityplot(c; legend=legend, args...) :
  ptype == :trace        ? traceplot(c; legend=legend, args...) :
    throw(ArgumentError("unsupported plot type $ptype"))
end

#################### Plot Engines ####################

function autocorplot(c::AbstractChains;
                     maxlag::Integer=round(Int, 10 * log10(length(c.range))),
                     legend::Bool=false, na...)
  nrows, nvars, nchains = size(c.value)
  plots = Array{Plot}(undef, nvars)
  pos = legend ? :right : :none
  lags = 0:maxlag
  ac = autocor(c, lags=collect(lags))
  for i in 1:nvars
    plots[i] = plot(y=vec(ac.value[i, :, :]),
                    x=repeat(collect(lags * step(c)), outer=[nchains]),
                    Geom.line(),
                    color=repeat(c.chains, inner=[length(lags)]),
                    Scale.color_discrete(), Guide.colorkey(title="Chain"),
                    Guide.xlabel("Lag", orientation=:horizontal),
                    Guide.ylabel("Autocorrelation", orientation=:vertical),
                    Guide.title(c.names[i]), Theme(key_position=pos))
  end
  return plots
end

function barplot(c::AbstractChains; legend::Bool=false,
                 position::Symbol=:stack, na...)
  nrows, nvars, nchains = size(c.value)
  plots = Array{Plot}(undef, nvars)
  pos = legend ? :right : :none
  for i in 1:nvars
    S = unique(c.value[:, i, :])
    n = length(S)
    x = repeat(S, 1, nchains)
    y = zeros(n, nchains)
    for j in 1:nchains
      m = countmap(c.value[:, i, j])
      for k in 1:n
        if S[k] in keys(m)
          y[k, j] = m[S[k]] / nrows
        end
      end
    end
    ymax = maximum(position == :stack ? mapslices(sum, y, dims=2) : y)
    plots[i] = plot(x=vec(x), y=vec(y), Geom.bar(position=position),
                    color=repeat(c.chains, inner=[n]),
                    Scale.color_discrete(), Guide.colorkey(title="Chain"),
                    Guide.xlabel("Value", orientation=:horizontal),
                    Guide.ylabel("Density", orientation=:vertical),
                    Guide.title(c.names[i]), Theme(key_position=pos),
                    Scale.y_continuous(minvalue=0.0, maxvalue=ymax))
  end
  return plots
end

function contourplot(c::AbstractChains; bins::Integer=100, na...)
  nrows, nvars, nchains = size(c.value)
  plots = Plot[]
  offset = 1e4 * eps()
  n = nrows * nchains
  for i in 1:(nvars - 1)
    X = c.value[:, i, :]
    qx = range(minimum(X) - offset, stop=maximum(X) + offset, length=bins + 1)
    mx = map(k -> mean([qx[k], qx[k + 1]]), 1:bins)
    idx = Int[findfirst(k -> qx[k] <= x < qx[k + 1], 1:bins) for x in X]
    for j in (i + 1):nvars
      Y = c.value[:, j, :]
      qy = range(minimum(Y) - offset, stop=maximum(Y) + offset, length=bins + 1)
      my = map(k -> mean([qy[k], qy[k + 1]]), 1:bins)
      idy = Int[findfirst(k -> qy[k] <= y < qy[k + 1], 1:bins) for y in Y]
      density = zeros(bins, bins)
      for k in 1:n
        density[idx[k], idy[k]] += 1.0 / n
      end
      p = plot(x=mx, y=my, z=density, Geom.contour(),
               Guide.colorkey(title="Density"),
               Guide.xlabel(c.names[i], orientation=:horizontal),
               Guide.ylabel(c.names[j], orientation=:vertical))
      push!(plots, p)
    end
  end
  return plots
end

function densityplot(c::AbstractChains; legend::Bool=false,
                     trim::Tuple{Real, Real}=(0.025, 0.975), na...)
  nrows, nvars, nchains = size(c.value)
  plots = Array{Plot}(undef, nvars)
  pos = legend ? :right : :none
  for i in 1:nvars
    val = Array{Vector{Float64}}(undef, nchains)
    for j in 1:nchains
      qs = quantile(c.value[:, i, j], [trim[1], trim[2]])
      val[j] = c.value[qs[1] .<= c.value[:, i, j] .<= qs[2], i, j]
    end
    plots[i] = plot(x=[val...;], Geom.density(),
                    color=repeat(c.chains, inner=[length(c.range)]),
                    Scale.color_discrete(), Guide.colorkey(title="Chain"),
                    Guide.xlabel("Value", orientation=:horizontal),
                    Guide.ylabel("Density", orientation=:vertical),
                    Guide.title(c.names[i]), Theme(key_position=pos))
  end
  return plots
end

function meanplot(c::AbstractChains; legend::Bool=false, na...)
  nrows, nvars, nchains = size(c.value)
  plots = Array{Plot}(undef, nvars)
  pos = legend ? :right : :none
  val = cummean(c.value)
  for i in 1:nvars
    plots[i] = plot(y=vec(val[:, i, :]),
                    x=repeat(collect(c.range), outer=[nchains]),
                    Geom.line(),
                    color=repeat(c.chains, inner=[length(c.range)]),
                    Scale.color_discrete(), Guide.colorkey(title="Chain"),
                    Guide.xlabel("Iteration", orientation=:horizontal),
                    Guide.ylabel("Mean", orientation=:vertical),
                    Guide.title(c.names[i]), Theme(key_position=pos))
  end
  return plots
end

function mixeddensityplot(c::AbstractChains;
                          barbounds::Tuple{Real, Real}=(0, Inf), args...)
  plots = Array{Plot}(undef, size(c, 2))
  discrete = indiscretesupport(c, barbounds)
  plots[discrete] = plot(c[:, discrete, :], :bar; args...)
  plots[.!discrete] = plot(c[:, .!discrete, :], :density; args...)
  return plots
end

function traceplot(c::AbstractChains; legend::Bool=false, na...)
  nrows, nvars, nchains = size(c.value)
  plots = Array{Plot}(undef, nvars)
  pos = legend ? :right : :none
  for i in 1:nvars
    plots[i] = plot(y=vec(c.value[:, i, :]),
                    x=repeat(collect(c.range), outer=[nchains]),
                    Geom.line(),
                    color=repeat(c.chains, inner=[length(c.range)]),
                    Scale.color_discrete(), Guide.colorkey(title="Chain"),
                    Guide.xlabel("Iteration", orientation=:horizontal),
                    Guide.ylabel("Value", orientation=:vertical),
                    Guide.title(c.names[i]), Theme(key_position=pos))
  end
  return plots
end


text:
module StarProblem

import Random
import Base.copy, Base.size

ACTION_1 = 0
ACTION_2 = 1

PHI = [[1,2,0,0,0,0,0,0], # state 1
       [1,0,2,0,0,0,0,0], # state 2
       [1,0,0,2,0,0,0,0], # state 3
       [1,0,0,0,2,0,0,0], # state 4
       [1,0,0,0,0,2,0,0], # state 5
       [1,0,0,0,0,0,2,0], # state 6
       [2,0,0,0,0,0,0,1]] # state 7

INITIAL = [1,1,1,1,1,1,1,10]

mutable struct environment
    agent_pos::Int64
    environment() = new(rand(1:6))
end

is_terminal(env::environment) = false
reward(env::environment) = 0
get_state(env::environment) = PHI[env.agent_pos]
size(env::environment) = 8

function start!(env; rng=Random.GLOBAL_RNG)
    env.agent_pos = rand(rng, 1:6)
    return get_state(env)
end

function step(env::environment, action; rng=Random.GLOBAL_RNG)
    new_env = environment(env.agent_pos, env.size; rng=rng)
    return step!(new_env, action)
end

function step!(env::environment, action; rng=Random.GLOBAL_RNG)
    @assert !is_terminal(env)
    if action == ACTION_1
        # env.agent_pos -= 1
        env.agent_pos = 7
    elseif action == ACTION_2
        env.agent_pos = rand(rng, 1:6)
    else
        @error "Action invalid"
    end
    return get_state(env), reward(env), is_terminal(env)
end

TRUTH = [0]

end


text:
# ---
# title: 879. Profitable Schemes
# id: problem879
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/profitable-schemes/description/>
# hidden: true
# ---
#
# There is a group of `G` members, and a list of various crimes they could
# commit.
#
# The `ith` crime generates a `profit[i]` and requires `group[i]` members to
# participate in it.
#
# If a member participates in one crime, that member can't participate in
# another crime.
#
# Let's call a _profitable  scheme_ any subset of these crimes that generates at
# least `P` profit, and the total number of members participating in that subset
# of crimes is at most `G`.
#
# How many schemes can be chosen?  Since the answer may be very large, **return
# it modulo** `10^9 + 7`.
#
#
#
# **Example 1:**
#
#
#
#     Input: G = 5, P = 3, group = [2,2], profit = [2,3]
#     Output: 2
#     Explanation:
#     To make a profit of at least 3, the group could either commit crimes 0 and 1, or just crime 1.
#     In total, there are 2 schemes.
#
#
# **Example 2:**
#
#
#
#     Input: G = 10, P = 5, group = [2,3,5], profit = [6,7,8]
#     Output: 7
#     Explanation:
#     To make a profit of at least 5, the group could commit any crimes, as long as they commit one.
#     There are 7 possible schemes: (0), (1), (2), (0,1), (0,2), (1,2), and (0,1,2).
#
#
#
#
# **Note:**
#
#   1. `1 <= G <= 100`
#   2. `0 <= P <= 100`
#   3. `1 <= group[i] <= 100`
#   4. `0 <= profit[i] <= 100`
#   5. `1 <= group.length = profit.length <= 100`
#
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
"""
    grade_histogram(results, exam)
"""
function grade_histogram(results, exam)
    grades = get_grades(exam)

    passed = [g["passed"] ? :green : :red for g in exam["grades"]]

    counter_grades = counter(results[:, :grade])
    occ = [g ∈ keys(counter_grades) ? counter_grades[g] : 0 for g in grades]

    p = bar(occ, xticks=(1:length(grades), grades), color=passed, legend=false, dpi=300)
    ylabel!("Occurence")
    xlabel!("Grade")
    title!("Grade Distribution - " * exam["title"])

    savefig(p, exam["title"] * "_grade_histogram" * ".png")
end


text:
"""
    SeisGeometry(in;<keyword arguments>)

Update headers with geometry information. Offsets and azimuths are calculated from source and receivers coordinates.

# Arguments
* `in`: input filename

# Keyword arguments
* `ang=90`: inline direction measured in degrees CC from East
* `gamma=1`: vp/vs ratio for PS Asymptotic Conversion Point gathers (use gamma=1 for PP data)
* `osx=0`,`osy=0`,`ogx=0`,`ogy=0` : origin for source and receiver coordinate system
* `omx=0`,`omy=0`,`ohx=0`,`ohy=0`: origin for midpoint and offset coordinate system
* `oaz=0`,`oh=0` : origin for azimuth and offset coordinate system
* `dsx=1`,`dsy=1`,`dgx=1`,`dgy=1`: source and receiver step-size
* `dmx=1`,`dmy=1`,`dhx=1`,`dhy=1`: midpoint and offset step-size
* `dh=1`,`daz=1`: offset and azimuth step-size

# Outputs
the .seish file is updated with the following information:

* hx,hy,h,az,mx,my : calculated offset, azimuth and midpoint
* isx,isy,igx,igy,imx,imy,ihx,ihy,ih,iaz: calculated grid nodes for source and receiver position and midpoint, offset and azimuth.

*Credits: A. Stanton, 2017*

"""
function SeisGeometry(in; ang=90, gamma=1, osx=0, osy=0, ogx=0, ogy=0, omx=0,
                      omy=0, ohx=0, ohy=0, oh=0, oaz=0, dsx=1, dsy=1, dgx=1,
                      dgy=1, dmx=1, dmy=1, dhx=1, dhy=1, dh=1, daz=1)

    r2d = 180/pi
    d2r = pi/180
    gammainv = 1/gamma
    if (ang > 90)
	ang2=-d2r*(ang-90)
    else
	ang2=d2r*(90-ang)
    end

    filename = ParseHeaderName(in)
    stream = open(filename,"r+")
    nhead = 27
    nx = round(Int,filesize(stream)/(4*length(fieldnames(Header))))

    naz=convert(Int32,360/daz)

    for j=1:nx
	h = GrabHeader(stream,j)
	h.hx = h.gx - h.sx
	h.hy = h.gy - h.sy
	h.h = sqrt((h.hx^2) + (h.hy^2))
	h.az = r2d*atan((h.gy-h.sy),(h.gx-h.sx))
	if (h.az < 0)
	    h.az += 360.0
	end
	h.mx = h.sx + h.hx/(1 + gammainv);
	h.my = h.sy + h.hy/(1 + gammainv);
	sx_rot = (h.sx-osx)*cos(ang2) - (h.sy-osy)*sin(ang2) + osx;
	sy_rot = (h.sx-osx)*sin(ang2) + (h.sy-osy)*cos(ang2) + osy;
	gx_rot = (h.gx-ogx)*cos(ang2) - (h.gy-ogy)*sin(ang2) + ogx;
	gy_rot = (h.gx-ogx)*sin(ang2) + (h.gy-ogy)*cos(ang2) + ogy;
	mx_rot = (h.mx-omx)*cos(ang2) - (h.my-omy)*sin(ang2) + omx;
	my_rot = (h.mx-omx)*sin(ang2) + (h.my-omy)*cos(ang2) + omy;
	hx_rot = (h.hx-ohx)*cos(ang2) - (h.hy-ohy)*sin(ang2) + ohx;
	hy_rot = (h.hx-ohx)*sin(ang2) + (h.hy-ohy)*cos(ang2) + ohy;
	h.isx = convert(Int32,round((sx_rot-osx)/dsx))
	h.isy = convert(Int32,round((sy_rot-osy)/dsy))
	h.igx = convert(Int32,round((gx_rot-ogx)/dgx))
	h.igy = convert(Int32,round((gy_rot-ogy)/dgy))
	h.imx = convert(Int32,round((mx_rot-omx)/dmx))
	h.imy = convert(Int32,round((my_rot-omy)/dmy))
	h.ihx = convert(Int32,round((hx_rot-ohx)/dhx))
	h.ihy = convert(Int32,round((hy_rot-ohy)/dhy))
	h.ih = convert(Int32,round((h.h-oh)/dh))

	h.iaz = convert(Int32,round((h.az-oaz)/daz))<naz ? convert(Int32,round((h.az-oaz)/daz)) : 0

	PutHeader(stream,h,j)
  end
  close(stream)

end


text:
const PRECHELT_REF = "[Prechelt, Lutz (1998): \"Early Stopping"*
    "- But When?\", in *Neural Networks: Tricks of the Trade*, "*
    "ed. G. Orr, Springer.](https://link.springer.com/chapter"*
    "/10.1007%2F3-540-49430-8_3)"

const STOPPING_DOC = "An early stopping criterion for loss-reporting "*
    "iterative algorithms. "

## NEVER

"""
    Never()

$STOPPING_DOC

Indicates early stopping is to be disabled.

See also [`NotANumber`](@ref), for stopping on encountering `NaN`.

"""
struct Never <: StoppingCriterion end

## NOT A NUMBER

"""
    NotANumber()

$STOPPING_DOC

Stop if a loss of `NaN` is encountered.

"""
struct NotANumber <: StoppingCriterion end

# state = `true` when NaN has been encountered

update(::NotANumber, loss) = isnan(loss)
update_training(::NotANumber, loss) = isnan(loss)

update(::NotANumber, loss, state) = state || isnan(loss)
update_training(::NotANumber, loss, state) = state || isnan(loss)

done(::NotANumber, state) = state

message(::NotANumber, state) = "Stopping early as NaN encountered. "

## TIME LIMIT

"""
    TimeLimit(; t=0.5)

$STOPPING_DOC

Stopping is triggered after `t` hours have elapsed since the stopping
criterion was initiated.

Any Julia built-in `Real` type can be used for `t`. Subtypes of
`Period` may also be used, as in `TimeLimit(t=Minute(30))`.

Internally, `t` is rounded to nearest millisecond.

"""
struct TimeLimit <: StoppingCriterion
    t::Millisecond
    function TimeLimit(t::Millisecond)
        t > Millisecond(0) ||
            throw(ArgumentError("Time limit `t` must be positive. "))
        return new(t)
    end
end
TimeLimit(t::T) where T <: Period = TimeLimit(convert(Millisecond, t))
# for t::T a "numeric" time in hours; assumes `round(Int, ::T)` implemented:
TimeLimit(t) = TimeLimit(round(Int, 3_600_000*t) |> Millisecond)
TimeLimit(; t =Minute(30)) = TimeLimit(t)

# state = time at initialization

update(::TimeLimit, loss) = now()
update_training(::TimeLimit, loss) = now()
update(::TimeLimit, loss, state) = state
done(criterion::TimeLimit, state) = begin
    criterion.t < now() - state
end

## GL

# helper:

generalization_loss(E, E_opt) =  100*(E/abs(E_opt) - one(E_opt))

"""
    GL(; alpha=2.0)

$STOPPING_DOC

A stop is triggered when the (rescaled) generalization loss exceeds
the threshold `alpha`.

**Terminology.** Suppose ``E_1, E_2, ..., E_t`` are a sequence of
losses, for example, out-of-sample estimates of the loss associated
with some iterative machine learing algorithm. Then the
*generalization loss* at time `t`, is given by

`` GL_t = 100 (E_t - E_opt) \\over |E_opt|``

where `E_opt` is the minimum value of the sequence.

Reference: $PRECHELT_REF.

"""
struct GL <: StoppingCriterion
    alpha::Float64
    function GL(alpha)
        alpha > 0 ||
            throw(ArgumentError("Threshold `alpha` must be positive. "))
        return new(alpha)
    end
end
GL(; alpha=2.0) = GL(alpha)

update(::GL, loss) = (loss=loss, min_loss=loss)
update(::GL, loss, state) = (loss=loss, min_loss=min(loss, state.min_loss))
done(criterion::GL, state) =
    generalization_loss(state.loss, state.min_loss) > criterion.alpha

## PQ

# helpers:

function prepend(v, x, max_length)
    if length(v) < max_length
        vcat([x, ], v)
    else
        vcat([x, ], v[1:end-1])
    end
end

progress(losses::AbstractVector{T}) where T =
    1000*(mean(losses)/minimum(losses) - one(T))

_min(x, ::Nothing) = x
_min(::Nothing, x) = x
_min(x, y) = min(x, y)

"""
    PQ(; alpha=0.75, k=5, tol=eps(Float64))

A stopping criterion for training iterative supervised learners.

A stop is triggered when Prechelt's progress-modified generalization
loss exceeds the threshold `alpha`, or if the training progress drops
below `tol`. Here `k` is the maximum number of training (in-sample)
losses to be used to estimate the training progress.

**Context and explanation of terminology.** The progress-modified loss
is defined in the following scenario: Estimates, ``E_1, E_2, ...,
E_t``, of the out-of-sample loss of an iterative supervised learner
are being computed, but not necessarily at every iteration. However,
training losses for every iteration *are* being made available
(usually as a by-product of training) which can be used to quantify
recent training progress, as follows.

Fix a time ``j``, corresponding to some out-of-sample loss ``E_j``,
and let ``F_1`` be the corresponding training loss, ``F_2`` the
training loss in the previous interation of the model, ``F_3``, the
training loss two iterations previously, and so on. Let ``K`` denote
the number of model iterations since the last out-of-sample loss
``E_{j-1}`` was computed, or `k`, whichever is the smaller.  Then
the *training progress* at time ``j`` is defined by

`` P_j = 1000 |(M - m) \\over m| ``

where `M` is the mean of the training losses ``F_1, F_2, \\ldots ,
F_K`` and `m` the minimum value of those losses.

The *progress-modified generalization loss* at time ``t`` is given by

`` PQ_t = GL_t \\over P_t``

where ``GL_t`` is the generalization loss at time ``t``; see
[`GL`](@ref).

Reference: $PRECHELT_REF.

"""
struct PQ <: StoppingCriterion
    alpha::Float64
    k::Union{Int,Float64} # could be Inf
    tol::Float64
    function PQ(alpha, k, tol)
        alpha > 0 ||
            throw(ArgumentError("Threshold `alpha` must be positive. "))
        if k < 2 || !(k isa Int || isinf(k))
            throw(ArgumentError("`k` must `Int` or `Inf` and `k > 1`. "))
        end
        tol > 0 ||
            throw(ArgumentError("`tol` must be positive. "))

        return new(alpha, k, tol)
    end
end
PQ(; alpha=0.75, k=5, tol=eps(Float64)) = PQ(alpha, k, tol)

struct PQState{T}
    training_losses::Vector{T}
    waiting_for_out_of_sample::Bool
    loss::Union{Nothing,T}
    min_loss::Union{Nothing,T}
end

function update_training(criterion::PQ, loss)
    training_losses = [loss, ]
    return PQState(training_losses,
                   true,
                   nothing,
                   nothing)
end

update(::PQ, loss) = error("First loss reported to the GL early stopping "*
                           "algorithm must be a training loss. ")

function update_training(criterion::PQ, loss, state)
    training_losses = if state.waiting_for_out_of_sample
        prepend(state.training_losses, loss, criterion.k)
    else
        [loss, ]
    end
    return PQState(training_losses,
                   true,
                   state.loss,
                   state.min_loss)
end

function update(::PQ, loss, state)
    length(state.training_losses) > 1 ||
        error("The PQ stopping criterion requires at least two training "*
              "losses between out-of-sample loss updates. ")
    return PQState(state.training_losses,
                   false,
                   loss,
                   _min(loss, state.min_loss))
end

function done(criterion::PQ, state)
    state.waiting_for_out_of_sample && return false
    GL = generalization_loss(state.loss, state.min_loss)
    P = progress(state.training_losses)
    P > criterion.tol || return true
    PQ = GL/P
    return  PQ > criterion.alpha
end

needs_in_and_out_of_sample(::Type{<:PQ}) = true

## PATIENCE

# This is UP_s in Prechelt 1998

"""
    Patience(; n=5)

$STOPPING_DOC

A stop is triggered by `n` consecutive increases in the loss.

Denoted "_UP_s" in $PRECHELT_REF.

"""
mutable struct Patience <: StoppingCriterion
    n::Int
    function Patience(n::Int)
        n > 0 ||
            throw(ArgumentError("The patience level `n` must be positive. "))
        return new(n)
    end
end
Patience(; n=5) = Patience(n)

# Prechelt alias:
const UP = Patience

update(criterion::Patience, loss) = (loss=loss, n_increases=0)
@inline function update(criterion::Patience, loss, state)
    old_loss, n = state
    if loss > old_loss
        n += 1
    else
        n = 0
    end
    return (loss=loss, n_increases=n)
end

done(criterion::Patience, state) = state.n_increases == criterion.n


text:
struct SAFTVRSWParam <: EoSParam
    Mw::SingleParam{Float64}
    segment::SingleParam{Float64}
    sigma::PairParam{Float64}
    lambda::PairParam{Float64}
    epsilon::PairParam{Float64}
    epsilon_assoc::AssocParam{Float64}
    bondvol::AssocParam{Float64}
end

abstract type SAFTVRSWModel <: SAFTModel end
@newmodel SAFTVRSW SAFTVRSWModel SAFTVRSWParam

export SAFTVRSW
function SAFTVRSW(components;
    idealmodel=BasicIdeal,
    userlocations=String[],
    ideal_userlocations=String[],
    verbose=false,
    assoc_options = AssocOptions())

    params,sites = getparams(components, ["SAFT/SAFTVRSW","properties/molarmass.csv"]; userlocations=userlocations, verbose=verbose)

    segment = params["m"]
    k = params["k"]
    Mw = params["Mw"]
    params["sigma"].values .*= 1E-10
    sigma = sigma_LorentzBerthelot(params["sigma"])
    epsilon = epsilon_LorentzBerthelot(params["epsilon"], k)
    lambda = lambda_squarewell(params["lambda"], sigma)

    epsilon_assoc = params["epsilon_assoc"]
    bondvol = params["bondvol"]

    packagedparams = SAFTVRSWParam(Mw, segment, sigma, lambda, epsilon, epsilon_assoc, bondvol)
    references = ["todo"]

    model = SAFTVRSW(packagedparams, sites, idealmodel; ideal_userlocations, references, verbose, assoc_options)
    return model
end

function a_res(model::SAFTVRSWModel, V, T, z)
    return @f(a_mono) + @f(a_chain) + @f(a_assoc)
end

function a_mono(model::SAFTVRSWModel, V, T, z)
    return @f(a_hs) + @f(a_disp)
end

function a_disp(model::SAFTVRSWModel, V, T, z)
    return @f(a_1) + @f(a_2)
end

function a_hs(model::SAFTVRSWModel, V, T, z)
    ζ0 = @f(ζn,0)
    ζ1 = @f(ζn,1)
    ζ2 = @f(ζn,2)
    ζ3 = @f(ζn,3)
    Σz = sum(z)
    m = model.params.segment.values
    m̄ = dot(z, m)/Σz
    return m̄*6/π/@f(ρ_S)*(3ζ1*ζ2/(1-ζ3) + ζ2^3/(ζ3*(1-ζ3)^2) + (ζ2^3/ζ3^2-ζ0)*log(1-ζ3))
end

function ζn(model::SAFTVRSWModel, V, T, z, n)
    σ = model.params.sigma.diagvalues
    return π/6*@f(ρ_S)*∑(@f(x_S,i)*σ[i]^n for i ∈ @comps)
end

function ρ_S(model::SAFTVRSWModel, V, T, z)
    Σz = sum(z)
    m = model.params.segment.values
    m̄ = dot(z, m)
    N = N_A
    return N/V*m̄
end

function x_S(model::SAFTVRSWModel, V, T, z, i)
    Σz = sum(z)
    m = model.params.segment.values
    m̄ = dot(z, m)
    return z[i]*m[i]/m̄
end

function ζ_X(model::SAFTVRSWModel, V, T, z)
    comps = @comps
    σ = model.params.sigma.values
    return π/6*@f(ρ_S)*∑(@f(x_S,i)*@f(x_S,j)*σ[i,j]^3 for i ∈ comps for j ∈ comps)
end

function a_1(model::SAFTVRSWModel, V, T, z)
    comps = @comps
    Σz = sum(z)
    m = model.params.segment.values
    m̄ = dot(z, m)/Σz
    return -m̄/T*@f(ρ_S)*∑(@f(x_S,i)*@f(x_S,j)*@f(a_1,i,j) for i ∈ comps for j ∈ comps)
end

function a_1(model::SAFTVRSWModel, V, T, z, i, j)
    ϵ = model.params.epsilon.values
    λ = model.params.lambda.values
    σ = model.params.sigma.values
    αVDWij = 2π*ϵ[i,j]*σ[i,j]^3*(λ[i,j]^3-1)/3
    return αVDWij * @f(gHS_0,i,j)
end

function ζeff_X(model::SAFTVRSWModel, V, T, z, λ)
    A = SAFTVRSWconsts.A
    ζ_X_ = @f(ζ_X)
    return A * [1; λ; λ^2] ⋅ [ζ_X_; ζ_X_^2; ζ_X_^3]
end

function gHS_0(model::SAFTVRSWModel,V, T, z, i, j)
    λ = model.params.lambda.values
    ζeff_X_ = @f(ζeff_X,λ[i,j])
    return (1-ζeff_X_/2)/(1-ζeff_X_)^3
end

function a_2(model::SAFTVRSWModel, V, T, z)
    comps = @comps
    Σz = sum(z)
    m = model.params.segment.values
    m̄ = dot(z, m)/Σz
    return m̄/T^2*∑(@f(x_S,i)*@f(x_S,j)*@f(a_2,i,j) for i ∈ comps for j ∈ comps)
end

function a_2(model::SAFTVRSWModel, V, T, z, i, j)
    ϵ = model.params.epsilon.values
    ζ0 = @f(ζn,0)
    ζ1 = @f(ζn,1)
    ζ2 = @f(ζn,2)
    ζ3 = @f(ζn,3)
    KHS = ζ0*(1-ζ3)^4/(ζ0*(1-ζ3)^2+6*ζ1*ζ2*(1-ζ3)+9*ζ2^3)
    return 1/2*KHS*ϵ[i,j]*@f(∂a_1╱∂ρ_S,i,j)
end

function ∂a_1╱∂ρ_S(model::SAFTVRSWModel, V, T, z, i, j)
    ϵ = model.params.epsilon.values
    λ = model.params.lambda.values
    σ = model.params.sigma.values
    αij = 2π*ϵ[i,j]*σ[i,j]^3*(λ[i,j]^3-1)/3
    ζ_X_ = @f(ζ_X)
    ζeff_X_ = @f(ζeff_X,λ[i,j])
    A = SAFTVRSWconsts.A
    ∂ζeff_X╱∂ζ_X = A * [1; λ[i,j]; λ[i,j]^2] ⋅ [ζ_X_; 2ζ_X_^2; 3ζ_X_^3]
    # ∂ζeff_X╱∂ζ_X = A * [1; λ[i,j]; λ[i,j]^2] ⋅ [1; 2ζ_X_; 3ζ_X_^2]
    return -αij*(@f(ρ_S)*@f(gHS_0,i,j)+(5/2-ζeff_X_)/(1-ζeff_X_)^4*∂ζeff_X╱∂ζ_X)
end

function a_chain(model::SAFTVRSWModel, V, T, z)
    Σz = sum(z)
    m = model.params.segment.values
    return -∑(z[i]*(log(@f(γSW,i))*(m[i]-1)) for i ∈ @comps)/Σz
end

function γSW(model::SAFTVRSWModel,V, T, z, i)
    ϵ = model.params.epsilon.diagvalues
    return @f(gSW,i,i)*exp(-ϵ[i]/T)
end

function gSW(model::SAFTVRSWModel,V, T, z, i, j)
    ϵ = model.params.epsilon.values
    return @f(gHS,i,j)+ϵ[i,j]/T*@f(g_1,i,j)
end

function gHS(model::SAFTVRSWModel,V, T, z, i, j)
    σ = model.params.sigma.values
    ζ3 = @f(ζn,3)
    D = σ[i]*σ[j]/(σ[i]+σ[j])*∑(@f(x_S,k)*σ[k]^2 for k ∈ @comps)/∑(@f(x_S,k)*σ[k]^3 for k ∈ @comps)
    return 1/(1-ζ3)+3*D*ζ3/(1-ζ3)^2+2*(D*ζ3)^2/(1-ζ3)^3
end

function g_1(model::SAFTVRSWModel,V, T, z, i, j)
    λ = model.params.lambda.values
    ζ_X_ = @f(ζ_X)
    ζeff_X_ = @f(ζeff_X,λ[i,j])
    A = SAFTVRSWconsts.A
    ∂ζeff_X╱∂ζ_X = A * [1; λ[i,j]; λ[i,j]^2] ⋅ [1; 2ζ_X_; 3ζ_X_^2]
    ∂ζeff_X╱∂λ = A * [0; 1; 2λ[i,j]] ⋅ [ζ_X_; ζ_X_^2; ζ_X_^3]
    return @f(gHS_0,i,j)+(λ[i,j]^3-1)*(5/2-ζeff_X_)/(1-ζeff_X_)^4*(λ[i,j]/3*∂ζeff_X╱∂λ-ζ_X_*∂ζeff_X╱∂ζ_X)
end

function Δ(model::SAFTVRSWModel, V, T, z, i, j, a, b)
    ϵ_assoc = model.params.epsilon_assoc.values
    κ = model.params.bondvol.values
    g = @f(gSW,i,j)
    return g*(exp(ϵ_assoc[i,j][a,b]/T)-1)*κ[i,j][a,b]
end

const SAFTVRSWconsts = (
    A = [2.25855   -1.50349  0.249434;
    -0.66927  1.40049   -0.827739;
    10.1576   -15.0427   5.30827],
)


text:
# extract items useful for plotting from the history generated by the
# fallback `result` method:
function plotting_report(fields, scales, history)

    n_models = length(history)
    n_parameters = length(fields)

    A = Array{Any}(undef, (n_models, n_parameters))
    measurements = Vector{Float64}(undef, n_models)

    for j in eachindex(history)
        m, r = history[j]
        A[j,:] = [recursive_getproperty(m, fld) for fld in fields]
        measurements[j] = last(r)[1]
    end

    return plotting=(parameter_names=string.(fields) |> collect,
                     parameter_scales=scales |> collect,
                     parameter_values = A,
                     measurements = measurements)

end