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) <: CallableFunction
# 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 end

mutable struct UnivariateZeroState{T,S} <: AbstractUnivariateZeroState where {T,S}
    xn1::T
    xn0::T
    xstar::T
    m::Vector{T}
    fxn1::S
    fxn0::S
    fxstar::S
    fm::Vector{S}
    steps::Int
    fnevals::Int
    stopped::Bool             # stopped, butmay not have converged
    x_converged::Bool         # converged via |x_n - x_{n-1}| < ϵ
    f_converged::Bool         # converged via |f(x_n)| < ϵ
    convergence_failed::Bool
    message::String
end

# Main constructor with no defaults
function _init_state(
    x₀::T, x₁::T, fx₀::S, fx₁::S;
    m=T[],
    fm=S[],
    xstar = nan(T)*oneunit(T),
    fxstar = nan(S)*oneunit(S),
    steps = 0,
    fnevals= 0,
    stopped=false,
    x_converged=false, f_converged=false,
    convergence_failed=false,
    message::String="") where {T,S}

    UnivariateZeroState(x₁, x₀, xstar, m,
                        fx₁, fx₀, fxstar, fm,
                        steps, fnevals,
                        stopped, x_converged, f_converged,
                        convergence_failed,
                        message)

end

# init_state(M, F, x; kwargs...)
# init_state(M, xs..., fs...; kwargs...);
# init_state!(state, M, F)
# init_state!(state, M, F, x0; kwargs...) for backward compatability
#
# The state holds
#
# * the values xₙ₋₁, xₙ and f(xₙ₋₁), f(xₙ) along with
# * fields for the converged value (xstar, fxstar),
# * storage space for some algorithms, m, fm
# * storage to count function evaluations, algorithm steps `steps`, `fnevals`
# * convergence flags
# * a message field
#
# A state is initialized with `init_state(M, F, x)` which sets up xₙ₋₁, xₙ, f(xₙ₋₁), f(xₙ)
# which then calls `init_state!(state, M, F, clear=false, )` to finish the initializatoin
# This should allow a state object for method M to be repurposed for method N by `init_state!(state,N,F)`
#
# !!! note:
# in [NonlinearSolve](https://github.com/SciML/NonlinearSolve.jl), the state object is immutable
# and the `Setfield` package modifes the immutable state
# this may have performance benefits
#
init_state(M::AbstractUnivariateZeroMethod, F, x;
           fnevals=initial_fnevals(M), kwargs...) = nothing
init_state(M::AbstractUnivariateZeroMethod, x₀::T, x₁::T, fx₀::S, fx₁::S; kwargs...) where {T,S} =
    _init_state(x₀::T, x₁::T, fx₀::S, fx₁::S; kwargs...)

function init_state!(state, M::AbstractUnivariateZeroMethod, F; clear=true, kwargs...)
    clear && clear_convergence_flags!(state)

    for (k,v) ∈ kwargs
        setproperty!(state, k, v)
    end

    nothing
end

## This function is used to reset the state to an initial value
## As initializing a state is somewhat costly, this can be useful when many
## function calls would be used.
function init_state!(state, M::AbstractUnivariateZeroMethod, F, x0; clear=true, kwargs...) where {T,S}
    x₀, x₁ = x₀x₁(x0)
    fx₀, fx₁  = first(F(x₀)), first(F(x₁))
    state.xn0, state.fxn0 = x₀, fx₀
    state.xn1, state.fxn1 = x₁, fx₁
    state.fnevals += 2*fn_argout(M)
    init_state!(state, M, F; clear=clear, kwargs...)

    nothing
end

# convert to a different state, e.g. with bracketing
function clear_convergence_flags!(state)
    state.x_converged = state.f_converged = state.stopped = state.convergence_failed = false
    nothing
end

# how many function evaluations in init_state
initial_fnevals(M::AbstractUnivariateZeroState) = @warn "initial_fnevals fix $M"
incfn(o::AbstractUnivariateZeroState, k=1)    = o.fnevals += k
incsteps(o::AbstractUnivariateZeroState, k=1) = o.steps += k
xType(::UnivariateZeroState{T,S}) where {T,S} = T
fxType(::UnivariateZeroState{T,S}) where {T,S} = S

### 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`, the number of
function evaluations is not capped.

"""
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::UnivariateZeroState{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]))
    options
end

# # reset options to default values
@deprecate init_options!(options, M) init_options(M)
# just use options = init_options(M, state) or init_options(M, xType(state), fxType(state); kwargs...)
# function init_options!(options::UnivariateZeroOptions{Q,R,S,T}, M::AbstractUnivariateZeroMethod) where {Q, R, S, T}

#     defs = default_tolerances(M, Q, S)
#     options.xabstol = defs[1]
#     options.xreltol = defs[2]
#     options.abstol = defs[3]
#     options.reltol = defs[4]
#     options.maxevals = defs[5]
#     options.maxfnevals = defs[6]
#     options.strict = defs[7]

#     nothing
# end

## Tracks (for logging actual steps)
## when no logging this should get optimized out to avoid a branch
abstract type AbstractTracks end
struct NullTracks <: AbstractTracks end
# api
log_step(s::NullTracks, M, x, init=false) = nothing
log_step(::Nothing, M, x, init=false) = nothing

mutable struct Tracks{T,S} <: AbstractTracks
xs::Vector{T}
fs::Vector{S}
end
Tracks(s::UnivariateZeroState{T,S}) where {T, S} = Tracks(T[],S[])
Tracks(verbose, tracks, state) = (verbose && isa(tracks, NullTracks)) ? Tracks(state) : tracks

log_step(s::Tracks, M::Any, o, ::Any) = log_step(s, M, o)

function log_step(s::Tracks, M::Any, o)
    push!(s.xs, o.xn1)
    push!(s.fs, o.fxn1)
    nothing
end
function show_tracks(s::Tracks, M::AbstractUnivariateZeroMethod)
    for (i, (xi, fxi)) in enumerate(zip(s.xs, s.fs))
        println(@sprintf("%s = % 18.16f,\t %s = % 18.16f", "x_$(i-1)", float(xi), "fx_$(i-1)", float(fxi)))
    end
    println("")
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 CallableFunction(M,F::Callable_Function, p=F.p)
    CallableFunction(M, F.f, p)
end

# Callable_Function(output_arity, input_arity, F, p)
(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 {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)

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

(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

(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[1])(x, F.p), (F.f[2])(x, F.p), (F.f[3])(x, F.p)
    (f, f/f′, f′/f′′)
end

_apply(f,x) = f(x)
_apply(f,x, p) = f(x, p)

(F::Callable_Function{S,T,𝑭,P})(x) where {𝐍, S <: Val{𝐍}, T <: Val{true}, 𝑭, P<:Nothing} = begin
    fs = _apply.(F.f, Ref(x))
    Tuple(iszero(i) ? fs[1] : fs[i]/fs[i+1] for i ∈ 0:length(fs)-1)
end

(F::Callable_Function{S,T,𝑭,P})(x) where {𝐍, S <: Val{𝐍}, T <: Val{true}, 𝑭, P} = begin
    fs = _apply.(F.f, Ref(x), Ref(p))
    Tuple(iszero(i) ? fs[1] : fs[i]/fs[i+1] for i ∈ 0:length(fs)-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.

If alogrithm hasn't converged returns `false`.

If algorithm has stopped or converged, return `true` and sets one of `state.stopped`, `state.x_converged`,  `state.f_converged`, or `state.convergence_failed`; as well, a message may be set.

* `state.x_converged = true` if `abs(xn1 - xn0) < max(xatol, max(abs(xn1), abs(xn0)) * xrtol)`

* `state.f_converged = true` if  `|f(xn1)| < max(atol, |xn1|*rtol)`

* `state.convergence_failed = true` if xn1 or fxn1 is `NaN` or an infinity

* `state.stopped = true` if the number of steps exceed `maxevals` or the number of function calls exceeds `maxfnevals`.

In `find_zero`, stopped values (and x_converged) are checked for convergence with a relaxed tolerance.

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

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

    if (state.x_converged || state.f_converged) #|| state.stopped)
        if isnan(state.xstar)
            state.xstar, state.fxstar =  xn1, fxn1
        end
        return true
    end

    if isnan(xn1) || isnan(fxn1)
        state.convergence_failed = true
        state.message *= "NaN produced by algorithm. "
        return true
    end

    if isinf(xn1) || isinf(fxn1)
        state.convergence_failed = true
        state.message *= "Inf produced by algorithm. "
        return true
    end

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

    # stop when xn1 ~ xn.
    # in find_zeros there is a check that f could be a zero with a relaxed tolerance
    if abs(xn1 - xn0) < max(options.xabstol, max(abs(xn1), abs(xn0)) * options.xreltol)
        state.xstar, state.fxstar = xn1, fxn1
        state.message *= "x_n ≈ x_{n-1}. "
        state.x_converged = true
        return true
    end

    check_steps_fnevals(state, options) && return true

    return false
end

## too many steps or function evaluations
function check_steps_fnevals(state, options)

    if state.steps > options.maxevals
        state.stopped = true
        state.message *= "Too many steps taken. "
        return true
    end

    if state.fnevals > options.maxfnevals
        state.stopped = true
        state.message *= "Too many function evaluations taken. "
        return true
    end

    return false
end

# state has stopped, this identifies if it has converged
function decide_convergence(M::AbstractUnivariateZeroMethod,  F, state::UnivariateZeroState{T,S}, options) where {T,S}
    xn1 = state.xstar
    fxn1 = state.fxstar

    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

    state.xstar

end

## ----

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

"""

    find_zero(fs, 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.

# 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) implement methods 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 methods.

For more detail, see the help page for each method (e.g., `?Order1`). Many methods are not exported, so much 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 `fs` 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
* `maxfnevals` - limit on maximum number of function evaluations
* `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 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 guaranteed to converge, non-bracketing methods may not.)
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 if the error aswill 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.1415926535897936
```

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> try
       find_zero(fn, x0, Order2(), atol=0.0, rtol=0.0)  # Too many steps taken.
       catch err
       typeof(err)
       end
Roots.ConvergenceFailed

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> try
       find_zero(fn, x0, Order2(), maxevals=3) # Roots.ConvergenceFailed: 26 iterations needed, not 3
       catch err; typeof(err); end
Roots.ConvergenceFailed
```

# Tracing

Passing `verbose=true` will show details on the steps of the algorithm:

```jldoctest find_zero
julia> find_zero(x->sin(x), 3.0, Order2(), verbose=true)   # 2 iterations
Results of univariate zero finding:

* Converged to: 3.1415926535897936
* Algorithm: Order2()
* iterations: 2
* function evaluations: 5
* stopped as |f(x_n)| ≤ max(δ, max(1,|x|)⋅ϵ) using δ = atol, ϵ = rtol

Trace:
x_0 =  3.0000000000000000,	 fx_0 =  0.1411200080598672
x_1 =  3.1425464815525403,	 fx_1 = -0.0009538278181169
x_2 =  3.1415926535897936,	 fx_2 = -0.0000000000000003

3.1415926535897936
```

For more detail on the algorithm, the underlying `state` contains the
number of steps and function evaluations; 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(fs, x0, M::AbstractUnivariateZeroMethod;
                   p = nothing,
                   verbose=false,
                   tracks::AbstractTracks=NullTracks(),
                   kwargs...)

    Z = ZeroProblem(fs, x0)
    ZPI = init(Z, M, p;
               verbose=verbose, tracks=tracks,
               kwargs...)

    xstar = solve!(ZPI)
    verbose && show_trace(M, nothing, ZPI.state, ZPI.logger)
    isnan(xstar) && throw(ConvergenceFailed("Stopped at: xn = $(ZPI.state.xn1). $(ZPI.state.message)"))

    return 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 passed to an iterator to be solved by `find_zero!` or `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

## 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=false,
              tracks = NullTracks(),
              kwargs...)

    F = Callable_Function(M, 𝑭𝑿.F, p === nothing ? p′ : p)  #⁺
    #state = init_state(M, F, 𝑭𝑿.x₀)
    state = init_state(M, F, 𝑭𝑿.x₀)
    options = init_options(M, state; kwargs...)
    l = Tracks(verbose, tracks, state)
    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

# Iteration interface to handle looping
function Base.iterate(P::ZeroProblemIterator, st=nothing)

    M, F, state, options, l = P.M, P.F, P.state, P.options, P.logger
    st == nothing && log_step(l, M, state, :init)

    if assess_convergence(M, state, options)
        decide_convergence(M,F,state,options)
        return  nothing
    end

    update_state(M, F, state, options)
    log_step(l, M, state)
    incsteps(state)
    (state.xn1, false)

end

log_step(P::ZeroProblemIterator) = log_step(P.logger, P.M, P.state)
decide_convergence(P::ZeroProblemIterator) =  decide_convergence(P.M, P.F, P.state, P.options)

"Get last element in the iteration, which is xstar, or throw a warning"
function Base.last(P::ZeroProblemIterator)
#    state.convergence_failed && @warn "The problem failed to converge"
    #    !(state.x_converged || state.f_converged) && @warn "The problem has not converged. Try calling `solve! first"
    P.state.xstar
end

"""
    tracks(P::ZeroProblemIterator)

Show trace of output when `verbose=true` is specified to the problem
"""
function tracks(P::ZeroProblemIterator{M,F,S,O,L}) where {M,F,S,O,L<:Tracks}
     show_trace(P.M, nothing, P.state, P.logger)
end
tracks(P::ZeroProblemIterator) = error("Set verbose=true when specifying the problem to see the tracks")
function show_trace(P::ZeroProblemIterator{M,F,S,O,L}) where {M,F,S,O,L<:Tracks}
     show_trace(P.M, nothing, P.state, P.logger)
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 specificy 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. (Also [`find_zero!`](@ref).)

Returns `NaN`, not an error, when the problem can not be solved.

## Examples:

```
fx = ZeroProblem(sin, 3)
solve(fx)
```

Or, if the iterator is required

```
fx = ZeroProblem(sin, 3)
problem = init(fx)
solve!(fx)
```

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.

```
fx = ZeroProblem(sin,3)
solve(fx, Order5(), atol=1/100)
```

The above is equivalent to:

```
fx = ZeroProblem(sin, 3)
problem = init(fx, Order5(), atol=1/100)
solve!(problem)
```

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

```
f(x,p) = exp(-x) - p # to solve p = exp(-x)
fx = ZeroProblem(f, 1)
solve(fx, 1/2)  # log(2)
```

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; these may be viewed with the method `Roots.tracks` for the iterator.

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.)

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

"""
function solve!(P::ZeroProblemIterator)
    for _ in P end         # iterate to completion
    last(P)
end

## -----
## deprecate this interface at some time.
@deprecate find_zero!(P::ZeroProblemIterator) solve!(P)
# """
#     find_zero!(P::ZeroProblemIterator)

# An alternate interface to `find_zero` whereby a problem is created with `ZeroProblemIterator` and solved
# with `find_zero!`. The generic [`solve!`](@ref) method is recommened for familiarity.

# ```jldoctest find_zero
# julia> using Roots

# julia> P = ZeroProblem(Order1(), sin, 3, verbose=true);

# julia> find_zero!(P)
# 3.141592653589793

# julia> last(P)
# 3.141592653589793

# julia> Roots.tracks(P) # possible when `verbose=true` is specified
# Results of univariate zero finding:

# * Converged to: 3.141592653589793
# * Algorithm: Roots.Secant()
# * iterations: 4
# * function evaluations: 6
# * stopped as |f(x_n)| ≤ max(δ, max(1,|x|)⋅ϵ) using δ = atol, ϵ = rtol

# Trace:
# x_0 =  3.0000000000000000,	 fx_0 =  0.1411200080598672
# x_1 =  3.1425464815525403,	 fx_1 = -0.0009538278181169
# x_2 =  3.1415894805773834,	 fx_2 =  0.0000031730124098
# x_3 =  3.1415926535902727,	 fx_3 = -0.0000000000004795
# x_4 =  3.1415926535897931,	 fx_4 =  0.0000000000000001
# ```
# """
#function find_zero!(P::ZeroProblemIterator)
#    solve!(P)
#end

@deprecate ZeroProblem(M::AbstractUnivariateZeroMethod,fs,x0;kwargs...) init(ZeroProblem(fs, x0), M;kwargs...)

# """
#     ZeroProblem(M, fs, x0; verbose=false, kwargs...)

# Setup an interator interface for the zero problem. Call `find_zero!` or solve! to solve.

# * `M`: Method. A non-hybrid method (not `Order0`).
# * `fs`: a function or for some methods a tuple of functions
# * `x0`: the initial guess
# * `verbose`: if `true`, then calling `Roots.tracks` on the output will show the steps on the algorithm
# * `kwargs`: passed to `Roots.init_options` to adjust tolerances

# """
# function ZeroProblem(M::AbstractUnivariateZeroMethod,
#                      fs,
#                      x0;
#                      verbose=false,
#                      kwargs...)
#     fx = ZeroProblem(fs, x0)
#     problem = init(fx, M; verbose=verbose, kwargs...)

# end


text:
using LinearAlgebra
using FFTW
using Pipe.Pipe

"""
Cooley-Tukey recursive Fast Walsh-Hadamard Transform.
"""
function fwht(x::Array{Float64,1}; recurses=false)
    N = length(x)
    if N == 1
        return x
    else
        X_even = fwht(x[1 : N÷2]; recurses=true)
        X_odd = fwht(x[N÷2+1 : N]; recurses=true)
        res = cat(X_even + X_odd, X_even - X_odd, dims=1)
        return res
    end
end

"""
Takes in a binary array and returns the corresponding integer.
"""
function bin_to_dec(x::AbstractArray{Bool})
    return 2 .^(length(x)-1 : -1 : 0) ⋅ x
end

"""
Takes in an integer and returns the corresponding binary array.
"""
function dec_to_bin(x::Int64, num_bits::Int64)
    @assert x < 2^num_bits "number of bits too small"
    return @pipe string(x, base=2) |> lpad(_, num_bits, "0") |> split(_, "") |> parse.(Bool, _)
end

"""
Returns a matrix where row 'i' is dec_to_bin(i, m),
for i from 0 to 2 ** m - 1.
"""
function binary_ints(m::Int64)
    return @pipe dec_to_bin.(0:2^m-1, m) |> hcat(_...)
end

"""
Sign function that matches the convention (footnote 2 on page 11):
returns 0 for positive reals and 1 for negative reals.
"""
function sign_spright(x)
    return Bool(x < 0)
end

function expected_bin(k::BitArray{1}, M::BitArray{2})
    return @pipe M'k .% 2 |> Bool.(_) |> bin_to_dec(_) + 1
end

ext: jl


text:
"""
    LatentGP(f<:GP, lik, Σy)

 - `f` is a `AbstractGP`.
 - `lik` is the log likelihood function which maps sample from f to corresposing
 conditional likelihood distributions.
 - `Σy` is the observation noise

"""
struct LatentGP{Tf<:AbstractGP,Tlik,TΣy}
    f::Tf
    lik::Tlik
    Σy::TΣy
end

"""
    LatentFiniteGP(fx<:FiniteGP, lik)

 - `fx` is a `FiniteGP`.
 - `lik` is the log likelihood function which maps sample from f to corresposing
 conditional likelihood distributions.

"""
struct LatentFiniteGP{Tfx<:FiniteGP,Tlik}
    fx::Tfx
    lik::Tlik
end

(lgp::LatentGP)(x) = LatentFiniteGP(lgp.f(x, lgp.Σy), lgp.lik)

function Distributions.rand(rng::AbstractRNG, lfgp::LatentFiniteGP)
    f = rand(rng, lfgp.fx)
    y = rand(rng, lfgp.lik(f))
    return (f=f, y=y)
end

"""
    logpdf(lfgp::LatentFiniteGP, y::NamedTuple{(:f, :y)})

```math
    log p(y, f; x)
```
Returns the joint log density of the gaussian process output `f` and real output `y`.
"""
function Distributions.logpdf(lfgp::LatentFiniteGP, y::NamedTuple{(:f, :y)})
    return logpdf(lfgp.fx, y.f) + logpdf(lfgp.lik(y.f), y.y)
end


text:
module EllipticCurves

##############################################################################
##
## TODO
##
##############################################################################

# - Should remove module and rely on package-level module as per convention

##############################################################################
##
## Exported functions and types
##
##############################################################################

export Curve, ConcretePoint, IdealPoint

##############################################################################
##
## Implementation
##
##############################################################################

immutable Curve
  a::Number
  b::Number
  discriminant

  function Curve(a::Number, b::Number)
    d = -16 * (4 * a * a * a + 27 * b * b)
    if d == 0
      error("Curve is singular.")
    end
    new(a, b, d)
  end
end

function show(io::IO, e::Curve)
    print(io, "y^2 = x^3 + ", e.a, "x + ", e.b, "\n")
end

function ==(e::Curve, f::Curve)
  return e.a == f.a && e.b == f.b
end

##############################################################################
##
## Point type definition
##
##############################################################################

abstract Point

type ConcretePoint <: Point
  x::Number
  y::Number
  c::Curve

  function ConcretePoint(x::Number, y::Number, c::Curve)
    y * y == x * x * x + c.a * x + c.b ? new(x, y, c) : error("Point is not on curve.")
  end
end

type IdealPoint <: Point
  c::Curve
  IdealPoint(c::Curve) = new(c)
end

function ==(p::IdealPoint, q::IdealPoint)
  return p.c == q.c
end

function ==(p::IdealPoint, q::ConcretePoint)
  return false
end

function ==(p::ConcretePoint, q::IdealPoint)
  return false
end

function ==(p::ConcretePoint, q::ConcretePoint)
  return p.x == q.x && p.y == q.y && p.c == q.c
end

function -(p::IdealPoint)
  return p
end

function -(p::ConcretePoint)
  return ConcretePoint(p.x,-p.y,p.c)
end

function +(p1::IdealPoint, p2::IdealPoint)
  @assert p1.c == p2.c
  return p1
end

function +(p1::IdealPoint, p2::ConcretePoint)
  @assert p1.c == p2.c
  return p2
end

function +(p1::ConcretePoint, p2::IdealPoint)
  @assert p1.c == p2.c
  return p1
end

function +(p1::ConcretePoint, p2::ConcretePoint)
  @assert p1.c == p2.c
  if p1 == p2
    if p1.y == 0
      return IdealPoint(p1.c)
    end
    m = (3*p1.x*p1.x+p1.c.a) / (2*p1.y)
  else
    if p1.x == p2.x
      return IdealPoint(p1.c)
    end
    m = (p2.y-p1.y) / (p2.x-p1.x)
  end

  new_x = m * m - p2.x - p1.x
  new_y = m * (new_x - p1.x) + p1.y

  return ConcretePoint(new_x, -new_y, p1.c)
end

function -(p1::Point, p2::Point)
  return p1 + (-p2)
end

function *(p::IdealPoint, n::Integer)
  return p
end

function *(n::Integer, p::IdealPoint)
  return p
end

function *(p::ConcretePoint, n::Integer)
  if n == 0
    return IdealPoint(p.c)
  elseif n < 0
    return -p * -n
  else
    q::Point = p
    r::Point = ( n & 1 == 1 ? p : IdealPoint(p.c))
    i = 2
    while i <= n
      q = q+q
      if n & i == i
        r = q+r
      end
      i = i << 1
    end
    return r
  end
end

function *(n::Integer, p::ConcretePoint)
  return p*n
end

end # module EllipticCurves


text:
module number

function gcd(a::T, b::T) where {T <: Integer}
    a, b = (a < b) ? (b, a) : (a, b)
    while b != 0
        c = b
        d = mod(a, b)
        a = c
        b = d
    end
    return a
end

function lcm(a, b)
    (a, b) = (a < b) ? (a, b) : (b, a)
    d = gcd(b, a)
    return div(b, d) * a
end

function sieve_of_eratosthenes(n)
    a = trues(n)
    max_n = ceil(Int, sqrt(Float64(n)))

    for ii = 2:max_n
        if a[ii]
            jj = ii^2
            while jj <= n
                a[jj] = false
                jj += ii
            end
        end
    end
    a[1] = false # for 1
    pr = findall(a)
    return pr
end

end

#=
module ntest
using ..number

a = 30 # 2 * 3 * 5
b = 21 # 3 * 7

d = number.gcd(a, b) # 3
l = number.lcm(a, b) # 7 * 2 * 5 * 3

a = 28 # 2*2 * 7
b = 21 # 3 * 7

d = number.gcd(a, b) # 3
l = number.lcm(a, b) # 2 * 2 * 3 * 7

pr = number.sieve_of_eratosthenes(1000)

end
=#

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    RegularDiscretization(n1, n2, ..., np)

A method to discretize primitive geometries using regular
samples along each parametric dimension. The number of
samples `n1`, `n2`, ..., `np` is passed to [`RegularSampling`](@ref).
"""
struct RegularDiscretization{N} <: DiscretizationMethod
  sizes::Dims{N}
end

RegularDiscretization(sizes::Vararg{Int,N}) where {N} =
  RegularDiscretization(sizes)

function discretize(sphere::Sphere{2,T},
                    method::RegularDiscretization) where {T}
  sz = fitdims(method.sizes, paramdim(sphere))
  nx = sz[1]

  # sample points regularly
  sampler = RegularSampling(nx)
  points  = collect(sample(sphere, sampler))

  # connect regular samples with segments
  topo   = GridTopology(nx-1)
  connec = collect(elements(topo))
  push!(connec, connect((nx,1)))

  SimpleMesh(points, connec)
end

function discretize(sphere::Sphere{3,T},
                    method::RegularDiscretization) where {T}
  nx, ny = fitdims(method.sizes, paramdim(sphere))

  # sample points regularly
  sampler = RegularSampling(nx, ny)
  points  = collect(sample(sphere, sampler))

  # connect regular samples with quadrangles
  topo   = GridTopology((nx-1, ny-1))
  middle = collect(elements(topo))
  offset = nx*ny - nx
  for i in 1:nx-1
    u = offset + i
    v = offset + i + 1
    w = i + 1
    z = i
    quad = connect((u, v, w, z))
    push!(middle, quad)
  end

  # add north and south poles
  c = center(sphere)
  r = radius(sphere)
  e⃗ = Vec{3,T}(0, 0, 1)
  push!(points, c + r*e⃗)
  push!(points, c - r*e⃗)

  # connect north pole with triangles
  north = map(1:ny-1) do j
    u = nx*ny + 1
    v = 1 + (j-1)*nx
    w = 1 + (j  )*nx
    connect((u, v, w))
  end
  u = nx*ny + 1
  v = 1 + (ny-1)*nx
  w = 1
  push!(north, connect((u, v, w)))

  # connect south pole with triangles
  south = map(1:ny-1) do j
    u = nx*ny + 2
    v = (j  )*nx
    w = (j+1)*nx
    connect((u, w, v))
  end
  u = nx*ny + 2
  v = ny*nx
  w = nx
  push!(south, connect((u, w, v)))

  connec = [middle; north; south]

  SimpleMesh(points, connec)
end

function discretize(ball::Ball{2,T},
                    method::RegularDiscretization) where {T}
  nx, ny = fitdims(method.sizes, paramdim(ball))

  # sample points regularly
  sampler = RegularSampling(nx, ny)
  points  = collect(sample(ball, sampler))

  # connect regular samples with quadrangles
  topo   = GridTopology((nx-1, ny-1))
  rings  = collect(elements(topo))
  for j in 1:ny-1
    u = (j  )*nx
    v = (j-1)*nx + 1
    w = (j  )*nx + 1
    z = (j+1)*nx
    quad = connect((u, v, w, z))
    push!(rings, quad)
  end

  # add point at center
  push!(points, center(ball))

  # connect center with triangles
  tris = map(1:nx-1) do i
    u = nx*ny + 1
    v = i + 1
    w = i
    connect((u, v, w))
  end
  u = nx*ny + 1
  v = 1
  w = nx
  push!(tris, connect((u, v, w)))

  connec = [rings; tris]

  SimpleMesh(points, connec)
end


text:
using Discord
using DotEnv

DotEnv.config()
token = ENV["TOKEN"]
prefix = ENV["PREFIX"]

client = Client(token; prefix=prefix)

function handler(c::Client, e::MessageCreate)
    if e.message.content == prefix * "about"
        embed = Embed(
            missing, missing,
            "Hyunjin Choi in the Land of Cute Bois.\nHyunjin Choi is from Camp Buddy: Scoutmasters' Season by [BLits Games](https://www.blitsgames.com/).\nHyunjin was made and developed by: \n**Arch#0226**, **Tetsuki Syu#1250**\nWritten with: \n[Julia](https://julialang.org/) and [Discord.jl](https://github.com/Xh4H/Discord.jl) library.",
            missing, missing, 0x43729F,
            EmbedFooter("Hyunjin Bot: Release 0.2 | 2021-04-06", missing, missing),
            missing,
            EmbedThumbnail("https://cdn.discordapp.com/attachments/811517007446671391/828962181915476008/743164.png", missing, missing, missing),
            missing, missing,
            EmbedAuthor("Hyunjin from Camp Buddy: Scoutmasters' Season", missing, "https://images-ext-2.discordapp.net/external/rOMIWJyMHkRZ76A-26UoPThECFQEIlmVtZph1p24KCI/https/cdn.discordapp.com/emojis/593518771554091011.png", missing),
            missing
        )
        create_message(c, e.message.channel_id; embed=embed)
    end

    if e.message.content == prefix * "ping"
        create_message(c, e.message.channel_id; content="Pong!")
    end
end

add_handler!(client, MessageCreate, handler)
open(client)
wait(client)


text:
using TensorOperations
using KrylovKit
using LinearAlgebra
using Random

function init_random_MPS(d,m,N)

    MPS=Dict()
    MPS[1] = im*rand(1,d,m)
    for i in 2:N-1
        MPS[i]= im*rand(m,d,m)
    end
    MPS[N] = im*rand(m,d,1)
    return MPS
end

function Construct_Ising_MPO(h,J,N)
    d=2
    D=3
    id = [ 1 0 ; 0 1 ]
    sp = [ 0 1 ; 0 0 ]
    sm = [ 0 0 ; 1 0 ]
    sz = [ 1 0 ; 0 -1 ]
    sx = [ 0 1 ; 1 0 ]
    id = [ 1 0 ; 0 1 ]
    W = im *  zeros(D,D,d,d)
    W1 = im *  zeros(1,D,d,d)
    W2 = im *  zeros(D,1,d,d)
    W[1,1,:,:]=id
    W[2,1,:,:]=sz
    W[3,1,:,:]=-h*sx
    W[3,2,:,:]=-J*sz
    W[3,3,:,:]=id;

    W1[1,1,:,:]=-h*sx
    W1[1,2,:,:]=-J*sz
    W1[1,3,:,:]=id

    W2[1,1,:,:]=id
    W2[2,1,:,:]=sz
    W2[3,1,:,:]=-h*sx

    MPO = Dict()
    MPO[1] = W1
    for i in 2:(N-1)
        MPO[i] = W
    end
    MPO[N] = W2

    return MPO
end

function Normalize(MPS)
    norm = MPS_dot(MPS,MPS)
    for i in 1:length(MPS)
        MPS[i] = MPS[i]/norm^(1/(2*length(MPS)))
    end
    return MPS
end

function MPS_dot(A,B)
    E=ones(1,1)
    for i in 1:length(A)
        E = contract_from_left_MPS( A[i] , E , B[i] )
    end
    return E[1]
end

## initial E and F matrices for the left and right vacuum states
function initial_L(W,N)
    sW=size(W)
    L = ones(sW[1],1,1)#im*zeros(sW[1],1,1)
    #L[1] = 1
    #println(size(L))
    return L
end

function initial_R(W,N)
    sW=size(W)
    R = ones(sW[2],1,1)#im*zeros(sW[2],1,1)
    #R[end] = 1
    #println(size(R))
    return R
end

function construct_R(MPS, MPO, Blist, N)
    R = Dict()
    R[N] = initial_R(MPO[N],N)
    for i in N:-1:2
        R[i-1] = contract_from_right(R[i], MPS[i], MPO[i])
    end
    return R
end

function construct_L(Alist, MPO, Blist,N)
    L = Dict()
    L[1] = initial_L(MPO[1],N)
    return L
end

max_stars_repo_name: aanniikkeettr/Tensor-Networks

max_issues_repo_name: aanniikkeettr/Tensor-Networks

max_forks_repo_name: aanniikkeettr/Tensor-Networks


text:
mvn = Sys.iswindows() ? "mvn.cmd" : "mvn"
which = Sys.iswindows() ? "where" : "which"

try
    run(`$which $mvn`)
catch
    error("Cannot find maven. Is it installed?")
end

SPARK_VERSION = get(ENV, "BUILD_SPARK_VERSION", "2.2.0")
SCALA_VERSION = get(ENV, "BUILD_SCALA_VERSION", "2.11.8")
SCALA_BINARY_VERSION = match(r"^\d+\.\d+", SCALA_VERSION).match

cd(joinpath(dirname(@__DIR__), "jvm/sparkjl")) do
    run(`$mvn clean package -Dspark.version=$SPARK_VERSION -Dscala.version=$SCALA_VERSION -Dscala.binary.version=$SCALA_BINARY_VERSION`)
end


text:
using DiffEqBase, StochasticDiffEq, Test

f1(u,p,t) = zero(u)
g1(u,p,t) = ones(size(u))
dt = 1//2^(4)
prob1 = SDEProblem{false}(f1,g1,0.,(0.0,1.0))
integrator = init(prob1,EM(),dt=dt,save_noise=true)

k = integrator.W.dW
@test integrator.W.dW != 0
solve!(integrator)
@test integrator.sol.W[2] == k

prob1 = SDEProblem{false}(f1,g1,zeros(4),(0.0,1.0))
sol = solve(prob1,EM(),dt=dt,save_noise=true)
@test sol.W[2] != zeros(4)

f1(du,u,p,t) = du.=0.
g1(du,u,p,t) = du.=1.
dt = 1//2^(4)
prob1 = SDEProblem(f1,g1,zeros(4),(0.0,1.0))
sol = solve(prob1,EM(),dt=dt,save_noise=true)
@test sol.W[2] != zeros(4)


text:
#=
    Scenario reduction methods
=#
abstract type AbstractScenariosReducer end
abstract type AbstractDimensionReducer end
abstract type AbstractClusteringMethod end

# Manual scenario reduction
mutable struct ManualReducer <: AbstractScenariosReducer
    h::Union{UnitRange{Int64}, Int64}
    y::Union{UnitRange{Int64}, Int64}
    s::Union{UnitRange{Int64}, Int64}
    ManualReducer(; h = 1:8760, y = 2:2, s = 1:1) = new(h, y, s)
end

function reduce(reducer::ManualReducer, ω::Scenarios)
    # Parameters
    h, y, s = reducer.h, reducer.y, reducer.s
    # Initialize
    demands, generations, storages, converters, grids = similar(ω.demands), similar(ω.generations), similar(ω.storages), similar(ω.converters), similar(ω.grids)
    # Demands
    for (k, a) in enumerate(ω.demands)
        demands[k] = (t = a.t[h, y, s], power = a.power[h, y, s])
    end
    # Generations
    for (k, a) in enumerate(ω.generations)
        generations[k] = (t = a.t[h, y, s], power = a.power[h, y, s], cost = a.cost[y, s])
    end
    # Storages
    for (k, a) in enumerate(ω.storages)
        storages[k] = (cost = a.cost[y, s],)
    end
    # Converters
    for (k, a) in enumerate(ω.converters)
        converters[k] = (cost = a.cost[y, s],)
    end
    # Grids
    for (k, a) in enumerate(ω.grids)
        grids[k] = (cost_in =  a.cost_in[h, y, s], cost_out =  a.cost_out[h, y, s])
    end

    return Scenarios(demands, generations, storages, converters, grids), ones(length(s)) / length(s)
end

# Sample Average Approximation scenario reduction
mutable struct SAAReducer <: AbstractScenariosReducer
    nsample::Int64

    SAAReducer(; nsample = 100) = new(nsample)
end

function reduce(reducer::SAAReducer, ω::Scenarios; y::Int64 = 1, s::Int64 = 1)
    # Parameters
    _, ny, ns = size(ω.demands[1].power)
    # Monte carlo indices
    idx = zip(rand(y:ny, reducer.nsample), rand(s:ns, reducer.nsample))
    # Initialize
    demands, generations, storages, converters, grids = similar(ω.demands), similar(ω.generations), similar(ω.storages), similar(ω.converters), similar(ω.grids)
    # Demands
    for (k, a) in enumerate(ω.demands)
        demands[k] = (t = reshape(hcat([a.t[:, y, s] for (y,s) in idx]...),:,1,reducer.nsample), power = reshape(hcat([a.power[:, y, s] for (y,s) in idx]...),:,1,reducer.nsample))
    end
    # Generations
    for (k, a) in enumerate(ω.generations)
        generations[k] = (t = reshape(hcat([a.t[:, y, s] for (y,s) in idx]...),:,1,reducer.nsample), power = reshape(hcat([a.power[:, y, s] for (y,s) in idx]...),:,1,reducer.nsample), cost = repeat(a.cost[y:y, s:s],1,reducer.nsample))
    end
    # Storages
    for (k, a) in enumerate(ω.storages)
        storages[k] = (cost =  repeat(a.cost[y:y, s:s],1,reducer.nsample),)
    end
    # Converters
    for (k, a) in enumerate(ω.converters)
        converters[k] = (cost =  repeat(a.cost[y:y, s:s],1,reducer.nsample),)
    end
    # Grids
    for (k, a) in enumerate(ω.grids)
        grids[k] = (cost_in = reshape(hcat([a.cost_in[:, y, s] for (y,s) in idx]...),:,1,reducer.nsample), cost_out = reshape(hcat([a.cost_out[:, y, s] for (y,s) in idx]...),:,1,reducer.nsample))
    end

    return Scenarios(demands, generations, storages, converters, grids), ones(reducer.nsample) / reducer.nsample
end

# Expected value scenario reduction
mutable struct MeanValueReducer <: AbstractScenariosReducer
    MeanValueReducer() = new()
end

function reduce(reducer::MeanValueReducer, ω::Scenarios; y::Int64 = 1, s::Int64 = 1)
    # Mean value
    # Initialize
    demands, generations, storages, converters, grids = similar(ω.demands), similar(ω.generations), similar(ω.storages), similar(ω.converters), similar(ω.grids)
    # Demands
    for (k, a) in enumerate(ω.demands)
        demands[k] = (t =a.t[:, y:y, s:s], power = mean(a.power, dims=[2,3]))
    end
    # Generations
    for (k, a) in enumerate(ω.generations)
        generations[k] = (t = a.t[:, y:y, s:s], power =  mean(a.power, dims=[2,3]), cost = a.cost[y:y, s:s])
    end
    # Storages
    for (k, a) in enumerate(ω.storages)
        storages[k] = (cost =  a.cost[y:y, s:s],)
    end
    # Converters
    for (k, a) in enumerate(ω.converters)
        converters[k] = (cost =  a.cost[y:y, s:s],)
    end
    # Grids
    for (k, a) in enumerate(ω.grids)
        grids[k] = (cost_in = mean(a.cost_in, dims=[2,3]), cost_out =  mean(a.cost_out, dims=[2,3]))
    end

    return Scenarios(demands, generations, storages, converters, grids), [1.]
end

# Clustering scenario reduction
mutable struct FeatureBasedReducer <: AbstractScenariosReducer
    transformation::Union{UnionAll, Nothing}
    reduction::Union{AbstractDimensionReducer, Nothing}
    clustering::AbstractClusteringMethod

    FeatureBasedReducer(; transformation = UnitRangeTransform,
                        reduction = StatsReduction(),
                        clustering = KmedoidsClustering()) = new(transformation, reduction, clustering)
end

function reduce(reducer::FeatureBasedReducer, ω::Scenarios; y::Int64 = 1, s::Int64 = 1)
    # Parameters
    nh, ny, ns = size(ω.demands[1].power)
    # Initialize
    demands, generations, storages, converters, grids = similar(ω.demands), similar(ω.generations), similar(ω.storages), similar(ω.converters), similar(ω.grids)
    # Formatting
    t_d, t_g = [reshape(a.t[:,2:end,:], nh, :) for a in ω.demands], [reshape(a.t[:,2:end,:], nh, :) for a in ω.generations]
    data_d, data_g, data_gd = [reshape(a.power[:,2:end,:], nh, :) for a in ω.demands], [reshape(a.power[:,2:end,:], nh, :) for a in ω.generations], [reshape(a.cost_in[:,2:end,:], nh, :) for a in ω.grids]
    # Transformation
    norm = replace!.([Genesys.StatsBase.standardize(reducer.transformation, d, dims = 1) for d in vcat(data_d, data_g, data_gd)], NaN => 0.)
    # Dimension reduction
    embedding = replace!(dimension_reduction(reducer.reduction, norm), NaN => 0.)
    # Clustering
    medoids, counts, assignments = clustering(reducer.clustering, embedding)
    # Building reduced scenario
    for (k, a) in enumerate(ω.demands)
        demands[k] = (t = reshape(t_d[k][:,medoids], nh, 1, :), power = reshape(data_d[k][:,medoids], nh, 1, :))
    end
    # Generations
    for (k, a) in enumerate(ω.generations)
        generations[k] = (t = reshape(t_g[k][:,medoids], nh, 1, :), power = reshape(data_g[k][:,medoids], nh, 1, :), cost =  repeat(a.cost[y:y, s:s], 1, length(medoids)))
    end
    # Storages
    for (k, a) in enumerate(ω.storages)
        storages[k] = (cost =  repeat(a.cost[y:y, s:s], 1, length(medoids)),)
    end
    # Converters
    for (k, a) in enumerate(ω.converters)
        converters[k] = (cost =  repeat(a.cost[y:y, s:s], 1, length(medoids)),)
    end
    # Grids
    for (k, a) in enumerate(ω.grids)
        grids[k] = (cost_in = reshape(data_gd[k][:,medoids], nh, 1, :), cost_out = reshape(reshape(a.cost_out[:,2:end,:], nh, :)[:,medoids], nh, 1, :))
    end

    return Scenarios(demands, generations, storages, converters, grids), counts / sum(counts), assignments
end

# Transformation
StatsBase.standardize(DT::Nothing, X; dims=nothing, kwargs...) = X

# Dimension reduction
# PCA
mutable struct PCAReduction <: AbstractDimensionReducer
    n_components::Int64

    PCAReduction(; n_components = 2) = new(n_components)
end

function dimension_reduction(reducer::PCAReduction, data::Array{Array{Float64,2}}; aggregated::Bool=false)
    # data is a vector of d x n matrix with d dimension and n observation
    if aggregated
        m = MultivariateStats.fit(PCA, vcat(data...), maxoutdim = reducer.n_components)
        return MultivariateStats.transform(m, data)
    else
        M = [MultivariateStats.fit(PCA, d, maxoutdim = reducer.n_components) for d in data]
        return vcat([MultivariateStats.transform(M[k], data[k]) for k in 1:length(data)]...)
    end
end

# Statistical moments
struct StatsReduction <: AbstractDimensionReducer end

function dimension_reduction(reducer::StatsReduction, data::Array{Array{Float64,2}}; aggregated::Bool=false)
    # data is a vector of d x n matrix with d dimension and n observation
    if aggregated
        d = vcat(data...)
        # Sum
        s = sum(d, dims = 1)
        # Max
        max = maximum(d, dims = 1)
        # 4 moments
        m = mean(d, dims = 1)
        v = var(d, dims = 1)
        kurt = permutedims([kurtosis(d[:,j]) for j in 1:size(d, 2)])
        skew = permutedims([skewness(d[:,j]) for j in 1:size(d, 2)])
    else
        # Sum
        s = vcat([sum(d, dims = 1) for d in data]...)
        # Max
        max = vcat([maximum(d, dims = 1) for d in data]...)
        # 4 moments
        m = vcat([mean(d, dims = 1) for d in data]...)
        v = vcat([var(d, dims = 1) for d in data]...)
        kurt = vcat([permutedims([kurtosis(d[:,j]) for j in 1:size(d, 2)]) for d in data]...)
        skew = vcat([permutedims([skewness(d[:,j]) for j in 1:size(d, 2)]) for d in data]...)
    end
    # Return aggregated values
    return vcat(s, max, m, v, kurt, skew)
    # return vcat(m, v, kurt, skew)
end

# No reduction
dimension_reduction(reducer::Nothing, data::Array{Array{Float64,2}}) = vcat(data...)

# Clustering methods
# K-medoids
mutable struct KmedoidsClustering <: AbstractClusteringMethod
    n_clusters::Int64
    distance::Distances.SemiMetric
    log::Bool

    KmedoidsClustering(; n_clusters = 20, distance = Distances.Euclidean(), log = true) = new(n_clusters, distance, log)
end

function clustering(method::KmedoidsClustering, embedding::AbstractArray{Float64,2})
    # data is a d x n matrix with d dimension and n observation
    # Distance matrix
    D = pairwise(method.distance, embedding, dims = 2)
    # Clustering
    results = kmedoids(D, method.n_clusters, display = method.log ? :iter : :none)

    return results.medoids, results.counts, results.assignments
end


text:
store = [KPCSupply("beef", 38//10, 36),
         KPCSupply("pork", 54//10, 43),
         KPCSupply("ham", 36//10, 90),
         KPCSupply("greaves", 24//10, 45),
         KPCSupply("flitch", 4//1, 30),
         KPCSupply("brawn", 25//10, 56),
         KPCSupply("welt", 37//10, 67),
         KPCSupply("salami", 3//1, 95),
         KPCSupply("sausage", 59//10, 98)]

sack = solve(store, 15)
println("The store contains:\n - ", join(store, "\n - "))
println("\nThe thief should take::\n - ", join(sack, "\n - "))
@printf("\nTotal value in the sack: %.2f €\n", sum(getfield.(sack, :value)))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# sde deals with the static data export,
# ideally it will download and parse the SDE

module SDE
using SQLite
using Requests
bitstype 32 Region

const sdemd5Addr = "https://www.fuzzwork.co.uk/dump/sqlite-latest.sqlite.bz2.md5"
const sdeAddr = "https://www.fuzzwork.co.uk/dump/sqlite-latest.sqlite.bz2"

global db = SQLite.DB(Pkg.dir("eve", "data", "sqlite-latest.sqlite"))
#increase database read speed at the cost of
#data integrety. We don't care since we never write
#so even without good transaction garentees we're fine
query(db, "PRAGMA synchronous = OFF")
query(db, "PRAGMA journal_mode = OFF")

function itemID(itemName :: AbstractString)
  result = SQLite.query(db, "SELECT typeID FROM invTypes WHERE typeName=?", [itemName])
  get(result.data[1][1])
end
function itemName(itemID :: Int)
  result = SQLite.query(db, "SELECT typeName FROM invTypes WHERE typeID=?", [itemID])
  get(result.data[1][1])
end
function solarSystemID(sysName :: AbstractString)
  result = query(db, "SELECT solarSystemID FROM mapSolarSystems WHERE solarSystemName=?", [sysName])
  get(result.data[1][1])
end

function regionID(regionName :: AbstractString)
  result = query(db, "SELECT regionID FROM mapRegions WHERE regionName=?", [regionName])
  get(result.data[1][1])
end

function regionName(regionID :: Int)
  result = query(db, "SELECT regionName FROM mapRegions WHERE regionID=?", [regionID])
  get(result.data[1][1])
end
"""
Get the market group ID of the group named `marketGroupName`.
Note that group names are not nessassarly unique, for example many group names
exist under both blueprints and somewhere else. This function will
simply return the first groupID (ordered by groupID). This should be
the "real" group and not the blueprint
"""
function marketGroupID(marketGroupName :: AbstractString)
  db = SQLite.DB("sqlite-latest.sqlite")
  result = query(db,
    "select marketGroupID from invMarketGroups where marketGroupName=? order by marketGroupID limit 1",
    [marketGroupName])
  get(result.data[1][1])
end

function getItemsFromGroup(groupID :: Int)
  db = SQLite.DB("sqlite-latest.sqlite")
  q = "WITH RECURSIVE child_groups(marketGroupID) AS "*
      "( Values(?) union select invMarketGroups.marketGroupID from "*
      "invMarketGroups, child_groups where invMarketGroups.parentGroupID=child_groups.marketGroupID) "*
      "select typeid from child_groups join invTypes on invTypes.marketGroupID=child_groups.marketGroupID "*
      "order by typeid"
  result = query(db, q, [groupID])
  result.data[1].values
end

"""
just gets a list of all the market items in the game
"""
function getAllMarketItems()
  db = SQLite.DB("sqlite-latest.sqlite")
  q = "WITH ids as (SELECT marketGroupID from invMarketGroups where hasTypes=1) "*
      "SELECT typeID from invTypes INNER JOIN ids ON invTypes.marketGroupID=ids.marketGroupID"
  result = query(db, q)
  result.data[1].values
end

"""
gets all the itemsIDs on the market in EVE, but leaves out DUST514 related items
because lol they won't sell.
"""
function getAllEveItems()
  itms = getAllMarketItems()
  dust = getItemsFromGroup(marketGroupID("Infantry Gear"))
  setdiff(itms, dust)
end

function Base.string(c :: Region)
  regionName(c)
end

export itemID, solarSystemID, regionID
export itemName
export marketGroupID, getItemsFromGroup
export getAllEveItems, getAllMarketItems
export Region
end


text:
struct PreRenderer <: AbstractRenderSystem end

function Overseer.update(::PreRenderer, m::AbstractLedger)
    glEnable(GL_STENCIL_TEST)
    glStencilOp(GL_KEEP, GL_KEEP, GL_REPLACE)
    glStencilMask(0xff)
    return glClearStencil(0)
end


text:
using LearnConvection
using Plots

D=32
N=4

function animate_error(ℳ, 𝒟)

    xlims = (-0.004, 0.004)
    predi = predict(ℳ, 𝒟; postprocessed=true)

    animation_set = 1:30:(𝒟.Nt-2)
    anim = @animate for i in animation_set
        day_string = string(floor(Int, 𝒟.t[i]/86400))
        plot(𝒟.vavg[i].-predi[i], 𝒟.zavg, legend = :topleft, label = "LES - GP", xlabel = "Temperature [C]", ylabel = "Depth [m]", title = "Error (LES-GP) with stencil size 3, day " * day_string, xlims=xlims)
    end

    return anim
end

train = ["general_strat_4_profiles.jld2", "general_strat_8_profiles.jld2"]
test = ["general_strat_4_profiles.jld2", "general_strat_8_profiles.jld2"]

problem  = Slack("TKE"; parameters=TKEMassFlux.TKEParameters())
# problem  = Sequential("dT")

train = ["general_strat_4_profiles.jld2", "general_strat_8_profiles.jld2", "general_strat_12_profiles.jld2", "general_strat_16_profiles.jld2"]
test = ["general_strat_28_profiles.jld2", "general_strat_32_profiles.jld2"]

k = 3
logγ = -1.7
distance = euclidean_distance
kernel   = get_kernel(k, logγ, 0.0, distance)

𝒟_train  = LearnConvection.Data.data(train, problem; D=D, N=N);
𝒟_test   = LearnConvection.Data.data(test, problem; D=D, N=N);

stencil_size=3
ℳ = LearnConvection.GaussianProcess.model(𝒟_train; kernel = kernel, stencil_size = stencil_size)

anim = animate_error(ℳ, 𝒟_train)

gif(anim, "all.gif")

max_stars_repo_name: adelinehillier/LearnConvection

max_issues_repo_name: adelinehillier/LearnConvection

max_forks_repo_name: adelinehillier/LearnConvection


text:
using KNITRO
using LogisticOptTools
const LOT = LogisticOptTools

function format!(tracer::LOT.Tracer; f♯=Inf)
    tracer.time .-= tracer.time[1]
    if ~isinf(f♯)
        tracer.obj .-= f♯
    end
end

function callbackEvalF(kc, cb, evalRequest, evalResult, dat)
    nfeat = LOT.nfeatures(dat) + dat.fit_intercept
    x = evalRequest.x[1:nfeat]
    evalResult.obj[1] = LOT.loss(x, dat)
    return 0
end

function callbackEvalG!(kc, cb, evalRequest, evalResult, dat)
    nfeat = LOT.nfeatures(dat) + dat.fit_intercept
    ω = evalRequest.x[1:nfeat]
    LOT.gradient!(evalResult.objGrad, ω, dat)
    #= if isa(dat.logger, LOT.Tracer) =#
    #=     push!(dat.logger.grad, norm(evalResult.objGrad, Inf)) =#
    #=     push!(dat.logger.obj, LOT.loss(evalRequest.x[1:n], dat)) =#
    #=     push!(dat.logger.time, time()) =#
    #= end =#
    return 0
end

function callbackEvalHv!(kc, cb, evalRequest, evalResult, dat)
    nfeat = LOT.nfeatures(dat) + dat.fit_intercept
    h = evalRequest.x[1:nfeat]
    vec = evalRequest.vec[1:nfeat]
    LOT.hessvec!(evalResult.hessVec, h, vec, dat)
    return 0
end

function callbackEvalH!(kc, cb, evalRequest, evalResult, dat)
    nfeat = LOT.nfeatures(dat) + dat.fit_intercept
    x = evalRequest.x[1:nfeat]
    LOT.hessian!(evalResult.hess, x, dat)
    return 0
end

function callbackNewPoint(kc, x, lambda_, m)
    dat = m
    # Query information about the current problem.
    dFeasError = KNITRO.KN_get_abs_opt_error(kc)
    if isdefined(dat.logger)
        push!(dat.logger.grad, dFeasError)
        push!(dat.logger.time, time())
    end
    return 0
end

knfit(X, y; x0=zeros(size(X, 2)), penalty=LOT.L2Penalty(0.0), options...) = knfit(LOT.LogitData(X, y), x0, penalty, options)
function knfit(logit, x0::AbstractVector; options...) where P <: Union{LOT.L2Penalty, LOT.L1Penalty}
    p = length(x0)
    kc = KNITRO.KN_new()
    KNITRO.KN_add_vars(kc, p)
    KNITRO.KN_set_var_primal_init_values(kc, x0)
    cb = KNITRO.KN_add_objective_callback(kc, callbackEvalF)
    KNITRO.KN_set_cb_grad(kc, cb, callbackEvalG!)
    # Parse options
    for val in options
        optval = string(val.first)
        if optval == "trace" && val.second == true
            KNITRO.KN_set_newpt_callback(kc, callbackNewPoint)
            continue
        end
        if optval == "hessopt" && val.second == 5
            KNITRO.KN_set_cb_hess(kc, cb, 0, callbackEvalHv!)
            KNITRO.KN_set_param(kc, "algorithm", 2)
        elseif optval == "hessopt" && val.second == 1
            KNITRO.KN_set_cb_hess(kc, cb, KNITRO.KN_DENSE_ROWMAJOR, callbackEvalH!)
        end
        KNITRO.KN_set_param(kc, optval, val.second)
    end

    KNITRO.KN_set_cb_user_params(kc, cb, logit)

    nStatus = @time KNITRO.KN_solve(kc)
    w_opt = KNITRO.get_solution(kc)
    KNITRO.KN_free(kc)
    return w_opt, logit.logger
end

function knfit(dat::LOT.AbstractDataset, x0, penalty::LOT.LinearizedL1Penalty; options...)
    n_features = LOT.dim(dat)
    kc = KNITRO.KN_new()
    KNITRO.KN_add_vars(kc, 2 * LOT.dim(dat))
    index_w = collect(Cint, 0:(n_features-1))
    index_z = collect(Cint, n_features:(2*n_features-1))
    w0 = x0
    z0 = -w0
    KNITRO.KN_set_var_primal_init_values(kc, index_w, w0)
    KNITRO.KN_set_var_primal_init_values(kc, index_z, z0)
    #= KNITRO.KN_set_var_primal_init_values(kc, randn(2*dim(dat))) =#
    cb = KNITRO.KN_add_objective_callback(kc, callbackEvalF)
    KNITRO.KN_set_cb_grad(kc, cb, callbackEvalG!)
    # Parse options
    for val in options
        optval = string(val.first)
        if optval == "trace" && val.second == true
            KNITRO.KN_set_newpt_callback(kc, callbackNewPoint)
            continue
        end
        if optval == "hessopt" && val.second == 5
            KNITRO.KN_set_cb_hess(kc, cb, 0, callbackEvalHv!)
            KNITRO.KN_set_param(kc, "algorithm", 2)
        elseif optval == "hessopt" && val.second == 1
            # Define Hessian sparsity pattern
            xHess1 = Cint[]
            xHess2 = Cint[]
            for i in 1:n_features
                for j in i:n_features
                    push!(xHess1, i-1)
                    push!(xHess2, j-1)
                end
            end
            nnzh = length(xHess2)
            KNITRO.KN_set_cb_hess(kc, cb, nnzh, callbackEvalH!,
                                  hessIndexVars1=xHess1,
                                  hessIndexVars2=xHess2)
        end
        KNITRO.KN_set_param(kc, optval, val.second)
    end

    # Constraints definition
    # 2n linear constrains for l1 norm :
    # 0 <= z_i + w_i,
    # -z_i + w_i <= 0,      for i = 0, ..., n-1
    KNITRO.KN_add_cons(kc, 2 * n_features)

    index_cons = collect(Cint, 0:(2*n_features-1))
    KNITRO.KN_set_con_lobnds(kc, index_cons[1:n_features], zeros(n_features))
    KNITRO.KN_set_con_upbnds(kc, index_cons[n_features+1:end], zeros(n_features))

    index_vars_z = collect(Cint, vcat(n_features:(2*n_features-1), n_features:(2*n_features-1)))
    index_vars_w = collect(Cint, vcat(0:(n_features-1), 0:(n_features-1)))
    coefs_w = ones(2 * n_features)
    coefs_z = vcat(ones(n_features), -ones(n_features))

    KNITRO.KN_add_con_linear_struct(kc, index_cons, index_vars_z, coefs_z)
    KNITRO.KN_add_con_linear_struct(kc, index_cons, index_vars_w, coefs_w)

    # Linear objective
    coefs = fill(penalty.constant, n_features)
    index_obj = collect(Cint, n_features:(2*n_features-1))
    KNITRO.KN_add_obj_linear_struct(kc, index_obj, coefs)

    logit = LOT.GeneralizedLinearModel(dat, penalty, LOT.Tracer{Float64}())
    KNITRO.KN_set_cb_user_params(kc, cb, logit)

    nStatus = @time KNITRO.KN_solve(kc)
    w♯ = KNITRO.get_solution(kc)
    KNITRO.KN_free(kc)
    return w♯, logit.logger
end

function knfit(dat::LOT.AbstractDataset, x0,
               penalty::LOT.L0Penalty;
               big_m=5.0, options...)
    # Formulation choices for (1 - z) w
    # * 0: qcqp
    # * 1: big-M
    # * 2: mpec
    formulation = 1
    # total number of features in dataset
    n_features = LOT.nfeatures(dat)
    # zero vector
    zero_n = zeros(n_features)
    # one vector
    one_n = ones(n_features)

    kc = KNITRO.KN_new()

    # Regularization parameters
    index_w = KNITRO.KN_add_vars(kc, n_features)
    KNITRO.KN_set_var_primal_init_values(kc, index_w, x0)
    KNITRO.KN_set_var_lobnds(kc, index_w, fill(-100.0, n_features))
    KNITRO.KN_set_var_upbnds(kc, index_w, fill(100.0, n_features))

    # Add callbacks
    cb = KNITRO.KN_add_objective_callback(kc, callbackEvalF)
    KNITRO.KN_set_cb_grad(kc, cb, callbackEvalG!)

    # Parse options
    for val in options
        optval = string(val.first)
        if optval == "trace" && val.second == true
            KNITRO.KN_set_newpt_callback(kc, callbackNewPoint)
            continue
        end
        if optval == "hessopt" && val.second == 5
            KNITRO.KN_set_cb_hess(kc, cb, 0, callbackEvalHv!)
            KNITRO.KN_set_param(kc, "algorithm", 2)
        elseif optval == "hessopt" && val.second == 1
            # Define Hessian sparsity pattern
            xHess1 = Cint[]
            xHess2 = Cint[]
            for i in 1:n_features
                for j in i:n_features
                    push!(xHess1, i-1)
                    push!(xHess2, j-1)
                end
            end
            nnzh = length(xHess2)
            KNITRO.KN_set_cb_hess(kc, cb, nnzh, callbackEvalH!,
                                  hessIndexVars1=xHess1,
                                  hessIndexVars2=xHess2)
        else
            KNITRO.KN_set_param(kc, optval, val.second)
        end
    end

    # Add constraint \|w\|_0 <= k
    # Formulate with binary variables z:
    index_z = KNITRO.KN_add_vars(kc, n_features)
    KNITRO.KN_set_var_upbnds(kc, index_z, fill(1.0, n_features))
    KNITRO.KN_set_var_lobnds(kc, index_z, fill(0.0, n_features))
    # Add constraint sum z_i <= k
    l0_cons = KNITRO.KN_add_cons(kc, 1)
    KNITRO.KN_set_con_upbnds(kc, l0_cons, [penalty.constant])
    KNITRO.KN_add_con_linear_struct(kc,
                                    repeat(l0_cons, n_features),
                                    index_z,
                                    one_n)

    if formulation == 0
        # Flag z as binary
        KNITRO.KN_set_var_type.(Ref(kc),
                                index_z,
                                fill(KNITRO.KN_VARTYPE_BINARY, n_features))
        # Add constraint w (1 - z) = 0
        comp_cons = KNITRO.KN_add_cons(kc, n_features)
        KNITRO.KN_set_con_eqbnds(kc, comp_cons, fill(0.0, n_features))
        KNITRO.KN_add_con_linear_struct(kc, comp_cons, index_w, fill(1.0, n_features))
        KNITRO.KN_add_con_quadratic_struct(kc, comp_cons,
                                    index_w,
                                    index_z,
                                    fill(-1.0, n_features))
    elseif formulation == 1
        # Flag z as binary
        KNITRO.KN_set_var_type.(Ref(kc),
                                index_z,
                                fill(KNITRO.KN_VARTYPE_BINARY, n_features))
        # Add Big-M formulation
        comp_cons = KNITRO.KN_add_cons(kc, 2 * n_features)
        KNITRO.KN_set_con_upbnds(kc, comp_cons, fill(0.0, 2* n_features))
        # w <= M z
        up_cons = comp_cons[1:n_features]
        KNITRO.KN_add_con_linear_struct(kc, up_cons, index_w, fill(1.0, n_features))
        KNITRO.KN_add_con_linear_struct(kc, up_cons, index_z, fill(-big_m, n_features))
        # - M z <= w
        lo_cons = comp_cons[n_features+1:end]
        KNITRO.KN_add_con_linear_struct(kc, lo_cons, index_w, fill(-1.0, n_features))
        KNITRO.KN_add_con_linear_struct(kc, lo_cons, index_z, fill(-big_m, n_features))
    elseif formulation == 2
        # Flag z as binary
        KNITRO.KN_set_var_type.(Ref(kc),
                                index_z,
                                fill(KNITRO.KN_VARTYPE_BINARY, n_features))
        # x+
        x_pos = KNITRO.KN_add_vars(kc, n_features)
        KNITRO.KN_set_var_lobnds(kc, x_pos, zero_n)
        # x-
        x_neg = KNITRO.KN_add_vars(kc, n_features)
        KNITRO.KN_set_var_lobnds(kc, x_neg, zero_n)
        # |x|
        x_abs = KNITRO.KN_add_vars(kc, n_features)
        KNITRO.KN_set_var_lobnds(kc, x_abs, zero_n)
        # x = x+ - x-
        cons_xeq = KNITRO.KN_add_cons(kc, n_features)
        KNITRO.KN_set_con_eqbnds(kc, cons_xeq, zero_n)
        KNITRO.KN_add_con_linear_struct(kc,
            repeat(cons_xeq, 3),
            vcat(index_w, x_pos, x_neg),
            vcat(one_n, -one_n, one_n))
        # |x| = x+ + x-
        cons_absx = KNITRO.KN_add_cons(kc, n_features)
        KNITRO.KN_set_con_eqbnds(kc, cons_absx, zero_n)
        KNITRO.KN_add_con_linear_struct(kc,
            repeat(cons_absx, 3),
            vcat(x_abs, x_pos, x_neg),
            vcat(one_n, -one_n, -one_n))
        # Specify complementarities
        # |x| comp z
        cc_types = fill(KNITRO.KN_CCTYPE_VARVAR, 2*n_features)
        KNITRO.KN_set_compcons(kc,
                               cc_types,
                               vcat(x_abs, x_pos),
                               vcat(index_z, x_neg))
    end

    logit = LOT.LogisticRegressor(dat, penalty)
    KNITRO.KN_set_cb_user_params(kc, cb, logit)

    # Resolve!
    nStatus = @time KNITRO.KN_solve(kc)
    # Get optimal solution
    w♯ = KNITRO.get_solution(kc)
    # Free Knitro instance
    KNITRO.KN_free(kc)

    return w♯, logit.logger
end


text:
using Gadfly, DataArrays, DataFrames

a = DataFrame(diff = PooledDataArray(Float64[1,2,3,3,3,4,3,2]))
plot(a, x="diff", Geom.histogram)


text:
# EXPERIMENTAL - Threaded version 3. Max number of tasks = 5 * Threads.nthreads()

function corkendallthreads_v3(X::Union{RealMatrix,RealOrMissingMatrix}, Y::Union{RealMatrix,RealOrMissingMatrix})
    nr = size(X, 2)
    nc = size(Y, 2)
    C = zeros(float(eltype(X)), nr, nc)

    numtasks = min(nc, 10 * Threads.nthreads())
    chunksize = Int(max(1, round(nc / numtasks)))
    numtasks = Int(round(nc / chunksize, RoundUp))

    tasks = Array{Task,1}(undef, numtasks)
    for j = 1:numtasks
        fromcol = (j - 1) * chunksize + 1
        tocol = min(fromcol + chunksize - 1, nc)
        tasks[j] = @spawn corkendall(X, Y[:,fromcol:tocol])
    end

    for j = 1:numtasks
        fromcol = (j - 1) * chunksize + 1
        tocol = min(fromcol + chunksize - 1, nc)
        C[:,fromcol:tocol] = fetch(tasks[j])
    end

    return C
end


text:
# Filesystem

function Base.open(ctx::AbstractContext, filename::AbstractString, mode::AbstractString; kws...)
    io = open(filename, mode; kws...)
    @defer ctx close(io)
    io
end

function Base.open(ctx::AbstractContext, filename::AbstractString; kws...)
    io = open(filename; kws...)
    @defer ctx close(io)
    io
end

function Base.mktemp(ctx::AbstractContext, parent=tempdir())
    (path, io) = mktemp(parent, cleanup=false)
    @defer ctx begin
        close(io)
        rm(path, force=true)
    end
    (path, io)
end

function Base.mktempdir(ctx::AbstractContext, parent=tempdir(); prefix="jl_")
    path = mktempdir(parent, cleanup=false)
    @defer ctx rm(path, recursive=true, force=true)
    path
end

function Base.cd(ctx::AbstractContext, path::AbstractString=homedir())
    oldpath = pwd()
    cd(path)
    @defer ctx cd(oldpath)
    # TODO: `Base.cd(::Function, path)` has some magic here for unix which we
    # could replicate.
end

# Locks

function Base.lock(ctx::AbstractContext, lk::Base.AbstractLock)
    lock(lk)
    @defer ctx unlock(lk)
end

# Standard streams

function Base.redirect_stdout(ctx::AbstractContext, stream)
    prev_stream = stdout
    x = redirect_stdout(stream)
    @defer ctx redirect_stdout(prev_stream)
    x
end

function Base.redirect_stderr(ctx::AbstractContext, stream)
    prev_stream = stderr
    x = redirect_stderr(stream)
    @defer ctx redirect_stderr(prev_stream)
    x
end

function Base.redirect_stdin(ctx::AbstractContext, stream)
    prev_stream = stdin
    x = redirect_stdin(stream)
    @defer ctx redirect_stdin(prev_stream)
    x
end


text:
module TensorTrainFEM

using TensorTrains, LinearAlgebra

export ttsin, ttcos,
    ttrhsdn, ttlaplacedd, ttlaplacedn, ttmassdn,
    ttBPXdnQ, ttBPXdnLambdaId, ttBPXdnLambdaD2, ttBPXdnQLaplace, ttBPXdnC

module Poisson
include("Poisson.jl")
end  # module Poisson

module PoissonAMEn
include("PoissonAMEn.jl")
end  # module PoissonAMEn

module Multiscale
include("Multiscale.jl")
end  # module Multiscale

module MultiscaleAMEn
include("MultiscaleAMEn.jl")
end # module MultiscaleAMEn

eye(n::Integer) = Matrix{Float64}(I, n, n)

# grid values of sine function
function ttsin(L::Int64, ω::Float64, s::Float64 = 0.0)::Tensor
  t = Tensor(2, L, 2)
  t[1][1,:,1] = [1., cos(ω*0.5)]
  t[1][1,:,2] = [0., sin(ω*0.5)]
  for i = 2:(L-1)
    x = 2.0^(-i) * ω
    t[i][1,:,1] = [1., cos(x)]
    t[i][1,:,2] = [0., sin(x)]
    t[i][2,:,1] = [0., -sin(x)]
    t[i][2,:,2] = [1., cos(x)]
  end
  x = 2.0^(-L) * ω
  t[L][1,:,1] = [sin(s*x), sin((1.0+s)*x)]
  t[L][2,:,1] = [cos(s*x), cos((1.0+s)*x)]
  return t
end

# grid values of cosine function
ttcos(L::Int64, ω::Float64, s::Float64 = 0.0)::Tensor = ttsin(L, ω, s + (.5*π)/ω*2.0^L)

# constant right-hand side
# (homogeneous Dirichlet-Neumann boundary conditions, with L^2 scaling)
function ttrhsdn(L::Int64)::Tensor
  f = ttones(L); d = ttdelta(L, [2 for i = 1:L]);
  add!(f, -0.5, d);
  scale!(f, 1. / 2.0^(L/2.))
end

# 1D Dirichlet-Dirichlet Laplacian
function ttlaplacedd(L::Int64)::TensorMatrix
  T = TensorMatrix(2, L, 3)
  T[1][1,:,:,1] = eye(2)
  T[1][1,:,:,2] = [ 0. 0.; 1. 0.]
  T[1][1,:,:,3] = [ 0. 1.; 0. 0.]
  T[1] *= 4.0
  for i = 2:(L-1)
    fill!(T[i], 0.)
    T[i][1,:,:,1] = eye(2)
    T[i][1,:,:,2] = [ 0. 0.; 1. 0.]
    T[i][1,:,:,3] = [ 0. 1.; 0. 0.]
    T[i][2,:,:,2] = [ 0. 1.; 0. 0.]
    T[i][3,:,:,3] = [ 0. 0.; 1. 0.]
    T[i] *= 4.0
  end
  T[L][1,:,:,1] = [ 2. -1.; -1. 2.]
  T[L][2,:,:,1] = [ 0. -1.; 0. 0.]
  T[L][3,:,:,1] = [ 0. 0.; -1. 0.]
  T[L] *= 4.0 * (1.0 + 2.0^(-L))^2
  return T
end

# 1D Dirichlet-Neumann Laplacian
function ttlaplacedn(L::Int64)::TensorMatrix
  T = TensorMatrix(2, L, 4)
  T[1][1,:,:,1] = eye(2)
  T[1][1,:,:,2] = [ 0. 0.; 1. 0.]
  T[1][1,:,:,3] = [ 0. 1.; 0. 0.]
  T[1][1,:,:,4] = [ 0. 0.; 0. 1.]
  T[1] *= 4.0
  for i = 2:(L-1)
    fill!(T[i], 0.)
    T[i][1,:,:,1] = eye(2)
    T[i][1,:,:,2] = [ 0. 0.; 1. 0.]
    T[i][1,:,:,3] = [ 0. 1.; 0. 0.]
    T[i][2,:,:,2] = [ 0. 1.; 0. 0.]
    T[i][3,:,:,3] = [ 0. 0.; 1. 0.]
    T[i][4,:,:,4] = [ 0. 0.; 0. 1.]
    T[i] *= 4.0
  end
  T[L][1,:,:,1] = [ 2. -1.; -1. 2.]
  T[L][2,:,:,1] = [ 0. -1.; 0. 0.]
  T[L][3,:,:,1] = [ 0. 0.; -1. 0.]
  T[L][4,:,:,1] = [ 0. 0.; 0. -1.]
  T[L] *= 4.0
  return T
end

function ttmassdn(L::Int64)::TensorMatrix
  T = TensorMatrix(2, L, 4)
  T[1][1,:,:,1] = eye(2)
  T[1][1,:,:,2] = [ 0. 0.; 1. 0.]
  T[1][1,:,:,3] = [ 0. 1.; 0. 0.]
  T[1][1,:,:,4] = [ 0. 0.; 0. 1.]
  for i = 2:(L-1)
    fill!(T[i], 0.)
    T[i][1,:,:,1] = eye(2)
    T[i][1,:,:,2] = [ 0. 0.; 1. 0.]
    T[i][1,:,:,3] = [ 0. 1.; 0. 0.]
    T[i][2,:,:,2] = [ 0. 1.; 0. 0.]
    T[i][3,:,:,3] = [ 0. 0.; 1. 0.]
    T[i][4,:,:,4] = [ 0. 0.; 0. 1.]
  end
  T[L][1,:,:,1] = [ 2.0/3.0  1.0/6.0; 1.0/6.0 2.0/3.0]
  T[L][2,:,:,1] = [ 0. 1.0/6.; 0. 0.]
  T[L][3,:,:,1] = [ 0. 0.; 1.0/6. 0.]
  T[L][4,:,:,1] = [ 0. 0.; 0. -1.0/3.]
  return T
end

function ttBPXdnC(L::Int64, D::Int64)::TensorMatrix
	c = 2.0.^(-collect(0:L))

	I = eye(2)
	J = [0 1; 0 0]
	I1 = [1 0; 0 0]
	I2 = [0 0; 0 1]
	O = zeros(2,2)

	U = [ I J  J' I2 ;
	      O J' O   O ;
	      O O  J   O ;
	      O O  O  I1 ]
	W = [ 1 2 0 1 0 0 0 0;
	 	    2 4 0 2 1 2 0 1;
				1 0 2 1 0 0 0 0;
				2 0 4 2 1 0 2 1;
				1 2 0 1 2 4 0 2;
				0 0 0 0 1 2 0 1;
				1 0 2 1 2 0 4 2;
				0 0 0 0 1 0 2 1 ] / 4/2
	UD = 1
	WD = 1
	for d ∈ 1:D
		UD = kron(U, UD)
		WD = kron(W, WD)
	end
	sz = repeat([2, 4], 2*D, 1)
	prm = collect(1:D)
	prm = [2*prm; 2*prm .- 1; 2*D .+ 2*prm .- 1; 2*D .+ 2*prm]
	UD = reshape(UD, sz...)
	UD = permutedims(UD, prm)
	UD = reshape(UD, 4^D, 2^D, 2^D, 4^D)
	WD = reshape(WD, sz...)
	WD = permutedims(WD, prm)
	WD = reshape(WD, 4^D, 2^D, 2^D, 4^D)

	W = zeros(2*4^D, 2^D, 2^D, 2*4^D);
	W[1:4^D,:,:,1:4^D] = UD;
	W[1:4^D,:,:,4^D+1:2*4^D] = UD;
	W[4^D+1:2*4^D,:,:,4^D+1:2*4^D] = WD;

	Q = [ copy(W) for l ∈ 1:L ]
	for l ∈ 1:L
		Q[l][1:4^D,:,:,4^D+1:2*4^D] *= c[l+1];
	end
  Q[1] = Q[1][1:1,:,:,:] .+ c[1]*Q[1][4^D+1:4^D+1,:,:,:]
	Q[L] = Q[L][:,:,:,4^D+1:4^D+1]

  return Q
end

function ttBPXdnLambdaD2(L::Int64, D::Int64, K::Int64)::TensorMatrix
	W = eye(2^D) / 2^D
	W = reshape(W, (1,2^D,2^D,1))

	S = [2] / 2
	M = [2 0; 0 2/3] / 2
	V = 1;
	for k=1:D
		if k == K
			V = kron(S,V);
		else
			V = kron(M,V);
		end
	end
	V = reshape(V,(1,2^(D-1),2^(D-1),1))
	Λ  = [ copy(W) for l ∈ 1:L+1 ]
	Λ[L+1] = V;

	return Λ
end

function ttBPXdnLambdaId(L::Int64, D::Int64)::TensorMatrix
	W = eye(2^D) / 2^D
	W = reshape(W, (1,2^D,2^D,1))

	M = [2 0; 0 2/3] / 2
	V = 1;
	for k=1:D
		V = kron(M,V);
	end
	V = reshape(V, (1,2^D,2^D,1))
	Λ  = [ copy(W) for l ∈ 1:L+1 ]
	Λ[L+1] = V;

	return Λ
end

function ttBPXdnQ(L::Int64, D::Int64, K::Int64)::TensorMatrix
		if K ∉ 0:D
			throw(ArgumentError("K should be in 0:D"))
		end
		if L < 2
			throw(ArgumentError("L should be larger than one"))
		end

		if K ∈ 1:D
			c = ones(L+1)
		else
			c = 2.0.^(-collect(0:L))
		end

    U0 = [  1  0  0  1  0  0  0  0 ;
            0  1  0  0  1  0  0  1 ;
            0  0  0  0  0  0  0  0 ;
            0  0  1  0  0  0  0  0 ;
            0  0  0  0  0  1  0  0 ;
            0  0  0  0  0  0  0  0 ;
            0  0  0  0  0  0  1  0 ;
            0  0  0  0  0  0  0  0 ]
    U0 = U0 * sqrt(2) * 1 * 1

    U1 = [  1  0  0  1  0  0  0  0 ;
            0  1  0  0  1  0  0  1 ;
            0  0  0  0  0  0  0  0 ;
            0  0  1  0  0  0  0  0 ;
            0  0  0  0  0  1  0  0 ;
            0  0  0  0  0  0  0  0 ;
            0  0  0  0  0  0  1  0 ;
            0  0  0  0  0  0  0  0 ]
    U1 = U1 * sqrt(2) * 2 * 1 / 2

    Y0 = [  2  4  0  2  0  0  0  0 ;
            2  4  0  2  0  0  0  0 ;
            2  0  4  2  0  0  0  0 ;
            2  0  4  2  0  0  0  0 ;
           -1 -2  0 -1  1  2  0  1 ;
            1  2  0  1  1  2  0  1 ;
           -1  0 -2 -1  1  0  2  1 ;
            1  0  2  1  1  0  2  1 ] / 2 / 2
    Y0 = Y0 * sqrt(2) * 1 / sqrt(2) / sqrt(2)

    Y1 = [  1  2  0  1 ;
            1  2  0  1 ;
            1  0  2  1 ;
            1  0  2  1 ] / 2 / 2
    Y1 = Y1 * sqrt(2) * 2 / sqrt(2) / sqrt(2)

    UT0 = [  1  0  0  1  1  0  0  1 ;
             1  1  0  1 -1  1  0 -1 ;
             0  0  0  0  0  0  0  0 ;
             0  0  1  0  0  0  1  0 ;
             0  1  0  0  0 -1  0  0 ;
             0  0  0  0  0  0  0  0 ;
             0  0  1  0  0  0 -1  0 ;
             0  0  0  0  0  0  0  0 ]
    UT0 = UT0 * sqrt(2) * 1 * 1

    UT1 = [  1  0  0  1 ;
            -1  1  0 -1 ;
             0  0  0  0 ;
             0  0  1  0 ;
             0 -1  0  0 ;
             0  0  0  0 ;
             0  0 -1  0 ;
             0  0  0  0 ]
    UT1 = UT1 * sqrt(2) * 2 * 1 / 2

    TY0 = [  1  2  0  1  1  2  0  1 ;
             3  6  0  3  1  2  0  1 ;
             1  0  2  1  1  0  2  1 ;
             3  0  6  3  1  0  2  1 ;
             3  6  0  3 -1 -2  0 -1 ;
             1  2  0  1 -1 -2  0 -1 ;
             3  0  6  3 -1  0 -2 -1 ;
             1  0  2  1 -1  0 -2 -1 ] / 2 / 2
    TY0 = TY0 * sqrt(2) * 1 / sqrt(2) / sqrt(2)

    TY1 = [  1  2  0  1 ;
             1  2  0  1 ;
             1  0  2  1 ;
             1  0  2  1 ;
            -1 -2  0 -1 ;
            -1 -2  0 -1 ;
            -1  0 -2 -1 ;
            -1  0 -2 -1 ] / 2 / 2
    TY1 = TY1 * sqrt(2) * 2 / sqrt(2) / sqrt(2)

    M0 = [ 1 ;
           0 ;
           0 ;
           0 ;
           0 ;
           1 ;
           0 ;
           0 ] / 2

    M1 = [ 1 ;
           0 ] / 2
    M1 = M1 * 2

	U  = 1
	UT = 1
	TY = 1
	Y  = 1
	M  = 1
	for k ∈ 1:D
		if k == K
			U  = kron(U1, U)
			UT = kron(UT1, UT)
			TY = kron(TY1, TY)
			Y  = kron(Y1, Y)
			M  = kron(M1, M)
		else
			U  = kron(U0, U)
			UT = kron(UT0, UT)
			TY = kron(TY0, TY)
			Y  = kron(Y0, Y)
			M  = kron(M0, M)
		end
	end

	prm = collect(1:D)
	prm = [2*prm; 2*prm .- 1; 2*D .+ 2*prm .- 1; 2*D .+ 2*prm]

	sz = repeat([2, 4], 2*D, 1)
  U = reshape(U, sz...)
	U = permutedims(U, prm)
	sz = sz[collect(prm)]
	sz = prod(reshape(sz, (D, 4)), dims=1)
	U = reshape(U, sz...)

	sz = repeat([2, 4], 2*D, 1)
	if K ∈ 1:D
		sz[2*D+2*K] = 2
	end
  UT = reshape(UT, sz...)
	UT = permutedims(UT, prm)
	sz = sz[collect(prm)]
	sz = prod(reshape(sz, (D, 4)), dims=1)
	UT = reshape(UT, sz...)

	sz = repeat([2, 4], 2*D, 1)
	if K ∈ 1:D
		sz[2*D+2*K] = 2
	end
  TY = reshape(TY, sz...)
	TY = permutedims(TY, prm)
	sz = sz[collect(prm)]
	sz = prod(reshape(sz, (D, 4)), dims=1)
	TY = reshape(TY, sz...)

	sz = repeat([2, 4], 2*D, 1)
	if K ∈ 1:D
		sz[2*K] = 2
		sz[2*D+2*K] = 2
	end
  Y = reshape(Y, sz...)
	Y = permutedims(Y, prm)
	sz = sz[collect(prm)]
	sz = prod(reshape(sz, (D, 4)), dims=1)
	Y = reshape(Y, sz...)

	sz = repeat([2, 4], D, 1)
	sz = [sz; ones(Int64, 2*D, 1)]
	if K ∈ 1:D
		sz[2*K-1] = 1
		sz[2*K] = 2
	end
  M = reshape(M, sz...)
	M = permutedims(M, prm)
	sz = sz[collect(prm)]
	sz = prod(reshape(sz, (D, 4)), dims=1)
	M = reshape(M, sz...)

	Q = TensorMatrix(L+1)
	Q[1] = cat(U[1:1,:,:,:], c[2]*UT[1:1,:,:,:] .+ c[1]*TY[1:1,:,:,:], dims = 4);
	O = zeros(size(Y, 1), 2^D, 2^D, size(U, 4))
	for ℓ ∈ 2:L-1
		Q[ℓ]= cat(cat(U, O, dims=1), cat(c[ℓ+1]*UT, Y, dims=1), dims=4)
	end
	Q[L] = cat(1, c[L+1]*UT, Y)
	Q[L+1] = M

  return Q
end

function ttBPXdnQLaplace(L::Int64, D::Int64, K::Int64)::TensorMatrix
		if K ∉ 0:D
			throw(ArgumentError("K should be in 0:D"))
		end
		if L < 2
			throw(ArgumentError("L should be larger than one"))
		end

		if K ∈ 1:D
			c = ones(L+1)
		else
			c = 2.0.^(-collect(0:L))
		end

    U0 = [  1  0  0  1  0  0  0  0 ;
            0  1  0  0  1  0  0  1 ;
            0  0  0  0  0  0  0  0 ;
            0  0  1  0  0  0  0  0 ;
            0  0  0  0  0  1  0  0 ;
            0  0  0  0  0  0  0  0 ;
            0  0  0  0  0  0  1  0 ;
            0  0  0  0  0  0  0  0 ]
    U0 = U0 * sqrt(2) * 1 * 1 / sqrt(2)

    U1 = [  1  0  0  1  0  0  0  0 ;
            0  1  0  0  1  0  0  1 ;
            0  0  0  0  0  0  0  0 ;
            0  0  1  0  0  0  0  0 ;
            0  0  0  0  0  1  0  0 ;
            0  0  0  0  0  0  0  0 ;
            0  0  0  0  0  0  1  0 ;
            0  0  0  0  0  0  0  0 ]
    U1 = U1 * sqrt(2) * 2 * 1 / 2 / sqrt(2)

    Y0 = [  2  4  0  2  0  0  0  0 ;
            2  4  0  2  0  0  0  0 ;
            2  0  4  2  0  0  0  0 ;
            2  0  4  2  0  0  0  0 ;
           -1 -2  0 -1  1  2  0  1 ;
            1  2  0  1  1  2  0  1 ;
           -1  0 -2 -1  1  0  2  1 ;
            1  0  2  1  1  0  2  1 ] / 2 / 2
    Y0 = Y0 * sqrt(2) * 1 / sqrt(2) / sqrt(2) / sqrt(2)

    Y1 = [  1  2  0  1 ;
            1  2  0  1 ;
            1  0  2  1 ;
            1  0  2  1 ] / 2 / 2
    Y1 = Y1 * sqrt(2) * 2 / sqrt(2) / sqrt(2) / sqrt(2)

    UT0 = [  1  0  0  1  1  0  0  1 ;
             1  1  0  1 -1  1  0 -1 ;
             0  0  0  0  0  0  0  0 ;
             0  0  1  0  0  0  1  0 ;
             0  1  0  0  0 -1  0  0 ;
             0  0  0  0  0  0  0  0 ;
             0  0  1  0  0  0 -1  0 ;
             0  0  0  0  0  0  0  0 ]
    UT0 = UT0 * sqrt(2) * 1 * 1 / sqrt(2)

    UT1 = [  1  0  0  1 ;
            -1  1  0 -1 ;
             0  0  0  0 ;
             0  0  1  0 ;
             0 -1  0  0 ;
             0  0  0  0 ;
             0  0 -1  0 ;
             0  0  0  0 ]
    UT1 = UT1 * sqrt(2) * 2 * 1 / 2 / sqrt(2)

    TY0 = [  1  2  0  1  1  2  0  1 ;
             3  6  0  3  1  2  0  1 ;
             1  0  2  1  1  0  2  1 ;
             3  0  6  3  1  0  2  1 ;
             3  6  0  3 -1 -2  0 -1 ;
             1  2  0  1 -1 -2  0 -1 ;
             3  0  6  3 -1  0 -2 -1 ;
             1  0  2  1 -1  0 -2 -1 ] / 2 / 2
    TY0 = TY0 * sqrt(2) * 1 / sqrt(2) / sqrt(2) / sqrt(2)

    TY1 = [  1  2  0  1 ;
             1  2  0  1 ;
             1  0  2  1 ;
             1  0  2  1 ;
            -1 -2  0 -1 ;
            -1 -2  0 -1 ;
            -1  0 -2 -1 ;
            -1  0 -2 -1 ] / 2 / 2
    TY1 = TY1 * sqrt(2) * 2 / sqrt(2) / sqrt(2) / sqrt(2)

    M0 = [ 1         ;
           0         ;
           0         ;
           0         ;
           0         ;
           1/sqrt(3) ;
           0         ;
           0         ] / 2

    M1 = [ 1 ;
           0 ] / 2
    M1 = M1 * 2

	U  = 1
	UT = 1
	TY = 1
	Y  = 1
	M  = 1
	for k ∈ 1:D
		if k == K
			U  = kron(U1, U)
			UT = kron(UT1, UT)
			TY = kron(TY1, TY)
			Y  = kron(Y1, Y)
			M  = kron(M1, M)
		else
			U  = kron(U0, U)
			UT = kron(UT0, UT)
			TY = kron(TY0, TY)
			Y  = kron(Y0, Y)
			M  = kron(M0, M)
		end
	end

	prm = collect(1:D)
	prm = [2*prm; 2*prm .- 1; 2*D .+ 2*prm .- 1; 2*D .+ 2*prm]

	sz = repeat([2, 4], 2*D, 1)
  U = reshape(U, sz...)
	U = permutedims(U, prm)
	sz = sz[collect(prm)]
	sz = prod(reshape(sz, (D, 4)), dims=1)
	U = reshape(U, sz...)

	sz = repeat([2, 4], 2*D, 1)
	if K ∈ 1:D
		sz[2*D+2*K] = 2
	end
  UT = reshape(UT, sz...)
	UT = permutedims(UT, prm)
	sz = sz[collect(prm)]
	sz = prod(reshape(sz, (D, 4)), dims=1)
	UT = reshape(UT, sz...)

	sz = repeat([2, 4], 2*D, 1)
	if K ∈ 1:D
		sz[2*D+2*K] = 2
	end
  TY = reshape(TY, sz...)
	TY = permutedims(TY, prm)
	sz = sz[collect(prm)]
	sz = prod(reshape(sz, (D, 4)), dims=1)
	TY = reshape(TY, sz...)

	sz = repeat([2, 4], 2*D, 1)
	if K ∈ 1:D
		sz[2*K] = 2
		sz[2*D+2*K] = 2
	end
  Y = reshape(Y, sz...)
	Y = permutedims(Y, prm)
	sz = sz[collect(prm)]
	sz = prod(reshape(sz, (D, 4)), dims=1)
	Y = reshape(Y, sz...)

	sz = repeat([2, 4], D, 1)
	sz = [sz; ones(Int64, 2*D, 1)]
	if K ∈ 1:D
		sz[2*K-1] = 1
		sz[2*K] = 2
	end
  M = reshape(M, sz...)
	M = permutedims(M, prm)
	sz = sz[collect(prm)]
	sz = prod(reshape(sz, (D, 4)), dims=1)
	M = reshape(M, sz...)

	Q = TensorMatrix(L+1)
	Q[1] = cat(U[1:1,:,:,:], c[2]*UT[1:1,:,:,:]+c[1]*TY[1:1,:,:,:], dims=4);
	O = zeros(size(Y, 1), 2^D, 2^D, size(U, 4))
	for ℓ ∈ 2:L-1
		Q[ℓ]= cat(cat(U, O, dims=1), cat(c[ℓ+1]*UT, Y, dims=1), dims=4)
	end
	Q[L] = cat(c[L+1]*UT, Y, dims=1)
	Q[L+1] = M

  return Q
end

end # module


text:
# Load Julia packages (libraries) needed  for the snippets in chapter 0

using DynamicHMCModels

ProjDir = @__DIR__
cd(ProjDir)

# Read in data

df = DataFrame(CSV.read(joinpath("..", "..", "data", "rugged.csv"), delim=';'))
df = filter(row -> !(ismissing(row[:rgdppc_2000])), df)
df[!, :log_gdp] = log.(df[:, :rgdppc_2000])
df[!, :cont_africa] = Array{Float64}(convert(Array{Int}, df[:, :cont_africa]))

Base.@kwdef mutable struct RuggedModel{Ty <: AbstractVector,
  Tx <: AbstractMatrix}
    "Observations."
    y::Ty
    "Covariates"
    x::Tx
end

# Write a function to return properly dimensioned transformation.

function make_transformation(model::RuggedModel)
  as((β = as(Array, size(model.x, 2)), σ = asℝ₊))
end

# Instantiate the model with data and inits.

x = hcat(ones(size(df, 1)), df[:, :rugged], df[:, :cont_africa],
  df[:, :rugged] .* df[:, :cont_africa]);
model = RuggedModel(;y=df[:, :log_gdp], x=x)

# Model callable with *a single argument*.

function (problem::RuggedModel)(θ)
    @unpack y, x = problem   # extract the data
    @unpack β, σ = θ            # works on the named tuple too
    ll = 0.0
    ll += logpdf(Normal(0, 100), x[1])
    ll += logpdf(Normal(0, 10), x[2])
    ll += logpdf(Normal(0, 10), x[3])
    ll += logpdf(Normal(0, 10), x[4])
    ll += logpdf(TDist(1.0), σ)
    ll += loglikelihood(Normal(0, σ), y .- x*β)
    ll
end

println()
model((β = [1.0, 2.0, 1.0, 2.0], σ = 1.0))
println()

# Wrap the problem with a transformation, then use Flux for the gradient.

P = TransformedLogDensity(make_transformation(model), model)
∇P = ADgradient(:ForwardDiff, P);

# Tune and sample.

results = mcmc_with_warmup(Random.GLOBAL_RNG, ∇P, 1000)
posterior = P.transformation.(results.chain)

println()
DynamicHMC.Diagnostics.EBFMI(results.tree_statistics) |> display
println()

DynamicHMC.Diagnostics.summarize_tree_statistics(results.tree_statistics) |> display
println()

a3d = Array{Float64, 3}(undef, 1000, 5, 1);
for j in 1:1
  for i in 1:1000
    a3d[i, 1:4, j] = values(posterior[i].β)
    a3d[i, 5, j] = values(posterior[i].σ)
  end
end

pnames = ["β[1]", "β[2]", "β[3]", "β[4]", "σ"]
sections =   Dict(
  :parameters => pnames,
)
chns = create_mcmcchains(a3d, pnames, sections, start=1);
chns = set_names(chns, Dict("β[1]" => "α", "β[2]" => "β[1]",
  "β[3]" => "β[2]", "β[4]" => "β[3]"))

# Result rethinking

rethinking = "
       mean   sd  5.5% 94.5% n_eff Rhat
a      9.22 0.14  9.00  9.46   282    1
bR    -0.21 0.08 -0.33 -0.08   275    1
bA    -1.94 0.24 -2.33 -1.59   268    1
bAR    0.40 0.14  0.18  0.62   271    1
sigma  0.96 0.05  0.87  1.04   339    1
"

# Summary

describe(chns)

# End of `08/m8.1s.jl`


text:
using AltDistributions, Test, LinearAlgebra, Distributions, Random
using ForwardDiff: jacobian

Random.seed!(1)

# utilities

"Elements below the diagonal as a vector (by row)."
lower_to_vec(L) = vcat((L[i, 1:(i-1)] for i in 1:size(L, 1))...)

"Elements below and including the diagonal as a vector (by row)."
lowerdiag_to_vec(L) = vcat((L[i, 1:i] for i in 1:size(L, 1))...)

"""
Reconstruct lower triangular matrix (including diagonal) from a vector of
elemenst (by row).
"""
function vec_to_lowerdiag(l::AbstractVector{T}, n::Int) where T
    A = zeros(T, n, n)
    cumulative_index = 0
    for i in 1:n
        A[i, 1:i] .= l[cumulative_index .+ (1:i)]
        cumulative_index += i
    end
    LowerTriangular(A)
end

@testset "vec lower utilities" begin
    l = Float64.(1:6)
    L = vec_to_lowerdiag(l, 3)
    @test size(L) == (3, 3)
    @test eltype(L) == eltype(l)
    @test L == [1.0 0.0 0.0;
                2.0 3.0 0.0;
                4.0 5.0 6.0]
    @test lowerdiag_to_vec(L) == l
    @test lower_to_vec(L) == [2.0, 4.0, 5.0]
end

"Random covariance matrix."
randΣ(n) = (A = randn(n, n); A'*A)

"""
Random correlation matrix `Ω` and its Cholesky factor `F` such that `FF'=Ω`. The standard
deviations `σ` are returned as the 3rd value.
"""
function randΩFσ(n)
    Σ = randΣ(n)
    σ = .√diag(Σ)
    d = Diagonal(σ)
    Ω = Symmetric(d \ Σ / d)
    F = cholesky(Ω).L
    Ω, F, σ
end

# tests

@testset "LKJL" begin
    n = 7
    Ω, F, _ = randΩFσ(n)
    η = abs2(randn())
    J_full = jacobian(lowerdiag_to_vec(F)) do z
        F = vec_to_lowerdiag(z, n)
        lowerdiag_to_vec(F*F')
    end
    @test logdet(J_full) + logdet(Ω)*(η-1) ≈ logpdf(LKJL(η), F)
end

@testset "AltMvNormal" begin
    ns = 3:7

    @testset "plain vanilla covariance matrix" begin
        for _ in 1:1000
            n = rand(ns)
            μ = randn(n)
            Σ = randΣ(n)
            x = randn(n)
            @test logpdf(MvNormal(μ, Σ), x) ≈ logpdf(AltMvNormal(μ, Symmetric(Σ)), x)
        end
    end

    @testset "Σ = I" begin
        for _ in 1:1000
            n = rand(ns)
            μ = randn(n)
            x = randn(n)
            @test logpdf(MvNormal(μ, Diagonal(ones(n))), x) ≈ logpdf(AltMvNormal(μ, I), x)
        end
    end

    @testset "Cholesky factorization" begin
        for _ in 1:1000
            n = rand(ns)
            μ = randn(n)
            x = randn(n)
            Σ = randΣ(n)
            @test logpdf(MvNormal(μ, Σ), x) ≈ logpdf(AltMvNormal(μ, cholesky(Σ)), x)
        end
    end

    @testset "Diagonal Σ" begin
        for _ in 1:1000
            n = rand(ns)
            μ = randn(n)
            x = randn(n)
            d = abs2.(randn(n))
            Σ = Diagonal(d)
            @test logpdf(MvNormal(μ, Σ), x) ≈ logpdf(AltMvNormal(μ, Σ), x)
        end
    end

    @testset "plain vanilla covariance matrix" begin
        for _ in 1:1000
            n = rand(ns)
            μ = randn(n)
            Ω, F, σ = randΩFσ(n)
            L = Diagonal(σ)*F
            Σ = L*L'
            G = StdCorrFactor(σ, F)
            x = randn(n)
            @test logpdf(MvNormal(μ, Σ), x) ≈ logpdf(AltMvNormal(Val(:L), μ, G), x) rtol = 1e-6
        end
    end

    @testset "errors" begin
        @test_throws ArgumentError AltMvNormal(ones(3), Diagonal(ones(2)))
        @test_throws DimensionMismatch logpdf(AltMvNormal(ones(3), Diagonal(ones(3))), ones(2))
    end

    @testset "random draws" begin
        N = 100000
        μ = [1.0, 2.0]
        Σ = Symmetric([1.0 0.5;
                       0.5 4.0])
        d = AltMvNormal(μ, Σ)
        zs = reduce(hcat, rand(d, N))
        @test vec(mean(zs, dims = 2)) ≈ μ atol = 0.01
        @test cov(zs, dims = 2) ≈ Σ atol = 0.01
    end
end

@testset "AltMultinomial" begin
    n = 10
    p = [0.1, 0.2]
    p′ = [0.25, 0.4]
    a = AltMultinomial(n, p)
    a′ = AltMultinomial(n, p′)
    m = Multinomial(n, vcat(p, [1-sum(p)]))
    m′ = Multinomial(n, vcat(p′, [1-sum(p′)]))
    for _ in 1:100
        x = rand(a)
        @test sum(x::Vector{Int}) == n
        @test length(x) == length(p) + 1
        @test logpdf(a, x) ≈ logpdf(m, x)
        @test logpdf(a, Fixed(x)) - logpdf(a′, Fixed(x)) ≈ logpdf(m, x) - logpdf(m′, x)
    end
end

@testset "AltBinomial" begin
    n = 10
    p = 0.1
    p′ = 0.25
    a = AltBinomial(n, p)
    a′ = AltBinomial(n, p′)
    b = Binomial(n, p)
    b′ = Binomial(n, p′)
    for _ in 1:100
        x = rand(a)
        @test 0 ≤ x::Int ≤ n
        @test logpdf(a, x) ≈ logpdf(b, x)
        @test logpdf(a, Fixed(x)) - logpdf(a′, Fixed(x)) ≈ logpdf(b, x) - logpdf(b′, x)
    end
end