text:
abstract type SimpleConeScaling{T <: AbstractFloat} end

@with_kw mutable struct NonnegConeScaling{T <: AbstractFloat} <: SimpleConeScaling{T}
    w::Vector{T}
    λ::Vector{T}
    @assert !any(isnan.(w))
    @assert !any(isnan.(λ))
end

NonnegConeScaling(C::NonnegCone) = NonnegConeScaling(get_id(C), get_id(C))

@with_kw mutable struct QuadConeScaling{T <: AbstractFloat} <: SimpleConeScaling{T}
    w::Vector{T}
    λ::Vector{T}
    η::T
    @assert !any(isnan.(w))
    @assert !any(isnan.(λ))
end

QuadConeScaling(C::QuadCone) = QuadConeScaling(get_id(C), get_id(C), 1.)

get_scaling(C::NonnegCone) = NonnegConeScaling(C)
get_scaling(C::QuadCone) = QuadConeScaling(C)

@with_kw mutable struct CompositeConeScaling{T <: AbstractFloat}
    scalings::Vector{SimpleConeScaling{T}}
end

CompositeConeScaling(C::CompositeCone{T}) where T = CompositeConeScaling{T}(get_scaling.(C))

iterate(W::CompositeConeScaling, args...) = iterate(W.scalings, args...)
getindex(W::CompositeConeScaling, args...) = getindex(W.scalings, args...)
length(W::CompositeConeScaling) = length(W.scalings)
size(W::CompositeConeScaling) = (length(W),)

dim(W::NonnegConeScaling) = length(W.w)
dim(W::QuadConeScaling)  = length(W.w)
dim(W::CompositeConeScaling) = sum(map(dim, W))

get_scaled_var(W::NonnegConeScaling) = W.λ
get_scaled_var(W::QuadConeScaling) = W.λ
function get_scaled_var(W::CompositeConeScaling{T}) where T
    return vcat(T[], get_scaled_var.(W)...)
end

*(W::NonnegConeScaling{T}, x::Vector{T}) where T <: AbstractFloat = W.w .* x
*(x::Vector{T}, W::NonnegConeScaling{T}) where T <: AbstractFloat = W * x
\(W::NonnegConeScaling{T}, x::Vector{T}) where T <: AbstractFloat = W.w .\ x
/(x::Vector{T}, W::NonnegConeScaling{T}) where T <: AbstractFloat = W \ x

function mult_or_div(W::QuadConeScaling{T}, x::Vector{T}; is_mult=true) where T <: AbstractFloat
    y = Vector{T}(undef, length(x))
    sign = is_mult ? 1 : -1
    w = W.w

    wTx = w[2:end]'*x[2:end]
    y[1] = w[1]*x[1] + sign*wTx
    y[2:end] = sign*w[2:end]*x[1] + x[2:end] + w[2:end]*wTx/(w[1] + 1)
    y *= W.η^sign
    return y
end

*(W::QuadConeScaling{T}, x::Vector{T}) where T <: AbstractFloat = mult_or_div(W, x, is_mult=true)
\(W::QuadConeScaling{T}, x::Vector{T}) where T <: AbstractFloat = mult_or_div(W, x, is_mult=false)

function *(W::CompositeConeScaling{T}, x::Vector{T}) where T <: AbstractFloat
    x_list = split_by_cone(x, W)
    return vcat(T[], (W .* x_list)...)
end

function \(W::CompositeConeScaling{T}, x::Vector{T}) where T <: AbstractFloat
    x_list = split_by_cone(x, W)
    y = Vector{T}(undef, length(x))
    idx = 1
    for Wi in W
        dim_Wi = dim(Wi)
        y[idx: idx+dim_Wi-1] = Wi \ x[idx: idx+dim_Wi-1]
        idx += dim_Wi
    end
    return y
end

function update!(W::NonnegConeScaling{T}, s::Vector{T}, z::Vector{T}) where T <: AbstractFloat
    W.w = sqrt.(s ./ z)
    W.λ = sqrt.(s) .* sqrt.(z)
end

function update!(W::QuadConeScaling{T}, s::Vector{T}, z::Vector{T}) where T <: AbstractFloat
    zbar = z / sqrt(pos(z'*(J*z)))
    sbar = s / sqrt(pos(s'*(J*s)))
    γ = sqrt((1 + zbar'*sbar)/2)
    w = .5*(sbar + J*zbar)/γ
    W.w = w
    W.η = pos((s'*(J*s)) / (z'*(J*z)))^(1/4)
    W.λ = W * z
end

function update!(W::CompositeConeScaling{T}, s::Vector{T}, z::Vector{T}) where T <: AbstractFloat
    s_list = split_by_cone(s, W)
    z_list = split_by_cone(z, W)
    update!.(W, s_list, z_list)
end

get_num_quad_cones(W::CompositeConeScaling) = sum([isa(Wi, QuadConeScaling) for Wi in W])

function quad_cone_idxs(W::CompositeConeScaling)
    n_qc = get_num_quad_cones(W)
    start_idxs = Vector{Int64}(undef, n_qc)
    wbars = Vector{Vector{Float64}}(undef, n_qc)
    idx = 1
    for Wi in W
        if isa(Wi, QuadConeScaling)
            wbars[i] = Wi.wbar
            start_idxs[i] = idx
            idx += length(Wi.wbar)
        end
    end
    return wbars, start_idxs
end

function get_lower_kkt_block(W::CompositeConeScaling{T}; identity=false) where T
    p = dim(W)
    q = get_num_quad_cones(W)
    Kl = spzeros(T, p + 2q, p + 2q)
    i = 1  # Cone dimension offset.
    k = 1  # Quad cone index.
    j = 1
    for Wi in W
        nW = dim(Wi)
        if isa(Wi, NonnegConeScaling)
            if identity
                Kl[i:i+nW-1, i:i+nW-1] = speye(dim(Wi))
            else
                Kl[i:i+nW-1, i:i+nW-1] = sparse(Diagonal(Wi.w.^2))
            end
        elseif isa(Wi, QuadConeScaling)
            if identity
                u = 1e-20*ones(T, dim(Wi))
                v = 1e-20*ones(T, dim(Wi))
                η2 = 1.
                D = speye(T, dim(Wi))
            else
                w = Wi.w
                η2 = Wi.η^2
                w0 = w[1]
                w1 = w[2:end]
                α = 1 + w0 + w1'*w1/(1 + w0)
                β = 1 + 2/(1 + w0) + w1'*w1/(1 + w0)^2
                d = .5*(w0^2 + w1'*w1*(1 - α^2/(1 + w1'*w1*β)))
                u0 = sqrt(w0^2 + w1'*w1 - d)
                u1 = α/u0
                v1 = sqrt(pos(u1^2 - β))
                u = [u0; u1*w1]
                v = [0; v1*w1]
                D = sparse(Diagonal([d; ones(T, nW-1)]))
            end

            full_block = true
            if full_block && !identity
                Wm = Wi.η*[w0  w1'; w1  speye(nW-1) + w1*w1'/(w0+1)]
                W2 = Wm'*Wm
                Kl[i:i+nW-1, i:i+nW-1] = W2
                Kl[p+k, p+k] = 1.
                Kl[p+k+1, p+k+1] = 1.
            else
                Kl[i:i+nW-1, i:i+nW-1] = η2*D
                Kl[p+k, i:i+nW-1] = η2*v'
                Kl[i:i+nW-1, p+k] = η2*v
                Kl[p+k+1, i:i+nW-1] = η2*u'
                Kl[i:i+nW-1, p+k+1] = η2*u
                Kl[p+k, p+k] = η2
                Kl[p+k+1, p+k+1] = -η2
            end
            k += 2
            j += nW
        else
            error("Unrecognized cone.")
        end
        i += nW
    end
    return Kl
end

function split_by_cone(x::Vector{T}, C::Union{CompositeCone{T}, CompositeConeScaling{T}}) where T
    j = 0
    x_list = []
    for Ci in C
        push!(x_list, x[j+1:j+dim(Ci)])
        j += dim(Ci)
    end
    return x_list
end


text:
# TODO: only store last N trees
_id = 0

treeid() = (global _id += 1)

const trees = Dict{Int,LazyTree}()

function render(i::Inline, tree::LazyTree)
  id = treeid()
  trees[id] = tree
  d(:type => :lazy,
    :head => render(i, tree.head),
    :id => id)
end

handle("getlazy") do id
  haskey(trees, id) || return [render(Inline(), fade("[out of date result]"))]
  return [render(Inline(), x) for x in pop!(trees, id).children()]
end

handle("clearLazy") do ids
  for id in ids
    delete!(trees, id)
  end
end

ext: jl


text:
# create a new "build_series_args" which converts all inputs into xs = Any[xitems], ys = Any[yitems].
# Special handling for: no args, xmin/xmax, parametric, dataframes
# Then once inputs have been converted, build the series args, map functions, etc.
# This should cut down on boilerplate code and allow more focused dispatch on type
# note: returns meta information... mainly for use with automatic labeling from DataFrames for now

const FuncOrFuncs{F} = Union{F, Vector{F}, Matrix{F}}
const MaybeNumber = Union{Number, Missing}
const MaybeString = Union{AbstractString, Missing}
const DataPoint = Union{MaybeNumber, MaybeString}

prepareSeriesData(x) = error("Cannot convert $(typeof(x)) to series data for plotting")
prepareSeriesData(::Nothing) = nothing
prepareSeriesData(t::Tuple{T, T}) where {T<:Number} = t
prepareSeriesData(f::Function) = f
prepareSeriesData(ar::AbstractRange{<:Number}) = ar
function prepareSeriesData(a::AbstractArray{<:MaybeNumber})
    f = isimmutable(a) ? replace : replace!
    a = f(x -> ismissing(x) || isinf(x) ? NaN : x, map(float, a))
end
prepareSeriesData(a::AbstractArray{<:Missing}) = fill(NaN, axes(a))
prepareSeriesData(a::AbstractArray{<:MaybeString}) = replace(x -> ismissing(x) ? "" : x, a)
prepareSeriesData(s::Surface{<:AMat{<:MaybeNumber}}) = Surface(prepareSeriesData(s.surf))
prepareSeriesData(s::Surface) = s  # non-numeric Surface, such as an image
prepareSeriesData(v::Volume) = Volume(prepareSeriesData(v.v), v.x_extents, v.y_extents, v.z_extents)

# default: assume x represents a single series
series_vector(x, plotattributes) = [prepareSeriesData(x)]

# fixed number of blank series
series_vector(n::Integer, plotattributes) = [zeros(0) for i in 1:n]

# vector of data points is a single series
series_vector(v::AVec{<:DataPoint}, plotattributes) = [prepareSeriesData(v)]

# list of things (maybe other vectors, functions, or something else)
function series_vector(v::AVec, plotattributes)
    if all(x -> x isa MaybeNumber, v)
        series_vector(Vector{MaybeNumber}(v), plotattributes)
    elseif all(x -> x isa MaybeString, v)
        series_vector(Vector{MaybeString}(v), plotattributes)
    else
        vcat((series_vector(vi, plotattributes) for vi in v)...)
    end
end

# Matrix is split into columns
function series_vector(v::AMat{<:DataPoint}, plotattributes)
    if all3D(plotattributes)
        [prepareSeriesData(Surface(v))]
    else
        [prepareSeriesData(v[:, i]) for i in axes(v, 2)]
    end
end

# --------------------------------------------------------------------
# Fillranges & ribbons

process_fillrange(range::Number, plotattributes) = [range]
process_fillrange(range, plotattributes) = series_vector(range, plotattributes)

process_ribbon(ribbon::Number, plotattributes) = [ribbon]
process_ribbon(ribbon, plotattributes) = series_vector(ribbon, plotattributes)
# ribbon as a tuple: (lower_ribbons, upper_ribbons)
process_ribbon(ribbon::Tuple{S, T}, plotattributes) where {S, T} = collect(zip(
    series_vector(ribbon[1], plotattributes),
    series_vector(ribbon[2], plotattributes),
))

# --------------------------------------------------------------------

compute_x(x::Nothing, y::Nothing, z)    = axes(z,1)
compute_x(x::Nothing, y, z)             = axes(y,1)
compute_x(x::Function, y, z)            = map(x, y)
compute_x(x, y, z)                      = x

compute_y(x::Nothing, y::Nothing, z)    = axes(z,2)
compute_y(x, y::Function, z)            = map(y, x)
compute_y(x, y, z)                      = y

compute_z(x, y, z::Function)            = map(z, x, y)
compute_z(x, y, z::AbstractMatrix)      = Surface(z)
compute_z(x, y, z::Nothing)             = nothing
compute_z(x, y, z)                      = z

nobigs(v::AVec{BigFloat}) = map(Float64, v)
nobigs(v::AVec{BigInt}) = map(Int64, v)
nobigs(v) = v

@noinline function compute_xyz(x, y, z)
    x = compute_x(x,y,z)
    y = compute_y(x,y,z)
    z = compute_z(x,y,z)
    nobigs(x), nobigs(y), nobigs(z)
end

# not allowed
compute_xyz(x::Nothing, y::FuncOrFuncs{F}, z) where {F<:Function}       = error("If you want to plot the function `$y`, you need to define the x values!")
compute_xyz(x::Nothing, y::Nothing, z::FuncOrFuncs{F}) where {F<:Function} = error("If you want to plot the function `$z`, you need to define x and y values!")
compute_xyz(x::Nothing, y::Nothing, z::Nothing)        = error("x/y/z are all nothing!")

# --------------------------------------------------------------------

# we are going to build recipes to do the processing and splitting of the args

# --------------------------------------------------------------------
# The catch-all SliceIt recipe
# --------------------------------------------------------------------

# ensure we dispatch to the slicer
struct SliceIt end

# The `SliceIt` recipe finishes user and type recipe processing.
# It splits processed data into individual series data, stores in copied `plotattributes`
# for each series and returns no arguments.
@recipe function f(::Type{SliceIt}, x, y, z)

    # handle data with formatting attached
    if typeof(x) <: Formatted
        xformatter := x.formatter
        x = x.data
    end
    if typeof(y) <: Formatted
        yformatter := y.formatter
        y = y.data
    end
    if typeof(z) <: Formatted
        zformatter := z.formatter
        z = z.data
    end

    xs = series_vector(x, plotattributes)
    ys = series_vector(y, plotattributes)
    zs = series_vector(z, plotattributes)

    fr = pop!(plotattributes, :fillrange, nothing)
    fillranges = process_fillrange(fr, plotattributes)
    mf = length(fillranges)

    rib = pop!(plotattributes, :ribbon, nothing)
    ribbons = process_ribbon(rib, plotattributes)
    mr = length(ribbons)

    mx = length(xs)
    my = length(ys)
    mz = length(zs)
    if mx > 0 && my > 0 && mz > 0
        for i in 1:max(mx, my, mz)
            # add a new series
            di = copy(plotattributes)
            xi, yi, zi = xs[mod1(i,mx)], ys[mod1(i,my)], zs[mod1(i,mz)]
            di[:x], di[:y], di[:z] = compute_xyz(xi, yi, zi)

            # handle fillrange
            fr = fillranges[mod1(i,mf)]
            di[:fillrange] = isa(fr, Function) ? map(fr, di[:x]) : fr

            # handle ribbons
            rib = ribbons[mod1(i,mr)]
            di[:ribbon] = isa(rib, Function) ? map(rib, di[:x]) : rib

            push!(series_list, RecipeData(di, ()))
        end
    end
    nothing  # don't add a series for the main block
end

# --------------------------------------------------------------------
# Apply type recipes
# --------------------------------------------------------------------

# this is the default "type recipe"... just pass the object through
@recipe f(::Type{T}, v::T) where T = v

# this should catch unhandled "series recipes" and error with a nice message
@recipe f(::Type{V}, x, y, z) where {V<:Val} = error("The backend must not support the series type $V, and there isn't a series recipe defined.")

function _apply_type_recipe(plotattributes, v, letter)
    _preprocess_axis_args!(plotattributes, letter)
    rdvec = RecipesBase.apply_recipe(plotattributes, typeof(v), v)
    warn_on_recipe_aliases!(plotattributes, :type, typeof(v))
    _postprocess_axis_args!(plotattributes, letter)
    return rdvec[1].args[1]
end

# Handle type recipes when the recipe is defined on the elements.
# This sort of recipe should return a pair of functions... one to convert to number,
# and one to format tick values.
function _apply_type_recipe(plotattributes, v::AbstractArray, letter)
    _preprocess_axis_args!(plotattributes, letter)
    # First we try to apply an array type recipe.
    w = RecipesBase.apply_recipe(plotattributes, typeof(v), v)[1].args[1]
    warn_on_recipe_aliases!(plotattributes, :type, typeof(v))
    # If the type did not change try it element-wise
    if typeof(v) == typeof(w)
        isempty(skipmissing(v)) && return Float64[]
        x = first(skipmissing(v))
        args = RecipesBase.apply_recipe(plotattributes, typeof(x), x)[1].args
        warn_on_recipe_aliases!(plotattributes, :type, typeof(x))
        _postprocess_axis_args!(plotattributes, letter)
        if length(args) == 2 && all(arg -> arg isa Function, args)
            numfunc, formatter = args
             return Formatted(map(numfunc, v), formatter)
         else
             return v
        end
    end
    _postprocess_axis_args!(plotattributes, letter)
    return w
end

# special handling for Surface... need to properly unwrap and re-wrap
_apply_type_recipe(
    plotattributes,
    v::Surface{<:AMat{<:Union{AbstractFloat, Integer, AbstractString, Missing}}},
) = v
function _apply_type_recipe(plotattributes, v::Surface)
    ret = _apply_type_recipe(plotattributes, v.surf)
    if typeof(ret) <: Formatted
        Formatted(Surface(ret.data), ret.formatter)
    else
        Surface(ret.data)
    end
end

# don't do anything vectors of datapoints and for nothing
_apply_type_recipe(plotattributes, v::Nothing, letter) = v
_apply_type_recipe(
    plotattributes,
    v::AbstractArray{<:Union{AbstractFloat, Integer, AbstractString, Missing}},
    letter,
) = v

# axis args before type recipes should still be mapped to all axes
function _preprocess_axis_args!(plotattributes)
    for (k, v) in plotattributes
        if is_axis_attr_noletter(k)
            pop!(plotattributes, k)
            for l in (:x, :y, :z)
                lk = Symbol(l, k)
                haskey(plotattributes, lk) || (plotattributes[lk] = v)
            end
        end
    end
end
function _preprocess_axis_args!(plotattributes, letter)
    plotattributes[:letter] = letter
    _preprocess_axis_args!(plotattributes)
end

# axis args in type recipes should only be applied to the current axis
function _postprocess_axis_args!(plotattributes, letter)
    pop!(plotattributes, :letter)
    if letter in (:x, :y, :z)
        for (k, v) in plotattributes
            if is_axis_attr_noletter(k)
                pop!(plotattributes, k)
                lk = Symbol(letter, k)
                haskey(plotattributes, lk) || (plotattributes[lk] = v)
            end
        end
    end
end

# --------------------------------------------------------------------
# Fallback user recipes calling type recipes
# --------------------------------------------------------------------

# handle "type recipes" by converting inputs, and then either re-calling or slicing
@recipe function f(x, y, z)
    wrap_surfaces!(plotattributes, x, y, z)
    did_replace = false
    newx = _apply_type_recipe(plotattributes, x, :x)
    x === newx || (did_replace = true)
    newy = _apply_type_recipe(plotattributes, y, :y)
    y === newy || (did_replace = true)
    newz = _apply_type_recipe(plotattributes, z, :z)
    z === newz || (did_replace = true)
    if did_replace
        newx, newy, newz
    else
        SliceIt, x, y, z
    end
end
@recipe function f(x, y)
    wrap_surfaces!(plotattributes, x, y)
    did_replace = false
    newx = _apply_type_recipe(plotattributes, x, :x)
    x === newx || (did_replace = true)
    newy = _apply_type_recipe(plotattributes, y, :y)
    y === newy || (did_replace = true)
    if did_replace
        newx, newy
    else
        SliceIt, x, y, nothing
    end
end
@recipe function f(y)
    wrap_surfaces!(plotattributes, y)
    newy = _apply_type_recipe(plotattributes, y, :y)
    if y !== newy
        newy
    else
        SliceIt, nothing, y, nothing
    end
end

# if there's more than 3 inputs, it can't be passed directly to SliceIt
# so we'll apply_type_recipe to all of them
@recipe function f(v1, v2, v3, v4, vrest...)
    did_replace = false
    newargs = map(v -> begin
        newv = _apply_type_recipe(plotattributes, v, :unknown)
        if newv !== v
            did_replace = true
        end
        newv
    end, (v1, v2, v3, v4, vrest...))
    if !did_replace
        error("Couldn't process recipe args: $(map(typeof, (v1, v2, v3, v4, vrest...)))")
    end
    newargs
end

# helper function to ensure relevant attributes are wrapped by Surface
function wrap_surfaces!(plotattributes, args...) end
wrap_surfaces!(plotattributes, x::AMat, y::AMat, z::AMat) = wrap_surfaces!(plotattributes)
wrap_surfaces!(plotattributes, x::AVec, y::AVec, z::AMat) = wrap_surfaces!(plotattributes)
function wrap_surfaces!(plotattributes, x::AVec, y::AVec, z::Surface)
    wrap_surfaces!(plotattributes)
end
function wrap_surfaces!(plotattributes)
    if haskey(plotattributes, :fill_z)
        v = plotattributes[:fill_z]
        if !isa(v, Surface)
            plotattributes[:fill_z] = Surface(v)
        end
    end
end

# --------------------------------------------------------------------
# 1 argument
# --------------------------------------------------------------------

@recipe f(n::Integer) = is3d(get(plotattributes,:seriestype,:path)) ? (SliceIt, n, n, n) : (SliceIt, n, n, nothing)

all3D(plotattributes) = trueOrAllTrue(
    st -> st in (
        :contour,
        :contourf,
        :heatmap,
        :surface,
        :wireframe,
        :contour3d,
        :image,
        :plots_heatmap,
    ),
    get(plotattributes, :seriestype, :none),
)

# return a surface if this is a 3d plot, otherwise let it be sliced up
@recipe function f(mat::AMat)
    if all3D(plotattributes)
        n, m = axes(mat)
        m, n, Surface(mat)
    else
        nothing, mat, nothing
    end
end

# if a matrix is wrapped by Formatted, do similar logic, but wrap data with Surface
@recipe function f(fmt::Formatted{<:AMat})
    if all3D(plotattributes)
        mat = fmt.data
        n, m = axes(mat)
        m, n, Formatted(Surface(mat), fmt.formatter)
    else
        nothing, fmt, nothing
    end
end

# assume this is a Volume, so construct one
@recipe function f(vol::AbstractArray{<:MaybeNumber, 3}, args...)
    seriestype := :volume
    SliceIt, nothing, Volume(vol, args...), nothing
end

# images - grays
function clamp_greys!(mat::AMat{<:Gray})
    for i in eachindex(mat)
        mat[i].val < 0 && (mat[i] = Gray(0))
        mat[i].val > 1 && (mat[i] = Gray(1))
    end
    mat
end

@recipe function f(mat::AMat{<:Gray})
    n, m = axes(mat)
    if is_seriestype_supported(:image)
        seriestype := :image
        yflip --> true
        SliceIt, m, n, Surface(clamp_greys!(mat))
    else
        seriestype := :heatmap
        yflip --> true
        cbar --> false
        fillcolor --> ColorGradient([:black, :white])
        SliceIt, m, n, Surface(clamp!(convert(Matrix{Float64}, mat), 0., 1.))
    end
end

# images - colors
@recipe function f(mat::AMat{T}) where T<:Colorant
	n, m = axes(mat)

    if is_seriestype_supported(:image)
        seriestype := :image
        yflip --> true
        SliceIt, m, n, Surface(mat)
    else
        seriestype := :heatmap
        yflip --> true
        cbar --> false
        aspect_ratio --> :equal
        z, plotattributes[:fillcolor] = replace_image_with_heatmap(mat)
        SliceIt, m, n, Surface(z)
    end
end

# plotting arbitrary shapes/polygons

@recipe function f(shape::Shape)
    seriestype --> :shape
    coords(shape)
end

@recipe function f(shapes::AVec{Shape})
    seriestype --> :shape
    coords(shapes)
end

@recipe function f(shapes::AMat{Shape})
    seriestype --> :shape
    for j in axes(shapes,2)
        @series coords(vec(shapes[:,j]))
    end
end

# Dicts: each entry is a data point (x,y)=(key,value)

@recipe f(d::AbstractDict) = collect(keys(d)), collect(values(d))

# function without range... use the current range of the x-axis

@recipe function f(f::FuncOrFuncs{F}) where F<:Function
    plt = plotattributes[:plot_object]
    xmin, xmax = if haskey(plotattributes, :xlims)
        plotattributes[:xlims]
    else
        try
            axis_limits(plt[1], :x)
        catch
            xinv = invscalefunc(get(plotattributes, :xscale, :identity))
            xm = PlotUtils.tryrange(f, xinv.([-5,-1,0,0.01]))
            xm, PlotUtils.tryrange(f, filter(x->x>xm, xinv.([5,1,0.99, 0, -0.01])))
        end
    end

    f, xmin, xmax
end

# --------------------------------------------------------------------
# 2 arguments
# --------------------------------------------------------------------

# if functions come first, just swap the order (not to be confused with parametric
# functions... as there would be more than one function passed in)

@recipe function f(f::FuncOrFuncs{F}, x) where F<:Function
    F2 = typeof(x)
    @assert !(F2 <: Function || (F2 <: AbstractArray && F2.parameters[1] <: Function))  # otherwise we'd hit infinite recursion here
    x, f
end

# --------------------------------------------------------------------
# 3 arguments
# --------------------------------------------------------------------

# surface-like... function
@recipe function f(x::AVec, y::AVec, zf::Function)
    x, y, Surface(zf, x, y)  # TODO: replace with SurfaceFunction when supported
end

# surface-like... matrix grid
@recipe function f(x::AVec, y::AVec, z::AMat)
    if !like_surface(get(plotattributes, :seriestype, :none))
        plotattributes[:seriestype] = :contour
    end
    x, y, Surface(z)
end

# images - grays
@recipe function f(x::AVec, y::AVec, mat::AMat{T}) where T<:Gray
    if is_seriestype_supported(:image)
        seriestype := :image
        yflip --> true
        SliceIt, x, y, Surface(mat)
    else
        seriestype := :heatmap
        yflip --> true
        cbar --> false
        fillcolor --> ColorGradient([:black, :white])
        SliceIt, x, y, Surface(convert(Matrix{Float64}, mat))
    end
end

# images - colors
@recipe function f(x::AVec, y::AVec, mat::AMat{T}) where T<:Colorant
    if is_seriestype_supported(:image)
        seriestype := :image
        yflip --> true
        SliceIt, x, y, Surface(mat)
    else
        seriestype := :heatmap
        yflip --> true
        cbar --> false
        z, plotattributes[:fillcolor] = replace_image_with_heatmap(mat)
        SliceIt, x, y, Surface(z)
    end
end

# --------------------------------------------------------------------
# Parametric functions
# --------------------------------------------------------------------

# special handling... xmin/xmax with parametric function(s)
@recipe function f(f::Function, xmin::Number, xmax::Number)
    xscale, yscale = [get(plotattributes, sym, :identity) for sym=(:xscale,:yscale)]
    _scaled_adapted_grid(f, xscale, yscale, xmin, xmax)
end
@recipe function f(fs::AbstractArray{F}, xmin::Number, xmax::Number) where F<:Function
    xscale, yscale = [get(plotattributes, sym, :identity) for sym=(:xscale,:yscale)]
    unzip(_scaled_adapted_grid.(fs, xscale, yscale, xmin, xmax))
end
@recipe f(fx::FuncOrFuncs{F}, fy::FuncOrFuncs{G}, u::AVec) where {F<:Function,G<:Function}  = mapFuncOrFuncs(fx, u), mapFuncOrFuncs(fy, u)
@recipe f(fx::FuncOrFuncs{F}, fy::FuncOrFuncs{G}, umin::Number, umax::Number, n = 200) where {F<:Function,G<:Function} = fx, fy, range(umin, stop = umax, length = n)

function _scaled_adapted_grid(f, xscale, yscale, xmin, xmax)
    (xf, xinv), (yf, yinv) =  ((scalefunc(s),invscalefunc(s)) for s in (xscale,yscale))
    xs, ys = adapted_grid(yf∘f∘xinv, xf.((xmin, xmax)))
    xinv.(xs), yinv.(ys)
end

# special handling... 3D parametric function(s)
@recipe function f(fx::FuncOrFuncs{F}, fy::FuncOrFuncs{G}, fz::FuncOrFuncs{H}, u::AVec) where {F<:Function,G<:Function,H<:Function}
    mapFuncOrFuncs(fx, u), mapFuncOrFuncs(fy, u), mapFuncOrFuncs(fz, u)
end
@recipe function f(fx::FuncOrFuncs{F}, fy::FuncOrFuncs{G}, fz::FuncOrFuncs{H}, umin::Number, umax::Number, numPoints = 200) where {F<:Function,G<:Function,H<:Function}
    fx, fy, fz, range(umin, stop = umax, length = numPoints)
end

# --------------------------------------------------------------------
# Lists of tuples and GeometryTypes.Points
# --------------------------------------------------------------------

@recipe f(v::AVec{<:Tuple})               = unzip(v)
@recipe f(v::AVec{<:GeometryTypes.Point}) = unzip(v)
@recipe f(tup::Tuple)             = [tup]
@recipe f(p::GeometryTypes.Point) = [p]

# Special case for 4-tuples in :ohlc series
@recipe f(xyuv::AVec{<:Tuple{R1,R2,R3,R4}}) where {R1,R2,R3,R4} = get(plotattributes,:seriestype,:path)==:ohlc ? OHLC[OHLC(t...) for t in xyuv] : unzip(xyuv)

# --------------------------------------------------------------------
# handle grouping
# --------------------------------------------------------------------

splittable_kw(key, val, lengthGroup) = false
splittable_kw(key, val::AbstractArray, lengthGroup) = !(key in (:group, :color_palette)) && length(axes(val,1)) == lengthGroup
splittable_kw(key, val::Tuple, lengthGroup) = all(splittable_kw.(key, val, lengthGroup))
splittable_kw(key, val::SeriesAnnotations, lengthGroup) = splittable_kw(key, val.strs, lengthGroup)

split_kw(key, val::AbstractArray, indices) = val[indices, fill(Colon(), ndims(val)-1)...]
split_kw(key, val::Tuple, indices) = Tuple(split_kw(key, v, indices) for v in val)
function split_kw(key, val::SeriesAnnotations, indices)
    split_strs = split_kw(key, val.strs, indices)
    return SeriesAnnotations(split_strs, val.font, val.baseshape, val.scalefactor)
end

function groupedvec2mat(x_ind, x, y::AbstractArray, groupby, def_val = y[1])
    y_mat = Array{promote_type(eltype(y), typeof(def_val))}(undef, length(keys(x_ind)), length(groupby.groupLabels))
    fill!(y_mat, def_val)
    for i in eachindex(groupby.groupLabels)
        xi = x[groupby.groupIds[i]]
        yi = y[groupby.groupIds[i]]
        y_mat[getindex.(Ref(x_ind), xi), i] = yi
    end
    return y_mat
end

groupedvec2mat(x_ind, x, y::Tuple, groupby) = Tuple(groupedvec2mat(x_ind, x, v, groupby) for v in y)

group_as_matrix(t) = false

# split the group into 1 series per group, and set the label and idxfilter for each
@recipe function f(groupby::GroupBy, args...)
    lengthGroup = maximum(union(groupby.groupIds...))
    if !(group_as_matrix(args[1]))
        for (i,glab) in enumerate(groupby.groupLabels)
            @series begin
                label     --> string(glab)
                idxfilter --> groupby.groupIds[i]
                for (key,val) in plotattributes
                    if splittable_kw(key, val, lengthGroup)
                        :($key) := split_kw(key, val, groupby.groupIds[i])
                    end
                end
                args
            end
        end
    else
        g = args[1]
        if length(g.args) == 1
            x = zeros(Int, lengthGroup)
            for indexes in groupby.groupIds
                x[indexes] = eachindex(indexes)
            end
            last_args = g.args
        else
            x = g.args[1]
            last_args = g.args[2:end]
        end
        x_u = unique(sort(x))
        x_ind = Dict(zip(x_u, eachindex(x_u)))
        for (key,val) in plotattributes
            if splittable_kw(key, val, lengthGroup)
                :($key) := groupedvec2mat(x_ind, x, val, groupby)
            end
        end
        label --> reshape(groupby.groupLabels, 1, :)
        typeof(g)((x_u, (groupedvec2mat(x_ind, x, arg, groupby, NaN) for arg in last_args)...))
    end
end


text: include("paddle.jl")

max_stars_repo_name: PangHua/InfiniteScience

max_issues_repo_name: PangHua/InfiniteScience

max_forks_repo_name: PangHua/InfiniteScience


text:
@testset "Vehicle tracking" begin
    n = 100
    H = [1 0
         0 1.0]
    Q = [1 0
        0 1.0]
    rho = 0.1
    model = VehicleTracking(rand(n, 2), rho, H, Q)

    # Not possible to fit
    @test !has_fit_methods(VehicleTracking)

    # Simply test if it runs
    initial_state = [0.0, 0, 0, 0]
    sim = simulate(model.system, initial_state, n)

    model = VehicleTracking(sim, 0.1, H, Q)
    kalman_filter(model)
    pos_pred = get_predictive_state(model)
    pos_filtered = get_filtered_state(model)
end


text:
# Co-occurence matrix
"""
    coo_matrix(::Type{T}, doc::Vector{AbstractString}, vocab::OrderedDict{AbstractString, Int}, window::Int, normalize::Bool)

Basic low-level function that calculates the co-occurence matrix of a document.
Returns a sparse co-occurence matrix sized `n × n` where `n = length(vocab)`
with elements of type `T`. The document `doc` is represented by a vector of its
terms (in order)`. The keywords `window` and `normalize` indicate the size of the
sliding word window in which co-occurrences are counted and whether to normalize
of not the counts by the distance between word positions.

# Examples
```
julia> using StringAnalysis
       doc = StringDocument("This is a text about an apple. There are many texts about apples.")
       docv = tokenize(text(doc))
       vocab = OrderedDict("This"=>1, "is"=>2, "apple."=>3)
       StringAnalysis.coo_matrix(Float16, docv, vocab, 5, true)
3×3 SparseArrays.SparseMatrixCSC{Float16,Int64} with 4 stored entries:
  [2, 1]  =  2.0
  [1, 2]  =  2.0
  [3, 2]  =  0.3999
  [2, 3]  =  0.3999
```
"""
function coo_matrix(::Type{T},
                    doc::Vector{<:AbstractString},
                    vocab::OrderedDict{<:AbstractString, Int},
                    window::Int,
                    normalize::Bool=true) where T<:AbstractFloat
    n = length(vocab)
    m = length(doc)
    coom = spzeros(T, n, n)
    # Count co-occurrences
    for (i, token) in enumerate(doc)
        @inbounds for j in max(1, i-window):min(m, i+window)
            wtoken = doc[j]
            nm = T(ifelse(normalize, abs(i-j), 1))
            row = get(vocab, token, nothing)
            col = get(vocab, wtoken, nothing)
            if i!=j && row != nothing && col != nothing
                coom[row, col] += one(T)/nm
                coom[col, row] = coom[row, col]
            end
        end
    end
    return coom
end

"""
Basic Co-occurrence Matrix (COOM) type.

# Fields
  * `coomm::SparseMatriCSC{T,Int}` the actual COOM; elements represent
co-occurrences of two terms within a given window
  * `terms::Vector{String}` a list of terms that represent the lexicon of
the document or corpus
  * `column_indices::OrderedDict{String, Int}` a map between the `terms` and the
columns of the co-occurrence matrix
"""
struct CooMatrix{T}
    coom::SparseMatrixCSC{T, Int}
    terms::Vector{String}
    column_indices::OrderedDict{String, Int}
end

"""
    CooMatrix{T}(crps::Corpus [,terms] [;window=5, normalize=true])

Auxiliary constructor(s) of the `CooMatrix` type. The type `T` has to be
a subtype of `AbstractFloat`. The constructor(s) requires a corpus `crps` and
a `terms` structure representing the lexicon of the corpus. The latter
can be a `Vector{String}`, an `AbstractDict` where the keys are the lexicon,
or can be omitted, in which case the `lexicon` field of the corpus is used.
"""
function CooMatrix{T}(crps::Corpus,
                      terms::Vector{String};
                      window::Int=5,
                      normalize::Bool=true) where T<:AbstractFloat
    column_indices = columnindices(terms)
    n = length(terms)
    coom = spzeros(T, n, n)
    for doc in crps
        docv = tokens(doc)
        coom .+= coo_matrix(T, docv, column_indices, window, normalize)
    end
    return CooMatrix{T}(coom, terms, column_indices)
end

CooMatrix(crps::Corpus, terms::Vector{String}; window::Int=5, normalize::Bool=true) =
    CooMatrix{DEFAULT_FLOAT_TYPE}(crps, terms, window=window, normalize=normalize)

CooMatrix{T}(crps::Corpus, lex::AbstractDict; window::Int=5, normalize::Bool=true
            ) where T<:AbstractFloat =
    CooMatrix{T}(crps, collect(keys(lex)), window=window, normalize=normalize)

CooMatrix(crps::Corpus, lex::AbstractDict; window::Int=5, normalize::Bool=true) =
    CooMatrix{DEFAULT_FLOAT_TYPE}(crps, lex, window=window, normalize=normalize)

CooMatrix{T}(crps::Corpus; window::Int=5, normalize::Bool=true) where T<:AbstractFloat = begin
    isempty(lexicon(crps)) && update_lexicon!(crps)
    CooMatrix{T}(crps, lexicon(crps), window=window, normalize=normalize)
end

CooMatrix(crps::Corpus; window::Int=5, normalize::Bool=true) = begin
    isempty(lexicon(crps)) && update_lexicon!(crps)
    CooMatrix{DEFAULT_FLOAT_TYPE}(crps, lexicon(crps), window=window, normalize=normalize)
end

# AbstractString, AbstractDocument methods
function CooMatrix{T}(doc,
                      terms::Vector{String};
                      window::Int=5,
                      normalize::Bool=true) where T<:AbstractFloat
    # Initializations
    column_indices = columnindices(terms)
    docv = tokens(doc)
    coom = coo_matrix(T, docv, column_indices, window, normalize)
    return CooMatrix{T}(coom, terms, column_indices)
end

CooMatrix(doc, terms::Vector{String}; window::Int=5, normalize::Bool=true
         ) where T<:AbstractFloat =
    CooMatrix{DEFAULT_FLOAT_TYPE}(doc, terms, window=window, normalize=normalize)

function CooMatrix{T}(doc; window::Int=5, normalize::Bool=true) where T<:AbstractFloat
    terms = unique(String.(tokens(doc)))
    CooMatrix{T}(doc, terms, window=window, normalize=normalize)
end

CooMatrix(doc; window::Int=5, normalize::Bool=true) where T<:AbstractFloat =
    CooMatrix{DEFAULT_FLOAT_TYPE}(doc, window=window, normalize=normalize)

"""
    coom(c::CooMatrix)

Access the co-occurrence matrix field `coom` of a `CooMatrix` `c`.
"""
coom(c::CooMatrix) = c.coom

"""
    coom(entity, eltype=DEFAULT_FLOAT_TYPE [;window=5, normalize=true])

Access the co-occurrence matrix of the `CooMatrix` associated
with the `entity`. The `CooMatrix{T}` will first have to
be created in order for the actual matrix to be accessed.
"""
coom(entity, eltype::Type{T}=DEFAULT_FLOAT_TYPE;
        window::Int=5, normalize::Bool=true) where T<:AbstractFloat =
    coom(CooMatrix{T}(entity, window=window, normalize=normalize))


text:
### ============================================================================
### A type related to Proximal Gradient and Lasso
### ============================================================================

mutable struct ProximalGradient
    """
        Dynamical, functional, keeps references to detailed implementation
        functions.
    """

    g::Function
    gradient::Function
    prox::Function
    h::Function

    β::Float64  # For beta convex function.
    tol::Float64
    maxItr::Int64
    solutionDim::Tuple

    function ProximalGradient(tol::Float64, maxItr::Int64)
        this = new()
        this.tol = tol
        this.maxItr = maxItr
        return this
    end

end

function OptimizeProximalGradient(
        this::ProximalGradient,
        warm_start::Union{Array, Nothing}=nothing
    )::Matrix
    """
        Implement FISTA, Accelerated Proximal Gradient, copied from my HW.
    """
    tol = this.tol
    max_itr = this.maxItr

    if warm_start === nothing
        warm_start = zeros(this.solutionDim)
    end

    x  = warm_start
    Δx = Inf
    y  = x
    t  = 1
    ∇f = this.gradient(y)
    δ  = 1/this.β # stepsize
    xNew = similar(x)
    yNew = similar(y)

    while Δx >= tol && max_itr >= 1
        xNew = x - δ*∇f
        this.prox(xNew, δ)
        tNew = (1 + sqrt(1 + 4t^2))/2
        yNew .= xNew
        yNew += ((t - 1)/tNew)*(xNew - x)

        ∇f .= this.gradient(yNew)
        Δx = norm(xNew - x, 1)/norm(x, 1)

        t = tNew
        x = xNew
        y = yNew

        max_itr -= 1
    end
    if max_itr == 0
        println()
        Warn("Maximal iteration $(this.maxItr) reached for Proximal gradient. ")
        @assert norm(x, Inf) != Inf "Solution is infinite, something blowed up."
        Warn("It will be bumped up by 200% more and then we try again.")
        this.maxItr = 2*this.maxItr
        return OptimizeProximalGradient(this, warm_start)
    end

    return x

end

### ----------------------------------------------------------------------------

### ============================================================================
### Lasso Proximal Operator
### ============================================================================

function L1LassoProximal!(y::Union{Vector, Matrix}, t::Float64, λ::Float64)
    """
        Textbook definition of the L1 Lasso proximal operator.
    """

    return map!(y, y) do yi
        if yi > λ*t
            return yi - λ*t
        elseif yi < -λ*t
            return yi + λ*t
        else
            return 0
        end
    end

end

### ----------------------------------------------------------------------------

### ============================================================================
### The 2 Norm loss function with L1 regularization
### ============================================================================

function BuildPG2NormLasso(
        A::Matrix,
        b::Matrix,
        λ::Float64,
        tol::Float64=1e-5,
        max_itr::Int64=1000
    )::ProximalGradient
    """
        Given matrix A, vector b, and λ the regularization perameter,
        this will build an instance of ProximalGradient.

        * Use Lambdas to capture variables references
        * A factory methods
        ! Pass copy of matrices will be better.
    """

    @assert size(b, 2) == 1 "Expect vector b in the"*
    " shape of (n, 1) but it turns out to be $(size(b))"
    @assert size(b, 1) == size(A, 1) "Expect the size of the matrix to match "*
    "the vector but it turns out to be A is in $(size(A)), and b is $(size(b))"
    ATA = A'*A
    ATb = A'*b
    β = 4*opnorm(ATA)  # convexity from spectral norm.
    t = 1/β
    f(x) = norm(A*x - b)^2
    dg(x) = 2(ATA*x - ATb)

    # ======= build ==================================
    proxVectorized(y, t) = L1LassoProximal!(y, t, λ)
    this = ProximalGradient(tol, max_itr)
    this.g = f
    this.gradient = dg
    this.prox = proxVectorized
    this.β = β
    this.solutionDim = (size(A, 2), size(b, 2))
    this.h = (x) -> norm(x, 1)

    return this
end

function ChangeProximalGradientLassoλ!(this::ProximalGradient, λ::Float64)
    """
        Change the proximal operator for a different λ
    """
    this.prox = (y, t) -> L1LassoProximal!(y, t, λ)
end

# ------------------------------------------------------------------------------

### ============================================================================
### Lasso huber loss function and l1 norm regularization
### ============================================================================

function BuildPGHuberLasso(
    A::Matrix,
    B::Array,
    δ::Float64=1,
    λ::Float64=0
)
    # TODO: implement this

end

### ============================================================================
### Lasso Cross Entropy loss function and l1 norm regularization
### ============================================================================

function BuildPGEntropyLasso()

end


text:
import RegisterMismatch  # we use qualified names to support simultaneous testing with RegisterMismatchCuda, which exports the same API
using ImageCore, ImageFiltering
using ImageCore.OffsetArrays
using CenterIndexedArrays, RegisterCore, RegisterMismatchCommon
using Test, Libdl

havecuda = isdefined(Main, :use_cuda) ? Main.use_cuda : !isempty(find_library(["libcudart", "cudart"], ["/usr/local/cuda"]))
if havecuda
    println("Running both CPU and CUDA versions")
    using CUDAdrv
    import RegisterMismatchCuda
    RMlist = (RegisterMismatch, RegisterMismatchCuda)
    devlist = devices() # dev->capability(dev)[1] >= 2, nmax=1)
else
    RMlist = (RegisterMismatch,)
end
accuracy = 1e-5 # new isapprox accuracy = 1e-6

# Iterators
@testset "Iterators" begin
    agrid = [(1,1)  (1,3)  (1,5)  (1,7);
             (4,1)  (4,3)  (4,5)  (4,7);
             (7,1)  (7,3)  (7,5)  (7,7);
             (10,1) (10,3) (10,5) (10,7)]
    @test RegisterMismatch.aperture_grid((10,7), (4,4)) == agrid
    for (i,pt) in enumerate(RegisterMismatch.each_point(agrid))
        @test pt == agrid[i]
    end
    agrida = reshape(reinterpret(Int, vec(agrid)), (2,4,4))
    for (i,pt) in enumerate(RegisterMismatch.each_point(agrida))
        @test pt == [agrid[i]...]
    end

    @test RegisterMismatch.aperture_range((15,), (7,)) == (12:18,)
end

# Bias correction
@testset "Bias correction" begin
    val = NumDenom(2.0,0.5)
    # 2d
    mma = fill(val, 3, 3)
    for i = 1:3
        mma[2,i] = NumDenom(rand(),rand())
        mma[i,2] = NumDenom(rand(),rand())
    end
    mm = CenterIndexedArray(mma)
    RegisterMismatch.correctbias!(mm)
    for I in eachindex(mm)
        mmI = mm[I]
        @test mmI.num ≈ val.num
        @test mmI.denom ≈ val.denom
    end
    # 3d
    mma = fill(val, 3, 3, 3)
    for i = 1:3, j = 1:3
        mma[2,i,j] = NumDenom(rand(),rand())
        mma[i,2,j] = NumDenom(rand(),rand())
    end
    mm = CenterIndexedArray(mma)
    RegisterMismatch.correctbias!(mm)
    for I in eachindex(mm)
        mmI = mm[I]
        @test mmI.num ≈ val.num
        @test mmI.denom ≈ val.denom
    end
end

@testset "Mismatch" begin
    mdutils = nothing
    devlist = nothing

    for imsz in ((7,10), (6,5))
        for maxshift in ((4,3), (3,2))
            Apad = parent(ImageFiltering.padarray(reshape(1:prod(imsz), imsz[1], imsz[2]), Fill(0, maxshift, maxshift)))
            Bpad = parent(ImageFiltering.padarray(rand(1:20, imsz[1], imsz[2]), Fill(0, maxshift, maxshift)))
            for RM in RMlist
                # intensity normalization
                mm = RM.mismatch(Apad, Bpad, maxshift, normalization=:intensity)
                num, denom = RegisterCore.separate(mm)
                mmref = CenterIndexedArray{Float64}(undef, (2 .* maxshift .+ 1)...)
                for j = -maxshift[2]:maxshift[2], i = -maxshift[1]:maxshift[1]
                    Bshift = circshift(Bpad,-[i,j])
                    df = Apad-Bshift
                    mmref[i,j] = sum(df.^2)
                end
                nrm = sum(Apad.^2)+sum(Bpad.^2)
                @test ≈(mmref.data, num.data, atol=accuracy*nrm)
                @test ≈(fill(nrm,size(denom)), denom.data, atol=accuracy*nrm)
                # pixel normalization
                mm = RM.mismatch(Apad, Bpad, maxshift, normalization=:pixels)
                _, denom = RegisterCore.separate(mm)
                n = Vector{Int}[size(Apad,i).-abs.(-maxshift[i]:maxshift[i]) for i = 1:2]
                @test ≈(denom.data, n[1].*n[2]', atol=accuracy*maximum(denom))
            end
        end
    end
end

# Test for denom overflow with mismatch0
@testset "Denom overflow" begin
    C16 = N0f16[0.6 0.1; 0.7 0.1]
    D16 = N0f16[0.7 0.1; 0.6 0.1]
    C = Float64.(C16)
    D = Float64.(D16)
    nd = RegisterMismatch.mismatch0(C16, D16; normalization=:intensity)
    @test nd.denom ≈ sum((C.^2).+(D.^2))
    @test nd.num ≈ sum((C.-D).^2)
end

# Compare direct, global, and apertured
@testset "Compare direct, global and apertured" begin
    C = rand(7,9)
    D = rand(7,9)
    mm = RegisterMismatch.mismatch(C, D, (3,3))
    nd = RegisterMismatch.mismatch0(C, D)
    @test mm[0,0].num ≈ nd.num
    @test mm[0,0].denom ≈ nd.denom
    mm = RegisterMismatch.mismatch(C, D, (3,3), normalization=:pixels)
    nd = RegisterMismatch.mismatch0(C, D, normalization=:pixels)
    @test mm[0,0].num ≈ nd.num
    @test mm[0,0].denom ≈ nd.denom

    mms = RegisterMismatch.mismatch_apertures(C, D, (2,2), (3,2), normalization=:intensity)
    nd0 = RegisterMismatch.mismatch0(C, D)
    nd1 = RegisterMismatch.mismatch0(mms)
    @test nd0.num ≈ nd1.num
    @test nd0.denom ≈ nd1.denom
end

# Now do it for block mismatch
# A key property we're testing here is that
#     sum(nums) == num
# where num is equivalent to what would be computed globally (using mismatch)
@testset "Block mismatch" begin
    for imsz in ((15,16), (14,17))
        for maxshift in ((4,3), (3,2))
            for gridsize in ((2,1), (2,3), (2,2), (1,3))
                Apad = parent(ImageFiltering.padarray(reshape(1:prod(imsz), imsz[1], imsz[2]), Fill(0, maxshift, maxshift)))
                Bpad = parent(ImageFiltering.padarray(rand(1:20, imsz[1], imsz[2]), Fill(0, maxshift, maxshift)))
                for RM in RMlist
                    # intensity normalization
                    mms = RM.mismatch_apertures(Float64, Apad, Bpad, gridsize, maxshift, normalization=:intensity, display=false)
                    nums, denoms = RegisterCore.separate(mms)
                    num = sum(nums)
                    denom = sum(denoms)
                    mm = CenterIndexedArray{Float64}(undef, (2 .* maxshift .+ 1)...)
                    for j = -maxshift[2]:maxshift[2], i = -maxshift[1]:maxshift[1]
                        Bshift = circshift(Bpad,-[i,j])
                        df = Apad-Bshift
                        mm[i,j] = sum(df.^2)
                    end
                    nrm = sum(Apad.^2)+sum(Bpad.^2)
                    @test ≈(mm.data, num.data, atol=accuracy*nrm)
                    @test ≈(fill(nrm,size(denom)), denom.data, atol=accuracy*nrm)
                    # pixel normalization
                    mms = RM.mismatch_apertures(Float64, Apad, Bpad, gridsize, maxshift, normalization=:pixels, display=false)
                    _, denoms = RegisterCore.separate(mms)
                    denom = sum(denoms)
                    n = Vector{Int}[size(Apad,i).-abs.(-maxshift[i]:maxshift[i]) for i = 1:2]
                    @test ≈(denom.data, n[1].*n[2]', atol=accuracy*maximum(denom))
                end
            end
        end
    end

    for RM in RMlist
        # Test 3d similarly
        Apad = parent(ImageFiltering.padarray(reshape(1:80*6, 10, 8, 6), Fill(0, (4,3,2))))
        Bpad = parent(ImageFiltering.padarray(rand(1:80*6, 10, 8, 6), Fill(0, (4,3,2))))
        mm = RM.mismatch(Apad, Bpad, (4,3,2))
        num, denom = RegisterCore.separate(mm)
        mmref = CenterIndexedArray{Float64}(undef, 9, 7, 5)
        for k=-2:2, j = -3:3, i = -4:4
            Bshift = circshift(Bpad,-[i,j,k])
            df = Apad-Bshift
            mmref[i,j,k] = sum(df.^2)
        end
        nrm = sum(Apad.^2)+sum(Bpad.^2)
        @test ≈(mmref.data, num.data, atol=accuracy*nrm)
        @test ≈(fill(nrm,size(denom)), denom.data, atol=accuracy*nrm)

        mms = RM.mismatch_apertures(Apad, Bpad, (2,3,2), (4,3,2), normalization=:intensity, display=false)
        nums, denoms = RegisterCore.separate(mms)
        num = sum(nums)
        denom = sum(denoms)
        @test ≈(mmref.data, num.data, atol=accuracy*nrm)
        @test ≈(fill(sum(Apad.^2)+sum(Bpad.^2),size(denom)), denom.data, atol=accuracy*nrm)
    end
end

### Test the utility functions that aren't implicitly tested by the above
# These only need to be tested once, since they come from RegisterMismatchCommon
@testset "Utility functions" begin
    RM = RegisterMismatch
    A = fill(7, 8, 4)
    Ahp = RM.highpass(Float32, A, (1,1))
    @test eltype(Ahp) == Float32
    @test maximum(abs.(Ahp)) < 100*eps(Float32)
    Ahp = RM.highpass(A, (1.2,Inf))
    @test A == Ahp
    A = Float64.(A)
    A[1,1] = NaN
    Ahp = RM.highpass(Float32, A, (1,1))
    @test sum(x->!isfinite(x), Ahp) == 1
    @test isnan(Ahp[1,1])
    Ahp[1,1] = 0
    @test maximum(abs, Ahp) < 3e-6

    num = float([1 2; 3 4])
    denom = [1 1e-6; 2 1]
    nd = RegisterCore._packnd!(similar(num, RegisterCore.NumDenom{Float64}), num, denom)
    RM.truncatenoise!(nd, 1e-5)
    numt, denomt = RegisterCore.separate(nd)
    @test numt == [1 0; 3 4]
    @test denomt == [1 0; 2 1]

    outer = reshape(1:120, 12, 10)
    A = outer[1:10,2:9]
    B = outer[3:12,1:8]
    maxshift = (3,11)
    @test register_translate(A, B, maxshift) == CartesianIndex((-2,1))
end

# Mismatched types
@testset "Mismatched types" begin
   A = rand(Float32, 5, 5)
    B = rand(5, 5)
    mm = RegisterMismatch.mismatch0(A, B)
    @test eltype(mm) == Float64
end

# Mismatched axes in the inputs (issue #7)
fullimg = collect(reshape(1:20, 4, 5))
fixed = fullimg[1:4, 2:5]
mxshift = (2, 2)
mm = RegisterMismatch.mismatch(fixed, fullimg, mxshift)
num, denom = RegisterCore.separate(mm)
mm = CenterIndexedArray{Float64}(undef, (2 .* mxshift .+ 1)...)
for j = -mxshift[2]:mxshift[2], i = -mxshift[1]:mxshift[1]
    mshift = OffsetArray(fullimg, (-i, -j))
    idx = intersect.(axes(fixed), axes(mshift))
    fsnip, msnip = fixed[idx...], mshift[idx...]
    mm[i,j] = sum((fsnip .- msnip).^2) / (sum(fsnip.^2) + sum(msnip.^2))
end
@test num ./ denom ≈ mm


text:
using Symbolics
using StaticArrays

"""
The potential and number of dimentions of a bead
"""
struct BeadDefinition{Ndims}
    chain_pe
    chain_force
    perbead_param_keys
    global_param_keys
    VNdims::Val{Ndims}
    chain_pe_expression
end

function BeadDefinition(chain_pe, perbead_param_keys::NTuple{Nbparams,Symbol}, global_param_keys::NTuple{Ngparams,Symbol} , VNdims::Val{Ndims}) where {Ndims, Nbparams, Ngparams}
    #symbolically diff potential energy to get force
    @variables p[1:Ndims]
    @variables b[1:Ndims]
    @variables n[1:Ndims]
    @variables bpars[1:Nbparams]
    @variables gpars[1:Ngparams]
    bsymparams= (; zip(perbead_param_keys,(bpars...,))...)
    gsymparams= (; zip(global_param_keys,(gpars...,))...)
    exprpe= chain_pe(p,b,n,bsymparams,gsymparams)
    g= Symbolics.gradient(-exprpe, [p;b;n]; simplify=true)
    @debug g
    gexpr= build_function(g, [p;b;n], bpars, gpars ,expression=Val{true})
    gfun= eval(gexpr[1])
    function chain_force(p,b,n, bpars, gpars)
        f= gfun([p;b;n], bpars, gpars)
        pf = ntuple((i -> @inbounds(f[i])), VNdims)
        bf = ntuple((i -> @inbounds(f[i+Ndims])), VNdims)
        nf = ntuple((i -> @inbounds(f[i+2Ndims])), VNdims)
        (SVector(pf),SVector(bf),SVector(nf))
    end
    BeadDefinition(chain_pe, chain_force, perbead_param_keys, global_param_keys, VNdims, exprpe)
end

struct Chain
    beaddef::BeadDefinition
    perbead_params
    global_params
    init_pos
    init_vel
    function Chain(beaddef::BeadDefinition, perbead_params, global_params, init_pos, init_vel)
        #error checking
        size(init_pos) == size(init_vel) || throw(DimensionMismatch("init_pos $(size(init_pos)) and init_vel $(size(init_vel)) don't match"))
        length(size(init_pos)) == 2 || throw(DimensionMismatch("init_pos $(length(size(init_pos))) must be 2"))
        Val(size(init_pos)[2]) == beaddef.VNdims || throw(DimensionMismatch("init_pos $(size(init_pos)) must match bead dimension $(beaddef.VNdims)"))
        Nbeads= size(init_pos)[1]
        Nbeads == length(perbead_params) || throw(DimensionMismatch("number of coords $(Nbeads) must match number of params $(length(perbead_params))"))
        #make sure forces are finite
        for beadi in 1:Nbeads
            p= init_pos[mod(beadi-2,Nbeads)+1,:]
            b= init_pos[beadi,:]
            n= init_pos[mod(beadi,Nbeads)+1,:]
            pf, bf, nf = beaddef.chain_force(p,b,n,perbead_params[beadi], global_params)
            e= beaddef.chain_pe(p,b,n,perbead_params[beadi], global_params)
            all(isfinite.([pf;bf;nf;e])) || @warn "some forces or energies are not finite on Bead number $(beadi)"
            all(abs.([pf;bf;nf;e]) .< 2.0^30) || @warn "some forces or energies are over 2^30 on Bead number $(beadi)"
        end
        new(beaddef::BeadDefinition, copy(perbead_params), global_params, copy(init_pos), copy(init_vel))
    end
end

"""
zero velocities by default
"""
function Chain(beaddef::BeadDefinition, perbead_params, global_params, init_pos)
    Chain(beaddef::BeadDefinition, perbead_params, global_params, init_pos, zero(init_pos))
end

"""
A force that acts on a specific list of beads
"""
struct SpecificForce{Ndims,Nways}
    pe
    force
    params
    interactions::AbstractArray{NTuple{Nways,NTuple{2,T}},1} where T<:Integer
    VNdims::Val{Ndims}
    pe_expression
end

"""
Create a specific force
pe

"""
function SpecificForce(pe, params,
    interactions::AbstractArray{NTuple{Nways,NTuple{2,T}},1} where T<:Integer,
    VNdims::Val{Ndims}) where {Ndims, Nways}
    #symbolically diff potential energy to get force
    a = Num(SymbolicUtils.Sym{Real}(:a))
    @variables r[1:Nways,1:Ndims]
    @variables pars[1:length(params)]
    @variables time
    symparams= (; zip(keys(params),(pars...,))...)
    exprpe= pe(ntuple(i->r[i,:],Nways), time ,symparams) + 0time#This 0time is to force exprpe into a Num type
    rssymb= vcat(ntuple(i->r[i,:],Nways)...)
    g= Symbolics.gradient(-exprpe, rssymb; simplify=true)
    @debug g
    gexpr= build_function(g, rssymb, time, pars, expression=Val{true})
    gfun= eval(gexpr[1])
    function force(rs, time, pars)
        f= gfun(vcat(rs...), time, pars)
        ntuple(j->SVector(ntuple(i->f[i+Ndims*(j-1)],Ndims)),Nways)
    end
    SpecificForce(pe,force,params,copy(interactions),VNdims,exprpe)
end

"""
A collection of Chains and SpecificForces
"""
struct ChainSystem
    starttime
    sp_pes
    sp_forces
    sp_interactions
    sp_params
    beaddef::BeadDefinition
    perbead_params
    global_params
    init_pos
    init_vel
    """
    chainbounds is a tuple of at least 1 Int,
        chainbounds[1] is 1
        chainbounds[end] is the total number of beads + 1
        the other element in chainbounds are the bead ids of the start of a new chain
    """
    chainbounds::Tuple{T,Vararg{T}} where T <: Integer
end

"""
Create an empty System
"""
function ChainSystem(starttime, beaddef::BeadDefinition, global_params)
    ChainSystem(starttime,(),(),(),(),
                    beaddef,
                    nothing,global_params,
                    nothing,nothing,(1,))
end

"""
return a new ChainSystem with added chain
"""
function append(s::ChainSystem, c::Chain)
    #error checking
    s.beaddef == c.beaddef || error("system and chain must have matching beaddef")
    s.global_params == c.global_params || error("system and chain must have matching global_params")
    Nbeads= size(c.init_pos)[1]
    if isnothing(s.init_pos)
        init_pos= c.init_pos
        init_vel= c.init_vel
        perbead_params= c.perbead_params
    else
        init_pos= [s.init_pos; c.init_pos]
        init_vel= [s.init_vel; c.init_vel]
        perbead_params= [s.perbead_params; c.perbead_params]
    end
    chainbounds= (s.chainbounds...,s.chainbounds[end]+Nbeads)
    news= ChainSystem(s.starttime,
                    s.sp_pes,s.sp_forces,s.sp_interactions,s.sp_params,
                    s.beaddef,
                    perbead_params,s.global_params,init_pos, init_vel,chainbounds)
    #check energy and forces for issues
    force_pe(news,news.init_pos,news.starttime)
    news
end

"""
Return a new ChainSystem with added force
"""
function append(s::ChainSystem, f::SpecificForce)
    #error checking
    s.beaddef.VNdims == f.VNdims || error("system and force must have matching beaddef")
    pes= (s.sp_pes..., f.pe)
    forces= (s.sp_forces..., f.force)
    #convert interactions to full bead id, instead of tuple.
    beadids= map(i->map(ids->tobeadid(ids...,s.chainbounds), i), f.interactions)
    interactions= (s.sp_interactions..., beadids)
    params= (s.sp_params..., f.params)
    news= ChainSystem(s.starttime,pes,forces,interactions,params,
                s.beaddef,
                s.perbead_params,s.global_params,s.init_pos, s.init_vel,s.chainbounds)
    #check energy and forces for issues
    force_pe(news,news.init_pos,news.starttime)
    news
end

"""
Return force, pe of the beads
"""
function force_pe(s::ChainSystem,pos,time)
    @assert size(s.init_pos) == size(pos)
    f= zero(pos)
    pe= 0.0
    for fi in 1:length(s.sp_forces)
        for interaction in s.sp_interactions[fi]
            rs= map(id->pos[id,:], interaction)
            fs= s.sp_forces[fi](rs,time,s.sp_params[fi])
            spe= s.sp_pes[fi](rs,time,s.sp_params[fi])
            all(isfinite.(vcat(fs..., spe))) || @warn "some forces or energies are not finite on specific force $(sf), interaction $(interaction)"
            all(abs.(vcat(fs..., spe)) .< 2.0^30) || @warn "some forces or energies are over 2^30 on specific force $(sf), interaction $(interaction))"
            for i in 1:length(interaction)
                f[interaction[i],:] .+= fs[i]
            end
            pe += spe
        end
    end
    Nbeads= size(pos)[1]
    #chain forces
    for beadid in 1:Nbeads
        previd = prev_beadid(beadid,s.chainbounds)
        nextid = next_beadid(beadid,s.chainbounds)
        p= pos[previd,:]
        b= pos[beadid,:]
        n= pos[nextid,:]
        pf, bf, nf = s.beaddef.chain_force(p,b,n,s.perbead_params[beadid], s.global_params)
        e= s.beaddef.chain_pe(p,b,n,s.perbead_params[beadid], s.global_params)
        all(isfinite.([pf;bf;nf;e])) || @warn "chain forces or energies are not finite on Bead number $(beadid)"
        all(abs.([pf;bf;nf;e]) .< 2.0^30) || @warn "chain forces or energies are over 2^30 on Bead number $(beadid)"
        pe +=e
        f[previd,:] .+= pf
        f[beadid,:] .+= bf
        f[nextid,:] .+= nf
    end
    f, pe
end

# """
# Return system stepper
# """
# function makeintegrator(kernel,s::ChainSystem)

# end


text:
using Plots
using LaTeXStrings
using Printf
using Polynomials
using LinearAlgebra
using OffsetArrays

using Oceananigans
using Oceananigans.Grids: min_Δx
using Oceananigans.Models: ShallowWaterModel

rate_of_convergence(::UpwindBiasedFirstOrder) = 1
rate_of_convergence(::CenteredSecondOrder)    = 2
rate_of_convergence(::UpwindBiasedThirdOrder) = 3
rate_of_convergence(::CenteredFourthOrder)    = 4
rate_of_convergence(::UpwindBiasedFifthOrder) = 5
rate_of_convergence(::WENO5)                  = 5

labels(::UpwindBiasedFirstOrder) = "Upwind1ˢᵗ"
labels(::CenteredSecondOrder)    = "Center2ⁿᵈ"
labels(::UpwindBiasedThirdOrder) = "Upwind3ʳᵈ"
labels(::CenteredFourthOrder)    = "Center4ᵗʰ"
labels(::UpwindBiasedFifthOrder) = "Upwind5ᵗʰ"
labels(::WENO5)                  = "WENO5ᵗʰ "

shapes(::UpwindBiasedFirstOrder) = :square
shapes(::CenteredSecondOrder)    = :diamond
shapes(::UpwindBiasedThirdOrder) = :dtriangle
shapes(::CenteredFourthOrder)    = :rect
shapes(::UpwindBiasedFifthOrder) = :star5
shapes(::WENO5)                  = :star6

colors(::UpwindBiasedFirstOrder) = :black
colors(::CenteredSecondOrder)    = :green
colors(::UpwindBiasedThirdOrder) = :red
colors(::CenteredFourthOrder)    = :cyan
colors(::UpwindBiasedFifthOrder) = :magenta
colors(::WENO5)                  = :purple

halos(::UpwindBiasedFirstOrder) = 1
halos(::CenteredSecondOrder)    = 1
halos(::UpwindBiasedThirdOrder) = 2
halos(::CenteredFourthOrder)    = 2
halos(::UpwindBiasedFifthOrder) = 3
halos(::WENO5)                  = 3

L  = 2
U  = 1
W  = 0.1
Ns = 2 .^ (6:10)
pnorm = 1

c(x, y, z, t, U, W) = exp( - (x - U * t)^2 / W^2 )
   h(x, y, z) = 1
  uh(x, y, z) = U * h(x, y, z)

schemes = (
 UpwindBiasedFirstOrder(),
 CenteredSecondOrder(),
 UpwindBiasedThirdOrder(),
 CenteredFourthOrder(),
 UpwindBiasedFifthOrder(),
 WENO5()
);

error = Dict()
ROC   = Dict()

for N in Ns, (adv, scheme) in enumerate(schemes)

    function Δx_str2(i, N)
        if i < N/4
            return 1
        elseif i > 3N/4
            return 1
        elseif i < N/2
            return 1.2 * (i - N/4) + 1
        else
            return 1.2 * (3N/4 - i) + 1
        end
     end

     xF = zeros(Float64, N+1)

     for i = 2:N+1
         xF[i] = xF[i-1] + Δx_str2(i-1, N)
     end

     xF ./= xF[end]
     xF = xF .*2 .- 1

    grid = RectilinearGrid(Float64; size=N,
                                x=xF,
                                halo=(halos(scheme)),
                                topology=(Periodic, Flat, Flat))

    Δt = 0.1 * min_Δx(grid)

    if adv == 7
        scheme = WENO5( grid = grid )
    end

    model = ShallowWaterModel(architecture = CPU(), grid = grid,
                                advection = scheme,
                                coriolis = nothing,
                                gravitational_acceleration = 0,
                                tracers = (:c))

    set!(model, h = h, uh = uh, c = (x, y, z) -> c(x, y, z, 0, U, W))

    simulation = Simulation(model, Δt=Δt, stop_iteration=1)

    run!(simulation)

    c₁  = c.(grid.xᶜᵃᵃ[:,1,1], 0, 0, Δt, U, W);

    error[(N, adv)] = norm(abs.(model.tracers.c[1:N, 1, 1] .- c₁[1:N]), pnorm)/N^(1/pnorm)

end

println("\nResults are for the L"*string(pnorm)*"-norm:\n")

for (adv, scheme) in enumerate(schemes)

    local best_fit = fit(log10.(Ns[2:end]),
                          log10.([error[(N, adv)] for N in Ns][2:end]), 1)

    ROC[adv] = best_fit[1]

    @printf("Method = % 24s, Rate of Convergence = %.2f, Expected = %d \n",
    scheme, -ROC[adv], rate_of_convergence(scheme))
end

function plot_solutions!(error, Ns, schemes, rate_of_convergence, shapes, colors, labels, pnorm, ROC)

    plt = plot()

    for (adv, scheme) in enumerate(schemes)

        plot!(
            plt,
            log2.(Ns),
            [error[(N, adv)] for N in Ns],
            seriestype = :scatter,
            shape = shapes(scheme),
            markersize = 8,
            markercolor = colors(scheme),
            xscale = :log10,
            yscale = :log10,
            xlabel = "log₂N",
            ylabel = "L"*string(pnorm)*"-norm: |cₛᵢₘ - c₁|",
            xticks = (log2.(Ns), string.(Int.(log2.(Ns)))),
            label =  string(labels(scheme))*" slope = "*@sprintf("%.2f", ROC[adv]),
            legend = :bottomleft,
            title = "Rates of Convergence"
        )

    end

    for (adv, scheme) in enumerate(schemes)

        roc = rate_of_convergence(scheme)

        best_line = [error[(Ns[1], adv)]] .* (Ns[1] ./ Ns) .^ roc

        plot!(
            plt,
            log2.(Ns),
            best_line,
            linestyle = :solid,
            lw = 3,
            linecolor = colors(scheme),
            label = "Expected slope = "*@sprintf("%d",-roc)
        )

    end

    return plt
end

plt = plot_solutions!(
    error,
    Ns,
    schemes,
    rate_of_convergence,
    shapes,
    colors,
    labels,
    pnorm,
    ROC)
savefig(plt, "convergence_rates")


text:
lines = readlines(ARGS[1])

function solution1(lines)
    horizontal = 0
    depth = 0
    for line in lines
        command, distance = split(line)
        int_distance = parse(Int, distance)
        if command == "up"
            depth -= int_distance
        elseif command == "down"
            depth += int_distance
        elseif command == "forward"
            horizontal += int_distance
        else
            throw(Exception(command, "unknown input"))
        end
    end
    return depth * horizontal
end

function solution2(lines)
    horizontal = 0
    depth = 0
    aim = 0
    for line in lines
        command, distance = split(line)
        int_distance = parse(Int, distance)
        if command == "up"
            aim -= int_distance
        elseif command == "down"
            aim += int_distance
        elseif command == "forward"
            horizontal += int_distance
            depth += aim * int_distance
        else
            throw(Exception(command, "unknown input"))
        end
    end
    return depth * horizontal
end

println(solution1(lines))
println(solution2(lines))


text:
tol = 1e-3
#
# Test Earth-like values against the online calculator.
# https://depts.washington.edu/naivpl/sites/default/files/hz.shtml
#
Teff   = 5780 # Sun's effective temperature.
L_star = 1.0

limits = habitable_zone(Teff, L_star)
@test abs(limits[1] - 0.750) < tol  # Recent Venus.
@test abs(limits[2] - 0.950) < tol  # Runaway Greenhouse.
@test abs(limits[3] - 1.676) < tol  # Maximum Greenhouse.
@test abs(limits[4] - 1.765) < 1e-2 # Early Mars. ---  FIXME (very close.)
@test abs(limits[5] - 0.917) < tol  # Runaway Greenhouse for 5.0 ME.
@test abs(limits[6] - 1.005) < tol  # Runaway Greenhouse for 0.1 ME.

#
# Test an M-dwarf against the online calculator.
# https://depts.washington.edu/naivpl/sites/default/files/hz.shtml
#
Teff   = 3700  # Temperature of a 0.5 M_sun star.
L_star = 0.04  # Luminosity of a 0.5 M_sun star.

limits = habitable_zone(Teff, L_star)
@test abs(limits[1] - 0.163) < tol # Recent Venus.
@test abs(limits[2] - 0.207) < tol # Runaway Greenhouse.
@test abs(limits[3] - 0.398) < tol # Maximum Greenhouse.
@test abs(limits[4] - 0.419) < tol # Early Mars.
@test abs(limits[5] - 0.199) < tol # Runaway Greenhouse for 5.0 ME.
@test abs(limits[6] - 0.219) < tol # Runaway Greenhouse for 0.1 ME.


text:
# Copyright 2021 Google LLC
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

import DataFrames
import CSV
using Plots
using LaTeXStrings

# Auxiliary functions used in the experiments.

"""
Saves the stats in the output of the primal-dual hybrid gradient solver into a
csv. The csv has two columns corresponding to primal and dual objectives
througout the execution of the solver; i.e., the rows correspond to iterations.
"""
function save_stats_to_csv(solver_output, csv_output_path::String)
  df = DataFrames.DataFrame(
    primal_objectives = solver_output.iteration_stats.primal_objectives,
    dual_objectives = solver_output.iteration_stats.dual_objectives,
    primal_solution_norms = solver_output.iteration_stats.primal_solution_norms,
    dual_solution_norms = solver_output.iteration_stats.dual_solution_norms,
    primal_delta_norms = solver_output.iteration_stats.primal_delta_norms,
    dual_delta_norms = solver_output.iteration_stats.dual_delta_norms,
    current_subspace_coordinate1 = solver_output.iteration_stats.current_subspace_coordinate1,
    current_subspace_coordinate2 = solver_output.iteration_stats.current_subspace_coordinate2,
    average_subspace_coordinate1 = solver_output.iteration_stats.average_subspace_coordinate1,
    average_subspace_coordinate2 = solver_output.iteration_stats.average_subspace_coordinate2,
  )
  CSV.write(csv_output_path, df)
end

"""
Generates plots from a stats csv generated by save_stats_to_csv.
"""
function generate_plots_from_stats_csv(csv_path::String, results_dir::String)
  df = CSV.read(csv_path)
  Plots.plot!(
    title = "Primal objective",
    xlabel = "Iteration",
    ylabel = "Primal objective",
  )
  plt_primal = Plots.plot(df.primal_objectives)
  Plots.savefig(plt_primal, joinpath(results_dir, "primal_objectives.png"))
  plt_dual = Plots.plot(df.dual_objectives)
  Plots.savefig(plt_dual, joinpath(results_dir, "dual_objectives.png"))

  residual = df.primal_delta_norms + df.dual_delta_norms
  plt_residual = Plots.plot(
    residual,
    title = "Residual",
    xlabel = "Iteration",
    ylabel = "Residual",
    yaxis = :log,
  )
  Plots.savefig(plt_residual, joinpath(results_dir, "residual.png"))

  if all(.!isnan.(df.current_subspace_coordinate1)) &&
     all(.!isnan.(df.current_subspace_coordinate2))
    plt_subspace = Plots.plot(
      df.average_subspace_coordinate1,
      df.average_subspace_coordinate2,
      title = "Iterates plotted on two dimensional subspace",
      xlabel = "Coordinate 1",
      ylabel = "Coordinate 2",
      label = "average iterate",
      markershape = :circle,
      linealpha = 0.0,
    )
    plt_subspace = Plots.plot!(
      df.current_subspace_coordinate1,
      df.current_subspace_coordinate2,
      label = "current iterate",
      markershape = :circle,
      linealpha = 0.0,
    )
    Plots.savefig(plt_subspace, joinpath(results_dir, "subspace.png"))
  end
end

function save_fraction_of_infeasible_instances_detected_to_csv(
  outputs::Vector{PdhgOutput},
  instance_names::Vector{String},
  csv_instance_path::String,
  csv_output_path::String,
)
  iterations = Set{Int64}([0])
  detection_iteration_average = Vector{Int64}()
  detection_iteration_iterate = Vector{Int64}()
  detection_iteration_difference = Vector{Int64}()
  detection_iteration_all = Vector{Int64}()
  for output in outputs
    detected_at = Vector{Int64}()
    if haskey(
      output.first_infeasibility_detection,
      FirstOrderLp.POINT_TYPE_AVERAGE_ITERATE,
    )
      iteration_average =
        output.first_infeasibility_detection[FirstOrderLp.POINT_TYPE_AVERAGE_ITERATE]
      push!(detection_iteration_average, iteration_average)
      push!(iterations, iteration_average)
      push!(detected_at, iteration_average)
    else
      push!(detection_iteration_average, typemax(Int64))
    end

    if haskey(
      output.first_infeasibility_detection,
      FirstOrderLp.POINT_TYPE_CURRENT_ITERATE,
    )
      iteration_iterate =
        output.first_infeasibility_detection[FirstOrderLp.POINT_TYPE_CURRENT_ITERATE]
      push!(detection_iteration_iterate, iteration_iterate)
      push!(iterations, iteration_iterate)
      push!(detected_at, iteration_iterate)
    else
      push!(detection_iteration_iterate, typemax(Int64))
    end

    if haskey(
      output.first_infeasibility_detection,
      FirstOrderLp.POINT_TYPE_ITERATE_DIFFERENCE,
    )
      iteration_difference =
        output.first_infeasibility_detection[FirstOrderLp.POINT_TYPE_ITERATE_DIFFERENCE]
      push!(detection_iteration_difference, iteration_difference)
      push!(iterations, iteration_difference)
      push!(detected_at, iteration_difference)
    else
      push!(detection_iteration_difference, typemax(Int64))
    end
    if length(detected_at) > 0
      push!(detection_iteration_all, minimum(detected_at))
    else
      push!(detection_iteration_all, typemax(Int64))
    end
  end
  iterations = sort(collect(iterations))
  fraction_identified_average = Vector{Float64}()
  fraction_identified_iterate = Vector{Float64}()
  fraction_identified_difference = Vector{Float64}()
  fraction_identified_all = Vector{Float64}()
  for iteration in iterations
    push!(
      fraction_identified_average,
      length(
        detection_iteration_average[detection_iteration_average.<=iteration],
      ) / length(outputs),
    )
    push!(
      fraction_identified_iterate,
      length(
        detection_iteration_iterate[detection_iteration_iterate.<=iteration],
      ) / length(outputs),
    )
    push!(
      fraction_identified_difference,
      length(
        detection_iteration_difference[detection_iteration_difference.<=iteration],
      ) / length(outputs),
    )
    push!(
      fraction_identified_all,
      length(detection_iteration_all[detection_iteration_all.<=iteration]) /
      length(outputs),
    )
  end
  df_instances = DataFrames.DataFrame(
    intance = instance_names,
    average = detection_iteration_average,
    iterate = detection_iteration_iterate,
    difference = detection_iteration_difference,
  )
  df = DataFrames.DataFrame(
    iterations = iterations,
    fraction_identified_average = fraction_identified_average,
    fraction_identified_iterate = fraction_identified_iterate,
    fraction_identified_difference = fraction_identified_difference,
    fraction_identified_all = fraction_identified_all,
  )
  CSV.write(csv_instance_path, df_instances)
  CSV.write(csv_output_path, df)
end

"""
Generates a plot of the fraction of infeasible instances detected vs the number of iterations.
"""
function generate_plot_fraction_infeasible_instances_detected(
  csv_path::String,
  output_path::String,
)
  df = CSV.read(csv_path)
  Plots.plot!(
    xlabel = "Iteration count",
    ylabel = "Fraction of infeasible instances detected",
  )
  gr()
  fntsm = font("serif-roman", pointsize = 12)
  fntlg = font("serif-roman", pointsize = 16)
  default(
    titlefont = fntlg,
    guidefont = fntlg,
    tickfont = fntsm,
    legendfont = fntsm,
  )
  plot(
    df.iterations,
    df.fraction_identified_difference,
    line = (3, :solid),
    label = "Difference",
    color = 4,
  )
  plot!(
    df.iterations,
    df.fraction_identified_average,
    line = (3, :solid),
    label = "Normalized average",
    color = 1,
  )
  plot!(
    df.iterations,
    df.fraction_identified_iterate,
    line = (3, :solid),
    label = "Normalized iterate",
    color = 2,
  )
  plot!(
    df.iterations,
    df.fraction_identified_all,
    line = (3, :solid),
    label = "All",
    color = 5,
  )
  plot!(xlims = (0.0, 1000000.0))
  ylims!((0.0, 1.0))
  savefig(output_path)
end

"""
Generates a convergence plot with respect to the infimal displacement vector.
The plot displays the distance to the differences, normalized iterates, and
normalized average. A vertical line displays the iteration where the support
last changed.
"""
function generate_infeasibility_convergence_plot(
  solver_output,
  plot_path::String,
)
  gr()
  fntsm = font("serif-roman", pointsize = 12)
  fntlg = font("serif-roman", pointsize = 16)
  default(
    titlefont = fntlg,
    guidefont = fntlg,
    tickfont = fntsm,
    legendfont = fntsm,
  )
  plot(
    solver_output.iteration_stats.distance_to_difference,
    line = (2, :solid),
    label = "Difference",
    color = 4,
  )
  plot!(
    solver_output.iteration_stats.distance_to_normalized_average,
    line = (2, :solid),
    label = "Normalized average",
    color = 1,
  )
  plot!(
    solver_output.iteration_stats.distance_to_normalized_current,
    line = (2, :solid),
    label = "Normalized iterate",
    color = 2,
  )
  plot!(
    [solver_output.last_active_set_change],
    seriestype = :vline,
    line = (2, :dash),
    label = "Last active set change",
    color = 6,
  )
  xaxis!("Iteration count")
  y_pow_10_range = -9:2
  ylims!((10.0^y_pow_10_range[1], 10.0^y_pow_10_range[end]))
  yticks!(
    10.0 .^ y_pow_10_range[1:2:end],
    [L"10^{%$i}" for i in y_pow_10_range[1:2:end]],
  )
  yaxis!("Distance to displacement vector", :log10)
  savefig(joinpath(plot_path))
end

"""
Generates a convergence plot showing scaled maximum error. A vertical line
displays the iteration where the support last changed.
"""
function generate_infeasibility_error_plot(solver_output, plot_path::String)
  gr()
  fntsm = font("serif-roman", pointsize = 12)
  fntlg = font("serif-roman", pointsize = 16)
  default(
    titlefont = fntlg,
    guidefont = fntlg,
    tickfont = fntsm,
    legendfont = fntsm,
  )
  plot(
    solver_output.iteration_stats.max_scaled_primal_certificate_error_of_difference,
    line = (2, :solid),
    label = "Difference",
    color = 4,
  )
  plot!(
    solver_output.iteration_stats.max_scaled_primal_certificate_error_of_normalized_average,
    line = (2, :solid),
    label = "Normalized average",
    color = 1,
  )
  plot!(
    solver_output.iteration_stats.max_scaled_primal_certificate_error_of_normalized_current,
    line = (2, :solid),
    label = "Normalized iterate",
    color = 2,
  )
  first_detection = minimum(values(solver_output.first_infeasibility_detection))
  x_limit = min(
    2 * first_detection,
    length(
      solver_output.iteration_stats.max_scaled_primal_certificate_error_of_difference,
    ),
  )
  if solver_output.last_active_set_change < x_limit
    plot!(
      [solver_output.last_active_set_change],
      seriestype = :vline,
      line = (2, :dash),
      label = "Last active set change",
      color = 6,
    )
  end
  xaxis!("Iteration count")
  xlims!(0, x_limit)
  y_pow_10_range = -9:1
  ylims!((10.0^y_pow_10_range[1], 10.0^y_pow_10_range[end]))
  yticks!(
    10.0 .^ y_pow_10_range[1:2:end],
    [L"10^{%$i}" for i in y_pow_10_range[1:2:end]],
  )
  yaxis!("Scaled certificate error", :log10)
  savefig(joinpath(plot_path))
end


text:
import Base.transpose, Base.adjoint, Base.conj

export AdjointLinearOperator,
  TransposeLinearOperator, ConjugateLinearOperator, adjoint, transpose, conj

# From julialang:stdlib/LinearAlgebra/src/adjtrans.jl
struct AdjointLinearOperator{T, S} <: AbstractLinearOperator{T}
  parent::S
  function AdjointLinearOperator{T, S}(A::S) where {T, S}
    new(A)
  end
end

struct TransposeLinearOperator{T, S} <: AbstractLinearOperator{T}
  parent::S
  function TransposeLinearOperator{T, S}(A::S) where {T, S}
    new(A)
  end
end

struct ConjugateLinearOperator{T, S} <: AbstractLinearOperator{T}
  parent::S
  function ConjugateLinearOperator{T, S}(A::S) where {T, S}
    new(A)
  end
end

AdjointLinearOperator(A) = AdjointLinearOperator{eltype(A), typeof(A)}(A)
TransposeLinearOperator(A) = TransposeLinearOperator{eltype(A), typeof(A)}(A)
ConjugateLinearOperator(A) = ConjugateLinearOperator{eltype(A), typeof(A)}(A)

adjoint(A::AbstractLinearOperator) = AdjointLinearOperator(A)
adjoint(A::AdjointLinearOperator) = A.parent
transpose(A::AbstractLinearOperator) = TransposeLinearOperator(A)
transpose(A::TransposeLinearOperator) = A.parent
conj(A::AbstractLinearOperator) = ConjugateLinearOperator(A)
conj(A::ConjugateLinearOperator) = A.parent

adjoint(A::ConjugateLinearOperator) = transpose(A.parent)
adjoint(A::TransposeLinearOperator) = conj(A.parent)
conj(A::AdjointLinearOperator) = transpose(A.parent)
conj(A::TransposeLinearOperator) = adjoint(A.parent)
transpose(A::AdjointLinearOperator) = conj(A.parent)
transpose(A::ConjugateLinearOperator) = adjoint(A.parent)

nprod(A::AdjointLinearOperator) = nctprod(A.parent)
ntprod(A::AdjointLinearOperator) = nprod(A.parent)   # transpose(A') = conj(A)
nctprod(A::AdjointLinearOperator) = nprod(A.parent)  # (A')' == A

nprod(A::TransposeLinearOperator) = ntprod(A.parent)
ntprod(A::TransposeLinearOperator) = nprod(A.parent)
nctprod(A::TransposeLinearOperator) = nprod(A.parent)  # (transpose(A))' = conj(A)

for f in [:nprod, :ntprod, :nctprod, :increase_nprod, :increase_ntprod, :increase_nctprod]
  @eval begin
    $f(A::ConjugateLinearOperator) = $f(A.parent)
  end
end

const AdjTrans = Union{AdjointLinearOperator, TransposeLinearOperator}

size(A::AdjTrans) = size(A.parent)[[2; 1]]
size(A::AdjTrans, d::Int) = size(A.parent, 3 - d)
size(A::ConjugateLinearOperator) = size(A.parent)
size(A::ConjugateLinearOperator, d::Int) = size(A.parent, d)

for f in [:hermitian, :ishermitian, :symmetric, :issymmetric]
  @eval begin
    $f(A::AdjTrans) = $f(A.parent)
    $f(A::ConjugateLinearOperator) = $f(A.parent)
  end
end

function show(io::IO, op::AdjointLinearOperator)
  println(io, "Adjoint of the following LinearOperator:")
  show(io, op.parent)
end

function show(io::IO, op::TransposeLinearOperator)
  println(io, "Transpose of the following LinearOperator:")
  show(io, op.parent)
end

function show(io::IO, op::ConjugateLinearOperator)
  println(io, "Conjugate of the following LinearOperator:")
  show(io, op.parent)
end

function mul!(res::AbstractVector, op::AdjointLinearOperator{T, S}, v::AbstractVector, α, β) where {T, S}
  p = op.parent
  (length(v) == size(p, 1) && length(res) == size(p, 2)) || throw(LinearOperatorException("shape mismatch"))
  if ishermitian(p)
    mul!(res, p, v, α, β)
    return nothing
  end
  if p.ctprod! !== nothing
    increase_nctprod(p)
    p.ctprod!(res, v, α, β)
    return nothing
  end
  tprod! = p.tprod!
  increment_tprod = true
  if p.tprod! === nothing
    if issymmetric(p)
      increment_tprod = false
      tprod! = p.prod!
    else
      throw(LinearOperatorException("unable to infer conjugate transpose operator"))
    end
  end
  if increment_tprod
    increase_ntprod(p)
  else
    increase_nprod(p)
  end
  tprod!(res, v, α, β)
  res .= conj.(res)
  return nothing
end

function mul!(res::AbstractVector, op::TransposeLinearOperator{T, S}, v::AbstractVector, α, β) where {T, S}
  p = op.parent
  (length(v) == size(p, 1) && length(res) == size(p, 2)) || throw(LinearOperatorException("shape mismatch"))
  if issymmetric(p)
    mul!(res, p, v, α, β)
    return nothing
  end
  if p.tprod! !== nothing
    increase_ntprod(p)
    p.tprod!(res, v, α, β)
    return nothing
  end
  increment_ctprod = true
  ctprod! = p.ctprod!
  if p.ctprod! === nothing
    if ishermitian(p)
      increment_ctprod = false
      ctprod! = p.prod!
    else
      throw(LinearOperatorException("unable to infer transpose operator"))
    end
  end
  if increment_ctprod
    increase_nctprod(p)
  else
    increase_nprod(p)
  end
  ctprod!(res, conj.(v), α, β)
  res .= conj.(res)
  return nothing
end

function mul!(res::AbstractVector, op::ConjugateLinearOperator{T, S}, v::AbstractVector, α, β) where {T, S}
  p = op.parent
  mul!(res, p, conj.(v), α, β)
  res .= conj.(res)
  return nothing
end

-(op::AdjointLinearOperator) = adjoint(-op.parent)
-(op::TransposeLinearOperator) = transpose(-op.parent)
-(op::ConjugateLinearOperator) = conj(-op.parent)

*(op::AdjointLinearOperator, x::Number) = adjoint(op.parent * conj(x))
*(op::TransposeLinearOperator, x::Number) = transpose(op.parent * x)
*(op::ConjugateLinearOperator, x::Number) = conj(op.parent * conj(x))

*(x::Number, op::AdjointLinearOperator) = adjoint(conj(x) * op.parent)
*(x::Number, op::TransposeLinearOperator) = transpose(x * op.parent)
*(x::Number, op::ConjugateLinearOperator) = conj(conj(x) * op.parent)


text:
"""
Calculate methanol critical properties
T --- temperature in Kelvin
p --- pressure in bar
ro --- density in mol/l
"""
function methanolCriticalF()
	sub=newCritical()
	sub.T=512.60
	sub.p=8.1035*10	#MPa to bar
	sub.ro=275.56/32.04	#mol/l #275.56 kg/m-3
	return sub
end

"""
Calculate mathanol properties at triple point
T --- temperature in Kelvin
p --- pressure in bar
vRo --- vapor density in mol/l
lRo --- liquid density in mol/l
url = https://www.engineeringtoolbox.com/methanol-methyl-alcohol-properties-CH3OH-d_2031.html
"""
function methanolTripleF()
	sub=newTriple()
	sub.T=175.50
	sub.p=1.86*10^(-8)	#MPa to bar
	sub.lRo=-1.0	#mol/l
	sub.vRo=-1.0	#mol/l
	return sub
end

"""
calculate mathanol properties at saturation line

T --- temperature in Kelvin
p --- pressure in bar
vRo --- vapor density in mol/l
lRo --- liquid density in mol/l
"""
function methanolSaturateF(T)
A=[-8.8570247, 2.4072447, -2.6452501, -1.5044111]
pA=[1.0, 1.5, 2.0, 3.5]

B=[-1.131507, 7.463057, -16.52860, -16.93269, -77.01263]
pB=[0.2, 0.8, 2.5, 4.5, 8.0]

C=[-0.0770237, 1.686259, -0.5426160, -0.3579365]
pC=[1.0/6.0, 0.3, 1.0 , 2.4]

sub=newSaturate()
crit=methanolCritical()
if (T<crit.T)
	v=1.0-T/crit.T
	sub.p=crit.p*exp( 1.0/(1.0-v) * (A[1]*v^pA[1] + A[2]*v^pA[2] + A[3]*v^pA[3] + A[4]*v^pA[4]) )
	sub.vRo=crit.ro*exp(B[1]*v^pB[1] + B[2]*v^pB[2] + B[3]*v^pB[3] + B[4]*v^pB[4] + B[5]*v^pB[5])
	sub.lRo=crit.ro*exp(C[1]*v^pC[1] + C[2]*v^pC[2] + C[3]*v^pC[3] + C[4]*v^pC[4])
	sub.T=T
else
	sub.p=-1
	sub.lRo=-1
	sub.vRo=-1
	sub.T=-1
end
return sub
end

"""
calculate methanol properties at single phase

T --- temperature, Kelvin
p --- pressure, bar
ro --- density, mol/l
z ---
s --- entropy,
s0 --- ideal part of entropy,
ds --- nonideal part derived k_B, dimensionless
g --- Gibbs energy
g0 --- ideal part of Gibbs energy
dg --- nonideal part of Gibbs energy
f ---  Helmholtz free energy
f0 --- ideal part
df --- nonideal part, dimensionless
u --- internal energy
du --- nonideal part
h --- entalpy
dh --- nonideal part
w --- speed of sound,
cp ---
cv ---
"""
function methanolSinglephaseF(T,ro)

sub=newSingle()
crit=methanolCritical()
delta=ro/crit.ro
tau=crit.T/T
R=8.3144598
mm=0.03204

nOne=6
nTwo=10
nThree=0
nFour=0

c=[0.000, 0.000, 0.000, 0.000, 0.000,
0.000, 1.000, 1.000, 1.000, 1.000,
1.000, 1.000, 1.000, 1.000, 1.000,
2.000, 2.000, 2.000, 2.000, 2.000,
2.000, 2.000, 2.000, 2.000, 3.000,
3.000, 4.000]

d=[1.0, 1.0, 2.0, 2.0, 3.0,
3.0, 1.0, 2.0, 3.0, 3.0,
4.0, 5.0, 6.0, 7.0, 8.0,
1.0, 1.0, 2.0, 2.0, 3.0,
4.0, 5.0, 6.0, 7.0, 7.0,
6.0, 5.0]

t=[0.5, 0,750, 0,125, 1.5, 0,375,
1.750, 0.000, 0.500, 1.000, 3.750,
2.000, 2.500, 0.000, 2.500, 4.500,
3.000, 4.000, 1.000, 3.000, 5.000,
6.000, 2.500, 5.000, 6.000, 12.000,
10.000, 15.000]

n=[0.12622395e02, -0.83224516e01,
-0.14647501e01, -0.12954522e1,
0.22417697, 0.22830533,
-0.46549039e01, -0.41099957e01,
0.70421007, 0.81617251e-01,
-0.37777607, 0.19627811,
0.45571723e-02, -0.11777859e-01,
-0.17277890e-04, 0.19096856e01,
-0.29551319e01, -0.28958480,
0.18169967e01, -0.96254996,
-0.11885503, -0.10730710,
0.15487654e-01, 0.66239025e-03,
0.15286750e-01, 0.31218155e-02,
0.14740469e-01]

al=[0.000, 0.000, 0.000, 0.000, 0.000,
0.000, 0.000, 0.000, 0.000, 0.000,
0.000, 0.000, 0.000, 0.000, 0.000,
0.000, 1.075, 0.463, 0.876, 1.108,
0.741, 4.032, 2.453, 2.300, 3.143]

bt=[0.000, 0.000, 0.000, 0.000, 0.000,
0.000, 0.000, 0.000, 0.000, 0.000,
0.000, 0.000, 0.000, 0.000, 0.000,
0.000, 1.207, 0.090, 0.581, 0.947,
2.356, 27.010, 4.542, 1.287, 3.090]

gm=[0.000, 0.000, 0.000, 0.000, 0.000,
0.000, 0.000, 0.000, 0.000, 0.000,
0.000, 0.000, 0.000, 0.000, 0.000,
0.000, 1.194, 1.986, 1.583, 0.756,
0.495, 1.002, 1.077, 1.493, 1.542]

ep=[0.000, 0.000, 0.000, 0.000, 0.000,
0.000, 0.000, 0.000, 0.000, 0.000,
0.000, 0.000, 0.000, 0.000, 0.000,
0.000, 0.779, 0.805, 1.869, 0.694,
1.312, 2.054, 0.441, 0.793, 0.313]

phidelta=0.0
phir=0.0
phitau=0.0
phidelta2=0.0
phitau2=0.0
dphideltadtau=0.0

i=0
for ni=1:nOne
	i=i+1
	phidelta=phidelta+n[i]*d[i]*delta^(d[i]-1.0)*tau^t[i]
	phir=phir+n[i]*delta^d[i]*tau^t[i]
	phitau=phitau+n[i]*t[i]*delta^d[i]*tau^(t[i]-1)
	phidelta2=phidelta2+n[i]*d[i]*(d[i]-1.0)*delta^(d[i]-2)*tau^t[i]
	phitau2=phitau2+n[i]*t[i]*(t[i]-1.0)*delta^d[i]*tau^(t[i]-2.0)
	dphideltadtau=dphideltadtau+n[i]*d[i]*t[i]*delta^(d[i]-1.0)*tau^(t[i]-1.0)
end
for ni=1:nTwo
	i+=1
	phidelta=phidelta+n[i]*exp(-delta^c[i])*(delta^(d[i]-1.0)*tau^t[i]*(d[i]-c[i]*delta^c[i]))
	phir=phir+n[i]*delta^(d[i])*tau^t[i]*exp(-delta^c[i])
	phitau=phitau+n[i]*t[i]*delta^d[i]*tau^(t[i]-1.0)*exp(-delta^c[i]);
	phidelta2=phidelta2+n[i]*exp(-delta^c[i])*(delta^(d[i]-2)*tau^(t[i])*((d[i]-c[i]*delta^c[i])*(d[i]-1.0-c[i]*delta^c[i])-c[i]^2*delta^c[i]))
	phitau2=phitau2+n[i]*t[i]*(t[i]-1.0)*delta^d[i]*tau^(t[i]-2)*exp(-delta^c[i])
	dphideltadtau=dphideltadtau+n[i]*t[i]*delta^(d[i]-1.0)*tau^(t[i]-1.0)*(d[i]-c[i]*delta^c[i])*exp(-delta^c[i])
end
for pn=1:nThree
	i+=1
	phidelta=phidelta+n[i]*delta^d[i]*tau^t[i]*exp(-al[i]*(delta-ep[i])^2-bt[i]*(tau-gm[i])^2)*(d[i]/delta-2.0*al[i]*(delta-ep[i]))
	phir=phir+n[i]*delta^d[i]*tau^t[i]*exp(-al[i]*(delta-ep[i])^2-bt[i]*(tau-gm[i])^2)
	phitau=phitau+n[i]*delta^d[i]*tau^t[i]*exp(-al[i]*(delta-ep[i])^2-bt[i]*(tau-gm[i])^2)*(t[i]/tau-2*bt[i]*(tau-gm[i]))
	phidelta2=phidelta2+n[i]*tau^t[i]*exp(-al[i]*(delta-ep[i])^2-bt[i]*(tau-gm[i])^2)*(-2*al[i]*delta^d[i]+4*al[i]^2*delta^d[i]*(delta-ep[i])^2-4*d[i]*al[i]*delta^(d[i]-1)*(delta-ep[i])+d[i]*(d[i]-1)*delta^(d[i]-2))
	phitau2=phitau2+n[i]*delta^d[i]*tau^t[i]*exp(-al[i]*(delta-ep[i])^2-bt[i]*(tau-gm[i])^2)*((t[i]/tau-2*bt[i]*(tau-gm[i]))^2-t[i]/tau^2-2*bt[i])
	dphideltadtau=dphideltadtau+n[i]*delta^d[i]*tau^t[i]*exp(-al[i]*(delta-ep[i])^2-bt[i]*(tau-gm[i])^2)*(d[i]/delta -2*al[i]*(delta-ep[i]))*(t[i]/tau-2*bt[i]*(tau-gm[i]))
end

sub.T=T
sub.ro=ro
sub.p=(1.0+phidelta*delta)*ro*R*T/100.0	#bar
sub.z=1.0+phidelta*delta	#
#ideal part
n0=[2.496674887, 2.900791185, -62.57135350, 10.99267739, 18.33682995, -16.36600476, -6.223234762, 2.803536282, 1.077809894, 0.969656970]
gamma0=[0.0, 0.0, 0.0, 4.119785, 3.264999, 3.769463, 2.931493, 8.225557, 10.31627, 0.5324892]

phi0 = n0[1] + n0[2]*log(tau) + n0[3]*tau + n0[4]*log(exp(gamma0[4]*tau)-1.0) + n0[5]*log(exp(gamma0[5]*tau)-1.0) + n0[6]*log(exp(gamma0[6]*tau)-1.0) + n0[7]*log(exp(gamma0[7]*tau)-1.0) + n0[8]*log(exp(gamma0[8]*tau)-1.0) + n0[9]*log(exp(gamma0[9]*tau)-1.0) + n0[10]*log(exp(gamma0*tau)-1.0)
phi0tau = n0[2]/tau + n0[3] + n0[4]*gamma0[4]*exp(gamma0[4]*tau)/(exp(gamma0[4]*tau)-1.0)+ n0[5]*gamma0[5]*exp(gamma0[5]*tau)/(exp(gamma0[5]*tau)-1.0) + n0[6]*gamma0[6]*exp(gamma0[6]*tau)/(exp(gamma0[6]*tau)-1.0) + n0[7]*gamma0[7]*exp(gamma0[7]*tau)/(exp(gamma0[7]*tau)-1.0) + n0[8]*gamma0[8]*exp(gamma0[8]*tau)/(exp(gamma0[8]*tau)-1.0) + n0[9]*gamma0[9]*exp(gamma0[9]*tau)/(exp(gamma0[9]*tau)-1.0) + n0[10]*gamma0[10]*exp(gamma0[10]*tau)/(exp(gamma0[10]*tau)-1.0)
phi0tau2 = -n0[2]/tau^2 + gamma0[4]*n0[4]*exp(gamma0[4]*tau)/(exp(gamma0[4]*tau)-1.0) - gamma0[4]^2*n0[4]*exp(2.0*gamma0[4]*tau)/(exp(gamma0[4]*tau)-1.0)^2 + gamma0[5]*n0[5]*exp(gamma0[5]*tau)/(exp(gamma0[5]*tau)-1.0) - gamma0[5]^2*n0[5]*exp(2.0*gamma0[5]*tau)/(exp(gamma0[5]*tau)-1.0)^2 + gamma0[6]*n0[6]*exp(gamma0[6]*tau)/(exp(gamma0[6]*tau)-1.0) - gamma0[6]^2*n0[6]*exp(2.0*gamma0[6]*tau)/(exp(gamma0[6]*tau)-1.0)^2 + gamma0[7]*n0[7]*exp(gamma0[7]*tau)/(exp(gamma0[7]*tau)-1.0) - gamma0[7]^2*n0[7]*exp(2.0*gamma0[7]*tau)/(exp(gamma0[7]*tau)-1.0)^2 + gamma0[8]*n0[8]*exp(gamma0[8]*tau)/(exp(gamma0[8]*tau)-1.0) - gamma0[8]^2*n0[8]*exp(2.0*gamma0[8]*tau)/(exp(gamma0[8]*tau)-1.0)^2 + gamma0[9]*n0[9]*exp(gamma0[9]*tau)/(exp(gamma0[9]*tau)-1.0) - gamma0[9]^2*n0[9]*exp(2.0*gamma0[9]*tau)/(exp(gamma0[9]*tau)-1.0)^2 + gamma0[10]*n0[10]*exp(gamma0[10]*tau)/(exp(gamma0[10]*tau)-1.0) - gamma0[10]^2*n0[10]*exp(2.0*gamma0[10]*tau)/(exp(gamma0[10]*tau)-1.0)^2

phi0delta = 0.0
phi0delta2 = 0.0
phi0taudelta=0.0
#entropy
#phi0=log(delta)+n0[1]+n0[2]*tau+n0[3]*log(tau)+n0[4]*log(1.0-exp(-gamma0[4]*tau))+n0[5]*log(1.0-exp(-gamma0[5]*tau))+n0[6]*log(1.0-exp(-gamma0[6]*tau))+n0[7]*log(1.0-exp(-gamma0[7]*tau))
#phi0tau=n0[2]+n0[3]/tau+n0[4]*gamma0[4]*(1.0/(1.0-exp(-gamma0[4]*tau))-1.0)+n0[5]*gamma0[5]*(1.0/(1.0-exp(-gamma0[5]*tau))-1.0)+n0[6]*gamma0[6]*(1.0/(1.0-exp(-gamma0[6]*tau))-1.0)+n0[7]*gamma0[7]*(1.0/(1.0-exp(-gamma0[7]*tau))-1.0)
#phi0delta=1/delta
#phi0delta2=-1/delta^2
#phi0tau2=-n0[3]/tau^2-n0[4]*gamma0[4]^2*exp(-gamma0[4]*tau)*(1.0-exp(-gamma0[4]*tau))^(-2)-n0[5]*gamma0[5]^2*exp(-gamma0[5]*tau)*(1.0-exp(-gamma0[5]*tau))^(-2)-n0[6]*gamma0[6]^2*exp(-gamma0[6]*tau)*(1.0-exp(-gamma0[6]*tau))^(-2)-n0[7]*gamma0[7]^2*exp(-gamma0[7]*tau)*(1.0-exp(-gamma0[7]*tau))^(-2)
#phi0taudelta=0.0

#Entropy
sub.s=(tau*(phi0tau+phitau)-phi0-phir)*R
sub.s0=(tau*(phi0tau)-phi0)*R
sub.ds=(tau*(phitau)-phir)	#devidet to R

#out.entrfrac=(tau*(phi0tau+phitau)-phi0-phir)/(tau*(phi0tau)-phi0)

dtaudtau=-tau/T
dphidt=dtaudtau*dphideltadtau

#sub.dpdt=(delta*cur_temp*dphidt+delta*phidelta+1)*R*T/100;
#phi0
#phir
#phidelta
#phidelta*delta
#(1+phi0+phir+phidelta*delta)
sub.g=(1+phi0+phir+phidelta*delta)*R*T	#R*cur_temp*
sub.g0=(1+phi0)*R*T
sub.dg=(phir+phidelta*delta)	#deriver RT
#
sub.f=(phi0+phir)*R*T
sub.f0=phi0*R*T
sub.df=phir

#Energy
sub.u=tau*(phi0tau+phitau)*R*T
sub.du=tau*(phitau)
#
sub.h=(1+tau*(phi0tau+phitau)+delta*phidelta)*R*T
sub.dh=1+tau*(phitau)+delta*phidelta;
#speed of sound
sub.w=sqrt(R/mm*T*(1+2*delta*phidelta+delta^2*phidelta2-(1+delta*phidelta-delta*tau*dphideltadtau)^2/(tau^2*(phi0tau2+phitau2))))	#0.018 перевод в килограмы
#
sub.cv=-tau^2*(phi0tau2+phitau2)*R
sub.cp=(-tau^2*(phi0tau2+phitau2)+(1+delta*phidelta-delta*tau*dphideltadtau)^2/(1+2*delta*phidelta+delta^2*phidelta2))*R
sub.dpdt=R/100.0*T*(1.0+2*delta*phidelta+delta^2*phidelta2)
#debug
# println("phi0 $(phi0)")
# println("phi0delta $(phi0delta)")
# println("phi0delta2 $(phi0delta2)")
# println("phi0tau $(phi0tau)")
# println("phi0tau2 $(phi0tau2)")
# println("phi0taudelta $(phi0taudelta)")
#
# println("phir $(phir)")
# println("phidelta $(phidelta)")
# println("phidelta2 $(phidelta2)")
# println("phitau $(phitau)")
# println("phitau2 $(phitau2)")
# println("dphideltadtau $(dphideltadtau)")
#enddebug
return sub
end


text:
# This file is a part of LegendDataTypes.jl, licensed under the MIT License (MIT).

# TODO: Replace by true custom types to avoid type piracy on NamedTuple
# (even though the named tuple types used here are quite specific).

const MaybeWithUnits{T} = Union{T,Quantity{<:T}}
const RealQuantity = MaybeWithUnits{<:Real}


text:
######### StanSample Bernoulli example  ###########

using StanSample

bernoulli_model = "
data {
  int<lower=1> N;
  int<lower=0,upper=1> y[N];
}
parameters {
  real<lower=0,upper=1> theta;
}
model {
  theta ~ beta(1,1);
  y ~ bernoulli(theta);
}
";

bernoulli_data = Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1])

# Keep tmpdir across multiple runs to prevent re-compilation
tmpdir = joinpath(@__DIR__, "tmp")

stanmodel = SampleModel(
  "bernoulli", bernoulli_model; tmpdir = tmpdir,
  method = StanSample.Sample(save_warmup=true,
    adapt = StanSample.Adapt(delta = 0.85)))

(sample_file, log_file) = stan_sample(stanmodel; data=bernoulli_data)

if !(sample_file == Nothing)
  # Use StanSamples to read a chain in NamedTupla format
  nt = read_samples(stanmodel.sm; chain = 3)

  # Convert to an MCMCChains.Chains object
  chns = read_samples(stanmodel)

  # Describe the MCMCChains using MCMCChains statistics
  # By default, just show the `parameters` section.
  # Use `chns.name_map` to see all sections.
  cdf = describe(chns)
  display(cdf)

  # Describe the `internals` section statistics
  icdf = describe(chns, sections=[:internals])
  display(icdf)

  # Show the same output in DataFrame format
  sdf = StanSample.read_summary(stanmodel)
  display(sdf)
end


text:
# set up ODE, get the derivatives
function Balances(t,y,dataDictionary)
  # correct for negatives
  idx_small = findall(y.<0)
  y[idx_small] .= 0.0

  # load data from data dictionary
  rnx_species = dataDictionary["rnx_species_array"]
  mRNA_species = dataDictionary["mRNA_species_array"]
  protein_species = dataDictionary["protein_species_array"]
  stoichiometry = dataDictionary["stoichiometric_matrix"]
  all_species_dict = dataDictionary["all_species_dict"]  # String --> number
  specific_growth_rate = dataDictionary["specificGrowthRate"]
  transcription_correction = dataDictionary["transcriptionSpecificCorrectionFactorDict"]
  translation_correction = dataDictionary["translationSpecificCorrectionFactor"]
  degradation_constant_mRNA = dataDictionary["degradationConstantmRNA"]
  degradation_constant_protein = dataDictionary["degradationConstantProtein"]

  # kinetics
  (rnx_rate, transcription_rate, translation_rate) = calculate_kinetics(y, dataDictionary)

  # calculate the derivatives
  dydt = zeros(length(all_species_dict))

  # signaling reaction network
  X2 = zeros(length(rnx_species)-3)  # species inside the cell
  for (id, st) in enumerate(rnx_species[3+1:end])  # get X2 value from y
    X2[id] = y[all_species_dict[st]]
  end
  dX1 = stoichiometry[1:3, :]*rnx_rate  # extracellular metabolites
  dX2 = stoichiometry[3+1:end, :]*rnx_rate - specific_growth_rate*X2  # INTRA
  dX = [dX1; dX2]
  for (id, st) in enumerate(rnx_species)  # return derivatives back to dydt
    dydt[all_species_dict[st]] = dX[id]
  end
  dydt[1] = specific_growth_rate*y[1]  # for Biomass

  # TXTL network
  for i = 1:length(mRNA_species)
  # seems easy to add additional terms to account for the signaling effect
    mRNA_id = all_species_dict[mRNA_species[i]]
    p_id = all_species_dict[protein_species[i]]
    dydt[mRNA_id] = transcription_rate[mRNA_species[i]] - (specific_growth_rate + transcription_correction[mRNA_species[i]]*degradation_constant_mRNA) * y[mRNA_id]
    dydt[p_id] = translation_rate[protein_species[i]] - (specific_growth_rate + translation_correction[protein_species[i]]*degradation_constant_protein) * y[p_id]
  end

  return dydt
end