text:
## Framework for setting up an iterative problem for finding a zero
## TODO
## * a graphic of trace when verbose=true?

#
# In McNamee & Pan (DOI:10.1016/j.camwa.2011.11.015 there are a number of
# results on efficiencies of a solution, (1/d) log_10(q)
# Those implemented here are:
# quadratic cut (Muller) .265 (in a42)
# Newton() newton = .1505   or 1/2log(2)
# Order1() secant method .20 (1/1 * log(1.6)
# FalsePostion(12) Anderson-Bjork [.226, .233]
# FalsePostion(3) (King?) .264
# A42() 0.191 but convergence guaranteed
# Order8() 8th order 4 steps: .225 (log10(8)/4
# Order16() 16th order 5 steps .240
# Order5(): 5th order, 4 steps. 0.1747

# A zero is found by specifying:
# the method to use <: AbstractUnivariateZeroMethod
# the function(s)
# the initial state through a value for x either x, [a,b], or (a,b) <: AbstractUnivariateZeroState
# the options (e.g., tolerances) <: UnivariateZeroOptions

# The minimal amount needed to add a method, is to define a Method and an update_state method.

### Methods
abstract type AbstractUnivariateZeroMethod end
abstract type AbstractBracketing <: AbstractUnivariateZeroMethod end
abstract type AbstractNonBracketing <: AbstractUnivariateZeroMethod end
abstract type AbstractSecant <: AbstractNonBracketing end

# indicate if we expect f() to return one or multiple values (e.g. Newton)
fn_argout(::AbstractUnivariateZeroMethod) = 1

### States
abstract type  AbstractUnivariateZeroState{T,S} end

struct UnivariateZeroState{T,S} <: AbstractUnivariateZeroState{T,S}
    xn1::T
    xn0::T
    fxn1::S
    fxn0::S
end

# init_state(M, F, state) -- convert
# basic idea to convert from N to M:
# Fₘ = copy(M, fₙ)
# stateₘ = init_state(M, stateₙ, Fₘ)
function init_state(M::AbstractUnivariateZeroMethod, state::AbstractUnivariateZeroState, F)
    init_state(M, F, state.xn0, state.xn1, state.fxn0, state.fxn1)
end

# init_state(M,F,x) --> call init_state(M,F,x₀,x₁,fx₀, fx₁)
function init_state(M::AbstractUnivariateZeroMethod, F, x)
    error("no default method")
end

# initialize from xs, fxs
function init_state(::AbstractUnivariateZeroMethod, F, x₀, x₁, fx₀, fx₁)
    error("no default method")
end

# init_state(M, F, x; kwargs...)
# init_state(M, F x₀,x₁,fx₀,fx₁; kwargs...)
# init_state(M, state, F)
#
# A state holds at a minimum:
#
# * the values xₙ₋₁, xₙ and f(xₙ₋₁), f(xₙ) along with
# * some method-specific values
#
# A state is initialized with `init_state(M, F, x)` which sets up xₙ₋₁, xₙ, f(xₙ₋₁), f(xₙ)
# which then calls `init_state(M, F, xₙ₋₁, xₙ, f(xₙ₋₁), f(xₙ))` to finish the initialization
# to change to a new state use `init_state(M, state, F)`
#

# how many function evaluations in init_state
# this is a worst case estimate leading to the function call count being an upper bound only
initial_fncalls(M::AbstractUnivariateZeroState) = @warn "initial_fncalls fix $M"

## -----
### Options
struct UnivariateZeroOptions{Q,R,S,T}
    xabstol::Q
    xreltol::R
    abstol::S
    reltol::T
    maxevals::Int
    maxfnevals::Int
    strict::Bool
end

"""
    default_tolerances(M::AbstractUnivariateZeroMethod, [T], [S])

The default tolerances for most methods are `xatol=eps(T)`,
`xrtol=eps(T)`, `atol=4eps(S)`, and `rtol=4eps(S)`, with the proper
units (absolute tolerances have the units of `x` and `f(x)`; relative
tolerances are unitless). For `Complex{T}` values, `T` is used.

The number of iterations is limited by `maxevals=40`.

"""
default_tolerances(M::AbstractUnivariateZeroMethod) = default_tolerances(M, Float64, Float64)
function default_tolerances(::AbstractUnivariateZeroMethod, ::Type{T}, ::Type{S}) where {T, S}
    xatol = eps(real(T)) * oneunit(real(T))
    xrtol = eps(real(T))  # unitless
    atol = 4 * eps(real(float(S))) * oneunit(real(S))
    rtol = 4 * eps(real(float(S))) * one(real(S))
    maxevals = 40
    maxfnevals = typemax(Int)
    strict = false
    (xatol, xrtol, atol, rtol, maxevals, maxfnevals, strict)
end

init_options(M::AbstractUnivariateZeroMethod,
             state::AbstractUnivariateZeroState{T,S};
             kwargs...
             ) where {T, S} = init_options(M, T, S; kwargs...)

function init_options(M, T=Float64, S=Float64; kwargs...)
    d = kwargs

    defs = default_tolerances(M, T, S)
    options = UnivariateZeroOptions(get(d, :xatol, get(d, :xabstol, defs[1])),
                                    get(d, :xrtol, get(d, :xreltol, defs[2])),
                                    get(d, :atol,  get(d, :abstol,  defs[3])),
                                    get(d, :rtol,  get(d, :reltol,  defs[4])),
                                    get(d, :maxevals,   get(d, :maxsteps, defs[5])),
                                    get(d, :maxfnevals, defs[6]),
                                    get(d, :strict,     defs[7]))
    if haskey(d, :maxfnevals)
        @warn("maxfnevals is ignored. See the test for an example to implement this featrue")
    end
    options
end

# # reset options to default values
@deprecate init_options!(options, M) init_options(M)

## Tracks (for logging actual steps)
## when no logging this should get optimized out
abstract type AbstractTracks end
struct NullTracks <: AbstractTracks end
# api
log_step(s::NullTracks, M, x, init=false) = nothing
log_steps(::NullTracks, n=1) = nothing
incfn(::NullTracks, i=1) = nothing
log_message(::NullTracks, msg) = nothing
log_convergence(::NullTracks, msg) = nothing
log_state(::NullTracks, state) = nothing
log_method(::NullTracks, method) = nothing
log_nmethod(::NullTracks, method) = nothing

# a tracks object to record tracks
"""
    Tracks(T, S)

Construct a `Tracks` object used to record the progress of the algorithm.

By default, a null tracks object used, but if `verbose=true` is specified of a `Tracks` object passed to the keyword `tracks` then the steps of the algorithm are recorded in the `Tracks` object. If `verbose=true`, the `Tracks` object will be shown.
"""
mutable struct Tracks{T,S} <: AbstractTracks
    xs::Vector{T}
    fs::Vector{S}
    steps::Int
    fncalls::Int
    convergence_flag::Symbol
    message::String
    state
    method
    nmethod
end
Tracks(T,S) = Tracks(T[],S[], 0, 0, :not_converged, "", nothing, nothing,nothing)
Tracks(s::AbstractUnivariateZeroState{T,S}) where {T, S} = Tracks(T,S)
Tracks(verbose, tracks, state) = (verbose && isa(tracks, NullTracks)) ? Tracks(state) : tracks

function log_step(l::Tracks, M::Any, state, init=nothing)
    if init != nothing
        x₀, fx₀ = state.xn0, state.fxn0
        push!(l.xs, x₀)
        push!(l.fs, fx₀)
    end
    x₁, fx₁ = state.xn1, state.fxn1
    push!(l.xs, x₁)
    push!(l.fs, fx₁)

    init == nothing && log_steps(l, 1)
    nothing
end

incfn(l::Tracks, i=1) = (l.fncalls += i; nothing)
log_steps(l::Tracks, n=1) = (l.steps += n; nothing)
log_message(l::Tracks, msg) = (l.message *= msg; nothing)
log_convergence(l::Tracks, msg) = (l.convergence_flag = msg; nothing)
log_state(l::Tracks, state) = (l.state = state; nothing)
log_method(l::Tracks, method) = (l.method = method; nothing)
log_nmethod(l::Tracks, method) = (l.nmethod = method; nothing)

Base.show(io::IO, l::Tracks) = show_trace(io, l.method, l.nmethod, l.state, l)

function show_tracks(io::IO, s::Tracks, M::AbstractUnivariateZeroMethod)
    for (i, (xi, fxi)) in enumerate(zip(s.xs, s.fs))
        println(io,@sprintf("%s = % 18.16f,\t %s = % 18.16f", "x_$(i-1)", float(xi), "fx_$(i-1)", float(fxi)))
    end
    println(io,"")
end

function show_trace(io::IO, method, N, state, tracks)

    if state == nothing
        print(io, "Algorithm has not been run")
        return nothing
    end

    converged = !isnan(state.xn1)
    println(io, "Results of univariate zero finding:\n")
    if converged
        println(io, "* Converged to: $(state.xn1)")
        if N == nothing || isa(method, AbstractBracketing)
            println(io, "* Algorithm: $(method)")
        else
            println(io, "* Algorithm: $(method), with possible bracketing with $N")
        end
        println(io, "* iterations: $(tracks.steps)")
        println("* function evaluations ≈ $(tracks.fncalls)")
        tracks.convergence_flag == :x_converged && println(io, "* stopped as x_n ≈ x_{n-1} using atol=xatol, rtol=xrtol")
        tracks.convergence_flag == :f_converged && tracks.message == "" && println(io, "* stopped as |f(x_n)| ≤ max(δ, |x|⋅ϵ) using δ = atol, ϵ = rtol")
        tracks.message != "" && println(io, "* Note: $(tracks.message)")
    else
        println(io, "* Convergence failed: $(tracks.message)")
        println(io, "* Algorithm $(method)")
    end
    println(io, "")
    println(io, "Trace:")
    show_tracks(io, tracks, method)
end

### Functions
# A hacky means to call a function so that parameters can be passed as desired
# and the correct number of outputs are computed
struct Callable_Function{Single, Tup, F, P}
    f::F
    p::P
    function Callable_Function(M, f, p=nothing)
        Single = Val{fn_argout(M)}
        Tup = Val{isa(f, Tuple)}
        F = typeof(f)
        P = typeof(p)
        new{Single, Tup, F, P}(f, p)
    end
end

function Callable_Function(M,F::Callable_Function, p=F.p)
    Callable_Function(M, F.f, p)
end

# return f(x); (f(x), f(x)/f'(x)); *or* f(x), (f(x)/f'(x), f'(x)/f''(x), ...) # so N=1, 2 are special cased
# Callable_Function(output_arity, input_arity, F, p)
# First handle: x -> (f,f/f', f'/f'', ...)
(F::Callable_Function{S,T,𝑭,P})(x) where {S <: Val{1}, T <: Val{false}, 𝑭, P<:Nothing} =
    first(F.f(x))
(F::Callable_Function{S,T,𝑭,P})(x) where {S <: Val{1}, T <: Val{false}, 𝑭, P} =
    first(F.f(x, F.p))

(F::Callable_Function{S,T,𝑭,P})(x) where {N, S <: Val{2}, T <: Val{false}, 𝑭, P<:Nothing} =
    F.f(x)[1:2]
(F::Callable_Function{S,T,𝑭,P})(x) where {N, S <: Val{2}, T <: Val{false}, 𝑭, P} =
    F.f(x, F.p)[1:2]

# N ≥ 3 returns (f, (...))
(F::Callable_Function{S,T,𝑭,P})(x) where {N, S <: Val{N}, T <: Val{false}, 𝑭, P<:Nothing} = begin
    fs = F.f(x)
    fs[1], ntuple(i->fs[i+1], Val(N-1))
end
(F::Callable_Function{S,T,𝑭,P})(x) where {N, S <: Val{N}, T <: Val{false}, 𝑭, P} = begin
    fs = F.f(x, F.p)
    fs[1], ntuple(i->fs[i+1], Val(N-1))
end

## f is specified as a tuple (f,f',f'', ...)
## N =1  return f(x)
(F::Callable_Function{S,T,𝑭,P})(x) where {S <: Val{1}, T <: Val{true}, 𝑭, P<:Nothing} =
    first(F.f)(x)
(F::Callable_Function{S,T,𝑭,P})(x) where {S <: Val{1}, T <: Val{true}, 𝑭, P} =
    first(F.f)(x, F.p)
## N=2 return (f(x), f(x)/f'(x))
(F::Callable_Function{S,T,𝑭,P})(x) where {S <: Val{2}, T <: Val{true}, 𝑭, P<:Nothing} = begin
    f, f′ = (F.f[1])(x), (F.f[2])(x)
    f, f/f′
end
(F::Callable_Function{S,T,𝑭,P})(x) where {S <: Val{2}, T <: Val{true}, 𝑭, P} = begin
    f, f′ = (F.f[1])(x, F.p), (F.f[2])(x, F.p)
    f, f/f′
end

## For N ≥ 3 we return (f, (f/f', f'/f'', ...);
## Pay no attention to this code; we hand write a bunch, as the
## general formula later runs more slowly.
(F::Callable_Function{S,T,𝑭,P})(x) where {S <: Val{3}, T <: Val{true}, 𝑭, P<:Nothing} = begin
    f, f′, f′′ = (F.f[1])(x), (F.f[2])(x), (F.f[3])(x)
    f, (f/f′, f′/f′′)
end
(F::Callable_Function{S,T,𝑭,P})(x) where {S <: Val{3}, T <: Val{true}, 𝑭, P} = begin
    f, f′, f′′ = (F.f[1])(x, F.p), (F.f[2])(x, F.p), (F.f[3])(x, F.p)
    f, (f/f′, f′/f′′)
end

(F::Callable_Function{S,T,𝑭,P})(x) where {S <: Val{4}, T <: Val{true}, 𝑭, P<:Nothing} = begin
    f, f′, f′′, f′′′ = (F.f[1])(x), (F.f[2])(x), (F.f[3])(x), (F.f[4])(x)
    f, (f/f′, f′/f′′, f′′/f′′′)
end
(F::Callable_Function{S,T,𝑭,P})(x) where {S <: Val{4}, T <: Val{true}, 𝑭, P} = begin
    f, f′,f′′,f′′′ = (F.f[1])(x, F.p), (F.f[2])(x, F.p), (F.f[3])(x, F.p), (F.f[4])(x, F.p)
    𝐓 = eltype(f/f′)
    f, NTuple{3,𝐓}((f/f′, f′/f′′, f′′/f′′′))
end

(F::Callable_Function{S,T,𝑭,P})(x) where {S <: Val{5}, T <: Val{true}, 𝑭, P<:Nothing} = begin
    f, f′, f′′, f′′′, f′′′′ = (F.f[1])(x), (F.f[2])(x), (F.f[3])(x), (F.f[4])(x), (F.f[5])(x)
    f, (f/f′, f′/f′′, f′′/f′′′, f′′′/f′′′′)
end
(F::Callable_Function{S,T,𝑭,P})(x) where {S <: Val{5}, T <: Val{true}, 𝑭, P} = begin
    f, f′,f′′,f′′′,f′′′′ = (F.f[1])(x, F.p), (F.f[2])(x, F.p), (F.f[3])(x, F.p), (F.f[4])(x, F.p), (F.f[5])(x, F.p)
    f, (f/f′, f′/f′′, f′′/f′′′, f′′′/f′′′′)
end

(F::Callable_Function{S,T,𝑭,P})(x) where {S <: Val{6}, T <: Val{true}, 𝑭, P<:Nothing} = begin
    f, f′, f′′, f′′′, f′′′′, f′′′′′ = (F.f[1])(x), (F.f[2])(x), (F.f[3])(x), (F.f[4])(x), (F.f[5])(x),(F.f[6])(x)
    f, (f/f′, f′/f′′, f′′/f′′′, f′′′/f′′′′, f′′′′/f′′′′′)
end
(F::Callable_Function{S,T,𝑭,P})(x) where {S <: Val{6}, T <: Val{true}, 𝑭, P} = begin
    f, f′,f′′,f′′′,f′′′′, f′′′′′ = (F.f[1])(x, F.p), (F.f[2])(x, F.p), (F.f[3])(x, F.p), (F.f[4])(x, F.p), (F.f[5])(x, F.p), (F.f[6])(x, F.p)
    f, (f/f′, f′/f′′, f′′/f′′′, f′′′/f′′′′, f′′′′/f′′′′′)
end

# faster with the above written out, should generate them...
(F::Callable_Function{S,T,𝑭,P})(x) where {𝐍, S <: Val{𝐍}, T <: Val{true}, 𝑭, P<:Nothing} = begin
    fs = ntuple(i -> F.f[i](x), Val(𝐍))
    first(fs), ntuple(i -> fs[i]/fs[i+1], Val(𝐍-1))
end

(F::Callable_Function{S,T,𝑭,P})(x) where {𝐍, S <: Val{𝐍}, T <: Val{true}, 𝑭, P} = begin
    fs = ntuple(i -> F.f[i](x, p), Val(𝐍))
    first(fs), ntuple(i -> fs[i]/fs[i+1], Val(𝐍-1))
end

## Assess convergence
@inline function _is_f_approx_0(fa, a, atol, rtol, relaxed::Any)
    aa, afa = abs(a), abs(fa)
    tol = max(_unitless(atol), _unitless(aa) * rtol)
    tol = cbrt(abs(_unitless(tol)))  # relax test
    afa <= tol * oneunit(afa)
end
@inline function _is_f_approx_0(fa, a, atol, rtol)
    aa, afa = abs(a), abs(fa)
    tol = max(_unitless(atol), _unitless(aa) * rtol)
    afa <= tol * oneunit(afa)
end

"""
    Roots.assess_convergence(method, state, options)

Assess if algorithm has converged.

Return a convergence flag and a Boolean indicating if algorithm has terminated (converged or not converged)

If algrithm hasn't converged returns `(:not_converged, false)`.

If algorithm has stopped or converged, return flag and `true`. Flags are:

* `:x_converged` if `abs(xn1 - xn0) < max(xatol, max(abs(xn1), abs(xn0)) * xrtol)`

* `:f_converged` if  `|f(xn1)| < max(atol, |xn1|*rtol)`

* `:nan`, `:inf` if xn1 or fxn1 is `NaN` or an infinity

* `:not_converged` if algorithm should continue

Does not check number of steps taken nor number of function evaluations.

In `decide_convergence`, stopped values (and `:x_converged` when `strict=false`) are checked for convergence with a relaxed tolerance.

"""
function assess_convergence(method::Any, state::AbstractUnivariateZeroState{T,S}, options) where {T,S}

    # return convergence_flag, boolean

    xn0, xn1 = state.xn0, state.xn1
    fxn1 = state.fxn1

    if isnan(xn1) || isnan(fxn1)
        return (:nan, true)
    end

    if isinf(xn1) || isinf(fxn1)
        return (:inf, true)
    end

    # f(xstar) ≈ xstar * f'(xstar)*eps(), so we pass in lambda
    if _is_f_approx_0(fxn1, xn1, options.abstol, options.reltol)
        return (:f_converged, true)
    end

    # stop when xn1 ~ xn.
    # in find_zeros there is a check that f could be a zero with a relaxed tolerance
    if isapprox(xn1, xn0, atol=options.xabstol, rtol=options.xreltol)
        return (:x_converged, true)
    end

    return (:not_converged, false)
end

# state has stopped, this identifies if it has converged
"""
    decice_convergence(M,F,state,options, convergence_flag)

When the algorithm terminates, this function decides the stopped value or returns NaN
"""
function decide_convergence(M::AbstractUnivariateZeroMethod,  F, state::AbstractUnivariateZeroState{T,S}, options, val) where {T,S}

    xn0, xn1 = state.xn0, state.xn1
    fxn1 = state.fxn1

    val ∈ (:f_converged, :exact_zero, :converged) && return xn1

    ## stopping is a heuristic, x_converged can mask issues
    ## if strict=true or the tolerance for f is 0 this will return xn1 if x_converged
    ## if strict == false, this will also check f(xn) ~ - with a relaxed
    ## tolerance
    if options.strict || (iszero(options.abstol) && iszero(options.reltol))
        val == :x_converged && return xn1
        _is_f_approx_0(fxn1, xn1, options.abstol, options.reltol) && return xn1
    else
        if val ∈ (:x_converged, :not_converged)
            # |fxn| <= relaxed

            _is_f_approx_0(fxn1, xn1, options.abstol, options.reltol, true) && return xn1
        end
    end

    # if (state.stopped || state.x_converged) && !(state.f_converged)
    #     ## stopped is a heuristic, x_converged can mask issues
    #     ## if strict == false, this will also check f(xn) ~ - with a relaxed
    #     ## tolerance

    #     ## are we at a crossing values?
    #     ## seems worth a check for 2 fn evals.
    #     if T <: Real && S <: Real
    #         for u in (prevfloat(xn1), nextfloat(xn1))
    #             fu = first(F(u))
    #             incfn(state)
    #             if iszero(fu) || _unitless(fu * fxn1) < 0
    #                 state.message *= "Change of sign at xn identified. "
    #                 state.f_converged = true
    #             end
    #         end
    #     end

    #     δ = maximum(_unitless, (options.abstol, options.reltol))
    #     if options.strict || iszero(δ)
    #         if state.x_converged
    #             state.f_converged = true
    #         else
    #             state.convergence_failed = true
    #         end

    #     else
    #         xstar, fxstar = state.xn1, state.fxn1
    #         if _is_f_approx_0(fxstar, xstar, options.abstol, options.reltol, :relaxed)
    #             state.xstar, state.fxstar = xstar, fxstar
    #             msg = "Algorithm stopped early, but |f(xn)| < ϵ^(1/3), where ϵ depends on xn, rtol, and atol. "
    #             state.message = state.message == "" ? msg : state.message * "\n\t" * msg
    #             state.f_converged = true
    #         else
    #             state.convergence_failed = true
    #         end
    #     end
    # end

    # if !state.f_converged
    #     state.xstar, state.fxstar = NaN*xn1, NaN*fxn1
    # end

    NaN * xn1

end

## ----

"""

    find_zero(f, x0, M, [N::AbstractBracketing]; kwargs...)

Interface to one of several methods for finding zeros of a univariate function, e.g. solving ``f(x)=0``.

# Initial starting value

For most methods, `x0` is a scalar value indicating the initial value
in the iterative procedure. (Secant methods can have a tuple specify
their initial values.) Values must be a subtype of `Number` and have
methods for `float`, `real`, and `oneunit` defined.

For bracketing intervals, `x0` is specified using a tuple, a vector, or any iterable with `extrema` defined. A bracketing interval, ``[a,b]``, is one where f(a) and f(b) have different signs.

# Return value

If the algorithm suceeds, the approximate root identified is returned. A `ConvergenceFailed` error is thrown if the algorithm fails. The alternate form `solve(ZeroProblem(f,x0), M)` returns `NaN` in case of failure.

# Specifying a method

A method is specified to indicate which algorithm to employ:

* There are methods for bisection where a bracket is specified: [`Bisection`](@ref), [`Roots.A42`](@ref), [`Roots.AlefeldPotraShi`](@ref), [`Roots.Brent`](@ref), and [`FalsePosition`](@ref)

* There are several derivative-free methods: cf. [`Order0`](@ref), [`Order1`](@ref) (also [`Roots.Secant`](@ref)), [`Order2`](@ref) (also [`Roots.Steffensen`](@ref)), [`Order5`](@ref), [`Order8`](@ref), and [`Order16`](@ref), where the number indicates the order of the convergence. Methods [`Roots.Order1B`](@ref) and [`Roots.Order2B`](@ref) are useful when the desired zero has a multiplicity.

* There are some classical methods where derivatives are required: [`Roots.Newton`](@ref), [`Roots.Halley`](@ref), [`Roots.Schroder`](@ref).

* The family [`Roots.LithBoonkkampIJzerman{S,D}`](@ref) for different `S` and `D` uses a linear multistep method root finder. The `(2,0)` method is the secant method, `(1,1)` is Newton's method.

For more detail, see the help page for each method (e.g., `?Order1`). Non-exported methods must be qualified with module name, as in `?Roots.Schroder`.

If no method is specified, the default method depends on `x0`:

* If `x0` is a scalar, the default is the slower, but more robust `Order0` method.

* If `x0` is a tuple, vector, or iterable with `extrema` defined indicating a *bracketing* interval, the `Bisection` method is used. (The exact algorithm depends on the number type and the tolerances.)

# Specifying the function

The function(s) are passed as the first argument.

For the few methods that use one or more derivatives (`Newton`, `Halley`,
`Schroder`, `LithBoonkkampIJzerman(S,D)`, and  `Order5Derivative`) a
tuple of functions is used. For the classical algorithms, a function returning `(f(x), f(x)/f'(x), [f'(x)/f''(x)])` may be used.

# Optional arguments (tolerances, limit evaluations, tracing)

* `xatol` - absolute tolerance for `x` values. Passed to `isapprox(x_n, x_{n-1})`.
* `xrtol` - relative tolerance for `x` values. Passed to `isapprox(x_n, x_{n-1})`.
* `atol`  - absolute tolerance for `f(x)` values.
* `rtol`  - relative tolerance for `f(x)` values.
* `maxevals`   - limit on maximum number of iterations.
* `strict` - if `false` (the default), when the algorithm stops, possible zeros are checked with a relaxed tolerance.
* `verbose` - if `true` a trace of the algorithm will be shown on successful completion. See the internal `Tracks` object to save this trace.

See the help string for `Roots.assess_convergence` for details on
convergence. See the help page for `Roots.default_tolerances(method)`
for details on the default tolerances.

In general, with floating point numbers, convergence must be
understood as not an absolute statement. Even if mathematically `α` is
an answer and `xstar` the floating point realization, it may be that
`f(xstar) - f(α)  ≈ xstar ⋅  f'(α) ⋅ eps(α)`, so tolerances must be
appreciated, and at times specified.

For the `Bisection` methods, convergence is guaranteed, so the tolerances are set to be ``0`` by default.

If a bracketing method is passed in after the method specification,
then whenever a bracket is identified during the algorithm, the method
will switch to the bracketing method to identify the zero. (Bracketing
methods are mathematicall guaranteed to converge, non-bracketing methods may or may not converge.)
This is what `Order0` does by default, with an initial secant method switching
to the `AlefeldPotraShi` method should a bracket be encountered.

Note: The order of the method is hinted at in the naming scheme. A
scheme is order `r` if, with `eᵢ = xᵢ - α`, `eᵢ₊₁ = C⋅eᵢʳ`. If the
error `eᵢ` is small enough, then essentially the error will gain `r`
times as many leading zeros each step. However, if the error is not
small, this will not be the case. Without good initial guesses, a high
order method may still converge slowly, if at all. The `OrderN`
methods have some heuristics employed to ensure a wider range for
convergence at the cost of not faithfully implementing the method,
though those are available through unexported methods.

# Examples:

Default methods.

```jldoctest find_zero
julia> using Roots

julia> find_zero(sin, 3)  # use Order0()
3.141592653589793

julia> find_zero(sin, (3,4)) # use Bisection()
3.141592653589793
```

Specifying a method,

```jldoctest find_zero
julia> find_zero(sin, (3,4), Order1())            # can specify two starting points for secant method
3.141592653589793

julia> find_zero(sin, 3.0, Order2())              # Use Steffensen method
3.1415926535897936

julia> find_zero(sin, big(3.0), Order16())        # rapid convergence
3.141592653589793238462643383279502884197169399375105820974944592307816406286198

julia> find_zero(sin, (3, 4), Roots.A42())      # fewer function calls than Bisection(), in this case
3.141592653589793

julia> find_zero(sin, (3, 4), FalsePosition(8))   # 1 of 12 possible algorithms for false position
3.141592653589793

julia> find_zero((sin,cos), 3.0, Roots.Newton())  # use Newton's method
3.141592653589793

julia> find_zero((sin, cos, x->-sin(x)), 3.0, Roots.Halley())  # use Halley's method
3.141592653589793
```

Changing tolerances.

```jldoctest find_zero
julia> fn = x -> (2x*cos(x) + x^2 - 3)^10/(x^2 + 1);

julia> x0, xstar = 3.0,  2.9947567209477;

julia> fn(find_zero(fn, x0, Order2())) <= 1e-14  # f(xₙ) ≈ 0, but Δxₙ can be largish
true

julia> find_zero(fn, x0, Order2(), atol=0.0, rtol=0.0) # error: x_n ≉ x_{n-1}; just f(x_n) ≈ 0
ERROR: Roots.ConvergenceFailed("Algorithm failed to converge")
[...]

julia> fn = x -> (sin(x)*cos(x) - x^3 + 1)^9;

julia> x0, xstar = 1.0,  1.112243913023029;

julia> find_zero(fn, x0, Order2()) ≈ xstar
true

julia> find_zero(fn, x0, Order2(), maxevals=3)    # need more steps to converge
ERROR: Roots.ConvergenceFailed("Algorithm failed to converge")
[...]
```

# Tracing

Passing `verbose=true` will show details on the steps of the algorithm.
The `tracks` argument allows
the passing of storage to record the values of `x` and `f(x)` used in
the algorithm.

!!! note
    See [`solve!`](@ref) and [`ZeroProblem`](@ref) for an alternate interface.
"""
function find_zero(f, x0, M::AbstractUnivariateZeroMethod;
                   p = nothing,
                   verbose=false,
                   tracks::AbstractTracks=NullTracks(),
                   kwargs...)

    xstar = solve(ZeroProblem(f, x0), M;
                  p=p, verbose=verbose, tracks=tracks,
                  kwargs...)

    isnan(xstar) && throw(ConvergenceFailed("Algorithm failed to converge"))

    xstar

end

# defaults when method is not specified
# if a number, use Order0
# O/w use a bracketing method of an assumed iterable
find_zero(f, x0::T; kwargs...)  where {T <: Number} = find_zero(f, x0, Order0(); kwargs...)
find_zero(f, x0; kwargs...) = find_zero(f, x0, Bisection(); kwargs...)

"""
    find_zero(M, F, state, [options], [l])

Find zero using method `M`, function(s) `F`, and initial state
`state`. Returns an approximate zero or `NaN`. Useful when some part
of the processing pipeline is to be adjusted.

* `M::AbstractUnivariateZeroMethod` a method, such as `Secant()`
* `F`: A callable object (or tuple of callable objects for certain methods)
* `state`: An initial state, as created by `init_state` (or `_init_state`).
* `options::UnivariateZeroOptions`: specification of tolerances
* `l::AbstractTracks`: used to record steps in algorithm, when requested.
```

!!! note
    To be deprecated in favor of `solve!(init(...))`.
"""
function find_zero(M::AbstractUnivariateZeroMethod,
                   F,
                   state::AbstractUnivariateZeroState,
                   options::UnivariateZeroOptions=init_options(M, state),
                   l::AbstractTracks=NullTracks()
                   )
    solve!(init(M, F, state, options, l))
end

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

## Create an Iterator interface
# returns NaN, not an error, if there are issues

"""
    ZeroProblem{F,X}

A container for a function and initial guess to be used with `solve`.
"""
struct ZeroProblem{F,X}
    F::F
    x₀::X
end

## The actual iterating object
struct ZeroProblemIterator{M,F,S,O,L}
    M::M
    F::F
    state::S
    options::O
    logger::L
end

Base.show(io::IO, Z::ZeroProblemIterator) = print(io, "A problem object to pass to `solve!`")

## Initialize a Zero Problem Iterator
## init(Z,p)
## init(Z,M,p)
## init(M,F,state, [options], [logger])
## want p to be positional, not named⁺
function init(𝑭𝑿::ZeroProblem, M::AbstractUnivariateZeroMethod, p′ = nothing;
              p = nothing,
              verbose::Bool=false,
              tracks = NullTracks(),
              kwargs...)

    F = Callable_Function(M, 𝑭𝑿.F, p === nothing ? p′ : p)  #⁺
    state = init_state(M, F, 𝑭𝑿.x₀)
    options = init_options(M, state; kwargs...)
    l = Tracks(verbose, tracks, state)
    incfn(l, initial_fncalls(M))
    ZeroProblemIterator(M,F,state,options,l)

end

function init(𝑭𝑿::ZeroProblem, p′=nothing; kwargs...)
    M = length(𝑭𝑿.x₀) == 1 ? Secant() : Bisection()
    init(𝑭𝑿, M, p′; kwargs...)
end

function init(M::AbstractUnivariateZeroMethod, F,
              state::AbstractUnivariateZeroState,
              options::UnivariateZeroOptions=init_options(M, state),
              l::AbstractTracks=NullTracks())
    ZeroProblemIterator(M, Callable_Function(M,F), state, options, l)
end

# Optional iteration interface to handle looping
# errors on non-convergence, unlike solve!(init(...)) which returns NaN
# allocates...
# This should be deprecated
function Base.iterate(P::ZeroProblemIterator, st=nothing)
    ## st = (val, (state, ctr, val))
    if st == nothing
        state = P.state
        ctr = 1
    else
        (state, ctr, val) = st
        ctr += 1
    end
    M, options, l = P.M, P.options, P.logger
    val,stopped = :no_convergence, false
    while true
        val, stopped = assess_convergence(M, state, options)
        stopped && break
        if ctr > options.maxevals
            stopped = true
            break
        end
        state, stopped = update_state(M, P.F, state, options, l)
        log_step(l, M, state)
        break
    end

    if stopped
        xstar = decide_convergence(P.M, P.F, state, P.options, val)
        isnan(xstar) && throw(ConvergenceFailed("Stopped at $(state.xn1)"))
        nothing
    else
        return (state.xn1, (state, ctr, val))
    end
end

"""
    solve(fx::ZeroProblem, [p=nothing]; kwargs...)
    solve(fx::ZeroProblem, M, [p=nothing]; kwargs...)
    init(fx::ZeroProblem, [M], [p=nothing];
         verbose=false, tracks=NullTracks(), kwargs...)
    solve!(P::ZeroProblemIterator)

Solve for the zero of a function specified through a  `ZeroProblem` or `ZeroProblemIterator`

The methods involved with this interface are:

* `ZeroProblem`: used to specify a problem with a function (or functions) and an initial guess
* `solve`: to solve for a zero in a `ZeroProblem`

The latter calls the following, which can be useful independently:

* `init`: to initialize an iterator with a method for solution, any adjustments to the default tolerances, and a specification to log the steps or not.
* `solve!` to iterate to convergence.

Returns `NaN`, not an error, when the problem can not be solved. Tested for zero allocations.

## Examples:

```jldoctest find_zero
julia> using Roots

julia> fx = ZeroProblem(sin, 3)
ZeroProblem{typeof(sin), Int64}(sin, 3)

julia> solve(fx)
3.141592653589793
```

Or, if the iterable is required

```jldoctest find_zero
julia> problem = init(fx);

julia> solve!(problem)
3.141592653589793
```

The default method is `Order1()`, when  `x0` is a number, or `Bisection()` when `x0` is an iterable with 2 or more values.

A second position argument for `solve` or `init` is used to specify a different method; keyword arguments can be used to adjust the default tolerances.

```jldoctest find_zero
julia> solve(fx, Order5(), atol=1/100)
3.1425464815525403
```

The above is equivalent to:

```jldoctest find_zero
julia> problem = init(fx, Order5(), atol=1/100);

julia> solve!(problem)
3.1425464815525403
```

The  argument `p` may be used if the function(s) to be solved depend on a parameter in their second positional argument (e.g., `f(x,p)`). For example

```jldoctest find_zero
julia> f(x,p) = exp(-x) - p # to solve p = exp(-x)
f (generic function with 1 method)

julia> fx = ZeroProblem(f, 1)
ZeroProblem{typeof(f), Int64}(f, 1)

julia> solve(fx, 1/2)  # log(2)
0.6931471805599453
```

This would be recommended, as there is no recompilation due to the function changing.

The argument `verbose=true` for `init` instructs that steps to be logged;

The iterator interface allows for the creation of hybrid solutions, for example, this is essentially how `Order0` is constructed (`Order0` follows secant steps until a bracket is identified, after which is switches to a bracketing algorithm.)

```jldoctest find_zero
julia> function order0(f, x)
           fx = ZeroProblem(f, x)
           p = init(fx, Roots.Secant())
           xᵢ,st = ϕ = iterate(p)
           while ϕ != nothing
               xᵢ, st = ϕ
               state, ctr = st
               fᵢ₋₁, fᵢ = state.fxn0, state.fxn1
               if sign(fᵢ₋₁)*sign(fᵢ) < 0 # check for bracket
                   x0 = (state.xn0, state.xn1)
                   fx′ = ZeroProblem(f, x0)
                   p = init(fx′, Bisection())
                   xᵢ = solve!(p)
                   break
               end
               ϕ = iterate(p, st)
           end
           xᵢ
       end
order0 (generic function with 1 method)

julia> order0(sin, 3)
3.141592653589793
```

"""
function solve!(P::ZeroProblemIterator; verbose=false)

    M, F, state, options, l = P.M, P.F, P.state, P.options, P.logger

    val, stopped = :not_converged, false
    ctr = 1
    log_step(l, M, state, :init)

    while !stopped
        val, stopped = assess_convergence(M, state, options)
        stopped && break
        ctr > options.maxevals && break

        state, stopped = update_state(M, F, state, options, l)

        log_step(l, M, state)
        ctr += 1

    end

    if !(l isa NullTracks)
        log_convergence(l, val); log_state(l, state); log_method(l, M)
        verbose && display(l)
    end

    val, stopped = assess_convergence(M, state, options) # udpate val flag
    decide_convergence(M, F, state, options, val)

end

# thread verbose through
CommonSolve.solve(F::ZeroProblem, args...; verbose=false, kwargs...) = solve!(init(F, args...; verbose=verbose, kwargs...); verbose=verbose)

ext: jl


text:
using BlackBoxOptim
using DwsimOpt
include("c:\\Users\\lfsfr\\Desktop\\DwsimOpt.jl\\src\\SimOpt.jl")

op1 = OptProblemDef()

function fpen(x)
    return Float64(op1.f(x)[1]) .+ sum(max.(0, op1.g(x)))
end

bound = [Tuple(op1.searchSpace[:, i]') for i = 1:size(op1.searchSpace)[2]]

res = bboptimize(fpen; SearchRange = bound, NumDimensions = op1.dim, MaxTime = 10.0)


text:
using Documenter, PDKernels

makedocs(
    sitename="PDKernels Documentation",
    modules = [Kernels, Interpolation, Parameters],
    pages = [],
)


text:
using Transits
using Limbdark
using LaTeXStrings
using BenchmarkTools
using DataFrames
using Statistics
using StatsPlots, ColorSchemes
using ProgressLogging

N_pts = [100, 316, 1000, 3160, 10000, 31600, 100000]
N_us = [1 2 3 5 8 13 21 34 55 90 144]

function benchmark(N_u, N_pts; r=0.1)
    u = ones(N_u) ./ N_u
    b = [sqrt(((i - N_pts/2) * 2/N_pts * (1 + 2 * r))^2) for i in 1:N_pts]

    ld = PolynomialLimbDark(u)
    output = similar(b)
    bench = @benchmark begin
        @inbounds for i in eachindex($b)
            $output[i] = compute($ld, $b[i], $r)
        end
    end
    return median(bench).time * 1e-9
end
vals = zeros(length(N_pts), length(N_us))
@progress for (j, Nu) in enumerate(N_us)
    @progress for (i, N) in enumerate(N_pts)
        vals[i, j] = benchmark(Nu, N)
    end
end

labels = reduce(hcat, "uₙ: $N" for N in N_us)

plot(
    N_pts,
    vals,
    scale=:log10,
    ls=:solid,
    m=:o,
    palette=palette(:plasma, length(N_us) + 2),
    msw=0,
    lab=labels,
    xlabel="number of points",
    ylabel="timing [s]",
    leg=:topleft
)

function benchmark_limbdark(N_u, N_pts; r=0.1)
    u = ones(N_u) ./ N_u
    b = [sqrt(((i - N_pts/2) * 2/N_pts * (1 + 2 * r))^2) for i in 1:N_pts]

    trans = Limbdark.transit_init(r, b[1], u, true)
    output = similar(b)
    bench = @benchmark begin
        @inbounds for i in eachindex($b)
            $trans.b = $b[i]
            $output[i] = Limbdark.transit_poly!($trans)
        end
    end
    return median(bench).time * 1e-9
end
vals_limbdark = zeros(length(N_pts), length(N_us))
@progress for (j, Nu) in enumerate(N_us)
    @progress for (i, N) in enumerate(N_pts)
        vals_limbdark[i, j] = benchmark_limbdark(Nu, N)
    end
end

plot(
    N_pts,
    vals_limbdark ./ vals,
    xscale=:log10,
    ls=:solid,
    m=:o,
    msw=0,
    palette=palette(:plasma, length(N_us) + 2),
    lab=labels,
    xlabel="number of points",
    ylabel="relative speedup (times faster) - Limbdark.jl"
)

tmed = median(vals ./ vals[:, 1], dims=1)
plot(
    N_us |> vec,
    tmed |> vec,
    scale=:log10,
    ls=:solid,
    m=:o,
    msw=0,
    lab="Transits.jl",
    xlabel="number of limb-darkening coefficients",
    ylabel="relative timing",
    ylims=(0.8, 10),
    leg=:topleft
)
plot!(
    N_us |> vec,
    median(vals_limbdark ./ vals_limbdark[:, 1], dims=1) |> vec,
    lab="Limbdark.jl",
    c=2,
    ls=:solid,
    m=:o,
    msw=0
)
plot!(N_us |> vec, n -> n^(0.2), ls=:dash, c=:gray, lab=L"N^{0.2}")
plot!(N_us |> vec, n -> tmed[end] * n / N_us[1, end], ls=:dash, c=:gray, lab=L"N^1")


text:
"""
    AlgoAcoustModule

Module for linear acoustics algorithms.
"""
module AlgoAcoustModule

using FinEtools.FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import FinEtools.AlgoBaseModule: dcheck!
import FinEtools.FieldModule: ndofs, setebc!, numberdofs!, applyebc!, scattersysvec!
import FinEtools.NodalFieldModule: NodalField, nnodes
import FinEtools.FEMMBaseModule: associategeometry!, distribloads
import ..FEMMAcoustModule: acousticmass, acousticstiffness, nzebcloadsacousticmass, nzebcloadsacousticstiffness
import ..FEMMAcoustSurfModule: acousticABC
import FinEtools.ForceIntensityModule: ForceIntensity
import SparseArrays: spzeros
import LinearAlgebra: norm, lu, cross

"""
    steadystate(modeldata::FDataDict)

Steady-state acoustics solver.

`modeldata` = dictionary with string keys

- `"fens"`  = finite element node set
- `"regions"`  = array of region dictionaries
- `"essential_bcs"` = array of essential boundary condition dictionaries
- `"ABCs"` = array of absorbing boundary condition dictionaries
- `"flux_bcs"` = array of flux boundary condition dictionaries

For each region (connected piece of the domain made of a particular material),
mandatory, the  region dictionary  contains items:
- `"femm"` = finite element mmodel machine (mandatory);

For essential boundary conditions (optional) each dictionary
would hold
  + `"pressure"` = fixed (prescribed) pressure (scalar),  or
            a function with signature
                function T = f(x)
            If not given, zero pressure assumed.
  + `"node_list"` = list of nodes on the boundary to which the condition applies
            (mandatory)

For absorbing boundary conditions (optional) each dictionary
may hold
  + `"femm"` = finite element mmodel machine (mandatory).

For flux boundary conditions (optional) each dictionary
would hold
  + `"femm"` = finite element mmodel machine (mandatory);
  + `"normal_flux"` = normal component of the flux through the boundary (scalar),
        which is the normal derivative of the pressure.

# Output
`modeldata` = the dictionary is augmented with
- `"geom"` = the nodal field that is the geometry
- `"P"` = the nodal field that is the computed pressure (in the general a
            complex-number field)
"""
function steadystate(modeldata::FDataDict)
    # Lists of recognized keys for the data dictionaries:
    modeldata_recognized_keys = ["fens", "regions", "essential_bcs", "flux_bcs"]
    essential_bcs_recognized_keys = ["pressure", "node_list"]
    ABCs_recognized_keys = ["femm"]
    flux_bcs_recognized_keys = ["femm", "normal_flux"]
    regions_recognized_keys = ["femm"]

    # Check that keys defined  in the dictionary are matched
    dcheck!(modeldata, modeldata_recognized_keys)

    # Extract the nodes
    fens = get(()->error("Must get fens!"), modeldata, "fens")

    # Construct the geometry field
    geom = NodalField(fens.xyz)

    # Construct the acoustic pressure field
    P = NodalField(zeros(FCplxFlt,size(fens.xyz,1),1))

    # The following properties are the same for all regions and all boundaries (ABC)
    omega = get(()->error("Must get angular frequency!"), modeldata, "omega");

    # Apply the essential boundary conditions on the acoustic pressure field
    essential_bcs = get(modeldata, "essential_bcs", nothing);
    if (essential_bcs != nothing)
        for j = 1:length(essential_bcs)
            ess = essential_bcs[j]
            dcheck!(ess, essential_bcs_recognized_keys)
            fenids = get(()->error("Must get node list!"), ess, "node_list");
            pressure = get(ess, "pressure", nothing);
            Pfixed = zeros(FCplxFlt,length(fenids)); # default is zero pressure
            if (pressure != nothing)
                if (typeof(pressure) <: Function) # pressure supplied through function
                    for k = 1:length(fenids)
                        Pfixed[k] = pressure(geom.values[fenids[k],:])[1];
                    end
                else # pressure given as  a constant
                    fill!(Pfixed, pressure);
                end
            end
            setebc!(P, fenids[:], true, 1, Pfixed);
            applyebc!(P);
        end
    end

    # Number the equations
    numberdofs!(P)           #,Renumbering_options); # NOT DONE <<<<<<<<<<<<<<<<<

    # Initialize the acoustic load vector
    F = zeros(FCplxFlt, P.nfreedofs);

    # Construct the system acoustic mass and stiffness matrix
    # and the absorbing boundary condition (ABC) matrix
    C=  spzeros(P.nfreedofs,P.nfreedofs); # (all zeros, for the moment)
    S=  spzeros(P.nfreedofs,P.nfreedofs); # (all zeros, for the moment)
    D=  spzeros(P.nfreedofs,P.nfreedofs); # (all zeros, for the moment)
    regions = get(()->error("Must get regions!"), modeldata, "regions")
    for i = 1:length(regions)
        region = regions[i]
        dcheck!(region, regions_recognized_keys)
        femm = get(()->error("Must get femm for the region!"), region, "femm");
        # Add up all the acoustic mass matrices for all the regions
        C = C + acousticmass(femm, geom, P);
        S = S + acousticstiffness(femm, geom, P);
        # Loads due to the essential boundary conditions on the pressure field
        essential_bcs = get(modeldata, "essential_bcs", nothing);
        if (essential_bcs != nothing)
            F = F + nzebcloadsacousticmass(femm, geom, P);
            F = F + nzebcloadsacousticstiffness(femm, geom, P);
        end
    end

    # Compute the ABC matrices
    ABCs = get(modeldata, "ABCs", nothing);
    if (ABCs != nothing)
        for j = 1:length(ABCs)
            ABC = ABCs[j]
            dcheck!(ABC, ABCs_recognized_keys)
            femm = get(()->error("Must get femm for the ABC!"), ABC, "femm");
            D = D + acousticABC(femm, geom, P);
        end
    end

    # Process the flux boundary condition:
    # dP/dn=-rho*a ... The normal derivative of the pressure in terms of the acceleration
    flux_bcs = get(modeldata, "flux_bcs", nothing);
    if (flux_bcs != nothing)
        for j = 1:length(flux_bcs)
            fluxbc = flux_bcs[j]
            dcheck!(fluxbc, flux_bcs_recognized_keys)
            normal_flux = get(()->error("Must get normal flux value!"), fluxbc, "normal_flux");
            if (typeof(normal_flux) <: Function)
                fi = ForceIntensity(FFlt, 1, normal_flux);
            else
                if typeof(normal_flux) <: AbstractArray
                else
                    normal_flux = FCplxFlt[normal_flux]
                end
                fi = ForceIntensity(normal_flux);
            end
            femm = get(()->error("Must get femm for the flux BC!"), fluxbc, "femm");
            F = F + distribloads(femm, geom, P, fi, 2);
        end
    end

    # Solve for the pressures
    K = lu((-omega^2*S +omega*1.0im*D + C));
    vP = K\F;
    scattersysvec!(P, vP[:])

    # Update the model data
    setindex!(modeldata, geom, "geom");
    setindex!(modeldata, P, "P");
    return modeldata            # ... And return the updated model data
end

end


text:
using KernelAbstractions
using KernelAbstractions.Extras
using CUDA
using CUDAKernels
using Random123
using WormlikeChain
using BenchmarkTools
using FastClosures

#BAOAB integration
@kernel function simulate!(positions,
                          velocities,
                          externalforce!,
                          bondedforce,
                          @Const(beadparams),
                          invmass,
                          γ,
                          invβ,
                          VNbeads::Val{Nbeads},
                          step0::UInt64,
                          Vrngkey::Val{rngkey},
                          Δt,
                          VNdims::Val{Ndims},
                          VChainbounds::Val{Chainbounds},
                          VNsteps::Val{Nsteps}) where {Ndims, Nbeads, rngkey, Chainbounds, Nsteps}
    @uniform Nthreads = groupsize()[1]
    @uniform Nbeadsperthread= (Nbeads-1)÷Nthreads + 1
    @uniform Nrngsperbead= (Ndims-1)÷4 + 1
    tid = @index(Local, Linear)
    groupid = @index(Group, Linear)
    forces = @localmem Int64 (Nbeads, Ndims)
    #@inbounds sharedA[I] = 0
    #TODO set up force calc shared mem
    for step in 1:Nsteps
        #Reset forces
        #Kick
        @unroll for beadloopi in 1:Nbeadsperthread
            beadid= tid + (beadloopi-1)*Nthreads
            previd = WormlikeChain.prev_beadid(beadid,Chainbounds)
            nextid = WormlikeChain.next_beadid(beadid,Chainbounds)
            #load bead1
            bead = ntuple(@closure(i -> @inbounds(positions[beadid,i])), Ndims)
            #load bead2
            #load bead3
            #load beadparams
            #calc bonded forces on bead1,2,3
            #force1, force2, force3= bondedforce(beadparam, bead1, bead2, bead3, step0+step)
            #inc force3
            @synchronize
            #inc force2
            @synchronize

            force = externalforce(beadparam, bead1)
            #externalforce!(forces,positions,beadparams,beadid)
            #apply total forces to each bead
            @unroll for i in 1:Ndims
                @inbounds velocities[beadid,i] += (invmass*Δt*(2f0^-32))*forces[beadid,i]
            end
        end
        #Drift with O step half way
        @unroll for beadloopi in 1:Nbeadsperthread
            mybeadid= tid + (beadloopi-1)*Nthreads
            randns= WormlikeChain.randn_32(rngkey, step0+step, (mybeadid<<32)%UInt64,VNdims)
            @unroll for i in 1:Ndims
                @inbounds positions[mybeadid,i] += (Δt/2) * velocities[mybeadid,i]
                @inbounds velocities[mybeadid,i] *= exp(-Δt*γ)
                @inbounds velocities[mybeadid,i] += (√(1-exp(-2γ*Δt))*√(invβ*invmass)) * randns[i]
                @inbounds positions[mybeadid,i] += (Δt/2) * velocities[mybeadid,i]
            end
        end
        @synchronize
    end
end

N=1024

simkh = simulate!(CPU(), N)
simkd = simulate!(CUDADevice(), N)

pos_h = zeros(Float32,N,3)
pos_d = cu(pos_h)
vel_h = zeros(Float32,N,3)
vel_d = cu(vel_h)
params_h = fill((k=4.0f0,),N)
params_d = cu(params_h)

function externalforce!(forces,positions,beadparams,id)
    @unroll for i in 1:3
        @inbounds f = -positions[id,i] * beadparams[id].k
        @inbounds forces[id,i] += unsafe_trunc(Int64,f*(2f0^32))
    end
end

function externalforce(beadparams,id)
    @unroll for i in 1:3
        @inbounds f = -positions[id,i] * beadparams[id].k
        @inbounds forces[id,i] += unsafe_trunc(Int64,f*(2f0^32))
    end
end

@btime wait(simkh(pos_h,
                vel_h,
                externalforce!,
                params_h,
                1f0,
                1f0/10f0,
                1f0,
                Val(N),
                1%UInt64,
                Val(124432%UInt64),
                0.1f0,
                Val(3),
                Val((1,N+1)),
                Val(1000),
                 ndrange=N))

@btime wait(simkd(pos_d,
                vel_d,
                externalforce!,
                params_d,
                1f0,
                1f0/10f0,
                1f0,
                Val(N),
                1%UInt64,
                Val(124432%UInt64),
                0.1f0,
                Val(3),
                Val((1,N+1)),
                Val(1000),
                 ndrange=N))


text:
using BenchmarkTools
const SUITE = BenchmarkGroup()

using Random
using Catlab.CategoricalAlgebra

bench = SUITE["Limits"] = BenchmarkGroup()

function benchmark_pullback(suite::BenchmarkGroup, name, cospan)
  for alg in (NestedLoopJoin(), SortMergeJoin(), HashJoin())
    suite["$name:$(nameof(typeof(alg)))"] =
      @benchmarkable limit($cospan, alg=$alg)
  end
end

sizes = (100, 150)
f, g = (FinFunction(ones(Int, size), 1) for size in sizes)
benchmark_pullback(bench, "pullback_terminal", Cospan(f, g))

f = FinFunction(identity, FinSet(10000))
benchmark_pullback(bench, "pullback_identity", Cospan(f, f))

n = 10000
Random.seed!(1)
f, g = FinFunction(randperm(n), n), FinFunction(randperm(n), n)
benchmark_pullback(bench, "pullback_randperm", Cospan(f, g))

k = 1000
sizes = (9000, 11000)
Random.seed!(1)
f, g = (FinFunction(rand(1:k, size), k) for size in sizes)
benchmark_pullback(bench, "pullback_randsparse", Cospan(f, g))

k = 50
sizes = (400, 500, 600)
Random.seed!(1)
f, g, h = (FinFunction(rand(1:k, size), k) for size in sizes)
benchmark_pullback(bench, "ternary_pullback_randsparse",
                   SMulticospan{3}(f, g, h))


text:
# Simulation of coupled Lorenz systems.

using Causal
using Plots

# Construct the model
ε = 10.

@defmodel model begin
    @nodes begin
        ds1 = ForcedLorenzSystem()
        ds2 = ForcedLorenzSystem()
        coupler = Coupler(conmat = ε*[-1. 1; 1 -1], cplmat=[1. 0 0; 0 0 0; 0 0 0])
        writer = Writer(input=Inport(6))
    end
    @branches begin
        ds1[1:3]        =>      coupler[1:3]
        ds2[1:3]        =>      coupler[4:6]
        coupler[1:3]    =>      ds1[1:3]
        coupler[4:6]    =>      ds2[1:3]
        ds1[1:3]        =>      writer[1:3]
        ds2[1:3]        =>      writer[4:6]
    end
end

# Plot signal flow diagram of model
display(signalflow(model))

# Simulate the model
simulate!(model, 0., 0.01, 100)

# Plot signal flow diagram of model
display(signalflow(model))

# Read simulation data
t, x = read(getnode(model, :writer).component)

# Compute errors
err = x[:, 1] - x[:, 4]

# Plot the results.
p1 = plot(x[:, 1], x[:, 2])
p2 = plot(x[:, 4], x[:, 5])
p3 = plot(t, err)
display(plot(p1, p2, p3, layout=(3, 1)))


text:
using StaticArrays
function cspline(s, t1, t2, p1, p2, m1, m2)
    d = t2-t1
    t = (s-t1)/(t2-t1)
    t2 = t*t
    t3 = t2*t
    z1, z2, z3, z4 = (1 - 3*t2 + 2*t3,  3*t2 - 2*t3,  t - 2*t2 + t3, 0 - t2 + t3)
    z1*p1 + z2*p2 + z3*d*m1 + z4*d*m2
end
function intcspline(s, t1, t2, p1, p2, m1, m2)
        d = t2-t1
        t = (s-t1)/(t2-t1)
        t2 = t*t
        t3 = t2*t
        t4 = t2*t2
        t4, t3, t2 = t4/4, t3/3, t2/2
        z1, z2, z3, z4 = (t - 3*t3 + 2*t4,  3*t3 - 2*t4,  t2 - 2*t3 + t4, 0 - t3 + t4)
        (z1*p1 + z2*p2 + z3*d*m1 + z4*d*m2)*d
end
intcspline(s, T, t1, t2, p1, p2, m1, m2) = intcspline(T, t1, t2, p1, p2, m1, m2) - intcspline(s, t1, t2, p1, p2, m1, m2)

struct CSpline{T}
    s; t; x::T; y::T; mx; my
end

"""
    CSpline(s, t, x, y = x, m0 = (y-x)/(t-s), m1 = (y-x)/(t-s))

Cubic spline parametrized by ``f(s) = x`` and ``f(t) = y``, ``f'(s) = m_0``, ``f'(t) = m_1``.
"""
CSpline(s, t, x::T, y = x, m0 = (y-x)/(t-s), m1 =  (y-x)/(t-s)) where {T} = CSpline{T}(s, t, x, y, mx, my)
(cs::CSpline)(t) =  cspline(t, cs.s, cs.t, cs.x, cs.y, cs.mx, cs.my)

"""
    integrate(cs::CSpline, s, t)

Integrate the cubic spline from `s` to `t`.
"""
integrate(cs::CSpline, s, t) = intcspline(s,t, cs.s, cs.t, cs.x, cs.y, cs.mx, cs.my)


text:
#
# AT.jl --
#
# Implement Julia interface to Andor SDK.
#
#-------------------------------------------------------------------------------
#
# This file is part of "AndorCameras.jl" released under the MIT license.
#
# Copyright (C) 2017-2019, Éric Thiébaut.
#

"""
```julia
using AndorCameras.AT
```

makes all Andor cameras constants and low-level functions available (they are
all prefixed with `AT.`).

"""
module AT

export
    @L_str

isfile(joinpath(@__DIR__,"..","deps","deps.jl")) ||
    error("Tcl not properly installed.  Please run `Pkg.build(\"Tcl\")` to create file \"",joinpath(@__DIR__,"..","deps","deps.jl"),"\"")

include("../deps/deps.jl")

struct Status
    func::Symbol
    code::STATUS
end

"""
```julia
@_call(func, rtype, proto, args...)
```

yields code to call C function `func` in Andor SDK library assuming `rtype`
is the return type of the function, `proto` is a tuple of the argument types
and `args...` are the arguments.

The return type `rtype` must be `STATUS` and the produced code is wrapped so
that an instance of `AT.Status` is returned with the status code and the
symbolic name of the called SDK function.

"""
macro _call(func, rtype, args...)
    rtype == :STATUS || error("return type must be STATUS")
    qfunc = _quoted(func)
    expr = Expr(:call, :ccall, Expr(:tuple, qfunc, :_DLL), rtype, args...)
    return quote
        Status($qfunc, $(esc(expr)))
    end
end
_quoted(x::QuoteNode) = x
_quoted(x::Symbol) = QuoteNode(x)
_quoted(x::AbstractString) = _quoted(Symbol(x))

InitialiseLibrary() = @_call(:AT_InitialiseLibrary, STATUS, ())

FinaliseLibrary() = @_call(:AT_FinaliseLibrary, STATUS, ())

Open(index) =
    (ref = Ref{HANDLE}();
     (@_call(:AT_Open, STATUS, (INDEX, Ref{HANDLE}),
             index, ref), ref[]))

Close(handle) = @_call(:AT_Close, STATUS, (HANDLE,), handle)

# FIXME: not yet interfaced:
# typedef int (AT_EXP_CONV *FeatureCallback)(AT_H Hndl, const AT_WC* Feature, void* Context);
# int AT_EXP_CONV AT_RegisterFeatureCallback(AT_H Hndl, const AT_WC* Feature,
#                                            FeatureCallback EvCallback, void* Context);
# int AT_EXP_CONV AT_UnregisterFeatureCallback(AT_H Hndl, const AT_WC* Feature,
#                                              FeatureCallback EvCallback, void* Context);

IsImplemented(handle, feature) =
    (ref = Ref{BOOL}(FALSE);
     (@_call(:AT_IsImplemented, STATUS, (HANDLE, FEATURE, Ref{BOOL}),
             handle, feature, ref), ref[] != FALSE))

IsReadable(handle, feature) =
    (ref = Ref{BOOL}(FALSE);
     (@_call(:AT_IsReadable, STATUS, (HANDLE, FEATURE, Ref{BOOL}),
             handle, feature, ref), ref[] != FALSE))

IsWritable(handle, feature) =
    (ref = Ref{BOOL}(FALSE);
     (@_call(:AT_IsWritable, STATUS, (HANDLE, FEATURE, Ref{BOOL}),
             handle, feature, ref), ref[] != FALSE))

IsReadOnly(handle, feature) =
    (ref = Ref{BOOL}(FALSE);
     (@_call(:AT_IsReadOnly, STATUS, (HANDLE, FEATURE, Ref{BOOL}),
             handle, feature, ref), ref[] != FALSE))

SetInt(handle, feature, value::Integer) =
    @_call(:AT_SetInt, STATUS, (HANDLE, FEATURE, INT),
           handle, feature, value)

GetInt(handle, feature) =
    (ref = Ref{INT}();
     (@_call(:AT_GetInt, STATUS, (HANDLE, FEATURE, Ref{INT}),
             handle, feature, ref), ref[]))

GetIntMin(handle, feature) =
    (ref = Ref{INT}();
     (@_call(:AT_GetIntMin, STATUS, (HANDLE, FEATURE, Ref{INT}),
             handle, feature, ref), ref[]))

GetIntMax(handle, feature) =
    (ref = Ref{INT}();
     (@_call(:AT_GetIntMax, STATUS, (HANDLE, FEATURE, Ref{INT}),
             handle, feature, ref), ref[]))

SetFloat(handle, feature, value::Real) =
    @_call(:AT_SetFloat, STATUS, (HANDLE, FEATURE, FLOAT),
           handle, feature, value)

GetFloat(handle, feature) =
    (ref = Ref{FLOAT}();
     (@_call(:AT_GetFloat, STATUS, (HANDLE, FEATURE, Ref{FLOAT}),
             handle, feature, ref), ref[]))

GetFloatMin(handle, feature) =
    (ref = Ref{FLOAT}();
     (@_call(:AT_GetFloatMin, STATUS, (HANDLE, FEATURE, Ref{FLOAT}),
             handle, feature, ref), ref[]))

GetFloatMax(handle, feature) =
    (ref = Ref{FLOAT}();
     (@_call(:AT_GetFloatMax, STATUS, (HANDLE, FEATURE, Ref{FLOAT}),
             handle, feature, ref), ref[]))

SetBool(handle, feature, value::Bool) =
    @_call(:AT_SetBool, STATUS, (HANDLE, FEATURE, BOOL),
           handle, feature, (value ? TRUE : FALSE))

GetBool(handle, feature) =
    (ref = Ref{BOOL}();
     (@_call(:AT_GetBool, STATUS, (HANDLE, FEATURE, Ref{BOOL}),
             handle, feature, ref), ref[] != FALSE))

SetEnumIndex(handle, feature, enum::Integer) =
    @_call(:AT_SetEnumIndex, STATUS, (HANDLE, FEATURE, ENUM),
           handle, feature, enum)

SetEnumerated(handle, feature, enum::Integer) =
    @_call(:AT_SetEnumerated, STATUS, (HANDLE, FEATURE, ENUM),
           handle, feature, enum)

SetEnumString(handle, feature, value) =
    @_call(:AT_SetEnumString, STATUS, (HANDLE, FEATURE, Ptr{WCHAR}),
           handle, feature, value)

SetEnumeratedString(handle, feature, value) =
    @_call(:AT_SetEnumeratedString, STATUS, (HANDLE, FEATURE, Ptr{WCHAR}),
           handle, feature, value)

GetEnumIndex(handle, feature) =
    (ref = Ref{ENUM}();
     (@_call(:AT_GetEnumIndex, STATUS, (HANDLE, FEATURE, Ref{ENUM}),
             handle, feature, ref), Int(ref[])))

GetEnumerated(handle, feature) =
    (ref = Ref{ENUM}();
     (@_call(:AT_GetEnumerated, STATUS, (HANDLE, FEATURE, Ref{ENUM}),
             handle, feature, ref), Int(ref[])))

GetEnumCount(handle, feature) =
    (ref = Ref{LENGTH}();
     (@_call(:AT_GetEnumCount, STATUS, (HANDLE, FEATURE, Ref{LENGTH}),
             handle, feature, ref), Int(ref[])))

GetEnumeratedCount(handle, feature) =
    (ref = Ref{LENGTH}();
     (@_call(:AT_GetEnumeratedCount, STATUS, (HANDLE, FEATURE, Ref{LENGTH}),
             handle, feature, ref), Int(ref[])))

IsEnumIndexAvailable(handle, feature, enum::Integer) =
    (ref = Ref{BOOL}();
     (@_call(:AT_IsEnumIndexAvailable, STATUS,
             (HANDLE, FEATURE, ENUM, Ref{BOOL}),
             handle, feature, enum, ref), ref[] != FALSE))

IsEnumeratedIndexAvailable(handle, feature, enum::Integer) =
    (ref = Ref{BOOL}();
     (@_call(:AT_IsEnumeratedIndexAvailable, STATUS,
             (HANDLE, FEATURE, ENUM, Ref{BOOL}),
             handle, feature, enum, ref), ref[] != FALSE))

IsEnumIndexImplemented(handle, feature, enum::Integer) =
    (ref = Ref{BOOL}();
     (@_call(:AT_IsEnumIndexImplemented, STATUS,
             (HANDLE, FEATURE, ENUM, Ref{BOOL}),
             handle, feature, enum, ref), ref[] != FALSE))

IsEnumeratedIndexImplemented(handle, feature, enum::Integer) =
    (ref = Ref{BOOL}();
     (@_call(:AT_IsEnumeratedIndexImplemented, STATUS,
             (HANDLE, FEATURE, ENUM, Ref{BOOL}),
             handle, feature, enum, ref), ref[] != FALSE))

GetEnumStringByIndex(handle, feature, enum::Integer, buf::DenseVector{WCHAR}) =
    GetEnumStringByIndex(handle, feature, enum, pointer(buf), length(buf))

GetEnumStringByIndex(handle, feature, enum::Integer, ptr::Ptr{WCHAR}, len::Integer) =
    @_call(:AT_GetEnumStringByIndex, STATUS,
           (HANDLE, FEATURE, ENUM, Ptr{WCHAR}, LENGTH),
           handle, feature, enum, ptr, len)

GetEnumeratedString(handle, feature, enum::Integer, buf::DenseVector{WCHAR}) =
    GetEnumeratedString(handle, feature, enum, pointer(buf), length(buf))

GetEnumeratedString(handle, feature, enum::Integer, ptr::Ptr{WCHAR}, len::Integer) =
    @_call(:AT_GetEnumeratedString, STATUS,
           (HANDLE, FEATURE, ENUM, Ptr{WCHAR}, LENGTH),
           handle, feature, enum, ptr, len)

Command(handle, feature) =
    @_call(:AT_Command, STATUS, (HANDLE, FEATURE), handle, feature)

SetString(handle, feature, value) =
    @_call(:AT_SetString, STATUS, (HANDLE, FEATURE, Ptr{WCHAR}),
           handle, feature, value)

GetString(handle, feature, buf::DenseVector{WCHAR}) =
    GetString(handle, feature, pointer(buf), length(buf))

GetString(handle, feature, ptr::Ptr{WCHAR}, len::Integer) =
    @_call(:AT_GetString, STATUS, (HANDLE, FEATURE, Ptr{WCHAR}, LENGTH),
           handle, feature, ptr, len)

GetStringMaxLength(handle, feature) =
    (ref = Ref{LENGTH}();
     (@_call(:AT_GetStringMaxLength, STATUS, (HANDLE, FEATURE, Ref{LENGTH}),
             handle, feature, ref), Int(ref[])))

QueueBuffer(handle, buf::DenseVector{BYTE}) =
    QueueBuffer(handle, pointer(buf), sizeof(buf))

QueueBuffer(handle, ptr::Ptr{BYTE}, siz::Integer) =
    @_call(:AT_QueueBuffer, STATUS, (HANDLE, Ptr{BYTE}, LENGTH),
           handle, ptr, siz)

WaitBuffer(handle, timeout::Integer) =
    (refptr = Ref{Ptr{BYTE}}();
     refsiz = Ref{LENGTH}();
     (@_call(:AT_WaitBuffer, STATUS, (HANDLE, Ref{Ptr{BYTE}}, Ref{LENGTH}, MSEC),
             handle, refptr, refsiz, timeout), refptr[], Int(refsiz[])))

Flush(handle) =
    @_call(:AT_Flush, STATUS, (HANDLE,), handle)

"""
```julia
widestring(str, len = strlen(str))
```

yields a vector of wide characters (`Cwchar_t`) with the contents of the string
`str` and properly zero-terminated.  This buffer is independent from the input
string and its contents can be overwritten.  An error is thrown if `str`
contains any embedded NULL characters (which would cause the string to be
silently truncated if the C routine treats NULL as the terminator).

An alternative (without the checking of embedded NULL characters) is:

```julia
push!(transcode(Cwchar_t, str), convert(Cwchar_t, 0))
```

This method is used to implement the `@L_str` macro which converts a
literal string into a wide character string.  For instance:

```julia
L"EventSelector"
```

"""
function widestring(str::AbstractString,
                    len::Integer = length(str)) :: Array{WCHAR}
    buf = Array{WCHAR}(undef, len + 1)
    i = 0
    @inbounds for c in str
        if i ≥ len
            break
        end
        c != '\0' || error("strings must not have embedded NULL characters")
        i += 1
        buf[i] = c
    end
    @inbounds while i ≤ len
        i += 1
        buf[i] = zero(WCHAR)
    end
    return buf
end

widestring(sym::Symbol) = widestring(string(sym))

function widestringtostring(arr::Array{WCHAR}) :: String
    len = length(arr)
    @inbounds while len > 0 && arr[len] == zero(WCHAR)
        len -= 1
    end
    buf = Vector{Char}(undef, len)
    @inbounds for i in 1:len
        c = arr[i]
        c != zero(WCHAR) || error("strings must not have embedded NULL characters")
        buf[i] = c
    end
    return String(buf)
end

"""

```julia
L"text"
```

yields an array of wide characters.

"""
macro L_str(str)
    :(widestring($str))
end

end


text:
"""
    FockBasis(N,offset=0)

Basis for a Fock space where `N` specifies a cutoff, i.e. what the highest
included fock state is. Similarly, the `offset` defines the lowest included
fock state (default is 0). Note that the dimension of this basis is `N+1-offset`.
"""
struct FockBasis{T} <: Basis
    shape::Vector{T}
    N::T
    offset::T
    function FockBasis(N::T,offset::T=0) where T
        if N < 0 || offset < 0 || N <= offset
            throw(DimensionMismatch())
        end
        new{T}([N-offset+1], N, offset)
    end
end

==(b1::FockBasis, b2::FockBasis) = (b1.N==b2.N && b1.offset==b2.offset)

"""
    number([T=ComplexF64,] b::FockBasis)

Number operator for the specified Fock space with optional data type `T`.
"""
function number(::Type{T}, b::FockBasis) where T
    diag = T.(b.offset:b.N)
    data = spdiagm(0 => diag)
    SparseOperator(b, data)
end

number(b::FockBasis) = number(ComplexF64, b)

"""
    destroy([T=ComplexF64,] b::FockBasis)

Annihilation operator for the specified Fock space with optional data type `T`.
"""
function destroy(::Type{C}, b::FockBasis) where C
    T = real(C)
    ns = b.offset+1.:b.N
    diag = @. C(sqrt(T(ns)))
    data = spdiagm(1 => diag)
    SparseOperator(b, data)
end

destroy(b::FockBasis) = destroy(ComplexF64, b)

"""
    create([T=ComplexF64,] b::FockBasis)

Creation operator for the specified Fock space with optional data type `T`.
"""
function create(::Type{C}, b::FockBasis) where C
    T = real(C)
    ns = b.offset+1.:b.N
    diag = @. C(sqrt(T(ns)))
    data = spdiagm(-1 => diag)
    SparseOperator(b, data)
end

create(b::FockBasis) = create(ComplexF64, b)

"""
    displace([T=ComplexF64,] b::FockBasis, alpha)

Displacement operator ``D(α)`` for the specified Fock space with optional data type `T`.
"""
function displace(::Type{T}, b::FockBasis, alpha::Number) where T
    alpha = T(alpha)
    exp(dense(alpha * create(T, b) - conj(alpha) * destroy(T, b)))
end

displace(b::FockBasis, alpha::T) where {T <: Number} = displace(ComplexF64, b, alpha)

"""
    fockstate([T=ComplexF64,] b::FockBasis, n)

Fock state ``|n⟩`` for the specified Fock space.
"""
function fockstate(::Type{T}, b::FockBasis, n::Integer) where T
    @assert b.offset <= n <= b.N
    basisstate(T, b, n+1-b.offset)
end
fockstate(b, n) = fockstate(ComplexF64, b, n)

"""
    coherentstate([T=ComplexF64,] b::FockBasis, alpha)

Coherent state ``|α⟩`` for the specified Fock space.
"""
function coherentstate(::Type{T}, b::FockBasis, alpha::Number) where T
    result = Ket(T, b)
    coherentstate!(result, b, alpha)
    return result
end
coherentstate(b, alpha) = coherentstate(ComplexF64, b, alpha)

"""
    coherentstate!(ket::Ket, b::FockBasis, alpha)

Inplace creation of coherent state ``|α⟩`` for the specified Fock space.
"""
function coherentstate!(ket::Ket, b::FockBasis, alpha::Number)
    C = eltype(ket)
    T = real(C)
    alpha = C(alpha)
    data = ket.data
    data[1] = exp(-abs2(alpha)/2)

    # Compute coefficient up to offset
    offset = b.offset
    @inbounds for n=1:offset
        data[1] *= alpha/sqrt(T(n))
    end

    # Write coefficients to state
    @inbounds for n=1:b.N-offset
        data[n+1] = data[n]*alpha/sqrt(T(n+offset))
    end

    return ket
end


text:
function _line!(simple_plot::SimplePlot, cur_mode::AbstractString, varargs...; kwargs...)
  if length(varargs) == 1
    first_item = first(varargs)

    if isa(first_item, Array)
      cur_list = first_item
    else
      cur_list = collect(first_item)
    end
  else
    cur_list = collect(varargs)
  end

  init_dict = Dict{Any,Any}(
    "x0" => 0,
    "x1" => 1,
    "y0" => 0,
    "y1" => 1,
  )

  for cur_entry in cur_list
    cur_dict = deepcopy(init_dict)

    if cur_mode == "x"
      cur_dict["y0"] = string(cur_entry)
      cur_dict["y1"] = string(cur_entry)
      cur_dict["xref"] = "paper"
    else
      @assert cur_mode == "y"

      cur_dict["x0"] = string(cur_entry)
      cur_dict["x1"] = string(cur_entry)
      cur_dict["yref"] = "paper"
    end

    cur_dict["type"] = "line"

    cur_label = nothing
    haskey(kwargs, :name) && ( cur_label = kwargs[:name] )
    haskey(kwargs, :label) && ( cur_label = kwargs[:label] )

    if !isnothing(cur_label)
      if strip(cur_label) == ""
        cur_dict["showlegend"] = false
      else
        cur_dict["name"] = cur_label
      end
    end

    cur_alpha = nothing
    haskey(kwargs, :alpha) && ( cur_alpha = kwargs[:alpha] )
    haskey(kwargs, :opacity) && ( cur_alpha = kwargs[:opacity] )
    isnothing(cur_alpha) || ( cur_dict["opacity"] = cur_alpha )

    if haskey(kwargs, :color)
      cur_index = kwargs[:color]
      @assert isa(cur_index, Int)
    else
      simple_plot.index += 1
      cur_index = simple_plot.index
    end

    cur_index = 1 + ( (cur_index-1) % length(_palette) )
    cur_color = _palette[cur_index]

    sub_dict = Dict()
    sub_dict["color"] = cur_color

    cur_width = nothing
    haskey(kwargs, :width) && ( cur_width = kwargs[:width] )
    haskey(kwargs, :linewidth) && ( cur_width = kwargs[:linewidth] )
    isnothing(cur_width) || ( sub_dict["width"] = cur_width )

    if haskey(kwargs, :linestyle) && string(kwargs[:linestyle]) != "auto"
      sub_dict["dash"] = kwargs[:linestyle]
    end

    cur_dict["line"] = sub_dict

    push!(simple_plot.layout["shapes"], cur_dict)
  end

  validate!(simple_plot)
end


text:
############################################################################################
## This file is part of RLinearAlgebra.jl
##
## Overview: abstractions and methods for performing various randomized linear algebra
## computations
##
## Contents
## - Dependencies
## - Export Statements
##  + Linear Sampler
##  + Linear Solver Routine
##  + Linear Solver Log
##  + Linear Solver Stopping Criteria
## - Source File Inclusions
##
############################################################################################

module RLinearAlgebra

###########################################
# Dependencies
###########################################

using LinearAlgebra, Random, Distributions

###########################################
# Exports
###########################################

#*****************************************#
# Linear Sampler Exports
#*****************************************#

# Abstract Types
export LinSysSampler, LinSysSketch, LinSysSelect
export LinSysVecRowSampler, LinSysVecRowSketch, LinSysVecRowSelect
export LinSysVecColSampler, LinSysVecColSketch, LinSysVecColSelect
export LinSysBlkRowSampler, LinSysBlkRowSketch, LinSysBlkRowSelect
export LinSysBlkColSampler, LinSysBlkColSketch, LinSysBlkColSelect

# Vector Row Samplers
export LinSysVecRowDetermCyclic, LinSysVecRowHopRandCyclic, LinSysVecRowOneRandCyclic,
    LinSysVecRowPropToNormSampler, LinSysVecRowSVSampler, LinSysVecRowRandCyclic,
    LinSysVecRowUnidSampler

# Vector Column Samplers
export LinSysVecColDetermCyclic

#*****************************************#
# Linear Solver Routine Exports
#*****************************************#

# Abstract Types
export LinSysSolveRoutine, LinSysVecRowProjection, LinSysVecColProjection,
    LinSysBlkRowProjection, LinSysBlkColProjection, LinSysPreconKrylov

# Vector Row Projection
export LinSysVecRowProjStd, Kaczmarz, ART, LinSysVecRowProjPO, LinSysVecRowProjFO

# Vector Column Projection
export LinSysVecColProjStd, CoordinateDescent, GaussSeidel, LinSysVecColProjPO,
    LinSysVecColProjFO

#*****************************************#
# Linear Solver Log Exports
#*****************************************#
export LinSysSolverLog
export LSLogOracle, LSLogFull

#*****************************************#
# Linear Solver Stopping Criteria Exports
#*****************************************#
export LinSysStopCriterion
export LSStopMaxIterations

#*****************************************#
# Randomized Linear Solver Exports
#*****************************************#
export RLSSolver, rsolve, rsolve!

###########################################
# Source File Inclusions
###########################################

include("tools.jl")
include("linear_samplers.jl")
include("linear_solver_routines.jl")
include("linear_solver_logs.jl")
include("linear_solver_stops.jl")
include("linear_rsolve.jl")

end # module


text:
function summarize(d::AbstractDocument; ns=5)
    sentences = sentence_tokenize(language(d), text(d))
    s = StringDocument.(sentences)
    c = Corpus(s)
    prepare!(c, strip_case | strip_stopwords | stem_words )
    update_lexicon!(c)
    t = tf_idf(dtm(c))
    T = t * t'
    p=pagerank(T)
    return sentences[sort(sortperm(vec(p), rev=true)[1:ns])]
end

function pagerank( A; Niter=20, damping=.15)
         Nmax = size(A, 1)
         r = rand(1,Nmax);              # Generate a random starting rank.
         r = r ./ norm(r,1);            # Normalize
         a = (1-damping) ./ Nmax;       # Create damping vector

         for i=1:Niter
             s = r * A
             rmul!(s, damping)
             r = s .+ (a * sum(r, dims=2));   # Compute PageRank.
         end

         r = r./norm(r,1);

         return r
end


text:
#####
##### `frule`/`rrule`
#####

#=
In some weird ideal sense, the fallback for e.g. `frule` should actually be "get
the derivative via forward-mode AD". This is necessary to enable mixed-mode
rules, where e.g. `frule` is used within a `rrule` definition. For example,
broadcasted functions may not themselves be forward-mode *primitives*, but are
often forward-mode *differentiable*.

ChainRulesCore, by design, is decoupled from any specific AD implementation. How,
then, do we know which AD to fall back to when there isn't a primitive defined?

Well, if you're a greedy AD implementation, you can just overload `frule` and/or
`rrule` to use your AD directly. However, this won't play nice with other AD
packages doing the same thing, and thus could cause load-order-dependent
problems for downstream users.

It turns out, Cassette solves this problem nicely by allowing AD authors to
overload the fallbacks w.r.t. their own context. Example using ForwardDiff:

```
using ChainRulesCore, ForwardDiff, Cassette

Cassette.@context MyChainRuleCtx

# ForwardDiff, itself, can call `my_frule` instead of
# `frule` to utilize the ForwardDiff-injected ChainRulesCore
# infrastructure
my_frule(args...) = Cassette.overdub(MyChainRuleCtx(), frule, args...)

function Cassette.execute(::MyChainRuleCtx, ::typeof(frule), f, x::Number)
    r = frule(f, x)
    if isa(r, Nothing)
        fx, df = (f(x), Rule(Δx -> ForwardDiff.derivative(f, x) * Δx))
    else
        fx, df = r
    end
    return fx, df
end
```
=#

"""
    frule(f, x...)

Expressing `x` as the tuple `(x₁, x₂, ...)` and the output tuple of `f(x...)`
as `Ω`, return the tuple:

    (Ω, (rule_for_ΔΩ₁::AbstractRule, rule_for_ΔΩ₂::AbstractRule, ...))

where each returned propagation rule `rule_for_ΔΩᵢ` can be invoked as

    rule_for_ΔΩᵢ(Δx₁, Δx₂, ...)

to yield `Ωᵢ`'s corresponding differential `ΔΩᵢ`. To illustrate, if all involved
values are real-valued scalars, this differential can be written as:

    ΔΩᵢ = ∂Ωᵢ_∂x₁ * Δx₁ + ∂Ωᵢ_∂x₂ * Δx₂ + ...

If no method matching `frule(f, xs...)` has been defined, then return `nothing`.

Examples:

unary input, unary output scalar function:

```jldoctest
julia> x = rand();

julia> sinx, dsin = frule(sin, x);

julia> sinx == sin(x)
true

julia> dsin(1) == cos(x)
true
```

unary input, binary output scalar function:

```jldoctest
julia> x = rand();

julia> sincosx, (dsin, dcos) = frule(sincos, x);

julia> sincosx == sincos(x)
true

julia> dsin(1) == cos(x)
true

julia> dcos(1) == -sin(x)
true
```

See also: [`rrule`](@ref), [`AbstractRule`](@ref), [`@scalar_rule`](@ref)
"""
frule(::Any, ::Vararg{Any}; kwargs...) = nothing

"""
    rrule(f, x...)

Expressing `x` as the tuple `(x₁, x₂, ...)` and the output tuple of `f(x...)`
as `Ω`, return the tuple:

    (Ω, (rule_for_Δx₁::AbstractRule, rule_for_Δx₂::AbstractRule, ...))

where each returned propagation rule `rule_for_Δxᵢ` can be invoked as

    rule_for_Δxᵢ(ΔΩ₁, ΔΩ₂, ...)

to yield `xᵢ`'s corresponding differential `Δxᵢ`. To illustrate, if all involved
values are real-valued scalars, this differential can be written as:

    Δxᵢ = ∂Ω₁_∂xᵢ * ΔΩ₁ + ∂Ω₂_∂xᵢ * ΔΩ₂ + ...

If no method matching `rrule(f, xs...)` has been defined, then return `nothing`.

Examples:

unary input, unary output scalar function:

```jldoctest
julia> x = rand();

julia> sinx, dx = rrule(sin, x);

julia> sinx == sin(x)
true

julia> dx(1) == cos(x)
true
```

binary input, unary output scalar function:

```jldoctest
julia> x, y = rand(2);

julia> hypotxy, (dx, dy) = rrule(hypot, x, y);

julia> hypotxy == hypot(x, y)
true

julia> dx(1) == (x / hypot(x, y))
true

julia> dy(1) == (y / hypot(x, y))
true
```

See also: [`frule`](@ref), [`AbstractRule`](@ref), [`@scalar_rule`](@ref)
"""
rrule(::Any, ::Vararg{Any}; kwargs...) = nothing


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*--------------------------------------------------------------------
#* Author: Gerson Balbueno Bicca
#* $Id$
#*--------------------------------------------------------------------
type NTU_Basic
	NTU_Basic()=begin
		new(
			positive(Dict{Symbol,Any}(
				:Brief=>"Hot Stream Heat Capacity",
				:Lower=>1e-3,
				:Default=>1e3,
				:Unit=>"W/K",
				:Protected=>true
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Cold Stream Heat Capacity",
				:Lower=>1e-3,
				:Default=>1e3,
				:Unit=>"W/K",
				:Protected=>true
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Heat Capacity Ratio",
				:Default=>0.5,
				:Lower=>1e-6,
				:Protected=>true
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Minimum Heat Capacity",
				:Lower=>1e-10,
				:Default=>1e3,
				:Unit=>"W/K",
				:Protected=>true
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Maximum Heat Capacity",
				:Lower=>1e-10,
				:Default=>1e3,
				:Unit=>"W/K",
				:Protected=>true
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Number of Units Transference",
				:Default=>0.05,
				:Lower=>1e-10,
				:Protected=>true
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Effectiveness",
				:Default=>0.5,
				:Lower=>1e-8,
				:Upper=>1,
				:Symbol=>"\\varepsilon",
				:Protected=>true
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Effectiveness Correction",
				:Lower=>1e-8,
				:Default=>0.5,
				:Symbol=>"\\hat {\\varepsilon}",
				:Protected=>true
			)),
			[:Ch,:Cc,:Cr,:Cmin,:Cmax,:NTU,:Eft,:Eft1,]
		)
	end
	Ch::positive
	Cc::positive
	Cr::positive
	Cmin::positive
	Cmax::positive
	NTU::positive
	Eft::positive
	Eft1::positive
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export NTU_Basic
function atributes(in::NTU_Basic,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Brief]="Number of Units Transference Method."
	fields[:Info]="to be documented"
	drive!(fields,_)
	return fields
end
NTU_Basic(_::Dict{Symbol,Any})=begin
	newModel=NTU_Basic()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
function kkrebook2(omega, imchi, alpha)
    if size(omega,1)>size(omega,2)
        omega = omega'
    end
    if size(imchi,1)>size(imchi,2)
        imchi=imchi'
    end
    g = size(omega,2)
    rechi = zeros(size(imchi))
    a = zeros(size(imchi))
    b = zeros(size(imchi))
    deltaomega = omega[2]-omega[1]
    j = 1
    beta1 = 0

    for k=2:g
        b[1] = beta1+imchi[k]*omega[k]^(2*alpha+1)/(omega[k]^2-omega[1]^2)
        beta1 = b[1]
    end
    rechi[1] = 2/pi*deltaomega*b[1]*omega[1]^(-2*alpha)

    j=g
    alpha1 = 0
    for k=1:g-1
        a[g] = alpha1+imchi[k]*omega[k]^(2*alpha+1)/(omega[k]^2-omega[g]^2)
        alpha1 = a[g]
    end
    rechi[g] = 2/pi*deltaomega*a[g]*omega[g]^(-2*alpha)

    for j=2:g-1
        alpha1 = 0
        beta1 = 0
        for k=1:j-1
            a[j] = alpha1+imchi[k]*omega[k]^(2*alpha+1)/(omega[k]^2-omega[j]^2)
            alpha1 = a[j]
        end
        for k=j+1:g
            b[j] = beta1+imchi[k]*omega[k]^(2*alpha+1)/(omega[k]^2-omega[j]^2)
            beta1 = b[j]
        end
        rechi[j] = 2/pi*deltaomega*(a[j]+b[j])*omega[j]^(-2*alpha)
        print("\r$j")
    end
    return rechi
end

imchi = readdlm(raw"C:\Users\ms246\Dropbox (Heriot-Watt University Team)\RES_EPS_Lupo\Projects\Mohammed\phd\ozone absorption spectrum\imchi_0.124nm.txt")
omega = readdlm(raw"C:\Users\ms246\Dropbox (Heriot-Watt University Team)\RES_EPS_Lupo\Projects\Mohammed\phd\ozone absorption spectrum\freqdata_0.124nm.txt")
kk = kkrebook2(omega[2:end], imchi, 0)

dir = raw"C:\Users\ms246\Dropbox (Heriot-Watt University Team)\RES_EPS_Lupo\Projects\Mohammed\phd\ozone absorption spectrum\\"

open(dir*"julia_kk.txt", "w") do io
    writedlm(io, kk)
end

open(dir*"freqJulia.txt", "w") do io
    writedlm(io, omega)
end

max_stars_repo_name: SabbahMohammed/Dev

max_issues_repo_name: SabbahMohammed/Dev

max_forks_repo_name: SabbahMohammed/Dev


text:
export projPCG

"""
	dm = projPCG(H,g,Active,Precond,cgTol,maxIter)

	Projected Preconditioned Conjugate Gradient method for solving

		H*dm = g    subject to    dm[!Active] == 0

	Input:

		H::Function       - computes action of Hessian
		g::Vector         - right hand side
		Active            - describes active cells
		Precond::Function - preconditioner
		cgTol             - tolerance
		maxIter           - maximum number of iterations
    out               - verbosity level
"""
function projPCG(H::Function,gc::Vector,Active::BitArray,Precond::Function,cgTol::Real,maxIter::Int;out::Int=0)
#  PCG over active cells
    his        = zeros(maxIter,3)
	if norm(gc)==0
    if out >= 2
      println("projPCG finished without iterating: trivial right hand side and solution")
    end
    return zeros(eltype(gc),length(gc)),his[1,:]
  end

	delm       = zeros(eltype(gc),size(gc))
	cgiter     = 0
	resid      = - .!Active.*gc
	normResid0 = norm(resid)
	rdlast = 0.0;
	pc = 0.0;
	while true
		cgiter += 1

		his[cgiter,2]= @elapsed dc = .!Active.*Precond(resid)

		rd = dot(resid,dc)

		#  Compute conjugate direction pc.
		if cgiter == 1
			pc = dc
		else
			betak = rd / rdlast
			pc = dc + betak * pc
		end

		#  Form product Hessian*pc.

		his[cgiter,3] = @elapsed Hp = H(pc)

		Hp = .!Active.*Hp
		#  Update delm and residual.
		alphak = rd / dot(pc,Hp)
		delm   = delm + alphak*pc
		resid  = resid - alphak*Hp
		rdlast = rd

	    his[cgiter,1] = norm(resid)/normResid0;
    if out >=3
      println(@sprintf("projPCG iteration %d: relres = %.3e", cgiter, his[cgiter,1]))
    end
		if (his[cgiter,1] <= cgTol) || (cgiter == maxIter )
	  	  break
		end
	end
  if out >= 2
    if his[cgiter,1] <= cgTol
      println(@sprintf("projPCG converged at iteration %d with relres = %.3e", cgiter, his[cgiter,1]))
    else
      println(@sprintf("projPCG stopped at iteration %d with relres = %.3e", cgiter, his[cgiter,1]))
    end
  end

	return delm,his[1:cgiter,:]
end

function projPCG(gc::Vector,pMis,pInv,sig::Vector,dsig,d2F,d2R,Active;out::Int=0)

		#  Set up Hessian and preconditioner.
		Hs(x) = dsig'*HessMatVec(dsig*x,pMis,sig,d2F) + d2R*x;
		#  build preconditioner
		pInv.HesPrec.setupPrec(Hs, d2R,pInv.HesPrec.param);
		PC(x) = pInv.HesPrec.applyPrec(Hs,d2R,x,pInv.HesPrec.param);
		# call projPCG and return
		return projPCG(Hs,gc,Active,PC,pInv.pcgTol,pInv.pcgMaxIter,out=out)
end


text:
# Comparing with Matt Taddy's gamlr.R
# To rebuild the test cases source(gammalasso.R)
using Lasso
using CSV, GLM, DataFrames, Random, SparseArrays

# often path length is different because of different stopping rules...
function issimilarhead(a::AbstractVector,b::AbstractVector;rtol=1e-4)
    n = min(size(a,1),size(b,1))
    isapprox(a[1:n],b[1:n];rtol=rtol) ? true : [a[1:n] b[1:n]]
end

function issimilarhead(a::AbstractMatrix,b::AbstractMatrix;rtol=1e-4)
    n = min(size(a,1),size(b,1))
    m = min(size(a,2),size(b,2))
    isapprox(a[1:n,1:m],b[1:n,1:m];rtol=rtol) ? true : [a[1:n,1:m] b[1:n,1:m]]
end

datapath = joinpath(dirname(@__FILE__), "data")

penaltyfactors = readcsvmat(joinpath(datapath,"penaltyfactors.csv"))

rtol=1e-2
Random.seed!(243214)
@testset "GammaLassoPath" begin
    @testset "$family" for (family, dist, link) in (("gaussian", Normal(), IdentityLink()), ("binomial", Binomial(), LogitLink()), ("poisson", Poisson(), LogLink()))
        data = readcsvmat(joinpath(datapath,"gamlr.$family.data.csv"))
        y = convert(Vector{Float64},data[:,1])
        X = convert(Matrix{Float64},data[:,2:end])
        (n,p) = size(X)
        @testset "penalty_factor=$pf" for pf = 1:size(penaltyfactors,2)
            penalty_factor = penaltyfactors[:,pf]
            if penalty_factor == ones(p)
                penalty_factor = nothing
            end
            @testset "γ=$γ" for γ in [0 2 10]
                fitname = "gamma$γ.pf$pf"

                # get gamlr.R prms and estimates
                prms = CSV.read(joinpath(datapath,"gamlr.$family.$fitname.params.csv"))
                fittable = CSV.read(joinpath(datapath,"gamlr.$family.$fitname.fit.csv"))
                gcoefs = readcsvmat(joinpath(datapath,"gamlr.$family.$fitname.coefs.csv");types=[Float64 for i=1:100])
                family = prms[1,Symbol("fit.family")]
                γ = prms[1,Symbol("fit.gamma")]
                λ = nothing #convert(Vector{Float64},fittable[!, Symbol("fit.lambda")]) # should be set to nothing evenatually

                # fit julia version
                glp = fit(GammaLassoPath, X, y, dist, link; γ=γ, stopearly=false,
                    λminratio=0.001, penalty_factor=penalty_factor, λ=λ,
                    standardize=false, standardizeω=false)

                # compare
                @test true==issimilarhead(glp.λ,fittable[!, Symbol("fit.lambda")];rtol=rtol)
                @test true==issimilarhead(glp.b0,fittable[!, Symbol("fit.alpha")];rtol=rtol)
                @test true==issimilarhead(convert(Matrix{Float64},glp.coefs'),gcoefs';rtol=rtol)
                # we follow GLM.jl convention where deviance is scaled by nobs, while in gamlr it is not
                @test true==issimilarhead(deviance(glp),fittable[!, Symbol("fit.deviance")]/nobs(glp);rtol=rtol)
                @test true==issimilarhead(deviance(glp,X,y),fittable[!, Symbol("fit.deviance")]/nobs(glp);rtol=rtol)
                # @test true==issimilarhead(round(df(glp)[2:end]),round(fittable[2:end,Symbol("fit.df")]))
                @test true==issimilarhead(loglikelihood(glp),fittable[!, Symbol("fit.logLik")];rtol=rtol)
                @test true==issimilarhead(aicc(glp),fittable[!, Symbol("fit.AICc")];rtol=rtol)

                # TODO: figure out why these are so off, maybe because most are corner solutions
                # and stopping rules for lambda are different
                # # what we really need all these stats for is that the AICc identifies the same minima:
                # if argmin(aicc(glp)) != lastindex(aicc(glp)) && argmin(fittable[!, Symbol("fit.AICc")]) != lastindex(fittable[!, Symbol("fit.AICc")])
                #     # interior minima
                #     println("comparing intereior AICc")
                #     @test argmin(aicc(glp)) == argmin(fittable[!, Symbol("fit.AICc")])
                # end

                # comparse CV, NOTE: this involves a random choice of train subsamples
                gcoefs_CVmin = vec(readcsvmat(joinpath(datapath,"gamlr.$family.$fitname.coefs.CVmin.csv")))
                gcoefs_CV1se = vec(readcsvmat(joinpath(datapath,"gamlr.$family.$fitname.coefs.CV1se.csv")))

                glp_CVmin = coef(glp,MinCVmse(glp, 10))
                glp_CV1se = coef(glp,MinCV1se(glp, 10))

                @test glp_CVmin ≈ gcoefs_CVmin rtol=0.3
                @test glp_CV1se ≈ gcoefs_CV1se rtol=0.3

                if γ==0
                    # Compare with LassoPath
                    lp = fit(LassoPath, X, y, dist, link; stopearly=false,
                        λminratio=0.001, penalty_factor=penalty_factor, λ=λ,
                        standardize=false, standardizeω=false)
                    @test glp.λ == lp.λ
                    @test glp.b0 ≈ lp.b0
                    @test glp.coefs ≈ lp.coefs
                end
            end
        end
    end
end
# ## the following code is useful for understanding comparison failures
#
# rdist(x::Number,y::Number) = abs(x-y)/max(abs(x),abs(y))
# rdist(x::AbstractArray{T}, y::AbstractArray{S}; norm::Function=norm) where {T<:Number,S<:Number} = norm(x - y) / max(norm(x), norm(y))

# (family, dist, link) = (("gaussian", Normal(), IdentityLink()), ("binomial", Binomial(), LogitLink()), ("poisson", Poisson(), LogLink()))[2]
# data = readcsvmat(joinpath(datapath,"gamlr.$family.data.csv"))
# y = convert(Vector{Float64},data[:,1])
# X = convert(Matrix{Float64},data[:,2:end])
# (n,p) = size(X)
# pf = (1:size(penaltyfactors,2))[1]
# penalty_factor = penaltyfactors[:,pf]
# if penalty_factor == ones(p)
#     penalty_factor = nothing
# end
#
# γ = [0 2 10][1]
# fitname = "gamma$γ.pf$pf"
#
# # get gamlr.R prms and estimates
# prms = CSV.read(joinpath(datapath,"gamlr.$family.$fitname.params.csv"))
# fittable = CSV.read(joinpath(datapath,"gamlr.$family.$fitname.fit.csv"))
# gcoefs = readcsvmat(joinpath(datapath,"gamlr.$family.$fitname.coefs.csv");types=[Float64 for i=1:100])
# family = prms[1,Symbol("fit.family")]
# γ = prms[1,Symbol("fit.gamma")]
# λ = nothing #convert(Vector{Float64},fittable[!, Symbol("fit.lambda")]) # should be set to nothing evenatually
#
# # fit julia version
# glp = fit(GammaLassoPath, X, y, dist, link; γ=γ, stopearly=false,
#     λminratio=0.001, penalty_factor=penalty_factor, λ=λ,
#     standardize=false, standardizeω=false)
#
# using InteractiveUtils, BenchmarkTools
# @code_warntype coef(glp; select=:all)
# @code_warntype minAICc(glp)
# @code_warntype segselect(glp, Lasso.MinAIC())
#
# @code_warntype coef(glp)
# @code_warntype coef(glp, AllSeg())
# @code_warntype coef(glp; select=:AICc)
# @code_warntype coef(glp, MinAICc())
#
# @btime coef(glp)
# s = AllSeg()
# @btime coef(glp, s)
# coef(glp; select=:AICc)
# @btime coef(glp; select=:AICc)
# @btime (s = MinAICc())
# s = MinAICc()
# @btime coef(glp, s)
#
# m = fit(LassoModel, X, y, dist, link; stopearly=false,
#     λminratio=0.001, penalty_factor=penalty_factor, λ=λ,
#     standardize=false, standardizeω=false)
# show(coef(m))
# predict(m)
#
# γ = 1
# m = fit(GammaLassoModel, X, y, dist, link; γ=γ, stopearly=false,
#     λminratio=0.001, penalty_factor=penalty_factor, λ=λ,
#     standardize=false, standardizeω=false)
# show(coef(m))
# predict(m)
#
# # compare
# @test true==issimilarhead(glp.λ,fittable[!, Symbol("fit.lambda")];rtol=rtol)
# @test true==issimilarhead(glp.b0,fittable[!, Symbol("fit.alpha")];rtol=rtol)
# @test true==issimilarhead(convert(Matrix{Float64},glp.coefs'),gcoefs';rtol=rtol)
# # we follow GLM.jl convention where deviance is scaled by nobs, while in gamlr it is not
# @test true==issimilarhead(deviance(glp),fittable[!, Symbol("fit.deviance")]/nobs(glp);rtol=rtol)
# @test true==issimilarhead(deviance(glp,X,y),fittable[!, Symbol("fit.deviance")]/nobs(glp);rtol=rtol)
# # @test true==issimilarhead(round(df(glp)[2:end]),round(fittable[2:end,Symbol("fit.df")]))
# @test true==issimilarhead(loglikelihood(glp),fittable[!, Symbol("fit.logLik")];rtol=rtol)
# @test true==issimilarhead(aicc(glp),fittable[!, Symbol("fit.AICc")];rtol=rtol)
#
# # TODO: figure out why these are so off, maybe because most are corner solutions
# # and stopping rules for lambda are different
# # # what we really need all these stats for is that the AICc identifies the same minima:
# # if argmin(aicc(glp)) != lastindex(aicc(glp)) && argmin(fittable[!, Symbol("fit.AICc")]) != lastindex(fittable[!, Symbol("fit.AICc")])
# #     # interior minima
# #     println("comparing intereior AICc")
# #     @test argmin(aicc(glp)) == argmin(fittable[!, Symbol("fit.AICc")])
# # end
#
# # comparse CV, NOTE: this involves a random choice of train subsamples
# gcoefs_CVmin = vec(readcsvmat(joinpath(datapath,"gamlr.$family.$fitname.coefs.CVmin.csv")))
# gcoefs_CV1se = vec(readcsvmat(joinpath(datapath,"gamlr.$family.$fitname.coefs.CV1se.csv")))
#
# @btime coef(glp,select=:CVmin,nCVfolds=10)
# @btime coef(glp, MinCVmse(glp, 10))
# @btime coef(glp, MinCV1se(glp, 10))
# glp_CVmin = coef(glp,select=:CVmin,nCVfolds=10)
#
# glp_CV1se = coef(glp,select=:CV1se,nCVfolds=10)
#
# @test glp_CVmin ≈ gcoefs_CVmin rtol=0.3
# @test glp_CV1se ≈ gcoefs_CV1se rtol=0.3
#
# if γ==0
#     # Compare with LassoPath
#     lp = fit(LassoPath, X, y, dist, link; stopearly=false,
#         λminratio=0.001, penalty_factor=penalty_factor, λ=λ,
#         standardize=false, standardizeω=false)
#     @test glp.λ == lp.λ
#     @test glp.b0 ≈ lp.b0
#     @test glp.coefs ≈ lp.coefs
# end


text:
include("analysis.jl")
include("indicator.jl")
include("rule.jl")
include("signal.jl")
include("strategy.jl")
include("utilities.jl")


text:
using LibSerialPort
using COBS, Statistics
using DataStructures
using AbstractPlotting, WGLMakie, JSServe, Markdown
using AbstractPlotting.MakieLayout
using JSServe: Slider
markdown_css = JSServe.Asset(JSServe.dependency_path("markdown.css"))

top_rpm = 11_500 + 1_150
t4 = 60e6/4
shortest_t = t4/1.1top_rpm
longest_t = 100_000
fps = 30
baudrate = 115200
history = 50
t₀ = time()

gettime() = time() - t₀

function toint(msg)
    y = zero(UInt32)
    for c in msg
        y <<= 8
        y += c
    end
    return y
end

function getrpm(t)
    if iszero(t)
        0.0
    elseif 0 < t ≤ shortest_t
        nothing
    elseif shortest_t < t ≤ longest_t
        t4/t
    else
        0.0
    end
end

last2 = last ∘ last

gety(line, ::Nothing) = last2(line)
gety(_, y) = y

function updatelinesegments!(linesegments, lines)
    empty!(linesegments)
    for line in lines
        for (x,y) in zip(line[1:end - 1], line[2:end])
            push!(linesegments, x => y)
        end
    end
end

struct Arduino
    c::ReentrantLock
    sp::SerialPort
    lines::Vector{CircularBuffer{Point2f0}}
    linesegments::Node{Vector{Pair{Point2f0, Point2f0}}}
    function Arduino(sp)
        lines = Vector{CircularBuffer{Point2f0}}(undef, 3)
        for i in 1:3
            lines[i] = CircularBuffer{Point2f0}(history)
            for j in 1:history
                push!(lines[i], Point2f0(j/fps, 0.0))
            end
        end
        linesegments = Node(Pair{Point2f0, Point2f0}[])
        updatelinesegments!(linesegments[], lines)
        c = ReentrantLock()
        new(c, sp, lines, linesegments)
    end
end

function sample!(a::Arduino)
    tss = lock(a.c) do
        sp_flush(a.sp, SP_BUF_INPUT)
        encode(a.sp, UInt8(0))
        decode(a.sp)
    end
    for (line, ts) in zip(a.lines, Iterators.partition(tss, 4))
        t = toint(ts)
        rpm = getrpm(t)
        y = gety(line, rpm)
        xy = Point2f0(gettime(), y)
        push!(line, xy)
    end
    updatelinesegments!(a.linesegments[], a.lines)
    a.linesegments[] = a.linesegments[]
end

function get_correct_port()
    for port in get_port_list()
        try
            sp = LibSerialPort.open(port, baudrate)
            manufacturer = LibSerialPort.sp_get_port_usb_manufacturer(sp.ref)
            close(sp)
            occursin(r"Arduino", manufacturer) && return port
        catch ex
            @warn ex
        end
    end
    @error "no Arduino connected"
end

port = get_correct_port()
sp = LibSerialPort.open(port, baudrate)
a = Arduino(sp)

reading = @async while isopen(sp)
    cond = Condition()
    Timer(x->notify(cond), 1/fps)
    t = @async sample!(a)
    wait(cond)
    wait(t)
end

pwm = Node(1)
on(pwm) do i
    lock(a.c) do
        encode(a.sp, UInt8(i))
    end
end

function handler(session, request)
    empty!(a.linesegments.listeners)
    scene, layout = layoutscene(0)
    ax = layout[1,1] = LAxis(scene, xlabel = "Time (s)", ylabel = "RPM", title = "Fans")
    linesegments!(ax, a.linesegments, color =  repeat([:red, :green, :blue], inner = history - 1))
    on(a.linesegments) do ps
        xmin = minimum(first ∘ first ∘ first, Iterators.partition(ps, history - 1))
        xmax = maximum(first ∘ last ∘ last, Iterators.partition(ps, history - 1))
        ax.targetlimits[] = FRect2D(xmin, 0, xmax - xmin, top_rpm)
    end
    rpm = lift(a.linesegments) do ps
        round(Int, mean(last ∘ last ∘ last, Iterators.partition(ps, history - 1)))
    end
    slider_s = Slider(1:255, pwm)
    dom = md"""
    $scene
    Speed setting: $pwm $slider_s

    RPM (μ): $rpm
    """
    return JSServe.DOM.div(markdown_css, dom)
end

app = JSServe.Application(handler, "0.0.0.0", 8082)


text:
module MolenaarBasic
using JuMP
using KNITRO

function molenaar_basic(sub_out_con)
    molenaar = Model(
        optimizer_with_attributes(
            KNITRO.Optimizer,
            "ms_enable" => 1,
            "opttol" => 1E-16,
            "opttolabs" => 1e-16,
            "honorbnds" => 1,
            "ms_maxsolves" => 10,
        ),
    )
    JuMP.set_silent(molenaar)

    michaelis_menten(metabolite, protein, kcat, Km) =
        (kcat * protein * metabolite) / (Km + metabolite)
    register(molenaar, :michaelis_menten, 4, michaelis_menten; autodiff = true)

    # Constraints (used the same for now)
    UB = 100
    LB = 1e-8

    @variables molenaar begin
        # growth rate
        LB <= mu <= UB

        # Concentrations
        LB <= sub_out <= UB # substarte outside concentration
        LB <= sub_in <= UB # subtrate inside concentration
        LB <= ribo <= UB # ribosome concentration
        LB <= lipid <= UB # lipid concentration
        LB <= prec <= UB # precursor concentration
        LB <= transp <= UB # transporter concentration
        LB <= metab <= UB # metabolic enzyme concentration
        LB <= lip_syn <= UB # lipid biosynthesis enzyme concentration

        # Rates
        LB <= v_transp <= UB # transporter rate
        LB <= v_metab <= UB  # matabolic rate
        LB <= v_lip <= UB    # lipid biosynthesis rate
        LB <= v_ribo <= UB   # ribosomale rate

        # Cellular constants
        LB <= beta <= UB         # surface to area ratio
        0 <= used_ribo[1:4] <= 1 # ribosomes used for synthesis
    end

    @NLparameters molenaar begin
        kcat_ribo == 3
        Km_ribo == 1
        kcat_transp == 7
        Km_transp == 1
        kcat_metab == 5
        Km_metab == 1
        kcat_lip == 5
        Km_lip == 1
        sub_out == sub_out_con
    end

    @NLconstraints molenaar begin
        v_ribo == michaelis_menten(prec, ribo, kcat_ribo, Km_ribo)
        v_transp == michaelis_menten(sub_out, transp, kcat_transp, Km_transp)
        v_metab == michaelis_menten(sub_in, metab, kcat_metab, Km_metab)
        v_lip == michaelis_menten(prec, lip_syn, kcat_lip, Km_lip)

        beta * (lipid + transp) == 1     # membrane proportions
        sum(used_ribo[i] for i = 1:4) == 1   # ribosome proportions
        metab + ribo + lip_syn <= 1.0      # intracellular density
        transp <= lipid                    # transporter density

        # ribosome activity producing enzymes
        used_ribo[1] * v_ribo - mu * transp == 0
        used_ribo[2] * v_ribo - mu * metab == 0
        used_ribo[3] * v_ribo - mu * ribo == 0
        used_ribo[4] * v_ribo - mu * lip_syn == 0

        # metabolic enzyme balances
        (v_transp - v_metab) - mu * sub_in == 0
        (v_metab - v_lip - v_ribo) - mu * prec == 0
        v_lip - mu * lipid == 0
    end

    @objective(molenaar, Max, mu)
    optimize!(molenaar)

    res = Dict(
        "mu" => value(mu),
        "used_ribo" => value.(used_ribo),
        "sub_out" => value(sub_out),
        "sub_in" => value(sub_in),
        "v_ribo" => value(v_ribo),
    )
    return res
end

end #module


text:
export calculate_sfc_fluxes_energy
export DryBulkFormulaWithRadiation, DryBulkFormula, LinearRelaxation
#surface_fluxes
"""
dummy for SurfaceFluxes.jl
"""
abstract type SurfaceFluxType end

struct LinearRelaxation <: SurfaceFluxType end
struct DryBulkFormula <: SurfaceFluxType end
struct DryBulkFormulaWithRadiation <: SurfaceFluxType end

calculate_sfc_fluxes_energy(forumation::LinearRelaxation, p, θ_sfc, θ1) = p.λ .* (θ_sfc .- θ1)

calculate_sfc_fluxes_energy(forumation::DryBulkFormula, p, θ_sfc, θ1, u_1, v_1, ρ_1) =
    p.Ch * p.C_p * ρ_1 * sqrt(u_1^2 + v_1^2) * (θ_sfc .- θ1)

function calculate_sfc_fluxes_energy(forumation::DryBulkFormulaWithRadiation, parameters, θ_sfc, θ1, u_1, v_1, ρ_1, t)

    p = parameters

    R_SW = (1 - p.α) * p.τ * p.F_sol * (1 .+ sin(t * 2π / p.τ_d))
    R_LW = p.ϵ * (p.σ * θ_sfc .^ 4 - p.σ * θ1 .^ 4)
    SH = p.Ch * p.C_p * ρ_1 * sqrt(u_1^2 + v_1^2) * (θ_sfc - θ1)

    F_tot = -(R_SW - R_LW - SH)
end

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine