text:
##
###  BinaryThresholdPredictor
##
"""
    BinaryThresholdPredictor(wrapped_model=ConstantClassifier(), threshold=0.5)

Wraps a `Probabilistic` model(`wrapped_model`) supporting binary
classification in a `BinaryThresholdPredictor` model. Training the
`BinaryThresholdPredictor` model results in training the underlying
`wrapped_model`.

`BinaryThresholdPredictor` model predictions are deterministic and
depend on the specified probability `threshold`(which lies in [0, 1)
). The positive class is predicted if the probability value exceeds
the threshold, where by convention the positive class is the second
class returned by `levels(y)`, where `y` is the target.

Note:
- The default `threshold`(0.5) predicts the modal class.
"""
mutable struct BinaryThresholdPredictor{M <: Probabilistic} <: Deterministic
    wrapped_model::M
    threshold::Float64
end

function clean!(model::BinaryThresholdPredictor)
    if !(AbstractVector{Multiclass{2}} <: target_scitype(model.wrapped_model) ||
        AbstractVector{OrderedFactor{2}} <: target_scitype(model.wrapped_model))
        throw(ArgumentError("`model` has unsupported target_scitype "*
              "`$(target_scitype(model.wrapped_model))`. "))
    end
    message = ""
    if model.threshold >= 1 || model.threshold < 0
        message = message*"`threshold` should be "*
        "in the range [0, 1). Resetting to 0.5. "
        model.threshold = 0.5
    end
    return message
end

function BinaryThresholdPredictor(;wrapped_model=ConstantClassifier(), threshold=0.5)
    model = BinaryThresholdPredictor(wrapped_model, Float64(threshold))
    message = clean!(model)
    isempty(message) || @warn message
    return model
end

function MMI.fit(model::BinaryThresholdPredictor, verbosity::Int, args...)
    scitype(args[2]) <: AbstractVector{Multiclass{2}} && begin
        first_class, second_class = levels(args[2])
        @warn "Taking positive class as `$(second_class)` and negative class as
        `$(first_class)`.
        Coerce target to `OrderedFactor{2}` to suppress this warning."
    end
    model_fitresult, model_cache, model_report = MMI.fit(
        model.wrapped_model, verbosity-1, args...
    )
    cache = (wrapped_model_cache = model_cache,)
    report = (wrapped_model_report = model_report,)
    fitresult = (model_fitresult, model.threshold)
    return fitresult, cache, report
end

function MMI.update(
    model::BinaryThresholdPredictor, verbosity::Int, old_fitresult, old_cache, args...
)
    model_fitresult, model_cache, model_report = MMI.update(
        model.wrapped_model, verbosity-1, old_fitresult[1], old_cache[1], args...
    )
    cache = (wrapped_model_cache = model_cache,)
    report = (wrapped_model_report = model_report,)
    fitresult = (model_fitresult, model.threshold)
    return fitresult, cache, report
end

function MMI.fitted_params(model::BinaryThresholdPredictor, fitresult)
    return (
        threshold= fitresult[2],
        wrapped_model_fitted_params = MMI.fitted_params(
            model.wrapped_model, fitresult[1]
        )
    )
end

function MMI.predict(model::BinaryThresholdPredictor, fitresult, X)
   yhat = MMI.predict(model.wrapped_model, fitresult[1], X)
   length(yhat.prob_given_ref) == 2 || begin
   # Due to resampling it's possible for Predicted `AbstractVector{<:UnivariateFinite}`
   # to contain one class. Hence the need for the following warning
   @warn "Predicted `AbstractVector{<:UnivariateFinite}`"*
       " contains only 1 class. Hence predictions will only contain this class "*
       "irrrespective of the set `threshold` "
   return mode.(yhat)
   end
   threshold = (1 - fitresult[2], fitresult[2])
   return _predict_threshold(yhat, threshold)
end

# `_div` and `_predict_threshold` methods are defined to help generalize to
# `MulticlassThresholdPredictor` (TODO)
_div(x,y) = ifelse( iszero(x) && x==y, Inf, x/y)

function _predict_threshold(yhat::UnivariateFinite, threshold)
    dict = yhat.prob_given_ref
    length(threshold) == length(dict) || throw(
        ArgumentError(
        "`length(threshold)` has to equal number of classes in specified "*
        "`UnivariateFinite` distribution."
        )
    )
    max_prob, max_class = findmax([_div(dict[ref], threshold[ref]) for ref in keys(dict)])
    return yhat.decoder(max_class)
end

function _predict_threshold(yhat::AbstractArray{<:UnivariateFinite}, threshold)
    return _predict_threshold.(yhat, (threshold,))
end

function _predict_threshold(yhat::UnivariateFiniteArray{S,V,R,P,N}, threshold) where{S,V,R,P,N}
    dict = yhat.prob_given_ref
    length(threshold) == length(dict) || throw(
        ArgumentError(
        "`length(threshold)` has to equal number of classes in specified "*
        "`UnivariateFiniteArray`."
        )
    )
    d = yhat.decoder(1)
    levs = levels(d)
    ord = isordered(d)
    # Array to house the predicted classes
    ret = CategoricalArray{V, N, R}(undef, size(yhat), levels=levs, ordered=ord)
    #ret = Array{CategoricalValue{V, R}, N}(undef, size(yhat))
    # `temp` vector allocted once to be used for calculations in each loop
    temp = Vector{Float64}(undef, length(dict))
    # allocate vectors of keys to enable use of @simd in innermost loop
    #kv = keys(dict) |> collect
    for i in eachindex(ret)
        #@simd for ind in eachindex(kv)
         @simd for ref in eachindex(temp)
           #@inbounds ref = kv[ind]
           #@inbounds temp[ind] = _div(dict[ref][i], threshold[ref])
           @inbounds temp[ref] = _div(dict[ref][i], threshold[ref])
        end
        max_prob, max_class = findmax(temp)
        @inbounds ret[i] = yhat.decoder(max_class)
    end
    return ret
end

## METADATA

# Note: input traits are inherited from the wrapped model

MMI.supports_weights(::Type{<:BinaryThresholdPredictor{M}}) where M =
    MMI.supports_weights(M)
MMI.load_path(::Type{<:BinaryThresholdPredictor}) =
    "MLJModels.BinaryThresholdPredictor"
MMI.package_name(::Type{<:BinaryThresholdPredictor}) = "MLJModels"
MMI.package_uuid(::Type{<:BinaryThresholdPredictor}) = ""
MMI.is_wrapper(::Type{<:BinaryThresholdPredictor}) = true
MMI.package_url(::Type{<:BinaryThresholdPredictor}) =
    "https://github.com/alan-turing-institute/MLJModels.jl"
MMI.is_pure_julia(::Type{<:BinaryThresholdPredictor{M}}) where M =
    MMI.is_pure_julia(M)
MMI.input_scitype(::Type{<:BinaryThresholdPredictor{M}}) where M =
    MMI.input_scitype(M)
MMI.target_scitype(::Type{<:BinaryThresholdPredictor}) =
    AbstractVector{<:Binary}


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ConnectionSharedKey <: SwaggerModel
    id::Any # spec type: Union{ Nothing, String } # spec name: id

    function ConnectionSharedKey(;id=nothing)
        o = new()
        validate_property(ConnectionSharedKey, Symbol("id"), id)
        setfield!(o, Symbol("id"), id)
        o
    end
end # type ConnectionSharedKey

const _property_map_ConnectionSharedKey = Dict{Symbol,Symbol}(Symbol("id")=>Symbol("id"))
const _property_types_ConnectionSharedKey = Dict{Symbol,String}(Symbol("id")=>"String")
Base.propertynames(::Type{ ConnectionSharedKey }) = collect(keys(_property_map_ConnectionSharedKey))
Swagger.property_type(::Type{ ConnectionSharedKey }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ConnectionSharedKey[name]))}
Swagger.field_name(::Type{ ConnectionSharedKey }, property_name::Symbol) =  _property_map_ConnectionSharedKey[property_name]

function check_required(o::ConnectionSharedKey)
    true
end

function validate_property(::Type{ ConnectionSharedKey }, name::Symbol, val)
end

ext: jl


text:
using IRTools: @dynamo, IR, xcall, arguments, insertafter!, recurse!

function hook(f, args...)
  print("Called ", f, "(")
  join(stdout, args, ", ")
  println(")")
end

@dynamo function logcalls(m...)
  ir = IR(m...)
  ir == nothing && return
  recurse!(ir)
  pushfirst!(ir, xcall(Main, :hook, arguments(ir)...))
  return ir
end

# @code_ir logcalls 2+3.0
logcalls(+, 2, 3.0)


text:
"""
    TestFunctions.jl

A collection of test functions for evaluating optimization algorithms.  Excerpted from Appendix B of:
Kochenderfer and Wheeler, "Algorithms for Optimization", MIT Press, 2019 (forthcoming)
"""
module TestFunctions

export
    ackley,
    booth,
    branin,
    flower,
    michalewicz,
    rosenbrock,
    wheeler,
    circle

"""
    ackley(x, a=20, b=0.2, c=2π)

Ackley's function.  N-D.
"""
function ackley(x, a=20, b=0.2, c=2π)
    d = length(x)
    return -a*exp(-b*sqrt(sum(x.^2)/d)) -
        exp(sum(cos.(c*xi) for xi in x)/d) + a + e
end

"""
    booth(x)

Booth's function. 2D.
"""
booth(x) = (x[1]+2x[2]-7)^2 + (2x[1]+x[2]-5)^2

"""
    branin(x; a=1, b=5.1/(4π^2), c=5/π, r=6, s=10, t=1/(8π))

Branin Function. 2D.
"""
function branin(x; a=1, b=5.1/(4π^2), c=5/π, r=6, s=10, t=1/(8π))
    return a*(x[2]-b*x[1]^2+c*x[1]-r)^2 + s*(1-t)*cos(x[1]) + s
end

"""
    flower(x; a=1, b=1, c=4)

Flower function. 2D.
"""
function flower(x; a=1, b=1, c=4)
    return a*norm(x) + b*sin(c*atan2(x[2],x[1]))
end

"""
    michalewicz(x; m=10)

Michaelewicz function.  N-D.
"""
function michalewicz(x; m=10)
    return -sum(sin(v)*sin(i*v^2/π)^(2m) for
                (i,v) in enumerate(x))
end

"""
    rosenbrock(x; a=1, b=5)

Rosenbrock's banana function. 2-D.
"""
rosenbrock(x; a=1, b=5) = (a-x[1])^2 + b*(x[2] - x[1]^2)^2

"""
    wheeler(x, a=1.5)

Wheeler's ridge. 2-D.
"""
wheeler(x, a=1.5) = -exp(-(x[1]*x[2] - a)^2 - (x[2]-a)^2)

"""
    circle(x)

Circle function. Multiobjective.  2-D.
"""
function circle(x)
    θ = x[1]
    r = 0.5 + 0.5*(2x[2]/(1+x[2]^2))
    y1 = 1 - r*cos(θ)
    y2 = 1 - r*sin(θ)
    return [y1, y2]
end

end # module


text:
function solve_004()
    # Get all products of three digit numbers
    prods = zeros(Int64, 900, 900)
    for i=100:999, j=100:999
        prods[i-99, j-99] = i*j
    end
    palindromes = filter(x -> x==reverse_int(x), prods)
    return maximum(palindromes)
end


text:
import Base: length
EmptyMat3 = Array{Float64}(undef, 0,0,0)
EmptyMat2 = Array{Float64}(undef, 0,0)
emptyMat3(P) = Array{P}(undef, 0,0,0)
emptyMat2(P) = Array{P}(undef, 0,0)
mutable struct Trace
    iter::Int64
    λ::Float64
    dλ::Float64
    cost::Float64
    α::Float64
    grad_norm::Float64
    improvement::Float64
    reduce_ratio::Float64
    time_derivs::Float64
    time_forward::Float64
    time_backward::Float64
    divergence::Float64
    η::Float64
    Trace() = new(0,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.)
end

(t::MVHistory)(args...) = increment!(t, args...)

"""
    `GaussianPolicy{P}`

# Fileds:
```
T::Int          # number of time steps
n::Int          # State dimension
m::Int          # Number of control inputs
K::Array{P,3}   # Time-varying feedback gain ∈ R(n,m,T)
k::Array{P,2}   # Open loop control signal  ∈ R(m,T)
Σ::Array{P,3}   # Time-varying controller covariance  ∈ R(m,m,T)
Σi::Array{P,3}  # The inverses of Σ
```
"""
mutable struct GaussianPolicy{P}
    T::Int
    n::Int
    m::Int
    K::Array{P,3}
    k::Array{P,2}
    Σ::Array{P,3}
    Σi::Array{P,3}
end

eye(P,n) = Matrix{P}(I,n,n)
GaussianPolicy(P) = GaussianPolicy(0,0,0,emptyMat3(P),emptyMat2(P),emptyMat3(P),emptyMat3(P))
GaussianPolicy(P,T,n,m) = GaussianPolicy(T,n,m,zeros(P,m,n,T),zeros(P,m,T),cat([eye(P,m) for t=1:T]..., dims=3),cat([eye(P,m) for t=1:T]..., dims=3))
Base.isempty(gp::GaussianPolicy) = gp.T == gp.n == gp.m == 0
Base.length(gp::GaussianPolicy) = gp.T

include("klutils.jl")

"""
iLQG - solve the deterministic finite-horizon optimal control problem.

minimize sum_i cost(x[:,i],u[:,i]) + cost(x[:,end])
s.t.  x[:,i+1] = f(x[:,i],u[:,i])

Inputs
======
`f, costfun, df`

1) step:
`xnew = f(x,u,i)` is called during the forward pass.
Here the state x and control u are vectors: size(x)==(n,),
size(u)==(m,). The time index `i` is a scalar.

2) cost:
`cost = costfun(x,u)` is called in the forward pass to compute
the cost per time-step along the trajectory `x,u`.

3) derivatives:
`fx,fu,fxx,fxu,fuu,cx,cu,cxx,cxu,cuu = df(x,u)` computes the
derivatives along a trajectory. In this case size(x)==(n, N) where N
is the trajectory length. size(u)==(m, N). The time indexes are I=(1:N).
Dimensions match the variable names e.g. size(fxu)==(n, n, m, N)
If cost function or system is time invariant, the dimension of the corresponding
derivatives can be reduced by dropping the time dimension

`x0` - The initial state from which to solve the control problem.
Should be a column vector. If a pre-rolled trajectory is available
then size(x0)==(n, N) can be provided and cost set accordingly.

`u0` - The initial control sequence. A matrix of size(u0)==(m, N)
where m is the dimension of the control and N is the number of state
transitions.

Outputs
=======
`x` - the optimal state trajectory found by the algorithm.
size(x)==(n, N)

`u` - the optimal open-loop control sequence.
size(u)==(m, N)

`traj_new` - A new `GaussianPolicy` object containing feedforward control trajectory and feedback-gains, these gains multiply the
deviation of a simulated trajectory from the nominal trajectory x. See `?GaussianPolicy` for more help.

`Vx` - the gradient of the cost-to-go. size(Vx)==(n, N)

`Vxx` - the Hessian of the cost-to-go. size(Vxx)==(n, n N)

`cost` - the costs along the trajectory. size(cost)==(1, N)
the cost-to-go is V = fliplr(cumsum(fliplr(cost)))

`trace` - a trace of various convergence-related values. One row for each
iteration, the columns of trace are
`[iter λ α g_norm Δcost z sum(cost) dλ]`
see below for details.

# Keyword arguments
`lims`,           [],            control limits\n
`α`,              logspace(0,-3,11), backtracking coefficients\n
`tol_fun`,         1e-7,          reduction exit criterion\n
`tol_grad`,        1e-4,          gradient exit criterion\n
`max_iter`,        500,           maximum iterations\n
`λ`,         1,             initial value for λ\n
`dλ`,        1,             initial value for dλ\n
`λfactor`,   1.6,           λ scaling factor\n
`λmax`,      1e10,          λ maximum value\n
`λmin`,      1e-6,          below this value λ = 0\n
`regType`,        1,             regularization type 1: q_uu+λ*I 2: V_xx+λ*I\n
`reduce_ratio_min`,           0,             minimal accepted reduction ratio\n
`diff_fun`,         -,             user-defined diff for sub-space optimization\n
`plot`,           1,             0: no  k>0: every k iters k<0: every k iters, with derivs window\n
`verbosity`,      2,             0: no  1: final 2: iter 3: iter, detailed\n
`plot_fun`,         x->0,          user-defined graphics callback\n
`cost`,           [],            initial cost for pre-rolled trajectory

This code consists of a port and extension of a MATLAB library provided by the autors of
`   INPROCEEDINGS{author={Tassa, Y. and Mansard, N. and Todorov, E.},
booktitle={Robotics and Automation (ICRA), 2014 IEEE International Conference on},
title={Control-Limited Differential Dynamic Programming},
year={2014}, month={May}, doi={10.1109/ICRA.2014.6907001}}`
"""
function iLQG(f,costfun,df, x0, u0;
    lims             = [],
    α                = exp10.(range(0, stop=-3, length=11)),
    tol_fun          = 1e-7,
    tol_grad         = 1e-4,
    max_iter         = 500,
    λ                = 1.,
    dλ               = 1.,
    λfactor          = 1.6,
    λmax             = 1e10,
    λmin             = 1e-6,
    regType          = 1,
    reduce_ratio_min = 0,
    diff_fun         = -,
    plot             = 1,
    verbosity        = 2,
    plot_fun         = x->0,
    cost             = [],
    traj_prev        = 0
    )
    debug("Entering iLQG")
    local fx,fu,fxx,fxu,fuu,cx,cu,cxx,cxu,cuu,xnew,unew,costnew,g_norm,Vx,Vxx,dV,αi
    # --- initial sizes and controls
    n   = size(x0, 1)          # dimension of state vector
    m   = size(u0, 1)          # dimension of control vector
    N   = size(u0, 2)          # number of state transitions
    u   = u0                   # initial control sequence
    traj_new  = GaussianPolicy(Float64)
    # traj_prev = GaussianDist(Float64)

    # --- initialize trace data structure
    trace = MVHistory()
    trace(:λ, 0, λ)
    trace(:dλ, 0, dλ)

    # --- initial trajectory
    debug("Setting up initial trajectory")
    if size(x0,2) == 1 # only initial state provided
        diverge = true
        for outer αi ∈ α
            debug("# test different backtracing parameters α and break loop when first succeeds")
            x,un,cost, = forward_pass(traj_new,x0[:,1],αi*u,[],1,f,costfun, lims,diff_fun)
            debug("# simplistic divergence test")
            if all(abs.(x) .< 1e8)
                u = un
                diverge = false
                break
            end
        end
    elseif size(x0,2) == N
        debug("# pre-rolled initial forward pass, initial traj provided")
        x        = x0
        diverge  = false
        isempty(cost) && error("Initial trajectory supplied, initial cost must also be supplied")
    else
        error("pre-rolled initial trajectory must be of correct length (size(x0,2) == N)")
    end

    trace(:cost, 0, sum(cost))
    #     plot_fun(x) # user plotting

    if diverge
        if verbosity > 0
            @printf("\nEXIT: Initial control sequence caused divergence\n")
        end
        return
    end

    # constants, timers, counters
    flg_change         = true
    Δcost              = 0.
    expected_reduction = 0.
    print_head         = 10 # print headings every print_head lines
    last_head          = print_head
    t_start            = time()
    verbosity > 0 && @printf("\n---------- begin iLQG ----------\n")
    satisfied          = true

    iter = accepted_iter = 1
    while accepted_iter <= max_iter
        reduce_ratio     = 0.
        # ====== STEP 1: differentiate dynamics and cost along new trajectory
        if flg_change
            _t = @elapsed fx,fu,fxx,fxu,fuu,cx,cu,cxx,cxu,cuu = df(x, u)
            trace(:time_derivs, iter, _t)
            flg_change   = false
        end
        # Determine what kind of system we are dealing with
        linearsys = isempty(fxx) && isempty(fxu) && isempty(fuu); debug("linear system: $linearsys")

        # ====== STEP 2: backward pass, compute optimal control law and cost-to-go
        back_pass_done = false
        while !back_pass_done
            _t = @elapsed diverge, traj_new,Vx, Vxx,dV = if linearsys
                back_pass(cx,cu,cxx,cxu,cuu,fx,fu,λ, regType, lims,x,u)
            else
                back_pass(cx,cu,cxx,cxu,cuu,fx,fu,fxx,fxu,fuu,λ, regType, lims,x,u)
            end
            increment!(trace, :time_backward, iter, _t)
            iter == 1 && (traj_prev = traj_new) # TODO: set k μu to zero fir traj_prev

            if diverge > 0
                verbosity > 2 && @printf("Cholesky failed at timestep %d.\n",diverge)
                dλ,λ = max(dλ*λfactor, λfactor), max(λ*dλ, λmin)
                if λ >  λmax; break; end
                continue
            end
            back_pass_done = true
        end

        k, K = traj_new.k, traj_new.K
        #  check for termination due to small gradient
        g_norm = mean(maximum(abs.(k) ./ (abs.(u) .+ 1), dims=1))
        trace(:grad_norm, iter, g_norm)
        if g_norm <  tol_grad && λ < 1e-5 && satisfied
            verbosity > 0 && @printf("\nSUCCESS: gradient norm < tol_grad\n")
            break
        end

        # ====== STEP 3: line-search to find new control sequence, trajectory, cost
        fwd_pass_done  = false
        if back_pass_done
            debug("#  serial backtracking line-search")
            @elapsed(for outer αi = α
                xnew,unew,costnew = forward_pass(traj_new, x0[:,1] ,u, x,αi,f,costfun, lims, diff_fun)
                Δcost    = sum(cost) - sum(costnew)
                expected_reduction = -αi*(dV[1] + αi*dV[2])
                reduce_ratio = if expected_reduction > 0
                    Δcost/expected_reduction
                else
                    @warn("negative expected reduction: should not occur")
                    sign(Δcost)
                end
                if reduce_ratio > reduce_ratio_min
                    fwd_pass_done = true
                    break
                end
            end) |> x -> trace(:time_forward, iter, x)

        end

        # ====== STEP 4: accept step (or not), print status

        #  print headings
        if verbosity > 1 && last_head == print_head
            last_head = 0
            @printf("%-12s", "iteration     cost    reduction     expected    gradient    log10(λ)    η    divergence\n")
        end

        if fwd_pass_done && satisfied # TODO: I added satisfied here, verify if this is reasonable
            if verbosity > 1
                @printf("%-12d%-12.6g%-12.3g%-12.3g%-12.3g%-12.1f\n",
                iter, sum(cost), Δcost, expected_reduction, g_norm, log10(λ))
                last_head += 1
            end
            dλ = min(dλ / λfactor, 1/ λfactor)
            λ *= dλ
            #  accept changes
            x,u,cost  = copy(xnew),copy(unew),copy(costnew)
            traj_new.k = copy(u)
            flg_change = true
            plot_fun(x)
            if Δcost < tol_fun
                verbosity > 0 &&  @printf("\nSUCCESS: cost change < tol_fun\n")
                break
            end
            accepted_iter += 1
        else #  no cost improvement
            αi =  NaN
            dλ,λ  = max(dλ * λfactor,  λfactor), max(λ * dλ,  λmin)#  increase λ
            if verbosity > 1
                @printf("%-12d%-12s%-12.3g%-12.3g%-12.3g%-12.1f\n",
                iter,"NO STEP", Δcost, expected_reduction, g_norm, log10(λ))
                last_head = last_head+1
            end
            if λ > λmax #  terminate ?
                verbosity > 0 && @printf("\nEXIT: λ > λmax\n")
                break
            end
        end
        #  update trace
        trace(:λ, iter, λ)
        trace(:dλ, iter, dλ)
        trace(:α, iter, αi)
        trace(:improvement, iter, Δcost)
        trace(:cost, iter, sum(cost))
        trace(:reduce_ratio, iter, reduce_ratio)
        iter += 1
    end

    iter ==  max_iter &&  verbosity > 0 && @printf("\nEXIT: Maximum iterations reached.\n")
    iter == 1 && error("Failure: no iterations completed, something is wrong. Try enabling the debug flag in DifferentialDynamicProgramming.jl for verbose printing.")

    verbosity > 0 && print_timing(trace,iter,t_start,cost,g_norm,λ)

    return x, u, traj_new, Vx, Vxx, cost, trace
end

function print_timing(trace,iter,t_start,cost,g_norm,λ)
    diff_t  = get(trace, :time_derivs)[2]
    diff_t  = sum(diff_t[.!isnan.(diff_t)])
    back_t  = get(trace, :time_backward)[2]
    back_t  = sum(back_t[.!isnan.(back_t)])
    fwd_t   = get(trace, :time_forward)[2]
    fwd_t   = sum(fwd_t[.!isnan.(fwd_t)])
    total_t = time()-t_start
    info = 100/total_t*[diff_t, back_t, fwd_t, (total_t-diff_t-back_t-fwd_t)]
    try
        @printf("\n iterations:   %-3d\n
        final cost:   %-12.7g\n
        final grad:   %-12.7g\n
        final λ: %-12.7e\n
        time / iter:  %-5.0f ms\n
        total time:   %-5.2f seconds, of which\n
        derivs:     %-4.1f%%\n
        back pass:  %-4.1f%%\n
        fwd pass:   %-4.1f%%\n
        other:      %-4.1f%% (graphics etc.)\n =========== end iLQG ===========\n",iter,sum(cost),g_norm,λ,1e3*total_t/iter,total_t,info[1],info[2],info[3],info[4])
    catch
        @show g_norm
    end
end


text:
# -------- Neural Model ------------------------------------------------------ #

"""
Trainable model that uses a neural network to generate the value and policy
for a game state. Optionally, the network can also predict features by
applying a dedicated dense layer on the representation used for value and
policy prediction.
"""
struct NeuralModel{G, GPU} <: Model{G, GPU}

  trunk :: Layer{GPU}           # Takes input and returns layer before logits
  features :: Vector{Feature}   # Features that the network must predict

  vhead :: Layer{GPU}                  # value head
  phead :: Layer{GPU}                  # policy head
  fhead :: Union{Nothing, Layer{GPU}}  # feature head

  vconv                         # Converts value-logit to value
  pconv                         # Converts policy-logits to policy

end

"""
    NeuralModel(G, trunk [, features; vhead, phead, fhead, vconv, pconv])

Construct a model for gametype `G` based on the neural network `trunk`,
optionally with `features` enabled. The heads `vhead`, `phead`, and `fhead` are
optional neural network layers that produce "logits" for the value, policy and
feature prediction. The functions `vconv` and `pconv` are used to map the
"logits" to values and policies. The respective feature converters are contained
in `features`.
"""
function NeuralModel( :: Type{G}
                    , trunk :: Layer{GPU}
                    , features = Feature[]
                    ; vhead :: Union{Nothing, Layer{GPU}} = nothing
                    , phead :: Union{Nothing, Layer{GPU}} = nothing
                    , fhead :: Union{Nothing, Layer{GPU}} = nothing
                    , vconv = Knet.tanh
                    , pconv = Knet.softmax
                    ) where {G, GPU}

  @assert valid_insize(trunk, size(G)) "Trunk incompatible with $G"

  features = check_features(Feature[features...])

  pl = policy_length(G)
  fl = feature_length(features, G)
  os = outsize(trunk, size(G))

  # Check the given heads and create linear default heads if not specified
  vhead = prepare_head(vhead, os, 1, GPU)
  phead = prepare_head(phead, os, pl, GPU)
  fhead = fl > 0 ? prepare_head(fhead, os, fl, GPU) : nothing

  NeuralModel{G, GPU}(trunk, features, vhead, phead, fhead, vconv, pconv)

end

# Low level access to neural models
function (m :: NeuralModel{G})(data, use_features = false) where {G <: Game}

  # Get the trunk output to calculate policy, value, features
  out = m.trunk(data)

  bs = size(out)[end]                 # batchsize
  pl = policy_length(G)               # policy length
  fl = feature_length(m.features, G)  # feature length

  # Apply the converters for value and policy on suitable reshapes
  v = m.vconv.(reshape(m.vhead(out), bs))
  p = m.pconv(reshape(m.phead(out), pl, bs), dims=1)

  # Apply the feature head if features are to be calculated
  if use_features && !isnothing(m.fhead)

    f = reshape(m.fhead(out), fl, bs)

    fs = map(features(l), feature_indices(features(l), G)) do feat, sel
      feature_conv(feat, f[sel,:])
    end

    # Collect the converted features
    f = isempty(fs) ? similar(p, 0, length(v)) : vcat(fs...)

  else

    f = similar(p, 0, length(v))

  end

  (v, p, f)

end

# Higher level access
function (m :: NeuralModel{G, GPU})( games :: Vector{G}
                                   , use_features = false
                                   ) where {G <: Game, GPU}

  at = atype(GPU)
  data = convert(at, representation(games))

  m(data, use_features)

end

function (m :: NeuralModel{G, GPU})( game :: G
                                   , use_features = false
                                   ) where {G <: Game, GPU}

  v, p, f = m([game], use_features)

  v[1], reshape(p, :), reshape(f, :)

end

function swap(m :: NeuralModel{G, GPU}) where {G, GPU}

  NeuralModel{G, !GPU}( swap(m.trunk)
                      , m.features
                      , swap(m.vhead)
                      , swap(m.phead)
                      , isnothing(m.fhead) ? nothing : swap(m.fhead)
                      , m.vconv
                      , m.pconv )

end

function Base.copy(m :: NeuralModel{G, GPU}) where {G, GPU}

  NeuralModel{G, GPU}( copy(m.trunk)
                     , m.features
                     , copy(m.vhead)
                     , copy(m.phead)
                     , isnothing(m.fhead) ? nothing : copy(m.fhead)
                     , m.vconv
                     , m.pconv )

end

features(m :: NeuralModel) = m.features

training_model(m :: NeuralModel) = m

# -------- Linear Neural Model ----------------------------------------------- #

function Shallow(:: Type{G}; kwargs...) where {G <: Game}
  NeuralModel(G, Pointwise(); kwargs...)
end

# -------- Multilayer Perceptron --------------------------------------------- #

function MLP(:: Type{G}, hidden, f = Knet.relu; kwargs...) where {G <: Game}

  widths = [ prod(size(G)), hidden..., policy_length(G) + 1 ]
  layers = [ Dense(widths[j], widths[j+1], f) for j in 1:length(widths) - 2 ]
  push!(layers, Dense(widths[end-1], widths[end], identity))

  NeuralModel(G, Chain(layers...); kwargs...)

end

# -------- Shallow Convolutional Network ------------------------------------ #

function ShallowConv( :: Type{G}
                    , filters
                    , f = Knet.relu
                    ; kwargs...
                    ) where {G <: Game}

  NeuralModel(G, Conv(size(G)[3], filters, f); kwargs...)

end

max_stars_repo_name: roSievers/jtac


text:
using LinearAlgebra

function load_sepsis_model()
    return []
end

function get_sepsis_score(data, model)
    x_mean = [
        83.8996, 97.0520,  36.8055,  126.2240, 86.2907,
        66.2070, 18.7280,  33.7373,  -3.1923,  22.5352,
        0.4597,  7.3889,   39.5049,  96.8883,  103.4265,
        22.4952, 87.5214,  7.7210,   106.1982, 1.5961,
        0.6943,  131.5327, 2.0262,   2.0509,   3.5130,
        4.0541,  1.3423,   5.2734,   32.1134,  10.5383,
        38.9974, 10.5585,  286.5404, 198.6777]
    x_std = [
        17.6494, 3.0163,  0.6895,   24.2988, 16.6459,
        14.0771, 4.7035,  11.0158,  3.7845,  3.1567,
        6.2684,  0.0710,  9.1087,   3.3971,  430.3638,
        19.0690, 81.7152, 2.3992,   4.9761,  2.0648,
        1.9926,  45.4816, 1.6008,   0.3793,  1.3092,
        0.5844,  2.5511,  20.4142,  6.4362,  2.2302,
        29.8928, 7.0606,  137.3886, 96.8997]
    c_mean = [60.8711, 0.5435, 0.0615, 0.0727, -59.6769, 28.4551]
    c_std = [16.1887, 0.4981, 0.7968, 0.8029, 160.8846, 29.5367]

    x = data[end, 1:34]
    c = data[end, 35:40]
    x_norm = nan_to_zero((x - x_mean) ./ x_std)
    c_norm = nan_to_zero((c - c_mean) ./ c_std)

    beta = [
        0.1806,  0.0249, 0.2120,  -0.0495, 0.0084,
        -0.0980, 0.0774, -0.0350, -0.0948, 0.1169,
        0.7476,  0.0323, 0.0305,  -0.0251, 0.0330,
        0.1424,  0.0324, -0.1450, -0.0594, 0.0085,
        -0.0501, 0.0265, 0.0794,  -0.0107, 0.0225,
        0.0040,  0.0799, -0.0287, 0.0531,  -0.0728,
        0.0243,  0.1017, 0.0662,  -0.0074, 0.0281,
        0.0078,  0.0593, -0.2046, -0.0167, 0.1239]
    rho = 7.8521
    nu = 1.0389

    xstar = vcat(x_norm, c_norm)
    exp_bx = exp(dot(xstar, beta))
    l_exp_bx = (4 / rho) ^ nu * exp_bx

    score = 1 - exp(-l_exp_bx)
    label = score > 0.45

    return score, label
end

function nan_to_zero(x)
    return map(y -> isnan(y) ? zero(y) : y, x)
end


text:
using ProfileSVG
using Test

# For these tests to work you need `rsvg-convert` installed.
# On Ubuntu this is `sudo apt install librsvg2-bin`.

function profile_test(n)
    for i = 1:n
        A = randn(100,100,20)
        m = maximum(A)
        Am = mapslices(sum, A; dims=2)
        B = A[:,:,5]
        Bsort = mapslices(sort, B; dims=1)
        b = rand(100)
        C = B.*b
    end
end

@testset "ProfileSVG.jl" begin
    profile_test(1)   # to compile
    @profview profile_test(10)
    mktemp() do path, io
        ProfileSVG.save(io)
        flush(io)
        # Validate the file by converting to PNG
        str = read(`rsvg-convert $path`, String)
        @test codeunits(str)[1:8] == UInt8[0x89, 0x50, 0x4e, 0x47, 0x0d, 0x0a, 0x1a, 0x0a]
    end
end


text:
using Colors, FixedPointNumbers, JLD2
using Test
using ColorTypes: eltype_default, parametric3

@testset "Conversion" begin
    r8(x) = reinterpret(N0f8, x)

    # Promotions
    a, b = promote(RGB(1,0,0), Gray(0.8))
    @test isa(a, RGB{Float64}) && isa(b, RGB{Float64})
    a, b = promote(RGBA(1,0,0), Gray(0.8))
    @test isa(a, RGBA{Float64}) && isa(b, RGBA{Float64})
    a, b = promote(RGBA(1,0,0), GrayA(0.8))
    @test isa(a, RGBA{Float64}) && isa(b, RGBA{Float64})
    a, b = promote(RGB(1,0,0), GrayA(0.8))
    @test isa(a, RGBA{Float64}) && isa(b, RGBA{Float64})
    a, b = promote(RGB(1,0,0), AGray(0.8))
    @test isa(a, ARGB{Float64}) && isa(b, ARGB{Float64})

    # srgb_compand / invert_srgb_compand
    @test Colors.srgb_compand(0.5) ≈ 0.7353569830524494 atol=eps()
    @test Colors.invert_srgb_compand(0.7353569830524494) ≈ 0.5 atol=eps()
    # issue #351
    l_pow_x_y() = for i=1:1000; (i/1000)^(1/2.4) end
    l_exp_y_log_x() = for i=1:1000; exp(1/2.4*log(i/1000)) end
    l_pow_x_y(); t_pow_x_y = @elapsed l_pow_x_y()
    l_exp_y_log_x(); t_exp_y_log_x = @elapsed l_exp_y_log_x()
    if t_exp_y_log_x > t_pow_x_y
        @warn "Optimization in `[invert_]srgb_compand()` may have the opposite effect."
    end

    fractional_types = (RGB, BGR, RGB1, RGB4)  # types that support Fractional

    redF64 = RGB{Float64}(1,0,0)
    redF32 = RGB{Float32}(1,0,0)
    red24 = reinterpret(RGB24, 0x00ff0000)
    red32 = reinterpret(ARGB32, 0xffff0000)
    @testset "type check: RGB, RGB{$T}" for T in (Float64, Float32, N0f8)
        c = RGB(one(T), zero(T), zero(T))
        @test eltype(c) == T
        @test convert(RGB{Float64}, c) === redF64
        @test convert(RGB{T}, redF64) == c
    end
    @test RGB(1,0,0) == redF64
    @test RGB(UInt8(1), 0, 0) == redF64
    @test RGB(UInt8(1), UInt8(0), UInt8(0)) == redF64 # != colorant"#010000"
    @test convert(RGB, red24) == redF64

    @testset "type check: RGB24-->$Cto" for Cto in parametric3
        @test typeof(convert(Cto, red24)) == Cto{eltype_default(Cto)}
        @test typeof(convert(Cto{Float64}, red24)) == Cto{Float64}
    end

    @testset "type check: C{Float}-->C{Float}" for Cfrom in parametric3, Tfrom in (Float32, Float64)
        c = convert(Cfrom{Tfrom}, redF64)
        @test typeof(c) == Cfrom{Tfrom}
        @testset "$Cfrom{$Tfrom}-->$Cto" for Cto in parametric3
            c1 = convert(Cto, c)
            @test eltype(c1) == Tfrom
        end
        @testset "$Cfrom{$Tfrom}-->$Cto{$Tto}" for Cto in parametric3, Tto in (Float32, Float64)
            c2 = convert(Cto{Tto}, c)
            @test typeof(c2) == Cto{Tto}
        end
    end

    normed_types = (N0f8, N6f10, N4f12, N2f14, N0f16)
    # Test conversion from Normed types
    @testset "type check: C{Normed}-->C{Float}" for Cfrom in fractional_types, Tfrom in normed_types
        c = convert(Cfrom{Tfrom}, redF64)
        @test typeof(c) == Cfrom{Tfrom}
        @testset "$Cfrom{$Tfrom}-->$Cto" for Cto in parametric3
            eltype_default(Cto) <: FixedPoint && continue
            c1 = convert(Cto, c)
            @test eltype(c1) == eltype_default(Cto)
        end
        @testset "$Cfrom{$Tfrom}-->$Cto{$Tto}" for Cto in parametric3, Tto in (Float32, Float64)
            c2 = convert(Cto{Tto}, c)
            @test typeof(c2) == Cto{Tto}
        end
    end

    # Test conversion to Normed types
    @testset "type check: C{Float}-->C{Normed}" for Cfrom in parametric3, Tfrom in (Float32, Float64)
        c = convert(Cfrom{Tfrom}, redF64)
        @test typeof(c) == Cfrom{Tfrom}
        @testset "$Cfrom{$Tfrom}-->$Cto{$Tto}" for Cto in fractional_types, Tto in normed_types
            c2 = convert(Cto{Tto}, c)
            @test typeof(c2) == Cto{Tto}
        end
    end

    ac = RGBA(redF64)

    @test convert(RGB, ac) == RGB(1,0,0)
    @test convert(RGB{N0f8}, ac) == RGB{N0f8}(1,0,0)
    @test convert(RGBA{N0f8}, ac) == RGBA{N0f8}(1,0,0,1)
    @test convert(HSVA, ac) == HSVA{Float64}(convert(HSV, redF64), 1.0)
    @test convert(HSVA{Float32}, ac) == HSVA{Float32}(convert(HSV{Float32}, redF64), 1.0f0)
    @test convert(RGBA, redF64) == ac

    @test convert(ARGB32, ac) == reinterpret(ARGB32, 0xffff0000)
    @test convert(ARGB32, RGB(1,0,0)) == reinterpret(ARGB32, 0xffff0000)
    ac32 = convert(ARGB32, RGB(1,0,0), 0.5)
    @test ac32 == reinterpret(ARGB32, 0x80ff0000)
    @test convert(ARGB32, ac32) === ac32
    @test convert(ARGB32, BGRA(1,0,0,0.5)) == reinterpret(ARGB32, 0x80ff0000)
    @test reinterpret(UInt32, convert(ARGB32, ac)) == 0xffff0000
    @test convert(RGB24, RGB(r8(0xff),r8(0x00),r8(0x00))) == reinterpret(RGB24, 0x00ff0000)
    @test reinterpret(UInt32, convert(RGB24, RGB(r8(0xff),r8(0x00),r8(0x00)))) == 0x00ff0000
    redhsv = convert(HSV, redF64)
    @test convert(RGB24, red24) === red24
    @test convert(RGB24, redhsv) == reinterpret(RGB24, 0x00ff0000)
    @test_throws ArgumentError convert(RGB24,  RGB(0, 1.1, 0))
    @test_throws ArgumentError convert(ARGB32, RGBA(0, 1.1, 0, 0.8))
    @test_throws ArgumentError convert(ARGB32, RGBA(0, 0.8, 0, 1.1))

    @test convert(RGB{N0f8}, red24) == RGB{N0f8}(1,0,0)
    @test convert(RGBA{N0f8}, red32) == RGBA{N0f8}(1,0,0,1)
    @test convert(HSVA{Float64}, red32) == HSVA{Float64}(360, 1, 1, 1)

    @test_throws MethodError AlphaColor(RGB(1,0,0), r8(0xff))

    # whitepoint conversions
    @test isa(convert(XYZ, convert(Lab, redF64), Colors.WP_DEFAULT), XYZ{Float64})
    @test isa(convert(XYZ{Float32}, convert(Lab, redF64), Colors.WP_DEFAULT), XYZ{Float32})
    @test isa(convert(XYZ, convert(Luv, redF64), Colors.WP_DEFAULT), XYZ{Float64})
    @test isa(convert(XYZ{Float32}, convert(Luv, redF64), Colors.WP_DEFAULT), XYZ{Float32})
    @test isa(convert(Lab, convert(XYZ, redF64), Colors.WP_DEFAULT), Lab{Float64})
    @test isa(convert(Lab{Float32}, convert(XYZ, redF64), Colors.WP_DEFAULT), Lab{Float32})
    @test isa(convert(Luv, convert(XYZ, redF64), Colors.WP_DEFAULT), Luv{Float64})
    @test isa(convert(Luv{Float32}, convert(XYZ, redF64), Colors.WP_DEFAULT), Luv{Float32})

    #59
    @test Colors.xyz_to_uv(XYZ{Float64}(0.0, 0.0, 0.0)) === (0.0, 0.0)
    @test Colors.xyz_to_uv(XYZ{Float64}(0.0, 1.0, 0.0)) === (0.0, 0.6)
    @test Colors.xyz_to_uv(XYZ{Float64}(0.0, 1.0, 1.0)) === (0.0, 0.5)
    @test Colors.xyz_to_uv(XYZ{Float64}(1.0, 0.0, 1.0)) === (1.0, 0.0)
    @test Colors.xyz_to_uv(XYZ{Float64}(1.0, 0.0, 0.0)) === (4.0, 0.0)

    # ColorTypes.jl issue #40
    @test convert(HSL, RGB{N0f8}(0.678, 0.847, 0.902)) ≈ HSL{Float32}(194.73685f0,0.5327105f0,0.7901961f0) atol=100eps(Float32)

    # YIQ
    @test convert(YIQ, RGB(1,0,0)) == YIQ{Float32}(0.299, 0.595716, 0.211456)
    @test convert(YIQ, RGB(0,1,0)) == YIQ{Float32}(0.587, -0.274453, -0.522591)
    @test convert(YIQ, RGB(0,0,1)) == YIQ{Float32}(0.114, -0.321263, 0.311135)
    @test convert(RGB, YIQ(1.0,0.0,0.0)) == RGB(1,1,1)
    v = 0.5957
    @test convert(RGB, YIQ(0.0,1.0,0.0)) == RGB(0.9563*v,0,0)
    v = -0.5226
    @test convert(RGB, YIQ(0.0,0.0,-1.0)) == RGB(0,-0.6474*v,0)

    # Gray
    @test convert(Gray{N0f8}, Gray{N0f8}(0.1)) == Gray{N0f8}(0.1)
    @test convert(Gray{N0f8}, Gray(0.1))     == Gray{N0f8}(0.1)
    @test convert(Gray{N0f8}, Gray24(0.1))   == Gray{N0f8}(0.1)
    @test convert(Gray24, Gray{N0f8}(0.1))   == Gray24(0.1)

    @test convert(RGB{N0f8}, Gray{N0f8}(0.1)) == RGB{N0f8}(0.1,0.1,0.1)
    @test convert(RGB{N0f8}, Gray24(0.1))   == RGB{N0f8}(0.1,0.1,0.1)

    grayN0f16 = Gray{N0f16}(0.8)
    @test convert(RGB, grayN0f16) == RGB{N0f16}(0.8,0.8,0.8)
    @test convert(RGB{Float32}, grayN0f16) == RGB{Float32}(0.8,0.8,0.8)

    @testset "$C-->Gray" for C in parametric3
        c = convert(C, RGB(1,1,1))
        g1 = convert(Gray, c)
        @test isa(g1, Gray)
        @test gray(g1) ≈ 1 atol=0.01
        g2 = convert(Gray{Float64}, c)
        @test typeof(g2) == Gray{Float64}
        @test gray(g2) ≈ 1 atol=0.01
    end

    # Images issue #382
    @test convert(Gray, RGBA(1,1,1,1)) == Gray(N0f8(1))

    # https://github.com/timholy/Images.jl/pull/445#issuecomment-189866806
    @test convert(Gray, RGB{N0f8}(0.145,0.145,0.145)) == Gray{N0f8}(0.145)

    # More AbstractRGB
    r4 = RGB4(1,0,0)
    @test convert(RGB, r4) == RGB(1,0,0)
    @test convert(RGB{N0f8}, r4) == RGB{N0f8}(1,0,0)
    @test convert(RGB4{N0f8}, r4) == RGB4{N0f8}(1,0,0)
    @test convert(RGB4{Float32}, r4) == RGB4{Float32}(1,0,0)
    @test convert(BGR{Float32}, r4) == BGR{Float32}(1,0,0)

    # Issue #257
    c = RGB{Float16}(0.9473,0.962,0.9766)
    hsi = convert(HSI, c)
    @test hsi.i > 0.96 && hsi.h ≈ 210

    # Test accuracy of conversion
    csconv = jldopen(joinpath(dirname(@__FILE__), "test_conversions.jld2")) do file
        read(file, "csconv")
    end

    # Since `colordiff`(e.g. `DE_2000`) involves a color space conversion, it is
    # not suitable for evaluating the conversion itself. On the other hand,
    # since the tolerance varies from component to component, a homogeneous
    # error evaluation function (e.g. a simple sum of differences) is also not
    # appropriate. Therefore, a series of `diffnorm`, which returns the
    # normalized Euclidian distance, is defined as follows. They are just for
    # testing purposes as the cyclicity of hue is ignored.
    sqd(a, b, s=1.0) = ((float(a) - float(b))/s)^2
    function diffnorm(a::T, b::T) where {T<:Color3} # RGB,XYZ,xyY,LMS
        sqrt(sqd(comp1(a), comp1(b)) + sqd(comp2(a), comp2(b)) + sqd(comp3(a),comp3(b)))/sqrt(3)
    end
    function diffnorm(a::T, b::T) where {T<:Union{HSV,HSL,HSI}}
        sqrt(sqd(a.h, b.h, 360) + sqd(a.s, b.s) + sqd(comp3(a), comp3(b)))/sqrt(3)
    end
    function diffnorm(a::T, b::T) where {T<:Union{Lab,Luv}}
        sqrt(sqd(a.l, b.l, 100) + sqd(comp2(a), comp2(b), 200) + sqd(comp3(a), comp3(b), 200))/sqrt(3)
    end
    function diffnorm(a::T, b::T) where {T<:Union{LCHab,LCHuv}}
        sqrt(sqd(a.l, b.l, 100) + sqd(a.c, b.c, 100) + sqd(a.h, b.h, 360))/sqrt(3)
    end
    function diffnorm(a::T, b::T) where {T<:Union{DIN99,DIN99d,DIN99o}} # csconv has no DIN99 case
        sqrt(sqd(a.l, b.l, 100) + sqd(a.a, b.a, 100) + sqd(a.b, b.b, 100))/sqrt(3)
    end
    function diffnorm(a::T, b::T) where {T<:YIQ}
        sqrt(sqd(a.y, b.y) + sqd(a.i, b.i, 1.2) + sqd(a.q, b.q, 1.2))/sqrt(3)
    end
    function diffnorm(a::T, b::T) where {T<:YCbCr}
        sqrt(sqd(a.y, b.y, 219) + sqd(a.cb, b.cb, 224) + sqd(a.cr, b.cr, 224))/sqrt(3)
    end

    for C in ColorTypes.parametric3
        y = convert(C, RGB(1.0,1.0,0.0))
        b = convert(C, RGB(0.1,0.1,0.2))
        diffnorm(y, b) < 0.5 && @warn("`diffnorm` for $C may be broken")
    end

    function convcompare(from, to, tol; showfailure::Bool=false)
        errmax = 0.0
        for i = 1:length(from)
            t = to[i]
            f = convert(typeof(t), from[i])
            diff = diffnorm(t, f)
            if showfailure && diff>tol
                original = from[i]
                @show original f t diff
            end
            errmax = max(errmax, diff)
        end
        errmax > tol && @warn("Error on conversions from $(eltype(from)) to $(eltype(to)), relative error = $errmax")
        errmax <= tol
    end

    base_t(i, from_to) = base_color_type(eltype(csconv[i][from_to]))
    @testset "accuracy test: $(base_t(i,1))-->$(base_t(i,2)) (index $i)" for i = 1:length(csconv)
        f, t = csconv[i]
        @test convcompare(f, t, 2e-3, showfailure=false)
    end
end


text:
using RecursiveArrayTools: ArrayPartition
using ProximalCore

struct ComposeAffine
    f
    A
    b
end

function (g::ComposeAffine)(x)
    res = g.A * x .+ g.b
    return g.f(res)
end

function compose_affine_gradient!(y, g::ComposeAffine, x)
    res = g.A * x .+ g.b
    gradres, v = gradient(g.f, res)
    mul!(y, adjoint(g.A), gradres)
    return v
end

ProximalCore.gradient!(y, g::ComposeAffine, x) = compose_affine_gradient!(y, g, x)
ProximalCore.gradient!(y::ArrayPartition, g::ComposeAffine, x::ArrayPartition) = compose_affine_gradient!(y, g, x)

function ProximalCore.gradient(h::ComposeAffine, x::ArrayPartition)
    grad_fx = similar(x)
    fx = ProximalCore.gradient!(grad_fx, h, x)
    return grad_fx, fx
end


text:
## This is until we have a good replacement for `info(io, msg)` on 0.7+
function info(io::IO, args...)
    printstyled(io, "Info: "; color=Base.info_color())
    printstyled(io, args...; color=Base.info_color())
    println(io)
end
function warn(io::IO, args...)
    printstyled(io, "Warning: "; color=Base.warn_color())
    printstyled(io, args...; color=Base.warn_color())
    println(io)
end

function extract_kwargs(kwargs, keys)
    return (k => v for (k, v) in kwargs if k in keys)
end


text:
function initialize_inner!(
    device_states,
    static::PSY.StaticInjection,
    dynamic_device::DynamicWrapper{
        PSY.DynamicInverter{C, O, PSY.VoltageModeControl, DC, P, F},
    },
    inner_vars::AbstractVector,
) where {
    C <: PSY.Converter,
    O <: PSY.OuterControl,
    DC <: PSY.DCSource,
    P <: PSY.FrequencyEstimator,
    F <: PSY.Filter,
}

    #Obtain external states inputs for component
    external_ix = get_input_port_ix(dynamic_device, PSY.VoltageModeControl)
    Ir_filter = device_states[external_ix[1]]
    Ii_filter = device_states[external_ix[2]]
    Ir_cnv = device_states[external_ix[3]]
    Ii_cnv = device_states[external_ix[4]]
    Vr_filter = device_states[external_ix[5]]
    Vi_filter = device_states[external_ix[6]]

    #Obtain inner variables for component
    ω_oc = get_ω_ref(dynamic_device)
    θ0_oc = inner_vars[θ_oc_var]
    Vdc = inner_vars[Vdc_var]

    #Obtain output of converter
    Vr_cnv0 = inner_vars[Vr_cnv_var]
    Vi_cnv0 = inner_vars[Vi_cnv_var]

    #Get Voltage Controller parameters
    inner_control = PSY.get_inner_control(dynamic_device)
    filter = PSY.get_filter(dynamic_device)
    kpv = PSY.get_kpv(inner_control)
    kiv = PSY.get_kiv(inner_control)
    kffi = PSY.get_kffi(inner_control)
    cf = PSY.get_cf(filter)
    rv = PSY.get_rv(inner_control)
    lv = PSY.get_lv(inner_control)

    #Get Current Controller parameters
    kpc = PSY.get_kpc(inner_control)
    kic = PSY.get_kic(inner_control)
    kffv = PSY.get_kffv(inner_control)
    lf = PSY.get_lf(filter)
    ωad = PSY.get_ωad(inner_control)
    kad = PSY.get_kad(inner_control)

    function f!(out, x)
        θ_oc = x[1]
        v_refr = x[2]
        ξ_d = x[3]
        ξ_q = x[4]
        γ_d = x[5]
        γ_q = x[6]
        ϕ_d = x[7]
        ϕ_q = x[8]

        #Reference Frame Transformations
        I_dq_filter = ri_dq(θ_oc + pi / 2) * [Ir_filter; Ii_filter]
        I_dq_cnv = ri_dq(θ_oc + pi / 2) * [Ir_cnv; Ii_cnv]
        V_dq_filter = ri_dq(θ_oc + pi / 2) * [Vr_filter; Vi_filter]
        V_dq_cnv0 = ri_dq(θ_oc + pi / 2) * [Vr_cnv0; Vi_cnv0]

        #Voltage controller references
        Vd_filter_ref = (v_refr - rv * I_dq_filter[d] + ω_oc * lv * I_dq_filter[q])
        Vq_filter_ref = (-rv * I_dq_filter[q] - ω_oc * lv * I_dq_filter[d])

        #Current controller references
        Id_cnv_ref = (
            kpv * (Vd_filter_ref - V_dq_filter[d]) + kiv * ξ_d -
            cf * ω_oc * V_dq_filter[q] + kffi * I_dq_filter[d]
        )
        Iq_cnv_ref = (
            kpv * (Vq_filter_ref - V_dq_filter[q]) +
            kiv * ξ_q +
            cf * ω_oc * V_dq_filter[d] +
            kffi * I_dq_filter[q]
        )

        #References for Converter Output Voltage
        Vd_cnv_ref = (
            kpc * (Id_cnv_ref - I_dq_cnv[d]) + kic * γ_d - ω_oc * lf * I_dq_cnv[q] +
            kffv * V_dq_filter[d] - kad * (V_dq_filter[d] - ϕ_d)
        )
        Vq_cnv_ref = (
            kpc * (Iq_cnv_ref - I_dq_cnv[q]) +
            kic * γ_q +
            ω_oc * lf * I_dq_cnv[d] +
            kffv * V_dq_filter[q] - kad * (V_dq_filter[q] - ϕ_q)
        )

        out[1] = Vd_filter_ref - V_dq_filter[d]
        out[2] = Vq_filter_ref - V_dq_filter[q]
        out[3] = Id_cnv_ref - I_dq_cnv[d]
        out[4] = Iq_cnv_ref - I_dq_cnv[q]
        out[5] = V_dq_filter[d] - ϕ_d
        out[6] = V_dq_filter[q] - ϕ_q
        out[7] = Vd_cnv_ref - V_dq_cnv0[d]
        out[8] = Vq_cnv_ref - V_dq_cnv0[q]
    end
    x0 = [θ0_oc, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
    sol = NLsolve.nlsolve(f!, x0)
    if !NLsolve.converged(sol)
        @warn("Initialization in Inner Control failed")
    else
        sol_x0 = sol.zero
        #Update angle:
        inner_vars[θ_oc_var] = sol_x0[1]
        outer_ix = get_local_state_ix(dynamic_device, O)
        outer_states = @view device_states[outer_ix]
        #Assumes that angle is in second position
        outer_states[1] = sol_x0[1]
        inner_vars[θ_oc_var] = sol_x0[1]
        set_V_ref(dynamic_device, sol_x0[2])
        PSY.set_V_ref!(
            PSY.get_reactive_power(PSY.get_outer_control(dynamic_device)),
            sol_x0[2],
        )
        inner_vars[V_oc_var] = sol_x0[2]
        #Update Converter modulation
        m0_dq = (ri_dq(sol_x0[1] + pi / 2) * [Vr_cnv0; Vi_cnv0]) ./ Vdc
        inner_vars[md_var] = m0_dq[d]
        inner_vars[mq_var] = m0_dq[q]
        #Update states
        inner_ix = get_local_state_ix(dynamic_device, PSY.VoltageModeControl)
        inner_states = @view device_states[inner_ix]
        inner_states[1] = sol_x0[3] #ξ_d
        inner_states[2] = sol_x0[4] #ξ_q
        inner_states[3] = sol_x0[5] #γ_d
        inner_states[4] = sol_x0[6] #γ_q
        inner_states[5] = sol_x0[7] #ϕ_d
        inner_states[6] = sol_x0[8] #ϕ_q
    end
    return
end

function initialize_inner!(
    device_states,
    static::PSY.StaticInjection,
    dynamic_device::DynamicWrapper{
        PSY.DynamicInverter{C, O, PSY.CurrentModeControl, DC, P, F},
    },
    inner_vars::AbstractVector,
) where {
    C <: PSY.Converter,
    O <: PSY.OuterControl,
    DC <: PSY.DCSource,
    P <: PSY.FrequencyEstimator,
    F <: PSY.Filter,
}

    #Obtain external states inputs for component
    external_ix = get_input_port_ix(dynamic_device, PSY.CurrentModeControl)

    Ir_cnv = device_states[external_ix[3]]
    Ii_cnv = device_states[external_ix[4]]
    Vr_filter = device_states[external_ix[5]]
    Vi_filter = device_states[external_ix[6]]

    #Obtain inner variables for component
    ω_oc = get_ω_ref(dynamic_device)
    θ0_oc = inner_vars[θ_freq_estimator_var]
    Vdc = inner_vars[Vdc_var]
    Id_cnv_ref = inner_vars[Id_oc_var]
    Iq_cnv_ref = inner_vars[Iq_oc_var]

    #Obtain output of converter
    Vr_cnv0 = inner_vars[Vr_cnv_var]
    Vi_cnv0 = inner_vars[Vi_cnv_var]

    #Get Current Controller parameters
    inner_control = PSY.get_inner_control(dynamic_device)
    filter = PSY.get_filter(dynamic_device)
    kpc = PSY.get_kpc(inner_control)
    kic = PSY.get_kic(inner_control)
    kffv = PSY.get_kffv(inner_control)
    lf = PSY.get_lf(filter)

    function f!(out, x)
        γ_d = x[1]
        γ_q = x[2]

        #Reference Frame Transformations
        I_dq_cnv = ri_dq(θ0_oc + pi / 2) * [Ir_cnv; Ii_cnv]
        V_dq_filter = ri_dq(θ0_oc + pi / 2) * [Vr_filter; Vi_filter]
        V_dq_cnv0 = ri_dq(θ0_oc + pi / 2) * [Vr_cnv0; Vi_cnv0]

        #References for Converter Output Voltage
        Vd_cnv_ref = (
            kpc * (Id_cnv_ref - I_dq_cnv[d]) + kic * γ_d - ω_oc * lf * I_dq_cnv[q] +
            kffv * V_dq_filter[d]
        )
        Vq_cnv_ref = (
            kpc * (Iq_cnv_ref - I_dq_cnv[q]) +
            kic * γ_q +
            ω_oc * lf * I_dq_cnv[d] +
            kffv * V_dq_filter[q]
        )

        out[1] = Vd_cnv_ref - V_dq_cnv0[d]
        out[2] = Vq_cnv_ref - V_dq_cnv0[q]
    end
    x0 = [0.0, 0.0]
    sol = NLsolve.nlsolve(f!, x0)
    if !NLsolve.converged(sol)
        @warn("Initialization in Inner Control failed")
    else
        sol_x0 = sol.zero
        #Update Converter modulation
        m0_dq = (ri_dq(θ0_oc + pi / 2) * [Vr_cnv0; Vi_cnv0]) ./ Vdc
        inner_vars[md_var] = m0_dq[d]
        inner_vars[mq_var] = m0_dq[q]
        #Update states
        inner_ix = get_local_state_ix(dynamic_device, PSY.CurrentModeControl)
        inner_states = @view device_states[inner_ix]
        inner_states[1] = sol_x0[1] #γ_d
        inner_states[2] = sol_x0[2] #γ_q
    end
    return
end

function initialize_inner!(
    device_states,
    static::PSY.StaticInjection,
    dynamic_device::DynamicWrapper{
        PSY.DynamicInverter{C, O, PSY.RECurrentControlB, DC, P, F},
    },
    inner_vars::AbstractVector,
) where {
    C <: PSY.Converter,
    O <: PSY.OuterControl,
    DC <: PSY.DCSource,
    P <: PSY.FrequencyEstimator,
    F <: PSY.Filter,
}
    # Obtain inner variables for component
    V_R = inner_vars[Vr_inv_var]
    V_I = inner_vars[Vi_inv_var]
    V_t = sqrt(V_R^2 + V_I^2)
    Iq_cmd = inner_vars[Iq_ic_var]
    Ip_oc = inner_vars[Id_oc_var]
    Iq_oc = inner_vars[Iq_oc_var]

    #Get Current Controller parameters
    inner_control = PSY.get_inner_control(dynamic_device)
    Q_Flag = PSY.get_Q_Flag(inner_control)
    PQ_Flag = PSY.get_PQ_Flag(inner_control)

    Ip_min, Ip_max, Iq_min, Iq_max =
        current_limit_logic(inner_control, Base.RefValue{PQ_Flag}, V_t, Ip_oc, Iq_cmd)

    if Ip_oc >= Ip_max + BOUNDS_TOLERANCE || Ip_min - BOUNDS_TOLERANCE >= Ip_oc
        @error(
            "Inverter $(PSY.get_name(static)) current $(Ip_oc) out of limits $(Ip_min) $(Ip_max). Check Power Flow or Parameters"
        )
    end

    if Iq_oc >= Iq_max + BOUNDS_TOLERANCE || Iq_min - BOUNDS_TOLERANCE >= Iq_oc
        @error(
            "Inverter $(PSY.get_name(static)) current $(Iq_oc) out of limits $(Iq_min) $(Iq_max). Check Power Flow or Parameters"
        )
    end

    if Q_Flag == 0
        local_ix = get_local_state_ix(dynamic_device, PSY.RECurrentControlB)
        #Define internal states for Inner Control
        internal_states = @view device_states[local_ix]
        internal_states[1] = V_t
        internal_states[2] = Iq_oc
    else
        local_ix = get_local_state_ix(dynamic_device, PSY.RECurrentControlB)
        K_vi = PSY.get_K_vi(inner_control)
        #Define internal states for Inner Control
        internal_states = @view device_states[local_ix]
        internal_states[1] = V_t
        internal_states[2] = Iq_cmd / K_vi
    end

    #Update additional variables
    # Based on PSS/E manual, if user does not provide V_ref0, then
    # V_ref0 is considered to be the output voltage of the PF solution
    if PSY.get_V_ref0(inner_control) == 0.0
        PSY.set_V_ref0!(inner_control, V_t)
    end
    return
end


text:
# Load ImplicitEquations

## meta these
preds = [(:Lt, :≪, :<), # \ll
         (:Le, :≦, :<=), # \leqq
         (:Eq, :⩵, :(==)), # \Equal
         (:Ge, :≧, :>=), # \gegg
         (:Gt, :≫, :>) # \gg
         ]

import SymPy: Lt, Le, Eq, Ge, Gt, ≪, ≦, ⩵, ≧, ≫
for (fn, uop, op) in preds
    fnname =  string(fn)

    @eval begin
        ($fn)(f::Function, x::Real) = Pred(f, $op, x)
        ($uop)(f::Function, x::Real) = ($fn)(f, x)
        ($fn)(f::Function, g::Function) = $(fn)((x,y) -> f(x,y) - g(x,y), 0)
        ($uop)(f::Function, g::Function) = ($fn)(f, g)
    end
    eval(Expr(:export, fn))
    eval(Expr(:export, uop))
end
import ImplicitEquations: Neq
export Neq

# some plotting utilities

"""
   trimplot(f, a, b, c=20; kwargs...)

Plot f over [a,b] but break graph if it exceeds c in absolute value.
"""
@userplot TrimPlot
@recipe function __(a::TrimPlot)
    if length(a.args) < 4
        f, x0, x1 = a.args
        c = 20
    else
        f, x0, x1, c = a.args
    end

    xs = range(x0, x1, length=251)
    ys = f.(xs)

    Xs = Any[]
    Ys = Any[]

    us, vs = Real[], Real[]
    for (x,y) in zip(xs, ys)
        if abs(y) <= c
            push!(us, x); push!(vs, y)
        else
            if length(us) > 0
                push!(Xs, copy(us))
                push!(Ys, copy(vs))
            end
            empty!(us); empty!(vs)
        end
    end
    if length(us) > 0
        push!(Xs, copy(us))
        push!(Ys, copy(vs))
    end

    seriestype := :path
    x := Xs
    y := Ys

    Xs, Ys
end

"""
    plotif(f, g, a, b)

Plot f colored depending on g < 0 or not.
"""
@userplot PlotIf
@recipe function __(a::PlotIf)
    f, g, x0, x1 = a.args
    xs = range(x0, x1, length=251)
    ys = f.(xs)

    seriestype := :path
    x := xs
    y := ys
    color := [gx < 0 ? :red : :blue for gx in g.(xs)]
    linewidth := 5
    legend := false

    xs, ys
end

@userplot SignChart
@recipe function __(a::SignChart)
    f, x0, x1 = a.args
    xs = range(x0, x1, length=251)
    ys = f.(xs)
    zs = zeros(length(xs))
    @series begin
        seriestype := :path
        x := xs
        y := zs
        linewidth := 1
        color := :black
        xs, zs
    end

    seriestype := :path
    x := xs
    y := ys
    color := [y < 0 ? :red : :blue for y in ys]

    linewidth --> 5
    legend --> false

    xs, ys
end

#### Multidimensional

"""
    plot_polar(r, a, b; [n=251]

Create a polar plot of `r(theta)` for `a <= theta <= b`.

An alternative to the three step approach:

```
ts = range(a, b, length=n)
rs = r.(ts)
plot(ts, rs, proj=:polar, legend=false)
```

"""
@userplot Plot_Polar
@recipe function __(rt::Plot_Polar; n=251)

    r, a, b = rt.args
    fn = t -> r(t) * [cos(t), sin(t)]
    ts = range(a, stop=b, length=n)
    xs = fn.(ts)

    xys = unzip(xs)
    seriestype := :path
    x := xys[1]
    y := xys[2]
    xys
end

"""
    plot_parametric_curve(r, a, b; [n=251])

For `r` taking `R -> R^n` make a plot over [a,b] using n points.

There are other ways to produce this plot:
* if we have two functions, `rx` and `ry` say, then passing as function arguments should work

```
plot(rx, ry, a, b)
```

* if we have a vector-valued function, then either of these work:

```
r(t) = [cos(t), sin(t), t]
plot(unzip(r, a, b, n)...)  # unzip reshapes the data
plot_parametric_curve(r, a, b)
```

We see that this function is a convenience for vector-valued descriptions of parameterized curves.
"""
@userplot Plot_Parametric_Curve
@recipe function __(r::Plot_Parametric_Curve; n=251)

    f, a, b = r.args
    ts = range(a, stop=b, length=n)
    xs = f.(ts)

    n = length(first(xs))
    if n == 2
        xys = unzip(xs)
        seriestype := :path
        x := xys[1]
        y := xys[2]
        xys
    else
        xyzs = unzip(xs)
        seriestype := :path3d
        x := xyzs[1]
        y := xyzs[2]
        z := xyzs[3]
        xyzs
    end
end

import SymPy: plot_parametric_surface

"""
    plot_parametric_surface(F, [xlim=(-5,5)], [ylim=(-5,5)], [nx=50], [ny=50])

Plot the paremeterized surface described by `F:R^2 -> R^3`. This works with `pyplot` and `plotly`, but not `gr`.

A surface `(x,y,f(x,y))` can be directly plotted with `surface` as:

```
surface(xs, ys, f)
```

where `xs` and `ys` provide a grid to plot over.

For parametrically described surfaces, the above doesn't work. This function provides an interface:

```
F(u,v) = [u*cos(v), u*sin(v), u] # a cone
plot_parametric_surface(F, xlim=(0, 1), ylim=(0,2pi))
```

The values of `xlim` are used as the range to plot over; `nx` specifies the number of points to use. Similarly for `ylim`.
"""
@userplot Plot_Parametric_Surface
@recipe function __(r::Plot_Parametric_Surface; nx=50, ny=50)

    F = first(r.args)

    xlim = get(plotattributes,:xlims, (-5,5))
    ylim = get(plotattributes,:ylims, (-5,5))

    us = range(xlim[1], xlim[2], length=nx)
    vs = range(ylim[1], ylim[2], length=ny)
    ws = unzip(F.(us, vs'))

    seriestype := :surface
    x := ws[1]
    y := ws[2]
    z := ws[3]
    xlims := extrema(ws[1])
    ylims := extrema(ws[2])
    ws

end

"""
   arrow(p, v)

Add vector, `v`, to plot anchored at point `p`.

Example
```
Fn = parametric(t -> [2cos(t), 3sin(t)])
Fnp = t -> ForwardDiff.derivative(Fn, t)
p = plot(Fn, 0, 2pi, legend=false)
for t in 0:pi/4:pi
   arrow!(Fn(t), Fnp(t))
end
p
```
"""
@userplot Arrow
@recipe function __(a::Arrow)
    v0, dv0 = a.args
    v = [v0...]; dv = [dv0...]
    if length(v) == 2
        seriestype := :quiver
        quiver := unzip([dv])
        x,y = unzip([v])
        x := x
        y := y
        ()
    else
        seriestype := :path3d
        M = hcat(v, v + dv)'
        (M[:,1],M[:,2], M[:,3])
    end

end

"""
    vectorfieldplot(F; [xlim=(-5,5)], [ylim=(-5,5)], [nx=8], [ny=8])

Create a vector field plot using a grid described by `xlim`, `ylim` with `nx` and `ny` grid points in each direction.

```
F(x,y) = [-y, x]
vectorfieldplot(F, xlim=(-4,4), ylim=(-4,4))
```
"""
@userplot VectorFieldPlot
@recipe function __(V::VectorFieldPlot; nx=8, ny=8)

    ars = V.args
    F = ars[1]
    xlim = get(plotattributes,:xlims, (-5,5))
    ylim = get(plotattributes,:ylims, (-5,5))

    dx, dy = (xlim[2]-xlim[1])/(nx-1), (ylim[2]-ylim[1])/(ny-1)
    xs, ys = xlim[1]:dx:xlim[2], ylim[1]:dy:ylim[2]

    ps = [[x,y] for x in xs for y in ys]
    vs = (xy -> F(xy[1], xy[2])).(ps)

    λ = 0.9 * minimum([u/maximum(getindex.(vs,i)) for (i,u) in enumerate((dx,dy))])

    seriestype := :quiver
    quiver := unzip(λ * vs)
    xs, ys = unzip(ps)
        x := xs
    y := ys
    ()
end

"""
    vectorfieldplot3d(F; [xlim=(-5,5)], [ylim=(-5,5)], [nx=5], [ny=5])

Create a 3 dimensional vector field plot using a grid described by `xlim`, `ylim`, `zlim` with `nx`, `ny`, and `nz` grid points in each direction.

Note: the vectors are represented with line, not arrow due to no implementation of `:quiver3d`.

```
F(x,y,z) = [-y, x,z]
vectorfieldplot3d(F, xlims=(-4,4), ylims=(-4,4), zlims=(0,3))
```
"""
@userplot VectorFieldPlot3d
@recipe function __(V::VectorFieldPlot3d;  nx=5, ny=5, nz=5)
    ars = V.args
    F = ars[1]
    xlim = get(plotattributes,:xlims, (-5,5))
    ylim = get(plotattributes,:ylims, (-5,5))
    zlim = get(plotattributes,:zlims, (-5,5))

    dx, dy, dz = (xlim[2]-xlim[1])/(nx-1), (ylim[2]-ylim[1])/(ny-1), (zlim[2]-zlim[1])/(nz-1)
    xs, ys, zs = xlim[1]:dx:xlim[2], ylim[1]:dy:ylim[2], zlim[1]:dz:zlim[2]

    ps = [[x,y,z] for x in xs for y in ys for z in zs]
    vs = (xyz -> F(xyz[1], xyz[2], xyz[3])).(ps)

    λ = 0.9 * minimum([u/maximum(getindex.(vs,i)) for (i,u) in enumerate((dx,dy,dx))])

    seriestype := :path3d
    dxs, dys, dzs = unzip(λ * vs)
    xs, ys, zs = unzip(ps)

    legend := false
    [[x,dx] for (x,dx) in zip(xs, xs+dxs)],    [[y,dy] for (y,dy) in zip(ys, ys+dys)],    [[z,dz] for (z,dz) in zip(zs, zs+dzs)]
end

import Contour # installed with the Plots package, so should be available
               # import -- not using -- to avoid name collision

"""
    Visualize `F(x,y,z) = c` by plotting assorted contour lines

This graphic makes slices in the `x`, `y`, and/or `z` direction of the 3-D level surface and plots them accordingly. Which slices (and their colors) are specified through a dictionary.

Examples:
```
F(x,y,z) = x^2 + y^2 + x^2
plot_implicit(F, 20)  # 20 slices in z direction
plot_implicit(F, 20, slices=Dict(:x=>:blue, :y=>:red, :z=>:green), nlevels=6) # all 3 shown

# A heart
a,b = 1,3
F(x,y,z) = (x^2+((1+b)*y)^2+z^2-1)^3-x^2*z^3-a*y^2*z^3
plot_implicit_surface(F, xlims=(-2,2),ylims=(-1,1),zlims=(-1,2))
```

Note: Idea [from](https://stackoverflow.com/questions/4680525/plotting-implicit-equations-in-3d).

Not exported?
"""
function plot_implicit_surface(F, c=0;
                       xlims=(-5,5), ylims=xlims, zlims=xlims,
                       nlevels=25,         # number of levels in a direction
                       slices=Dict(:z => :blue), # Dict(:x => :color, :y=>:color, :z=>:color)
                       kwargs...          # passed to initial `plot` call
                       )

    _linspace(rng, n=150) = range(rng[1], stop=rng[2], length=n)

    X1, Y1, Z1 = _linspace(xlims), _linspace(ylims), _linspace(zlims)

    p = Plots.plot(;legend=false,kwargs...)

    if :x ∈ keys(slices)
        for x in _linspace(xlims, nlevels)
            local X1 = [F(x,y,z) for y in Y1, z in Z1]
            cnt = Contour.contours(Y1,Z1,X1, [c])
            for line in Contour.lines(Contour.levels(cnt)[1])
                ys, zs = Contour.coordinates(line) # coordinates of this line segment
                plot!(p, x .+ 0 * ys, ys, zs, color=slices[:x])
          end
        end
    end

    if :y ∈ keys(slices)
        for y in _linspace(ylims, nlevels)
            local Y1 = [F(x,y,z) for x in X1, z in Z1]
            cnt = Contour.contours(Z1,X1,Y1, [c])
            for line in Contour.lines(Contour.levels(cnt)[1])
                xs, zs = Contour.coordinates(line) # coordinates of this line segment
                plot!(p, xs, y .+ 0 * xs, zs, color=slices[:y])
            end
        end
    end

    if :z ∈ keys(slices)
        for z in _linspace(zlims, nlevels)
            local Z1 = [F(x, y, z) for x in X1, y in Y1]
            cnt = Contour.contours(X1, Y1, Z1, [c])
            for line in Contour.lines(Contour.levels(cnt)[1])
                xs, ys = Contour.coordinates(line) # coordinates of this line segment
                plot!(p, xs, ys, z .+ 0 * xs, color=slices[:z])
            end
        end
    end

    p
end


text:
# DistDensity.jl
dat = [[0.0 1.0];[1.0 -1.0];[2.0 2.0];[2.0 3.0];[1.0 1.0]]
D = dist_matrix(dat)
knn_dists_out = knn_dists(D, 2, 0)
@test all(vec(knn_dists_out[4]) .== [5,5,4,3,1,2,1,5,5,3]) || all(vec(knn_dists_out[4]) .== [5,5,4,3,1,3,1,5,5,3])
# parameters
Q = zeros(Float64, 2,2); Q[1,1] = 1.0; Q[2,2] = 1.0
sigma = median(D)
K = kernel_matrix(D, sigma)
#####################
# DetectionAlgorithms.jl
# model
svdd_model = SVDD_train(K[1:4,1:4], 0.2);
knfst_model = KNFST_train(K[1:4,1:4])
# REC, checked manually
@test all(round.(REC(D, sigma, 0), digits = 1) .== round.(1 .- ([1,0,2,1,2] ./ 5), digits = 1))
# KNN Gamma, checked manually
@test all(round.(KNN_Gamma(knn_dists_out), digits = 1) .== [1.6, 2.1,1.2,1.6,1.2])
# KNN Delta, approximately
@test all(round.(KNN_Delta(knn_dists_out, dat), digits = 1) .== [1.4, 2.1, 0.5, 1.6, 0.5])
# KDE # results show exhibit similar ordering like REC
@test all(sortperm(KDE(K)) .== [5,3,1,4,2])
# Hotelling's T^2
# is also the quared mahalanobis distance to the data's mean
using Distances
@test all(round.(T2(dat, Q, mean(dat, dims = 1)), digits = 2) .== round.(pairwise(SqMahalanobis(Q), dat', mean(dat', dims = 2), dims = 2), digits = 2))

# SVDD
Ktest = exp.(-0.5 * pairwise(Euclidean(), dat[1:4,:]', dat') ./ sigma^2)
@test all(SVDD_predict(svdd_model, Ktest)[1] .== [0, 1, 0, 0, 1])
# KNFST, last data point (not seen in training) should differ, i.e. have largest values
@test sortperm(-KNFST_predict(knfst_model, Ktest)[1])[1] == 5

# high level functions
algorithms = ["T2", "REC", "KDE", "KNN_Gamma", "KNN_Delta"]
P = getParameters(algorithms, dat)
P.KNN_k = 2
P.Q = Q
@test P.K_sigma == sigma == P.REC_varepsilon
# detectAnomalies
detectAnomalies(dat, P)
# chekct detectAnomalies for self consistency
@test round.(P.REC, digits = 3) == round.(REC(D, sigma, 0), digits = 3)
@test round.(P.KDE, digits = 3) == round.(KDE(K), digits = 3)
@test round.(P.KNN_Gamma, digits = 3) ==  round.(KNN_Gamma(knn_dists_out), digits = 3)
@test round.(P.KNN_Delta[1], digits = 3) ==  round.(KNN_Delta(knn_dists_out, dat), digits = 3)
@test round.(P.T2[1], digits = 3) ==  round.(T2(dat, Q, mean(dat, dims = 1)), digits = 3)

algorithms = ["KNFST", "SVDD"]
P = getParameters(algorithms, dat[1:4,:])
P.K_sigma = sigma
P.SVDD_model = svdd_model
P.KNFST_model = knfst_model
detectAnomalies(dat, P)
(labels, decvalues) = SVDD_predict(svdd_model, Ktest)
@test round.(P.SVDD, digits = 3) == round.(decvalues, digits = 3) * -1
@test K[1:4,1:4] == P.D_train[1]
@test Ktest == P.D_test[1]
@test round.(P.KNFST[1], digits = 3) == round.(KNFST_predict(knfst_model, Ktest)[1], digits = 3)

data = rand(200,3)
methods = ["REC","KDE"]
@test all(detectAnomalies(data, getParameters(methods, data)) .== detectAnomalies(data, methods))

x = [[1.0,2.0,3.0] [-2.0,4.0,1.0]]
@test all(round.(UNIV(x), digits = 0) .== [0,1,1])

# test online algorithms with the non-online corresponding ones
using StatsBase

out = zeros(100)
x = randn(100, 3)
Q = StatsBase.cov(x)

@test all(round.(KDEonline!(out, x, 1.0), digits = 4) .== round.(KDE(kernel_matrix(dist_matrix(x), 1.0)), digits = 4))

@test all(round.(KDEonline!(out, x, Q, 1.0), digits = 4) .== round.(KDE(kernel_matrix(dist_matrix(x, dist = "Mahalanobis", Q = Q), 1.0)), digits = 4))

@test all(round.(REConline!(out, x, 1.0), digits = 4) .== round.(REC(dist_matrix(x), 1.0), digits = 4))

@test all(round.(REConline!(out, x, Q, 1.0), digits = 4) .== round.(REC(dist_matrix(x, dist = "Mahalanobis", Q = Q), 1.0), digits = 4))

@test all(round.(KNNonline!(out, x, 5), digits = 4) .== round.(KNN_Gamma(knn_dists(dist_matrix(x, dims = 2), 5)), digits = 4))

@test all(round.(KNNonline!(out, x, Q, 5), digits = 4) .== round.(KNN_Gamma(knn_dists(dist_matrix(x, dist = "Mahalanobis",Q = Q, dims = 2), 5)), digits = 4))

s = zeros(1)

SigmaOnline!(s, x)

@test floor.(Int, s)[1] == 2 ||  floor.(Int, s)[1] == 1

SigmaOnline!(s, x, Q)

@test floor.(Int, s)[1] == 2 ||  floor.(Int, s)[1] == 1


text:
function tellastory(http)
    story = """There once was a Julia programmer named Emmett.
    Emmett sought to create a great web-development framework inside of
    Julia. Although it was certainly a lot of work, and at times was
    surprisingly challenging, you are now viewing the results!"""
    for character in story
        write(http, string(character))
        sleep(.07)
    end
end
heading = html("""

Welcome to ToolipsApp

""")
sub_heading = html("
A collection of Toolips examples
")
summary = html("""
Toolips.jl is a web-development framework for Julia.
The framework features a full front-end templating and building engine as well
as a functional API. The following are some examples which utilize the various
capabilities of Toolips.jl.

""")
story_button = Toolips.Button("story", onAction = tellastory, label = "Tell a story")
text_editor = html("""<form action="/text-editor">
    <button value="Text Editor" />Text Editor</button>
</form>""")
style = Toolips.css("""div {
  background-image: url('https://i.postimg.cc/JzD6s6Vx/bgexamp.png');
  background-repeat: no-repeat;
  background-size: cover;
  }
  button {
	background-color:#db8cdb;
	border-radius:22px;
	display:inline-block;
	cursor:pointer;
	color:#ffffff;
	font-family:Impact;
	font-size:17px;
	font-weight:bold;
	padding:17px 39px;
	text-decoration:none;
}
button:hover {
	background-color:#80b5ea;
}
button:active {
	position:relative;
	top:1px;
}
  """)
spacing = html("""

This software, and its contents are MIT-licensed.
You can view more information on Toolips here.

Created by Emmy Boudreau :)""")
page_contents = [style, heading, sub_heading, summary, story_button, text_editor, spacing]

home = Page(page_contents, "ToolipsApp")

max_issues_repo_name: emmettgb/ToolipsApp

max_forks_repo_name: emmettgb/ToolipsApp


text:
using Slomo: mass, density, Halos, Constants
using Roots: ConvergenceFailed
using Test

rtol = 1e-3
maxevals = 100

# NFW
for Mvir in exp10.([8.0, 10.0, 12.0, 14.0])
    for mdef in ["vir", "200c", "200m", "500c"]
        for z in [0.0, 0.5, 1.0]
            if mdef == "vir" || mdef == "200c"
                cvir = Halos.hmcr(Mvir; mdef = mdef, z = z)
            else
                cvir = 10.0
            end
            halo = Halos.NFW_from_virial(Mvir, cvir; mdef = mdef, z = z)
            rs = Halos.scale_radius(halo)
            Rvir = Halos.virial_radius(halo; mdef = mdef, z = z, rtol = rtol, maxevals = maxevals)
            @test isapprox(Mvir, mass(halo, Rvir), rtol = rtol)
            @test isapprox(cvir, Rvir / rs, rtol = rtol)
            @test isapprox(Mvir, Halos.virial_mass(halo; mdef = mdef, z = z, rtol = rtol, maxevals = maxevals), rtol = rtol)
            @test isapprox(cvir, Halos.concentration(halo; mdef = mdef, z = z, rtol = rtol, maxevals = maxevals), rtol = rtol)
            @test isapprox(Mvir, Halos.Mvir_from_Rvir(Rvir; mdef = mdef, z = z), rtol = rtol)
            @test isapprox(Rvir, Halos.Rvir_from_Mvir(Mvir; mdef = mdef, z = z), rtol = rtol)
            for rs_fraction in [1, 1.0, [1], [1.0], [0.5, 1.0, 1.0]]
                r = rs .* rs_fraction
                @test size(mass(halo, r)) == size(r)
                @test size(density(halo, r)) == size(r)
            end
        end
    end
end

# GNFW
for gamma in [eps()^0.5, 0.5, 1.0, 1.5, 2.0 - eps()^0.5]
    for Mvir in exp10.([8.0, 10.0, 12.0, 14.0])
        for mdef in ["vir", "200c", "200m", "500c"]
            for z in [0.0, 0.5, 1.0]
                if mdef == "vir" || mdef == "200c"
                    cvir = Halos.hmcr(Mvir; mdef = mdef, z = z)
                else
                    cvir = 10.0
                end
                halo = Halos.GNFW_from_virial(Mvir, cvir, gamma; mdef = mdef, z = z)
                rs = Halos.scale_radius(halo)
                Rvir = Halos.virial_radius(halo; mdef = mdef, z = z, rtol = rtol, maxevals = maxevals)
                @test isapprox(Mvir, mass(halo, Rvir), rtol = rtol)
                @test isapprox(cvir, Rvir / rs, rtol = rtol)
                @test isapprox(Mvir, Halos.virial_mass(halo; mdef = mdef, z = z, rtol = rtol, maxevals = maxevals), rtol = rtol)
                @test isapprox(cvir, Halos.concentration(halo; mdef = mdef, z = z, rtol = rtol, maxevals = maxevals), rtol = rtol)
                @test isapprox(Mvir, Halos.Mvir_from_Rvir(Rvir; mdef = mdef, z = z), rtol = rtol)
                @test isapprox(Rvir, Halos.Rvir_from_Mvir(Mvir; mdef = mdef, z = z), rtol = rtol)
                for rs_fraction in [1, 1.0, [1], [1.0], [0.5, 1.0, 1.0]]
                    r = rs .* rs_fraction
                    @test size(mass(halo, r)) == size(r)
                    @test size(density(halo, r)) == size(r)
                end
            end
        end
    end
end

# CoreNFW
for Mvir in exp10.([8.0, 10.0, 12.0, 14.0])
    for mdef in ["vir", "200c", "200m", "500c"]
        for z in [0.0, 0.5, 1.0]
            if mdef == "vir" || mdef == "200c"
                cvir = Halos.hmcr(Mvir; mdef = mdef, z = z)
            else
                cvir = 10.0
            end
            Re = 10.0
            t_sf = 10.0
            halo = Halos.CoreNFW_from_virial(Mvir, cvir, Re, t_sf; mdef = mdef, z = z)
            rs = Halos.scale_radius(halo)
            Rvir = Halos.virial_radius(halo; mdef = mdef, z = z, rtol = rtol, maxevals = maxevals)
            @test isapprox(Mvir, mass(halo, Rvir), rtol = rtol)
            @test isapprox(cvir, Rvir / rs, rtol = rtol)
            @test isapprox(Mvir, Halos.virial_mass(halo; mdef = mdef, z = z, rtol = rtol, maxevals = maxevals), rtol = rtol)
            @test isapprox(cvir, Halos.concentration(halo; mdef = mdef, z = z, rtol = rtol, maxevals = maxevals), rtol = rtol)
            @test isapprox(Mvir, Halos.Mvir_from_Rvir(Rvir; mdef = mdef, z = z), rtol = rtol)
            @test isapprox(Rvir, Halos.Rvir_from_Mvir(Mvir; mdef = mdef, z = z), rtol = rtol)
            for rs_fraction in [1, 1.0, [1], [1.0], [0.5, 1.0, 1.0]]
                r = rs .* rs_fraction
                @test size(mass(halo, r)) == size(r)
                @test size(density(halo, r)) == size(r)
            end
        end
    end
end

# Einasto
# for Mvir in exp10.([8.0, 10.0, 12.0, 14.0])
#     for mdef in ["vir", "200c", "200m", "500c"]
#         for z in [0.0, 0.5, 1.0]
#             if mdef == "vir" || mdef == "200c"
#                 cvir = Halos.hmcr(Mvir; mdef = mdef, z = z)
#             else
#                 cvir = 10.0
#             end
#             alpha = 0.16
#             halo = Halos.Einasto_from_virial(Mvir, cvir, alpha; mdef = mdef, z = z)
#             rs = Halos.scale_radius(halo)
#             Rvir = cvir * rs
#             Mvir_calc = Mvir
#             cvir_calc = cvir
#             try
#                 Rvir = Halos.virial_radius(halo; mdef = mdef, z = z, rtol = rtol, maxevals = maxevals)
#             catch err
#                 @warn(err, attempt="Rvir", Mvir=Mvir, cvir=cvir, mdef=mdef, z=z)
#             end
#             try
#                 Mvir_calc = Halos.virial_mass(halo; mdef = mdef, z = z, rtol = rtol, maxevals = maxevals)
#             catch err
#                 @warn(err, attempt="Mvir", Mvir=Mvir, cvir=cvir, mdef=mdef, z=z)
#             end
#             try
#                 cvir_calc = Halos.concentration(halo; mdef = mdef, z = z, rtol = rtol, maxevals = maxevals)
#             catch err
#                 @warn(err, attempt="cvir", Mvir=Mvir, cvir=cvir, mdef=mdef, z=z)
#             end
#             @test isapprox(Mvir, mass(halo, Rvir), rtol = rtol)
#             @test isapprox(cvir, Rvir / rs, rtol = rtol)
#             @test isapprox(Mvir, Mvir_calc, rtol = rtol)
#             @test isapprox(cvir, cvir_calc, rtol = rtol)
#             @test isapprox(Mvir, Halos.Mvir_from_Rvir(Rvir; mdef = mdef, z = z), rtol = rtol)
#             @test isapprox(Rvir, Halos.Rvir_from_Mvir(Mvir; mdef = mdef, z = z), rtol = rtol)
#             for rs_fraction in [1, 1.0, [1], [1.0], [0.5, 1.0, 1.0]]
#                 r = rs .* rs_fraction
#                 @test size(mass(halo, r)) == size(r)
#                 @test size(density(halo, r)) == size(r)
#             end
#         end
#     end
# end

# SolNFW
# for Mvir in exp10.([8.0, 10.0, 12.0, 14.0])
#     for mdef in ["vir", "200c", "200m", "500c"]
#         for z in [0.0, 0.5, 1.0]
#             if mdef == "vir" || mdef == "200c"
#                 cvir = Halos.hmcr(Mvir; mdef = mdef, z = z)
#             else
#                 cvir = 10.0
#             end
#             rs = 0.0
#             Rvir = 0.0
#             m22 = 5.0
#             halo = nothing
#             try
#                 halo = Halos.SolNFW_from_virial(Mvir, cvir, m22; mdef = mdef, z = z)
#             catch err
#                 if isa(err, ConvergenceFailed)
#                     @warn(err, halo="SolNFW", Mvir=Mvir, cvir=cvir, mdef=mdef, z=z)
#                     continue
#                 end
#             end
#             rs = Halos.scale_radius(halo)
#             Rvir = Halos.virial_radius(halo; mdef = mdef, z = z, rtol = rtol, maxevals = maxevals)
#             @test isapprox(Mvir, mass(halo, Rvir), rtol = rtol)
#             @test isapprox(cvir, Rvir / rs, rtol = rtol)
#             @test isapprox(Mvir, Halos.virial_mass(halo; mdef = mdef, z = z, rtol = rtol, maxevals = maxevals), rtol = rtol)
#             @test isapprox(cvir, Halos.concentration(halo; mdef = mdef, z = z, rtol = rtol, maxevals = maxevals), rtol = rtol)
#             @test isapprox(Mvir, Halos.Mvir_from_Rvir(Rvir; mdef = mdef, z = z), rtol = rtol)
#             @test isapprox(Rvir, Halos.Rvir_from_Mvir(Mvir; mdef = mdef, z = z), rtol = rtol)
#             for rs_fraction in [1, 1.0, [1], [1.0], [0.5, 1.0, 1.0]]
#                 r = rs .* rs_fraction
#                 @test size(mass(halo, r)) == size(r)
#                 @test size(density(halo, r)) == size(r)
#             end
#         end
#     end
# end


text:
module OpenGLDisplayUtils
export  basicRender
using ModernGL
using GLFW

"""
As most functions will deal with just addind the quad to the screen
and swapping buffers
"""
function basicRender(window)
	glDrawElements(GL_TRIANGLES, 6, GL_UNSIGNED_INT, C_NULL)
	# Swap front and back buffers
	GLFW.SwapBuffers(window)

end

end #..OpenGLDisplayUtils


text:
using Base: Float64
using Dates: string
println("---> START <---")
# Set up environment
using Pkg
Pkg.activate(".")
using JSON
using Glob
using CSV
using Random
using DataFrames

using Statistics
using StatsBase

using Distributions
using DataStructures
using Parameters

using Plots
using Colors
using ColorSchemes
using Dates
using Plotly

include("post_simulation_plots_functions.jl")

# SETUP
outputs_path = "/Users/kingasiuta/Desktop/results/results_27_06"
vars = ["p0"] # what has to be constant

# reading possible parameter values
params_to_plot = CSV.read("$outputs_path/sim_params_comb_df.csv", DataFrame)

# file_list = glob("simul_res_parset_id_*.json", outputs_path)
#     for id in params_to_plot[:,:parset_id]
#         if string(outputs_path,"/simul_res_parset_id_",id,".json") ∉ file_list
#             filter!(r -> !in(r.parset_id, [id]), params_to_plot)
#         end
#     end
# CSV.write("$outputs_path/sim_params_comb_df.csv", params_to_plot)

values_to_loop = Dict()
for var in vars
    values_to_loop[var] = unique(params_to_plot[:,var])
end

# what to get
total_indic = "N_agents"
response_var = "total_infected"
variable = "TTC_car_freq"
div = 60

#what to plot
points_in_time = Int(3*24*3600)

df = Dict()
full_output = Dict()
for p0 in values_to_loop["p0"]
    full_output[p0], params_sets = read_results(outputs_path,
        control_vars = vars,
        control_val_equal = [p0],
        vars_to_check = ["TTC_car_freq"])

    println("$p0 read to dict")

    for (key, value) in full_output[p0]
        symb = value[variable]/div
        if !isnothing(total_indic) && length(value[total_indic])>1
            value[total_indic] = value[total_indic][1]
            if total_indic == "prct_of_agents_used_TTC"
                value[total_indic] = value[total_indic]*value["N_agents"]/100
            end
        end
        if p0 ∉ keys(df)
            df[p0] = DataFrame()
        end
        if response_var in ["total_infected", "TTC_infected", "street_infected"]
            df[p0][:,Symbol(symb)] = 100value[response_var][1]./(isnothing(total_indic) ? 1 : value[total_indic])
        else
            df[p0][:,Symbol(symb)] = value[response_var]./(isnothing(total_indic) ? 1 : value[total_indic])
        end
    end
    println("$p0 read to df")
end

out_df = DataFrame()
for p0 in values_to_loop["p0"]
    results = Dict()
    results[response_var] = df[p0][points_in_time, :]
    results[response_var] = DataFrame([[names(results[response_var])]; collect.(eachrow(DataFrame(results[response_var])))], [Symbol(variable); Symbol.("ITER_", points_in_time)])
    results[response_var][!,variable] = parse.(Float64, results[response_var][!,variable])
    sort!(results[response_var], Symbol(variable))

    min_row = DataFrame(sort!(results[response_var], Symbol("ITER_", points_in_time), rev=false)[1,:])
    p = DataFrame(p = p0)
    min_row = hcat(min_row,p)
    append!(out_df, min_row)
end

sort!(out_df, :p)
#out_df2 = out_df[1:9, :]
pyplot()
Plots.plot(out_df[!, :p], out_df[!, :TTC_car_freq], linewidth = 2.5,
legend = false,
linecolor = "purple")

xlabel!("base probability of infection")
ylabel!("optimal TTC frequency")
Plots.savefig("EOD_3_mins.eps")

max_stars_repo_path: codes/optimal freq plot.jl

max_stars_repo_name: NykPol/EpidemicInUrbanNetworkToronto

max_issues_repo_path: codes/optimal freq plot.jl

max_issues_repo_name: NykPol/EpidemicInUrbanNetworkToronto

max_forks_repo_path: codes/optimal freq plot.jl

max_forks_repo_name: NykPol/EpidemicInUrbanNetworkToronto


text:
import Theta: Lattice, svp, check_size_reduce, size_reduce!, check_lll, lll!

@testset "SVP correctness" begin
# Test data are generated using fplll
    @test norm(svp([116 303;331 963])[1]) == norm([-45, -30])
    @test norm(svp([116 456 99 566; 331 225 649 213; 303 827 395 694; 963 381 975 254])[1]) == norm([110 -12 -133 -127])
    @test norm(svp([116 827 566 597 912 611; 331 381 213 980 451 734; 303 99 694 462 926 668; 963 649 254 533 105 308; 456 395 629 231 162 101; 225 975 303 443 742 863])[1]) == norm([-301 283 -258 203 -61 121])
    @test norm(svp([116 99 629 912 668 438 621 371; 331 649 303 451 308 456 215 554; 303 395 597 926 101 887 6 434; 963 975 980 105 863 466 340 638; 456 566 462 162 279 634 142 939; 225 213 533 742 976 748 411 554; 827 694 231 611 127 307 250 840; 381 254 443 734 549 1002 873 633])[1]) == norm([-17 318 92 12 110 -12 -133 -127])
end

@testset "SVP length" begin
    @test begin
        S = svp(rand(1,1));
        (norm(S[1]))^2 ≈ S[2]
    end
    @test begin
        S = svp(rand(2,2));
        norm(S[1])^2 ≈ S[2]
    end
    @test begin
        S = svp(rand(5,5));
        norm(S[1])^2 ≈ S[2]
    end
    @test begin
        S = svp(rand(4,3));
        norm(S[1])^2 ≈ S[2]
    end
end

@testset "Size reduction" begin
    @test check_size_reduce(size_reduce!(Lattice(rand(1,1))))
    @test check_size_reduce(size_reduce!(Lattice(rand(2,2))))
    @test check_size_reduce(size_reduce!(Lattice(rand(3,3))))
    @test check_size_reduce(size_reduce!(Lattice(rand(4,4))))
    @test check_size_reduce(size_reduce!(Lattice(rand(5,5))))
    @test check_size_reduce(size_reduce!(Lattice(rand(10,10))))
end

@testset "LLL" begin
    @test check_lll(lll!(Lattice(rand(1,1))))
    @test check_lll(lll!(Lattice(rand(2,2))))
    @test check_lll(lll!(Lattice(rand(3,3))))
    @test check_lll(lll!(Lattice(rand(4,4))))
    @test check_lll(lll!(Lattice(rand(5,5))))
    @test check_lll(lll!(Lattice(rand(10,10))))
end


text:
# Copyright (c) 2021 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

export rm_eps!
# follows M. MOHRI GENERIC ∊-REMOVAL AND INPUT∊-NORMALIZATION ALGORITHMS FOR WEIGHTED TRANSDUCERS
# many optimizations are currently missing
"""
  `rm_eps!(fst::WFST)`

Same as `rm_eps` but operates in place.
"""
function rm_eps!(fst::WFST{W,A,D}) where {W,A,D}
  iseps_arc = arc -> iseps(get_ilabel(arc)) && iseps(get_olabel(arc))
  for (p,s) in enumerate(get_states(fst))
    arcs_to_rm = D[]
    new_arcs = Tuple{D,Arc{W,D}}[]
    for (i,arc) in enumerate(s)
      if iseps_arc(arc)
        push!(arcs_to_rm,i)
      end
      # get closure on Aϵ
      C = shortest_distance(fst,p; filter=iseps_arc)
      # TODO this allocate stuff... would make sense that C is also sparse?
      for (q,w) in enumerate(C)
        if w != zero(W)
          for a in fst[q]
            if !iseps_arc(a)
              new_arc=Arc{W,D}(
                               get_ilabel(a),
                               get_olabel(a),
                               get_weight(a) * w,
                               get_nextstate(a)
                              )
              push!(new_arcs, (p,new_arc))
            end
          end
          if isfinal(fst,q)
            if !isfinal(fst,p)
              final!(fst,p,zero(W))
            end
            fst.finals[p] += w * fst.finals[q]
          end
        end
      end
    end
    deleteat!(s,arcs_to_rm)
    for (p,a) in new_arcs
      push_arc!(fst,p,a)
    end
    unique!(fst[p]) # avoids multiple arcs with same i/o
  end
  connect!(fst)
  return fst
end

export rm_eps
"""
  `rm_eps(fst::WFST)`

Returns an equivalent WFST where arcs with input and output labels are removed.
"""
function rm_eps(fst::WFST)
  new_fst = deepcopy(fst)
  rm_eps!(new_fst)
end


text:
"""
    gather(x::PencilArray, [root::Integer=0])

Gather data from all MPI processes into one (big) array.

Data is received by the `root` process.

Returns the full array on the `root` process, and `nothing` on the other
processes.

This can be useful for testing, but it shouldn't be used with very large
datasets!
"""
function gather(x::PencilArray{T,N}, root::Integer=0) where {T, N}
    timer = Pencils.timer(pencil(x))

    @timeit_debug timer "gather" begin

    # TODO reduce allocations! see `transpose_impl!`
    comm = get_comm(x)
    rank = MPI.Comm_rank(comm)
    mpi_tag = 42
    pen = pencil(x)
    extra_dims = PencilArrays.extra_dims(x)

    # Each process sends its data to the root process.
    # If the local indices are permuted, the permutation is reverted before
    # sending the data.
    data = let perm = permutation(pen)
        if isidentity(perm)
            x.data
        else
            # Apply inverse permutation.
            p = append(inv(perm), Val(length(extra_dims)))
            permutedims(x.data, Tuple(p))  # creates copy!
        end
    end

    if rank != root
        # Wait for data to be sent, then return.
        # NOTE: When `data` is a ReshapedArray, I can't pass it directly to
        # MPI.Isend, because Base.cconvert(MPIPtr, ::ReshapedArray) is not
        # defined.
        # (Maybe it works in the current master of MPI.jl?)
        buf = data isa Base.ReshapedArray ? parent(data) : data
        send_req = MPI.Isend(buf, root, mpi_tag, comm)
        MPI.Wait!(send_req)
        return nothing
    end

    # Receive data (root only).
    topo = pen.topology
    Nproc = length(topo)
    recv = Vector{Array{T,N}}(undef, Nproc)
    recv_req = Vector{MPI.Request}(undef, Nproc)

    root_index = -1

    for n = 1:Nproc
        # Global data range that I will receive from process n.
        rrange = pen.axes_all[n]
        rdims = length.(rrange)

        src_rank = topo.ranks[n]  # actual rank of sending process
        if src_rank == root
            root_index = n
            recv_req[n] = MPI.REQUEST_NULL
        else
            # TODO avoid allocation?
            recv[n] = Array{T,N}(undef, rdims..., extra_dims...)
            recv_req[n] = MPI.Irecv!(recv[n], src_rank, mpi_tag, comm)
        end
    end

    # Unpack data.
    dest = Array{T,N}(undef, size_global(x))

    # Copy local data.
    colons_extra_dims = ntuple(n -> Colon(), Val(length(extra_dims)))
    dest[pen.axes_local..., colons_extra_dims...] .= data

    # Copy remote data.
    for m = 2:Nproc
        n, status = MPI.Waitany!(recv_req)
        rrange = pen.axes_all[n]
        dest[rrange..., colons_extra_dims...] .= recv[n]
    end

    end  # @timeit_debug

    dest
end