text:
using PSCOPF

using Test
using Dates
using DataStructures

@testset verbose=true "example_usecase_from_folder" begin

    @testset "example_2bus_mode1" begin
        data_path = joinpath(@__DIR__, "..", "..", "data", "2buses_usecase", "instance")
        out_path = joinpath(@__DIR__, "..", "..", "default_out", "example_2bus_mode1")
        rm(out_path, recursive=true, force=true)
        mode = PSCOPF.PSCOPF_MODE_1

        # load network
        network = PSCOPF.Data.pscopfdata2network(data_path)
        uncertainties = PSCOPF.PSCOPFio.read_uncertainties(data_path)
        gen_init_state = PSCOPF.PSCOPFio.read_initial_state(data_path)
        println("init_state: ", gen_init_state)

        ts1 = Dates.DateTime("2015-01-01T11:00:00")
        TS = PSCOPF.create_target_timepoints(ts1) #11h, 11h15, 11h30, 11h45

        ECH = PSCOPF.generate_ech(network, TS, mode)

        sequence = PSCOPF.generate_sequence(network, TS, ECH, mode)

        exec_context = PSCOPF.PSCOPFContext(network, TS, mode, gen_init_state, uncertainties,
                                            nothing,
                                            out_path)
        PSCOPF.run!(exec_context, sequence)
    end

end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
module TestAutoMLPipeline
using Test

# suppress warnings
@info "suppressing PyCall warnings"
using PyCall
warnings = pyimport("warnings")
warnings.filterwarnings("ignore")

# test modules
include("test_baseline.jl")
include("test_skpreprocessing.jl")
include("test_decisiontree.jl")
include("test_sklearner.jl")
include("test_ensemble.jl")
include("test_crossvalidator.jl")
include("test_skcrossvalidator.jl")
include("test_featureselector.jl")
include("test_pipeline.jl")
include("test_naremover.jl")

#include("test_valdate.jl")
#include("test_mlbase.jl")
#include("test_tsclassifier.jl")
#include("test_statifier.jl")
#include("test_monotonicer.jl")
#include("test_cliwrapper.jl")
#include("test_outliernicer.jl")
#include("test_plotter.jl")
#include("test_normalizer.jl")
#include("test_schemalizer.jl")

end

ext: jl


text:
setprecision(100)
T = BigFloat
PI = T(π)
πi = PI * 1im
πi2 = 2πi

function γ(z, ϵ)
    z = z - 1
    half = T(1) / T(2)
    a = floor(max(real(z), 2)) + half
    while true
        ϵ̂ = √(a) / (2 * PI)^(a + half) / (real(z) + a)
        ϵ̂ < ϵ && break
        a += 1
    end

    coef = 1
    k = 1
    fac = T(1)
    while k < a
        c_k = (a - k)^(k - half) * exp(a - k) * (-1)^(k - 1) / √(2 * PI) / fac
        fac *= k
        coef += c_k / (z + k)
        k += 1
    end
    coef * (z + a)^(z + half) * √(2 * PI) / exp(z + a)
end

function brentq(f, xa, xb, xtol, rtol, iter)
    xpre, xcur, xblk, fpre, fcur, fblk, spre, scur = xa, xb, 0, f(xa), f(xb), 0, 0, 0

    fpre * fcur > 0 && return 0
    fpre == 0 && return xpre
    fcur == 0 && return xcur

    for _ in 1:iter
        if fpre * fcur < 0
            xblk = xpre
            fblk = fpre
            spre = xcur - xpre
            scur = spre
        end

        if abs(fblk) < abs(fcur)
            xpre = xcur
            xcur = xblk
            xblk = xpre

            fpre = fcur
            fcur = fblk
            fblk = fpre
        end

        Δ = (xtol + rtol * abs(xcur)) / 2
        sbis = (xblk - xcur) / 2
        (fcur == 0 || abs(sbis) < Δ) && return xcur

        if abs(spre) < Δ && abs(fcur) < abs(fpre)
            if xpre == xblk
                stry = -fcur * (xcur - xpre) / (fcur - fpre)
            else
                dpre = (fpre - fcur) / (xpre - xcur);
                dblk = (fblk - fcur) / (xblk - xcur);
                stry = -fcur * (fblk * dblk - fpre * dpre) / (dblk * dpre * (fblk - fpre))
            end

            if 2 * abs(stry) < min(abs(spre), 3 * abs(sbis) - Δ)
                spre = scur
                scur = stry
            else
                spre, scur = sbis, sbis
            end
        else
            spre, scur = sbis, sbis
        end

        xpre = xcur
        fpre = fcur

        if abs(scur) > Δ
            xcur += scur
        elseif sbis > 0
            xcur += Δ
        else
            xcur -= Δ
        end

        fcur = f(xcur)
    end

    xcur
end

function I0(s, ϵ)
#     println("[I0] s = $s")
    g(z, s) = z^(-s) * exp(z * z * πi)
    H(w) = 1 / (1 - exp(πi2 * w))

    function f(z, s)
        a = exp(z * πi)
        g(z, s) / (a - 1 / a)
    end

    σ = real(s)
    z_s = √(s / πi2)
    n = max(floor(real(z_s) - imag(z_s)), 1)
#     println("z_s = $z_s")

    Δ = max(0, log(n^(-σ) / ϵ) / (2 * PI))
#     println("σ = $σ, n = $n, n^-σ = $(n^-σ)")
    direction = (1 + 1im) / √ T(2)

    criterion(α) = abs(f(n + 0.5 + direction * α * Δ, s)) - ϵ
    α_1 = brentq(criterion, 0,  2, ϵ, ϵ, 20)
    α_2 = brentq(criterion, 0, -2, ϵ, ϵ, 20)

    z_1 = n + 1/2 + direction * α_1 * Δ
    z_2 = n + 1/2 + direction * α_2 * Δ
#     println("α_1 = $α_1, α_2 = $α_1, $(Δ), $σ, $(n^(-σ))")

    z_l = 0
    z_r = 0
    h = 0

    m = 1
    found = false
    for i in 1:1000
        h = (z_2 - z_1) / m
        inv_2h = 1 / (2 * h)
        z_l = sqrt(inv_2h^2 + s / πi2) + inv_2h
        z_r = sqrt(inv_2h^2 + s / πi2) - inv_2h
        err_l = g(z_l, s) * exp((z_1 - z_l) / h * πi2)
        err_r = g(z_r, s) * exp((z_r - z_1) / h * πi2)
        if abs(err_l) ≤ ϵ && abs(err_r) ≤ ϵ
            found = true
            break
        end
        m = ceil(m * 1.1)
#         println("err_l = $err_l, err_r = $err_r, g(z_l, s) = $(g(z_l, s)), z_l = $z_l, z_r = $z_r")
    end
    @assert found

    n_l = ceil(real(z_l) - imag(z_l))
    n_r = floor(real(z_r) - imag(z_r))

    fs2(i) = i^(-s) * H((i - z_1) / h)
    fs3(i) = i^(-s) * H((z_1 - i) / h)
#     println("[ζ] n = $n, m = $m, n_l = $n_l, n_r = $n_r")

    s0 = sum(@. (1:n)^(-s))
    s1 = sum(@. f((0:m) * h + z_1, s))
    s2 = sum(@. fs2(n_l:n))
    s3 = sum(@. fs3((n+1):n_r))

    s0 + s1 * h - s2 + s3
end

χ(s, ϵ) = (2 * PI)^s / cos(PI * s / 2) / γ(s, ϵ) / 2
# ζ(s, ϵ) = I0(s, ϵ / 2) + χ(s, ϵ) * conj(I0(1 - conj(s), ϵ / 2))
function ζ(s, ϵ)
    ret = I0(s, ϵ / 2) + χ(s, ϵ) * conj(I0(1 - conj(s), ϵ / 2))
    println("[ζ] s = $s, ϵ = $ϵ, ret = $ret, $(χ(s, ϵ))")
    ret
end

function erfc(z, ϵ)
    if abs(z) > 2
        h = π / √(log(6 / ϵ))
        K = ceil(√(log(1 / ϵ)) / h)
        ret = 1 / z^2
        PI = T(π)
        for k in 1:K
            ret += 2 * exp(-h^2 * k^2) / (z^2 + h^2 * k^2)
        end
        return ret * exp(-z^2) * h * z / PI + 2 / (1 - exp(2 * PI * z / h))
    else
        if z ≥ 0
            s = 1
        else
            s = -1
            z = -z
        end

        ϵ0 = ϵ / 2
        K = ceil(log2(1 / ϵ) / 2)
        ϵ1 = ϵ0 / 2K
        ϵ1 = 2^floor(log2(ϵ1))

        t = z
        k = T(0)
        S = T(0)
        while abs(t) / (2k + 1) > ϵ0
            S = S + t / (2k + 1)
            t = -z^2 * t / (k + 1) # error in paper
            k += 1
        end
        S += t / (2k + 1)
        1 - s * S * 2 / √π
    end
end

using SpecialFunctions

function find_T(x, λ, σ, ϵ)
#     println("[find T] x = $x, λ = $λ, σ = $σ, ϵ = $ϵ")
    # only to determine bounds, so it's fine to use Float64 here.

    # σ is real and 1 ≤ σ ≤ 2

    E_Σ(T, x, λ, σ) = exp(λ^2 * σ^2 / 2) / 2π * x^σ * log(zeta(σ)) * expint(λ^2 * T^2 / 2)

    # now we solve min T s.t. E_Σ(T, x, λ, σ) < 0.75ϵ
    # T = √(2u) / λ where E_1(u) < 0.75 ϵ / exp(λ^2 * σ^2 / 2) / (2 * π) * log(gamma(σ))
    limit = 0.75ϵ / (exp(λ^2 * σ^2 / 2) / (2 * π) * log(zeta(σ))) / x^σ
#     println("[find T] limit = $limit, expint(0) = $(expint(0))")
    f(u) = expint(u) - limit

    u = brentq(f, λ^2 / 2, -log(limit), ϵ, ϵ, 100)
#     println("[debug] est u = $(-log(limit)), u = $u")
    T = √(2u) / λ
    println("[param] T = $T, residual = $(E_Σ(T, x, λ, σ)) < ϵ = $ϵ")

    T
end

function quad_π_star(x, σ, λ, ϵ)
    ϵ = ϵ / 4
    T = find_T(x, λ, σ, ϵ)

    h_L = 2 * π * (σ - 1) / (log(x / ϵ) + λ^2 / 2 + 1 / x)
    h_R = 2 * π / (log(x / 2) + σ * λ^2 + λ * √(2 * σ * log(x / 2) + σ^2 * λ^2 + 2 * log(3.4 / ϵ)))
    h = min(h_L, h_R)
    println("[param] h = $h, h_L = $h_L, h_R = $h_R")

    ϵ_0 = ϵ / x^σ / log(x)

    Ψ(s, x, λ, ϵ) = exp((λ * s)^2 / 2) * x^s * log(ζ(s, ϵ)) / s

    results = @. real(Ψ(σ + 1im * h * (1:ceil(T/h)+10), x, λ, ϵ_0))
#     println(results)
    s = sum(results)
    h / π * (Ψ(σ, x, λ, ϵ_0) / 2 + s)
end

function get_primes(l, r)
    a = ones(Bool, r - l + 1)
    for i in 2:Int64(floor(√r))
        li = max((l - 1) ÷ i + 1, 2)
        ri = r ÷ i
        for j in li : ri
            a[i * j - l + 1] = 0
        end
    end

    ret = []
    for i in max(l, 2):r
        if a[i - l + 1]
            push!(ret, i)
        end
    end
    ret
end

function Δ(x, λ, ϵ)
    Φ(ρ, ϵ) = erfc(ρ / √2, ϵ) / 2
    ϕ(u, x, λ, ϵ) = Φ(log(u / x) / λ, ϵ)
    E_plus(u, x, λ, ϵ) = x * exp(λ^2 / 2) * Φ(log(u / x) / λ - λ, ϵ / x) - u * Φ(log(u / x) / λ, ϵ / 8 / u)
    E_minus(u, x, λ, ϵ) = u * Φ(-log(u / x) / λ, ϵ / 8 / u) - x * exp(λ^2 / 2) * Φ(λ - log(u / x) / λ, ϵ / x)

    ϵ = ϵ / 2
    x1 = floor(x)
    if E_minus(x1, x, λ, ϵ) > ϵ
        g(u) = E_minus(u, x, λ, ϵ) - 0.75 * ϵ
        u = brentq(g, 2, x, 1e-8, 1e-8, 100)
        x1 = floor(u)
    end

    x2 = ceil(x)
    if E_plus(x2, x, λ, ϵ) > ϵ
        f(u) = E_plus(u, x, λ, ϵ) - 0.75 * ϵ
        u = brentq(f, x * 2, x, 1e-8, 1e-8, 100)
        x2 = ceil(u)
    end
    println("[Δ] x1 = $x1, x2 = $x2")
    println("$(E_minus(900, x, λ, ϵ))")
    println("$(E_plus(x2, x, λ, ϵ))")

    N = (x2 - x1 + 1) + ceil(√x2)
    ϵ_1 = ϵ * 2^ceil(-log2(N))

    ret = 0
    for p in get_primes(Int64(x1), Int64(x2))
        ret -= ϕ(T(p), x, λ, ϵ_1)
        if p ≤ x
            ret += 1
        end
    end
    for p in get_primes(2, Int64(floor(√x2)))
        m = 2
        while p^m < x2
            if p^m < x1
                ret -= T(1 / m)
            else
                ret -= ϕ(p^m, x, λ, ϵ_1) / m
            end
            m += 1
        end
    end

    ret
end

function analytic_π(x)
    μ = 1 + 0
    σ_LR = 2
    σ_0 = (μ + 1/2) / (μ - 1/2)
#     σ = min(σ_LR, σ_0)
    σ = 1.5
    println("[param] σ = $σ")

    λ = 1 / (√x * log(x)^0.25)
    println("[param] λ = $λ")

    ϵ = 1/4
    pi_star = quad_π_star(x, σ, λ, ϵ / 2)
    Delta = Δ(x, λ, ϵ / 2)
    println("pi_star = $pi_star, Delta = $Delta")
    round(Delta + pi_star)
end

print(analytic_π(10^4))

max_stars_repo_name: roosephu/apc

max_issues_repo_name: roosephu/apc

max_forks_repo_name: roosephu/apc


text:
import Rotations.params

params(::Type{<:UnitQuaternion}) = 4
params(::Type{<:RodriguesParam}) = 3
params(::Type{<:MRP}) = 3
params(::Type{<:RotMatrix3}) = 9
params(::Type{<:RotMatrix2}) = 4

"""
    LieState{R,P}

Specifies a state with rotational components mixed in with standard vector components. All
rotational components are assumed to be parameterizations of 3D rotations.

# Parameters
- `R <: Rotation` is the rotational representation used in the state vector. Must have
`params(::Type{R})` defined, which returns the number of parameters used by the rotation,
as well a constructor that takes each parameter as a separate scalar argument.

- `P <: Tuple{Vararg{Int}}` is a tuple of integers specifying the partitioning of the state
vector. Each element of `P` specifies the length of the vector component between the rotational
components, and `P[1]` and `P[end]` specify the number of vector states at the beginning and
end of the state vector.

# Examples
If we want to construct a state vector like the following: `[v3, q, v2, q, v3]` where `v2`
and `v3` and vector components of length 2 and 3, respectively, and `q` is a 4-dimensional
unit quaternion. The `LieState` for this state vector would be
`LieState{UnitQuaternion{Float64},3,2,3}`. The length should be (3+4+2+4+3) = 16, which can
be verified by `length(s::LieState)`.

# Constructors
    LieState(::Type{R}, P::Tuple{Vararg{Int}})
    LieState(::Type{R}, p1::Int, p2::Int, p3::Int...)
"""
struct LieState{R,P}
    function LieState{R,P}() where {R<:Rotation,P}
        new{R,P::Tuple{Vararg{Int}}}()
    end
end

@inline LieState(::Type{R}, P::Tuple{Vararg{Int}}) where R <: Rotation = LieState{R,P}()
@inline LieState(::Type{R}, P::Int...) where R <: Rotation = LieState{R,P}()

"""
    QuatState(n::Int, Q::StaticVector{<:Any,Int})
    QuatState(n::Int, Q::NTuple{<:Any,Int})

Create a `n`-dimensional `LieState` assuming `R = UnitQuaternion{Float64}` and `Q[i]` is the
first index of each quaternion in the state vector.

# Example
If we want to construct a state vector like the following: `[v3, q, v2, q, v3]` where `v2`
and `v3` and vector components of length 2 and 3, respectively, and `q` is a 4-dimensional
unit quaternion. Since the first quaternion starts at index 4, and the second starts at index
10, Q = [4,10]. The entire length of the vector is `n = 16 = 3 + 4 + 2 + 4 + 3`, so we would
call `QuatState(16, SA[4,10])`.
"""
@generated function QuatState(n::Int, Q::StaticVector{NR,Int}) where {NR}
    if NR == 1
        P0 = SVector{0,Int}()
    else
        diffs = [:(Q[$i] - Q[$(i-1)] - 4) for i = 2:NR]
        P0 = :(SVector{NR-1}(tuple($(diffs...))))
    end
    quote
        P = $P0
        P = push(P, n - (Q[end] + 4) + 1)
        P = pushfirst(P, Q[1] - 1)
        return LieState(UnitQuaternion{Float64}, Tuple(P))
    end
end
@inline QuatState(n::Int, Q::NTuple{<:Any,Int}) = QuatState(n, SVector(Q))

"number of rotations"
num_rotations(::LieState{<:Any,P}) where P = length(P) - 1

state_dim_vec(::LieState{<:Any,P}) where P = sum(P)
state_dim_rot(s::LieState{R}) where R = params(R)*num_rotations(s)

Base.length(s::LieState{R,P}) where {R,P} = params(R)*num_rotations(s) + sum(P)
Base.length(::Type{LieState{R,P}}) where {R,P} = params(R)*(length(P)-1) + sum(P)

state_diff_size(s::LieState{R,P}) where {R,P} = 3*num_rotations(s) + sum(P)

state_diff_size(model::LieGroupModel) = state_diff_size(LieState(model))

# Useful functions for meta-programming
rot_inds(R,P, i::Int) = (sum(P[1:i]) + (i-1)*params(R)) .+ (1:params(R))
vec_inds(R,P, i::Int) =
    ((i > 1 ? sum(P[1:i-1]) : 0) + (i-1)*params(R)) .+ (1:P[i])
inds(R,P, i::Int) = isodd(i) ? vec_inds(R,P, 1+i÷2) : rot_inds(R,P, i÷2)
@inline inds(s::LieState{R,P}, i::Int) where {R,P} = inds(R, P, i)
rot_state(R,P, i::Int, sym=:x) = [:($(sym)[$j]) for j in rot_inds(R,P,i)]

"""
    vec_states(model::LieGroupModel, x)
    vec_states(s::LieState, x)

Extracts the "vector" states out of the state vector `x` for a `LieGroupModel`. Returns
a tuple `v` of `SVector`s, where `length(v[i])` is equal to the length specified by the
`LieState`.
"""
@generated function vec_states(s::LieState{R,P}, x) where {R,P}
    T = eltype(x)
    inds = [vec_inds(R,P,i) for i = 1:length(P)]
    states = [[:(x[$i]) for i in ind] for ind in inds]
    vecs = [:(SVector{$(length(inds)),$T}($(inds...))) for inds in states]
    quote
        tuple($(vecs...))
    end
end
@inline vec_states(model::LieGroupModel, x) = vec_states(LieState(model), x)

"""
    vec_states(model::LieGroupModel, x)
    vec_states(s::LieState, x)

Extracts the rotations out of the state vector `x` for a `LieGroupModel`. Returns
a tuple rotations, whose type matches the rotation type specified in the `LieState`.
"""
@generated function rot_states(s::LieState{R,P}, x) where {R,P}
    T = eltype(x)
    inds = [rot_inds(R,P,i) for i = 1:length(P)-1]
    states = [[:(x[$i]) for i in ind] for ind in inds]
    if R <: UnitQuaternion
        vecs = [:(R($(inds...),false)) for inds in states]
    else
        vecs = [:(R($(inds...))) for inds in states]
    end
    quote
        tuple($(vecs...))
    end
end
@inline rot_states(model::LieGroupModel, x) = rot_states(LieState(model), x)

@generated function Base.rand(s::LieState{R,P}) where {R,P}
    nr = length(P) - 1   # number of rotations
    np = nr + length(P)  # number of partitions
    n = length(s)

    q = [:($(Symbol("q$i")) = Rotations.params(rand(R))) for i = 1:nr]

    x = Expr[]
    for i = 1:np
        if isodd(i)  # vector part
            vi = inds(R,P,i)
            for j in vi
                push!(x, :(rand()))
            end
        else
            r = i ÷ 2
            ri = inds(R,P,i)
            for j = 1:Rotations.params(R)
                push!(x, :($(Symbol("q$r"))[$j]))
            end
        end
    end
    quote
        $(Expr(:block, q...))
        $(:(SVector{$n}(tuple($(x...)))))
    end
end

function Base.rand(model::LieGroupModel)
    rand(LieState(model)), @SVector rand(control_dim(model))
end

@inline state_dim(model::LieGroupModel) = length(LieState(model))

@inline state_diff(model::LieGroupModel, x::AbstractVector, x0::AbstractVector, errmap=Rotations.CayleyMap()) =
    state_diff(LieState(model), x, x0, errmap)
@generated function state_diff(s::LieState{R,P}, x::AbstractVector, x0::AbstractVector,
        errmap=Rotations.CayleyMap()) where {R,P}
    nr = length(P) - 1   # number of rotations
    np = nr + length(P)  # number of partitions
    n̄ = 3*nr + sum(P)    # error state size

    # Generate a vector of δq = q0\q expressions for each q in the state
    dq = [:($(Symbol("q$i")) = Rotations.rotation_error(R($(rot_state(R,P,i)...)), R($(rot_state(R,P,i,:x0)...)), errmap)) for i = 1:nr]

    # Generate the vector of expressions for each element of the state differential
    dx = Expr[]
    for i = 1:np
        # for the vector parts, simply subtract the indices from the original states
        if isodd(i)
            vi = inds(R,P,i)
            for j in vi
                push!(dx, :(x[$j] - x0[$j]))
            end
        # for the rotational parts, use the elements of the rotational errors generated above
        else
            r = i ÷ 2  # rotation index
            for j = 1:3
                push!(dx, :($(Symbol("q$r"))[$j]))
            end
        end
    end
    quote
        $(Expr(:block, dq...))
        $(:(SVector{$n̄}(tuple($(dx...)))))
    end
end

@inline state_diff_jacobian!(G, model::LieGroupModel, x::StaticVector) =
    state_diff_jacobian!(G, LieState(model), x)
@generated function state_diff_jacobian!(G, s::LieState{R,P}, x::StaticVector) where {R,P}
    nr = length(P) - 1   # number of rotations
    np = nr + length(P)  # number of partitions
    nv = length(P)
    n̄ = 3*nr + sum(P)    # error state size
    n = params(R)*nr + sum(P)

    # Generate a vector of δq = q0\q expressions for each q in the state
    q = [:(R($(rot_state(R,P,i)...))) for i = 1:nr]

    # Generate a vector of expressions assigning 1s to all the vector state diagonals
    Gv = Expr[]
    r = 1
    c = 1
    for k = 1:nv
        for j = 1:P[k]
            push!(Gv, :(G[$(LinearIndices((n,n̄))[r,c])] = 1))
            r += 1
            c += 1
        end
        r += params(R)
        c += 3
    end

    # Generate a vector of expressions assigning the differential Jacobians for the rotations
    Gr = Expr[]
    for k = 1:nr
        rinds = rot_inds(R, P, k)
        cinds = (sum(P[1:k]) + 3*(k-1)) .+ (1:3)
        push!(Gr, :(G[$rinds,$cinds] .= Rotations.∇differential($(q[k]))))
    end
    quote
        $(Expr(:block, Gv...))
        $(Expr(:block, Gr...))
        return nothing
    end
end

@inline ∇²differential!(∇G, model::LieGroupModel, x::StaticVector, dx::AbstractVector) =
    ∇²differential!(∇G, LieState(model), x, dx)
@generated function ∇²differential!(∇G, s::LieState{R,P}, x::StaticVector, dx::AbstractVector) where {R,P}
    nr = length(P) - 1   # number of rotations
    np = nr + length(P)  # number of partitions
    nv = length(P)
    n̄ = 3*nr + sum(P)    # error state size
    n = params(R)*nr + sum(P)

    # Generate a vector of δq = q0\q expressions for each q in the state
    q  = [:(R($(rot_state(R,P,i,:x)...)))  for i = 1:nr]
    dq = [:(SVector{$(params(R))}($(rot_state(R,P,i,:dx)...))) for i = 1:nr]

    Gr = Expr[]
    for k = 1:nr
        cinds = (sum(P[1:k]) + 3*(k-1)) .+ (1:3)
        push!(Gr, :(view(∇G,$cinds,$cinds) .= SMatrix{3,3}(Rotations.∇²differential($(q[k]), $(dq[k])))))
    end
    quote
        $(Expr(:block, Gr...))
        return nothing
    end
end


text:
#=
  From the WebPPL model:https://github.com/probmods/ppaml2016/blob/gh-pages/chapters/5-data.md
  """
  Here, we explore the result of an experiment with 15 trials and binary outcomes (e.g., flipping a coin with an
  uncertain weight, asking people if they'll vote for Hillary Clinton or Donald Trump, ...)
  """
=#

using Turing
include("jl_utils.jl")

@model function bda(k=1,n=15)
    p ~ Uniform(0, 1)

    # Observed k number of successes, assuming a binomial
    k ~ Binomial(n, p)

    # sample from binomial with updated p
    posteriorPredictive ~ Binomial(n,p)

    # sample fresh p
    prior_p ~ Uniform(0,1)
    # sample from binomial with fresh p
    priorPredictive ~ Binomial(n,prior_p)

end

# data
k = 1  # number of successes
n = 15 # number of attempts

model = bda(k,n)

# chns = sample(model, MH(), 10_000)
# chns = sample(model, PG(15), 10_000)
# chns = sample(model, IS(), 10_000)
chns = sample(model, SMC(), 10_000)
# chns = sample(model, SMC(), MCMCThreads(), 10_000, 4)

display(chns)

show_var_dist_pct(chns, :priorPredictive)
show_var_dist_pct(chns, :posteriorPredictive)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
using LightGraphs: degree, vertices, rem_vertex!, AbstractGraph
using MetaGraphs
using DelimitedFiles
using DataFrames

"""
source: https://discourse.julialang.org/t/get-the-argument-names-of-an-function/32902/3
I modified this function to retrieve the keyword arguments.
"""
function getkwargs(f)
    ms = collect(methods(f))
    m = last(ms)
    argnames = ccall(:jl_uncompress_argnames, Vector{Symbol}, (Any,), m.slot_syms)
    isempty(argnames) && return argnames
    return argnames[m.nargs + 1:end]
end

function get_valid_temp_file(instance_name)
    index =  0
    file =  "$instance_name.temp"
    while ispath(file)
        file = "$(instance_name)_$(index).temp"
        index += 1
    end
    return file
end

function load_matpower(filename)
    instance_name = split(basename(filename), '.')[1]
    temp_filename = get_valid_temp_file(instance_name)
    touch(temp_filename)
    f = open(filename)
    out = open(temp_filename, "w")
    # removing all ';' at end of lines
    while !eof(f)
        line = readline(f)
        if length(line) > 0 && line[1] != '%' && line[1] != 'f'
            s = split(line, ";")
            println(out, s[1])
        end
    end
    close(f)
    close(out)
    data = DelimitedFiles.readdlm(temp_filename)
    rm(temp_filename)
    return data
end

function find_numarray(i_start, data)
    i_debut = i_start
    while !isa(data[i_debut, 1], Int)
        i_debut += 1
    end
    i_fin=i_debut
    while !isa(data[i_fin,1], SubString)
        i_fin += 1
    end
    i_debut, i_fin - 1
end

function checkfor(data, line_ind, name)
    (data[line_ind, 1] == name) || error("Expected $name at line $line_ind got $(data[line_ind,1]) instead.")
end

function get_names_buses(data, nb_bus)
    index = 1
    while data[index, 1] != "mpc.bus_name"
        index += 1
    end
    index += 1
    names = [join(data[i, :]) for i in index:index + nb_bus - 1]
    return names
end

function convertkeytosymbol(dict)
    symboldict = Dict()
    for (key, val) in dict
        if typeof(val) <: Dict
            symboldict[Symbol(key)] = Dict(Symbol(subkey) => subval for (subkey, subval) in dict[key])
        else
            symboldict[Symbol(key)] = val
        end
    end
    return symboldict
end

function parsecolumns!(type::Type, df::DataFrame, exclude::Vector{String}=["solver.nb_iter", "solver.solving_time"])
    for (index, name) in enumerate(names(df))
        name in exclude && continue
        try
            values = tryparse.(type, df[!,name])
            if !(typeof(values) in [Array{Union{Nothing,type},1}, Array{Nothing,1}])
                df[!,name] = values
            end
        catch e
            continue
        end
    end
end

function split_index(length, nb_split)
    if nb_split > length
        return [i:i for i in 1:length]
    end
    quotient = div(length, nb_split)
    remainder = length % nb_split
    indexes = []
    extra = remainder == 0 ? 0 : 1
    offset = 0
    for i in 1:nb_split
        start = (i - 1) * (quotient + extra) + 1 + offset
        stop = i * (quotient + extra) + offset
        if i == remainder
            offset = remainder
            extra = 0
        end
        push!(indexes, [start, stop])
    end
    return indexes
end

function format_options(options)
    return Dict(string(key) => Dict(string(key2) => val2 for (key2, val2) in val)
           for (key, val) in options)
end

function get_param_merge(path)
    open(path, "r") do io
        params = split(strip(read(io, String)), ',')
        return parse.(Float64, params)
    end
end

function get_nb_constraint(path)
    return length(readlines(path)) - 1
end

max_stars_repo_name: charlyalizadeh/OPFdecgen

max_issues_repo_name: charlyalizadeh/OPFdecgen

max_forks_repo_name: charlyalizadeh/OPFdecgen


text:
struct SIGPLAN end
textwidth(::Type{SIGPLAN}) = 7.00137inch
linewidth(::Type{SIGPLAN}) = 3.33461inch
defaultsm(::Type{SIGPLAN}) = 8

scale(::Type{Plots.PyPlotBackend}) = 4.629 * 0.6485500820772881
scale(::Type{Plots.PGFPlotsBackend}) = 4.629  * 0.6485500820772881

dpi(::Type{Plots.PGFPlotsBackend}) = 300
dpi(::Type{Plots.PyPlotBackend}) = 300


text:
using FPRS
using Test

@testset "FPRS.jl" begin
    include("test-utilities.jl")
    include("test-densities.jl")
    include("test-projections.jl")
    include("test-loadElser.jl")
end

max_stars_repo_name: tsipenyuk/FPRS

max_issues_repo_name: tsipenyuk/FPRS

max_forks_repo_name: tsipenyuk/FPRS


text:
module AA222Project2_2019

using JSON
using PyCall
using Random
using Statistics
using ForwardDiff

const project2 = PyNULL()
const err = Ref{String}("")

function append_to_python_search_path(str::AbstractString)
    pushfirst!(PyVector(pyimport("sys")."path"), str)
end

const COUNTERS = Dict{String, Int}()
const PROBS = Dict{String, Dict{Int, Any}}()

"""
Count the number of times a function was evaluated
"""
macro counted(f)
    name = f.args[1].args[1]
    name_str = String(name)
    body = f.args[2]
    update_counter = quote
        if !haskey(COUNTERS, $name_str)
            COUNTERS[$name_str] = 0
        end
        COUNTERS[$name_str] += 1
    end
    insert!(body.args, 1, update_counter)
    return f
end

function check_constraints(c, x::Vector)
    vals = c(x)
    COUNTERS[String(typeof(c).name.mt.name)] -= 1
    return all(vals .<= 0)
end

"""
    get_score(f, g, x_star_hat, n)

Arguments
    - `f`: function
    - `g`: gradient of function
    - `c`: constraint function
    - `x_star_hat`: Found optima
    - `n`: Number of evaluations allowed (f+2g)
Returns
    - `num_evals`
    - `score`
"""
function get_score(f, g, c, x_star_hat, n)
    f_name = String(typeof(f).name.mt.name)
    g_name = String(typeof(g).name.mt.name)
    c_name = String(typeof(c).name.mt.name)

    if haskey(COUNTERS, f_name)
        f_evals = COUNTERS[f_name]
    else
        f_evals = 0
    end

    if haskey(COUNTERS, g_name)
        g_evals = COUNTERS[g_name]
    else
        g_evals = 0
    end

    if haskey(COUNTERS, c_name)
        c_evals = COUNTERS[c_name]
    else
        c_evals = 0
    end

    score = typemin(Int32)
    num_evals = (f_evals + 2*g_evals)
    if check_constraints(c, x_star_hat)
        if num_evals <= n && c_evals <= n
            score = -f(x_star_hat)
            COUNTERS[f_name] -= 1
        end
    else
        err[] = "Constraints not satisfied!"
    end
    return num_evals, c_evals, score
end

if isfile(joinpath(@__DIR__, "simple.jl"))
    include(joinpath(@__DIR__, "simple.jl"))
    PROBS["simple"] = simple_problems
end

"""
    main(mode, pidx, repeat)

Arguments:
    - `mode`: simple or secret
    - `pidx`: problem index
    - `repeat`: Number of Monte Carlo evaluations
Returns:
    - nothing
"""
function main(mode::String, pidx::Int, repeat::Int)
    if mode ∉ ["simple", "secret"]
        err[] = "Invalid mode = $(mode)"
    end
    inseed = 42
    if mode == "simple"
        inseed = 1337
    end

    scores = Float64[]
    errors = String[]
    nevals = Int[]
    cevals = Int[]
    score = 0
    neval = 0
    ceval = 0
    if err[] == ""
        # Repeat the optimization with a different initialization
        for i in 1:repeat
            Random.seed!(inseed + i)
            prob = PROBS[mode][pidx]
            try
                if project2 != PyNULL()
                    opt_func = project2.optimize
                else
                    opt_func = optimize
                end
                if mode == "secret"
                    # Reset the affine_function
                    sz = size(prob.x0(), 1)
                    A[] = randn(sz, sz)
                    b[] = abs.(randn(sz))
                end

                x_star_hat = opt_func(prob.f, prob.g, prob.c, prob.x0(), prob.n, "$(mode)_$(pidx)")
                neval, ceval, score = get_score(prob.f, prob.g, prob.c, x_star_hat, prob.n)
                push!(errors, err[])
            catch e
                push!(errors, "$(e)")
                x_star_hat = typemin(Int32)
                score = typemin(Int32)
                neval = typemin(Int32)
                ceval = typemin(Int32)
            end

            push!(nevals, neval)
            push!(cevals, ceval)
            push!(scores, score)

            for k in keys(COUNTERS)
                delete!(COUNTERS, k)
            end
        end

        # So that we know if there was a problem with a particular seed
        enerrs = enumerate(errors)
        # Filter out empty
        enerrs = filter(x -> x[2] != "", collect(enerrs))
        err[] = join(enerrs, " \n")
    else
        push!(nevals, typemin(Int32))
        push!(cevals, typemin(Int32))
    end
    result = Dict("score" => mean(scores), "error"=> err[])

    open(".results_$(mode)_$(pidx).json", "w") do io
        write(io, JSON.json(result))
    end

    final_score = result["score"]
    max_evals = maximum(nevals)
    max_cevals = maximum(cevals)
    println("===================")
    println("Mode: $(mode)")
    println("Problem: $(pidx)")
    println("-------------------")
    if err[] != ""
        println("Error: $(err[])")
    end
    println("Max Num Function (and Gradient) Evals: $(max_evals)")
    println("Max Num Constraint Function Evals: $(max_cevals)")
    println("Avg Score: $(final_score)")
    println("===================")
end

function __init__()
    if length(ARGS) < 1
        filepath = "project2.jl"
    else
        filepath = ARGS[1]
    end

    if isfile(filepath)
        if filepath[end-1:end] == "jl"
            @info("Including Julia file")
            try
                include(filepath)
            catch e
                err[] = "$(e)"
            end
        elseif filepath[end-1:end] == "py"
            @info("Importing Python file")
            try
                dir = dirname(filepath)
                append_to_python_search_path(dir)
                copy!(project2, pyimport("project2"))
            catch e
                err[] = "$(e)"
            end
        else
            err[] = "Not a correct file type"
        end
    else
        err[] = "File: $(filepath) doesn't exist"
    end
end

end # module


text:
using BenchmarkTools, DynamicalSystemsBase

const SUITE = BenchmarkGroup(["DynamicalSystemsBase"])

include("reconstructions.jl")
include("integrators.jl")


text:
import Base: get, getindex, setindex!, haskey
import Pkg: TOML

mutable struct ConfigManager{DM <: DataManager}
    config_dict::Dict
    parsed_config::Vector
    total_combinations::Int
    spec::Dict
    rootdir::String
    dataManager::DM
end

function ConfigManager(config_file::String, rootdir::String, manager::DM = HDF5Manager()) where { DM<:DataManager }
    cfg = Dict()
    if splitext(config_file)[end] == ".toml"
        cfg = TOML.parsefile(config_file)
    else
        open(config_file, "r") do f
            txt = read(f,String)
            cfg = JSON.parse(txt)
        end
    end

    parsed = flattenToArray(cfg)
    total = getTotalCombinations(cfg)

    cm = ConfigManager(cfg, parsed, total, Dict(), rootdir, manager)

    # sometimes want this without having to parse!() the config
    cm.spec["save_path"] = getPermutation(cm.parsed_config, 1)["save_path"]

    _ensure([_data_root(cm), _output_root(cm), _checklists(cm) ])

    return cm
end

function Base.copy(cm::ConfigManager{DM}) where {DM<:DataManager}
    ConfigManager{DM}(
        deepcopy(cm.config_dict),
        copy(cm.parsed_config),
        cm.total_combinations,
        Dict("save_path"=>getPermutation(cm.parsed_config, 1)["save_path"]),
        cm.rootdir,
        typeof(cm.dataManager)()
    )
end

haskey(self::ConfigManager, k::String) = k ∈ keys(self.spec)
getindex(self::ConfigManager, param::String) = self.spec[param]
setindex!(self::ConfigManager, value::Any, param::String) = setindex!(self.spec, value, param)

_data_root(self::ConfigManager) = joinpath(self.rootdir, "data")
_output_root(self::ConfigManager) = joinpath(_data_root(self), "output")
_checklists(self::ConfigManager) = joinpath(_data_root(self), "checklists")

total_combinations(self::ConfigManager) = self.total_combinations

struct _CM_Iterator{CM<:ConfigManager}
    manager::CM
    num_runs::Int
end

Base.length(iter::_CM_Iterator) = total_combinations(iter.manager)*iter.num_runs

iterator(cm::ConfigManager, num_runs) = _CM_Iterator(cm, num_runs)

function Base.iterate(iter::_CM_Iterator, state=(1,1))

    if state[1] > total_combinations(iter.manager)
        return nothing
    end

    new_state = begin
        if state[2] == iter.num_runs
            (state[1] + 1, 1)
        else
            (state[1], state[2]+1)
        end
    end
    cm = copy(iter.manager)

    parse!(cm, state[1], state[2])
    return ((state[1]-1)*iter.num_runs + state[2], cm), new_state

end

"""
    parse!

    get spec while modifying config manager. repercusions.
"""
function parse!(self::ConfigManager, idx::Int, run::Int=1)
    self.spec = getPermutation(self.parsed_config, idx)

    self.spec["run"] = run
    self.spec["param_setting"] = idx

    _ensure([get_expdir(self), get_logdir(self)])

    return self
end

function param_values(self::ConfigManager, key::String)
    path = map(s->string(s), split(key,"/"))
    d = copy(self.config_dict)
    for ss in path
        d = d[ss]
    end
    return d
end

function params_with(self::ConfigManager, cfg_dict::Dict)
    indices = Int[]
    for idx in 1:self.total_combinations
        param_setting = getPermutation(self.parsed_config,idx)
        diff = getDiff(param_setting, cfg_dict)
        if isempty(diff)
            push!(indices, idx)
        end
    end
    return indices
end

function get_expdir(self::ConfigManager)
    return joinpath(_output_root(self), self["save_path"])
end

function get_logdir(self::ConfigManager)
    return joinpath(get_expdir(self), string(self["run"],"_run"), string(self["param_setting"], "_param_setting"))
end

function get_datafile(self::ConfigManager)
    path = joinpath(get_logdir(self),string("data", extension(self.dataManager)))
    if !isfile(path)
        touch(path)
    end
    return path
end

function log_config(self::ConfigManager, logger::Logging.SimpleLogger)
    function _log(keys, attrs)
        for param in keys
            value = attrs[param]
            if value == nothing
                continue
            end

            s = string(param, ": ", value)

            # log to stdout
            @info s

            # log to logfile
            Logging.with_logger(logger) do
                @info s
            end
        end
    end

    attrs = get_attrs(self)
    attr_lst = sort(collect(keys(attrs)), by = x->x[1])
    model_attrs = filter(k -> occursin("|", k), attr_lst)
    other_attrs = filter(k-> !occursin("|",k), attr_lst)

    _log(other_attrs, attrs)
    _log(model_attrs, attrs)
end

function get_attrs(self::ConfigManager)
    return copy(self.spec)
end

function dump_config(self::ConfigManager)
    # Dump the config to the data file so it's easier
    # to recreate a single specific experiment/setting

    path = joinpath(get_logdir(self), "config.json")
    attrs = get_attrs(self)
    delete!(attrs, "param_setting")
    delete!(attrs, "run")

    open(path,"w+") do f
        JSON.print(f, attrs, 4)
    end
end

function get_subconfig(self::ConfigManager, keys...)
    sub = self.spec[keys[1]]
    for k in keys[2:end]
        sub = sub[k]
    end
    return sub
end


text:
function differentiate(obj::UnivariateFunctionDifferentiator, arg0::UnivariateFunction)
    return jcall(obj, "differentiate", UnivariateDifferentiableFunction, (UnivariateFunction,), arg0)
end


text:
module ExtendedVariantModeSetJumpCountTrigger

using ..Ahorn, Maple

@mapdef Trigger "ExtendedVariantMode/SetJumpCountTrigger" SetJumpCountTrigger(x::Integer, y::Integer, width::Integer=16, height::Integer=16,
	jumpCount::Int=0, mode::String="Set")

const placements = Ahorn.PlacementDict(
	"Set Jump Count Trigger (Extended Variant Mode)" => Ahorn.EntityPlacement(
		SetJumpCountTrigger,
		"rectangle"
	)
)

Ahorn.editingOptions(trigger::SetJumpCountTrigger) = Dict{String, Any}(
	"mode" => String["Set", "Cap"]
)

end


text:
module Morse

using KnetLayers, ArgParse, Random
import KnetLayers: _pack_sequence

dir(path...) = joinpath(dirname(@__DIR__),path...)

include("parser.jl"); export Analysis, EncodedAnalysis, ParsedIO, EncodedIO,
                             IndexedDict, Vocabulary, Parser, parseDataLine,
                             parseFile, encode, specialTokens,
                             UDDataSet, TRDataSet, TRtoLower;

include("minibatch.jl"); export SentenceBatch, miniBatch;

include("code2lang.jl"); export CODE_TO_LANG;

include("models.jl"); export Sequential, Discriminative, Disambiguator,
                             MorseModel, S2S, S2SContext, Classifier, S2SContexDis,
                             WordEncoder,ContextEncoder, OutputEncoder, Decoder,
                             loss, predict, MorseDis, S2SContexDis;

include("util.jl"); export download, loadModel, saveModel, setoptim!, lrdecay!,
                           splitdata, createsplits, prepareData, StringAnalysis,
                           evaluate, F1average, F1update!, getLabels,
                           makeFormat, printFormat,trained;

include("train.jl"); export trainmodel!, trainepoch!;

end # module


text:
# ------------------------------------------------------------------------------------------------ #
#                                     Time series distributions                                    #
# ------------------------------------------------------------------------------------------------ #
@shorthands tsdensity
@recipe function f(::Type{Val{:tsdensity}}, plt::AbstractPlot; colordensity=false)
    x, y = plotattributes[:x], plotattributes[:y]
    i = isfinite.(x) .& isfinite.(y)
    x, y = x[i], y[i]
    xlims, ylims = extrema(x), widen(extrema(y), 0.1)
    ylims = (ylims[1] - 0.1*(ylims[1] + ylims[2])/2, ylims[2] + 0.1*(ylims[1] + ylims[2])/2)
    sz = (800, 200)
    size --> sz
    legend --> false
    seriescolor --> :black
    left_margin --> 5mm
    bottom_margin --> 5mm
    right_margin --> 20mm
    top_margin --> 5mm
    link := :both
    grid --> false
    #layout --> @layout [
    #    ts{0.9w,1.0h} density{0.1w,1.0h}
    #]
    @series begin
        if colordensity
            seriestype := :tscolordensity
        else
            seriestype := :path
        end
        framestyle --> :on
        xguide --> "t"
        xlims --> xlims
        ylims --> ylims
        x := x
        y := y
    end

    idxs = LinRange(ylims..., 1000)
    p = pdf(kde(y), idxs)
    @series begin
        seriestype := :shape
        ticks := nothing
        inset_subplots := (1, bbox(1.0, 0.0, 0.1, 1.0))
        subplot := 2
        background_color_inside := nothing
        foreground_color_border := nothing
        fillcolor --> plotattributes[:seriescolor]
        linecolor --> plotattributes[:seriescolor]
        framestyle --> :none
        yguide := ""
        xguide := ""
        xticks := nothing
        yticks := nothing
        ylims --> ylims
        xlims --> (0, max(p...))
        fillalpha := 0.25
        x := [0.0, p..., 0.0]
        y := [ylims[1], idxs..., ylims[2]]
    end
end

tsdensity(P::Process; vars=1, kwargs...) = tsdensity(timeseries(P, vars); kwargs...)
export tsdensity

@shorthands tscolordensity
@recipe function f(::Type{Val{:tscolordensity}}, plt::AbstractPlot)
    x, y = plotattributes[:x], plotattributes[:y]
    i = isfinite.(x) .& isfinite.(y)
    x, y = x[i], y[i]
    xlims, ylims = (extrema(x), extrema(y))
    ylims = (ylims[1] - 0.1*(ylims[1] + ylims[2])/2, ylims[2] + 0.1*(ylims[1] + ylims[2])/2)
    size := (800, 200)
    legend --> false
    seriescolor --> :black
    link := :both
    bottom_margin --> 5mm
    grid --> false

    #idxs = LinRange(ylims..., 1000)
    p = pdf(kde(y), y)

    @series begin
        seriestype := :path
        right_margin --> 0mm
        top_margin --> 0mm
        colorbar --> false
        framestyle --> :on
        xguide --> "t"
        subplot := 1
        xlims --> xlims
        ylims --> ylims
        line_z --> p
        x := x
        y := y
    end

end
tscolordensity(P::Process; vars=1, kwargs...) = tscolordensity(timeseries(P, vars); kwargs...)
export tscolordensity

# ------------------------------------------------------------------------------------------------ #
#                                Plot 2D trajectory with time series                               #
# ------------------------------------------------------------------------------------------------ #
@shorthands marginaltrajectory2
@recipe function f(::Type{Val{:marginaltrajectory2}}, plt::AbstractPlot; marginaldownsample=1, marginallength=2000, marginalcolor=:black, colormode=:velocity, mainalpha=0.5, trajectorycolor=:black)
    cmode = colormode
    x, y = plotattributes[:x], plotattributes[:y]
    i = isfinite.(x) .& isfinite.(y)
    x, y = x[i], y[i]
    xlims, ylims = (extrema(x), extrema(y))
    top_margin --> 20Plots.mm
    right_margin --> 20Plots.mm
    legend --> :none
    framestyle --> :box

    if cmode == :velocity
        velocity = sqrt.(sum([(r[2:end] .- collect(r[1:end-1])).^2 for r ∈ [x, y]], dims=1)[1])
    elseif !isnothing(cmode) && cmode != :none
        @error "Not a supported cmode"
    end

    inset_subplots := [(1, bbox(0.0, -0.1, 1.0, 0.1)), (1, bbox(1.0, 0.0, 0.1, 1.0))]
    @series begin
        seriestype := :path
        seriescolor := marginalcolor
        ticks := nothing
        linewidth --> 2
        subplot := 2
        framestyle := :none
        x := (1:marginaldownsample:length(y))[1:marginallength]
        y := (y[1:marginaldownsample:length(y)])[1:marginallength]
    end

    @series begin
        seriestype := :path
        seriescolor := marginalcolor
        ticks := nothing
        linewidth --> 2
        subplot := 3
        framestyle := :none
        y := (length(x):-marginaldownsample:1)[1:marginallength]
        x := (x[1:marginaldownsample:length(x)])[1:marginallength]
    end

    @series begin
        seriestype := :path
        if cmode == :velocity
            line_z := velocity
        end
        linecolor := trajectorycolor
        linealpha := mainalpha
        linewidth --> 0.5
        xlims --> xlims
        ylims --> ylims
        x := x
        y := y
    end
end

@shorthands marginaltrajectory3
@recipe function f(::Type{Val{:marginaltrajectory3}}, plt::AbstractPlot; buffer=0.5, linewidth=1.0, colormode=:velocity, marginalcolor=:black, mainalpha=0.75, trajectorycolor=:black)
    cmode = colormode
    x, y, z = plotattributes[:x], plotattributes[:y], plotattributes[:z]
    i = isfinite.(x) .& isfinite.(y) .& isfinite.(z)
    x, y, z = x[i], y[i], z[i]
    if all(haskey.((plotattributes,), (:xlims, :ylims, :zlims)))
        xlims, ylims, zlims = plotattributes[:xlims], plotattributes[:ylims], plotattributes[:zlims]
    else
        xlims, ylims , zlims = widen.([extrema(x), extrema(y), extrema(z)], (buffer,))
    end
    legend --> :none
    xlims --> xlims
    ylims --> ylims
    zlims --> zlims
    # TODO Add coloring by 4th variable
    if cmode == :velocity
        velocity = sqrt.(sum([(r[2:end] .- collect(r[1:end-1])).^2 for r ∈ [x, y, z]], dims=1)[1])
    elseif !isnothing(cmode) && cmode != :none
        @error "Not a supported cmode"
    end
    @series begin
        if isnothing(marginalcolor) && cmode == :velocity
            line_z := velocity
        end
        seriestype := :path
        linecolor := marginalcolor
        linealpha := 0.1
        linewidth := 0.5*linewidth
        x := fill(xlims[1], length(x))
        y := y
        z := z
    end

    @series begin
        if isnothing(marginalcolor) && cmode == :velocity
            line_z := velocity
        end
        seriestype := :path
        linecolor := marginalcolor
        linealpha := 0.1
        linewidth := 0.5*linewidth
        x := x
        y := fill(ylims[2], length(y))
        z := z
    end

    @series begin
        if isnothing(marginalcolor) && cmode == :velocity
            line_z := velocity
        end
        linecolor := marginalcolor
        seriestype := :path
        linealpha := 0.1
        linewidth := 0.5*linewidth
        x := x
        y := y
        z := fill(zlims[1], length(z))
    end

    @series begin
        if cmode == :velocity
            line_z := velocity
        end
        seriescolor := trajectorycolor
        seriestype := :path
        linealpha := mainalpha
        linewidth --> linewidth
        x := x
        y := y
        z := z
    end

end

function marginaltrajectory(S::Process; vars=1:length(getX0(S)), marginaldownsample=1, marginallength=2500, trajectorycolor=:turbo, marginalcolor=:black, size=(700, 600), kwargs...)
    if length(vars) == 1
        t = times(S)
        x = timeseries(S, vars)
        t = t[1:marginaldownsample:end][1:marginallength]
        x = x[1:marginaldownsample:end][1:marginallength]
        plot(t, x; xlabel="t", ylabel="x", st=:tsdensity, framestyle=nothing, xticks = optimize_ticks(extrema(t)...; k_min = 3, k_max = 5)[1], yticks = optimize_ticks(extrema(x)...; k_min = 3, k_max = 5)[1], kwargs...)
    elseif length(vars) == 2 # 2D
        marginaltrajectory2(S; vars, marginaldownsample, marginallength, trajectorycolor, marginalcolor, size,
        xticks=optimize_ticks(extrema(timeseries(S, 1))...; k_min = 3, k_max = 5)[1],
        yticks=optimize_ticks(extrema(timeseries(S, 2))...; k_min = 3, k_max = 5)[1],
        kwargs...)
    else # 3D
        marginaltrajectory3(S; vars, marginaldownsample, marginallength, mainalpha=0.25, linecolor=trajectorycolor, marginalcolor, size,
        xticks=optimize_ticks(widen(extrema(timeseries(S, 1)), 0.5)...; k_min = 3, k_max = 4)[1],
        yticks=optimize_ticks(widen(extrema(timeseries(S, 2)), 0.5)...; k_min = 3, k_max = 4)[1],
        zticks=optimize_ticks(widen(extrema(timeseries(S, 3)), 0.5)...; k_min = 3, k_max = 4)[1],
        kwargs...)
    end
end
export marginaltrajectory


text:
# TODO Check faster implementations in computational methods for local regression
# William S. Cleveland & E. Grosse https://link.springer.com/article/10.1007/BF01890836
function ghat(x::Float64;A,b,d=2,q,rho)

    xv = A[:,d]
    yv = b
    ## lambda q
    n = length(xv)
    q = min(q,n)
    xvx = @. abs(xv-x)
    qidx = sortperm(xvx)[1:q]
    qdist = abs(xv[last(qidx)]-x)*max(1,q/n)

    ## upsilon
    w = zeros(n)
    for wi in qidx
        aq = abs(xv[wi]-x)/qdist
        w[wi] = max((1-aq^3)^3,0)
    end

    A = @. A*(w*rho)
    b = @. b*(w*rho)

    lsq_x = A \ b

    d == 1 ? [x,1.0]'*lsq_x : [x^2.0,x,1.0]'*lsq_x

end

"""
Package: Forecast
    loess(xv,yv;
          d=2,
          q=Int64(round(3/4*length(xv))),
          rho=repeat([1.0],inner=length(xv)),
          predict = xv)
Smooth a vector of observations using locally weighted regressions.
Although loess can be used to smooth observations for any given number of independent variables, this implementation is univariate. The speed of loess can be greatly increased by using fast aproximations for the linear fitting calculations, however this implementation calculates only exact results.
The loess functionality and nomenclature follows the descriptions in:
"STL: A Seasonal, Trend Decomposition Procedure Based on Loess"
Robert B. Cleveland, William S. Cleveland, Jean E. McRae, and Irma Terpenning.
Journal of Official Statistics Vol. 6. No. 1, 1990, pp. 3-73 (c) Statistics Sweden.
# Arguments
- `xv`: Observations' support.
- `yv`: Observation values.
- `d`: Degree of the linear fit, it accepts values 1 or 2.
- `q`: As q increases loess becomes smoother, when q tends to infinity loess tends to an ordinary least square poynomial fit of degree `d`. It defaults to the rounding of 3/4 of xv's length.
- `rho`: Weights expressing the reliability of the observations (e.g. if yi had variances sigma^2*ki where ki where known, the rhoi could be 1/ki). It defaults to 1.0.
- `predict`: Vector containing the real values to be predicted, by default predicts xv.
# Returns
The loess values for the values contained in `predict`.
# Examples
```julia-repl
julia> loess(rand(5), rand(5); predict=rand(10))
10-element Array{Float64,1}:
[...]
```
"""
function loess(xv,yv;
               d=2,
               q=Int64(round(3/4*length(xv))),
               rho=repeat([1.0],inner=length(xv)),
               predict = xv)

    @assert (d==1) | (d==2) "Linear Regression must be of degree 1 or 2"
    @assert length(findall(x -> ismissing(x), xv)) == 0  "xv should not contain missing values"

    myi = findall(x -> !ismissing(x),yv)
    xv = xv[myi]
    yv = yv[myi]
    rho = rho[myi]

    res = zeros(length(predict))

    ## Ax = b
    A = hcat(xv,repeat([1.0],inner=length(xv)))
    b = yv
    d == 2 ? A = hcat(xv .^ 2.0, A) : nothing

    for (i,xi) in enumerate(predict)
        res[i] = ghat(xi;A=A,b=b,d=d,q=q,rho=rho)
    end
    res
end


text:
struct DmlType{T} end
struct DdlType{T} end

function make_url(::Type{DmlType{:_index}}, info::Esinfo, index::AbstractString, type::AbstractString,id::AbstractString)
	"$(info.url)/$index/$type/$id"
end

function make_url(::Type{DmlType{:_index}}, info::Esinfo, index::AbstractString, type::AbstractString)
	"$(info.url)/$index/$type"
end

function make_url(::Type{DmlType{:_create}}, info::Esinfo, index::AbstractString, type::AbstractString,id::AbstractString)
	"$(info.url)/$index/$type/$id/_create"
end

function make_url(::Type{DmlType{:_update}}, info::Esinfo, index::AbstractString, type::AbstractString,id::AbstractString)
	"$(info.url)/$index/$type/$id/_update"
end

function make_url(::Type{DmlType{:_delete}}, info::Esinfo, index::AbstractString, type::AbstractString,id::AbstractString)
	"$(info.url)/$index/$type/$id"
end

function make_url(::Type{DmlType{:_head}}, info::Esinfo, index::AbstractString, type::AbstractString,id::AbstractString)
	"$(info.url)/$index/$type/$id"
end

function make_url(::Type{DmlType{:_get}}, info::Esinfo, index::AbstractString, type::AbstractString,id::AbstractString)
	"$(info.url)/$index/$type/$id"
end

function make_url(::Type{DmlType{:_get_source}}, info::Esinfo, index::AbstractString, type::AbstractString,id::AbstractString)
	"$(info.url)/$index/$type/$id/_source"
end

function make_url(::Type{DdlType{:_explain}}, info::Esinfo, index::AbstractString, type::AbstractString,id::AbstractString)
	"$(info.url)/$index/$type/$id/_explain"
end

function make_url(::Type{DdlType{:_delete}}, info::Esinfo, index::AbstractString )
	"$(info.url)/$index/_delete_by_query"
end

function make_url(::Type{DdlType{:_delete}}, info::Esinfo, index::AbstractString , type::AbstractString)
	"$(info.url)/$index/$type/_delete_by_query"
end

function make_url(::Type{DdlType{:_update}}, info::Esinfo, index::AbstractString )
	"$(info.url)/$index/_update_by_query"
end

function make_url(::Type{DdlType{:_update}}, info::Esinfo, index::AbstractString , type::AbstractString)
	"$(info.url)/$index/$type/_update_by_query"
end

function make_url(::Type{DdlType{:_delete_script}}, info::Esinfo, id::AbstractString )
	"$(info.url)/_scripts/$id"
end

function make_url(::Type{DdlType{:_get_script}}, info::Esinfo, id::AbstractString )
	"$(info.url)/_scripts/$id"
end

function make_url(::Type{DdlType{:_field_caps}}, info::Esinfo, index::AbstractString )
	"$(info.url)/$index/_field_caps"
end

function make_url(::Type{DdlType{:_field_caps}}, info::Esinfo  )
	"$(info.url)/_field_caps"
end

function make_url(::Type{DdlType{:_delete_rethrottle}}, info::Esinfo, task_id::AbstractString )
	"$(info.url)/_delete_by_query/$task_id/_rethrottle"
end

function make_url(::Type{DdlType{:_update_rethrottle}}, info::Esinfo, task_id::AbstractString )
	"$(info.url)/_update_by_query/$task_id/_rethrottle"
end

function make_url(::Type{DdlType{:_mget}}, info::Esinfo  )
	"$(info.url)/_mget"
end

function make_url(::Type{DdlType{:_mget}}, info::Esinfo, index::AbstractString )
	"$(info.url)/$index/_mget"
end

function make_url(::Type{DdlType{:_mget}}, info::Esinfo, index::AbstractString , type::AbstractString)
	"$(info.url)/$index/$type/_mget"
end

function make_url(::Type{DdlType{:_mtermvectors}}, info::Esinfo  )
	"$(info.url)/_mtermvectors"
end

function make_url(::Type{DdlType{:_mtermvectors}}, info::Esinfo, index::AbstractString )
	"$(info.url)/$index/_mtermvectors"
end

function make_url(::Type{DdlType{:_mtermvectors}}, info::Esinfo, index::AbstractString , type::AbstractString)
	"$(info.url)/$index/$type/_mtermvectors"
end

function make_url(::Type{DdlType{:_ping}}, info::Esinfo)
	"$(info.url)/"
end

function make_url(::Type{DdlType{:_put_script}}, info::Esinfo, id::AbstractString )
	"$(info.url)/_scripts/$id"
end

function make_url(::Type{DdlType{:_put_script}}, info::Esinfo, id::AbstractString, context::AbstractString )
	"$(info.url)/_scripts/$id/$context"
end

function make_url(::Type{DdlType{:_rank_eval}}, info::Esinfo, index::AbstractString )
	"$(info.url)/_index/_rank_eval"
end

function make_url(::Type{DdlType{:_rank_eval}}, info::Esinfo )
	"$(info.url)/_rank_eval"
end

function make_url(::Type{DdlType{:_reindex}}, info::Esinfo)
	"$(info.url)/_reindex"
end

function make_url(::Type{DdlType{:_reindex_rethrottle}}, info::Esinfo, task_id::AbstractString )
	"$(info.url)/_reindex/$task_id/_rethrottle"
end

function make_url(::Type{DdlType{:_render_search_template}}, info::Esinfo, id::AbstractString )
	"$(info.url)/_render/template/$id"
end

function make_url(::Type{DdlType{:_render_search_template}}, info::Esinfo)
	"$(info.url)/_render/template"
end

function make_url(::Type{DdlType{:_scripts_painless_execute}}, info::Esinfo)
	"$(info.url)/_scripts/painless/_execute"
end

function make_url(::Type{DdlType{:_search_shards}}, info::Esinfo)
	"$(info.url)/_search_shards"
end

function make_url(::Type{DdlType{:_search_shards}}, info::Esinfo, index::AbstractString)
	"$(info.url)/$index/_search_shards"
end

function es_delete(info::Esinfo, index::AbstractString, type::AbstractString,id::AbstractString ; kw...)

	url   = make_url(DmlType{:_delete}, info, index ,type, id )
	@esdelete(info, url ,  Dict(kw...) )

end

function es_delete_by_query(info::Esinfo, scripts_id::AbstractString ; kw...)

	query = Dict(kw...)
	url   = make_url(DdlType{:_delete_script}, info, scripts_id  )
	@esdelete(info, url ,  Dict(kw...) )

end

function es_exists(info::Esinfo, index::AbstractString, type::AbstractString, id::AbstractString ; kw...)

	url   = make_url(DmlType{:_head}, info, index ,type, id )
	@eshead info  url   Dict(kw...)

end

function es_exists_source(info::Esinfo, index::AbstractString, type::AbstractString, id::AbstractString ; kw...)

	es_exists(info, index ,type, "$id/_source"  ; kw...)

end

function es_ping(info::Esinfo )

	url   = make_url(DdlType{:_ping}, info )
	@eshead info  url  Dict()

end

@genfunction "POST" es_create DmlType{:_create} 1 index type id
@genfunction "POST" es_delete_by_query_rethrottle DdlType{:_delete_rethrottle} 0 task_id
@genfunction "POST" es_update_by_query_rethrottle DdlType{:_update_rethrottle} 0 task_id
@genfunction "POST" es_delete_by_query DdlType{:_delete} 1 index type
@genfunction "POST" es_delete_by_query DdlType{:_delete} 1 index
@genfunction "POST" es_update_by_query DdlType{:_update} 1 index type
@genfunction "POST" es_update_by_query DdlType{:_update} 1 index
@genfunction "POST" es_explain DdlType{:_explain} 1 index type id
@genfunction "POST" es_field_caps DdlType{:_field_caps} 0 index
@genfunction "POST" es_field_caps DdlType{:_field_caps} 0
@genfunction "GET" es_get DmlType{:_get} 0 index type id
@genfunction "GET" es_getscript DdlType{:_get_script} 0 scripts_id
@genfunction "GET" es_getsource DmlType{:_get_source} 0 index type id
@genfunction "POST" es_index DmlType{:_index} 1 index type id
@genfunction "POST" es_index DmlType{:_index} 1 index type
@genfunction "POST" es_update DmlType{:_update} 1 index type id
@genfunction "POST" es_update DdlType{:_mget} 1 index type
@genfunction "POST" es_update DdlType{:_mget} 1 index
@genfunction "POST" es_update DdlType{:_mget} 1
@genfunction "POST" es_mtermvectors DdlType{:_mtermvectors} 1 index type
@genfunction "POST" es_mtermvectors DdlType{:_mtermvectors} 1 index
@genfunction "POST" es_mtermvectors DdlType{:_mtermvectors} 1
@genfunction "PUT" es_put_script DdlType{:_put_script} 1 scripts_id context
@genfunction "PUT" es_put_script DdlType{:_put_script} 1 scripts_id
@genfunction "PUT" es_rank_eval DdlType{:_rank_eval} 1 index
@genfunction "PUT" es_rank_eval DdlType{:_rank_eval} 1
@genfunction "POST" es_reindex DdlType{:_reindex} 1
@genfunction "POST" es_reindex_rethrottle DdlType{:_reindex_rethrottle} 0 task_id
@genfunction "POST" es_render_search_template DdlType{:_render_search_template} 1 id
@genfunction "POST" es_render_search_template DdlType{:_render_search_template} 1
@genfunction "POST" es_scripts_painless_execute DdlType{:_scripts_painless_execute} 1
@genfunction "POST" es_search_shards DdlType{:_search_shards} 0 index
@genfunction "POST" es_search_shards DdlType{:_search_shards} 0

max_stars_repo_name: frankill/ES

max_issues_repo_name: frankill/ES

max_forks_repo_name: frankill/ES


text:
using IntervalArithmetic
using StaticArrays

using Test

let X, A  # avoid problems with global variables

@testset "Operations on boxes" begin
    A = IntervalBox(1..2, 3..4)
    B = IntervalBox(0..2, 3..6)
    s = @SVector [1, 2]

    @test 2*A == A*2 == IntervalBox(2..4, 6..8)
    @test typeof(2*A) == IntervalBox{2, Float64}
    @test A + B == IntervalBox(1..4, 6..10)
    @test A + B.v == IntervalBox(1..4, 6..10)
    @test A.v + B == IntervalBox(1..4, 6..10)
    @test A + s == IntervalBox(2..3, 5..6)
    @test A - B == IntervalBox(-1..2, -3..1)
    @test A.v - B == IntervalBox(-1..2, -3..1)
    @test A - B.v == IntervalBox(-1..2, -3..1)
    @test A - s == IntervalBox(0..1, 1..2)
    @test 2 + A == IntervalBox(3..4,5..6)
    @test A + 2 == IntervalBox(3..4,5..6)
    @test -A == IntervalBox((-2)..(-1), (-4)..(-3))
    @test 2 - A == IntervalBox(0..1, (-2)..(-1))
    @test B - 2 == IntervalBox((-2)..0, 1..4)
    @test dot(A, B) == @interval(9, 28)
    @test dot(A, B.v) == @interval(9, 28)
    @test dot(A.v, B) == @interval(9, 28)
    @test A .* B == IntervalBox(0..4, 9..24)
    @test A ./ A == IntervalBox((0.5)..2, (0.75)..(4/3))
    @test 1 ./ B == IntervalBox((0.5)..Inf, (1/6)..(1/3))
    @test B ./ 1 == B
    @test A .^ 2 == IntervalBox(1..4, 9..16)
    @test B .^ 0.5 == IntervalBox(@interval(0,sqrt(2)), @interval(sqrt(3),sqrt(6)))

    @test A ⊆ A
    @test A ⊆ B
    @test A.v ⊆ B
    @test A ⊂ B.v
    @test !(A ⊂ A)
    @test A ⊂ B
    @test A.v ⊂ B
    @test A ⊆ B.v
    @test !(B ⊆ A)
    @test !(B ⊂ A)
    @test B ⊇ B
    @test B ⊇ A
    @test B.v ⊇ A
    @test B ⊇ A.v
    @test !(B ⊃ B)
    @test B ⊃ A
    @test B.v ⊃ A
    @test B ⊃ A.v
    @test !(A ⊇ B)
    @test !(A ⊃ B)

    @test A ∩ B == A
    @test A.v ∩ B == A
    @test A ∩ B.v == A

    @test A ∪ B == B
    @test A.v ∪ B == B
    @test A ∪ B.v == B

    X = IntervalBox(1..2, 3..4)
    Y = IntervalBox(3..4, 3..4)

    @test isempty(X ∩ Y)
    @test X ∪ Y == IntervalBox(1..4, 3..4)

    @test !contains_zero(X ∩ Y)
    @test contains_zero( (-1..1) × (-1..1) )
    @test !isinf( (-1..1) × (0..Inf) )
    @test isinf( entireinterval() × entireinterval() )

    X = IntervalBox(2..4, 3..5)
    Y = IntervalBox(3..5, 4..17)
    @test X ∩ Y == IntervalBox(3..4, 4..5)

    v = [@interval(i, i+1) for i in 1:10]
    V = IntervalBox(v...)
    @test length(V) == 10

    Y = IntervalBox(1..2)  # single interval
    @test isa(Y, IntervalBox)
    @test length(Y) == 1
    @test Y == IntervalBox( (Interval(1., 2.),) )
    @test typeof(Y) == IntervalBox{1, Float64}
end

@testset "Functions on boxes" begin
    A = IntervalBox(1..2, 3..4)

    @test exp.(A) == IntervalBox(exp(A[1]), exp(A[2]))
    @test typeof(exp.(A)) == IntervalBox{2,Float64}
    @test log.(A) == IntervalBox(log(A[1]), log(A[2]))
    @test sqrt.(A) == IntervalBox(sqrt(A[1]), sqrt(A[2]))
    @test inv.(A) == IntervalBox(inv(A[1]), inv(A[2]))
end

@testset "setdiff for IntervalBox" begin
    X = IntervalBox(2..4, 3..5)
    Y = IntervalBox(3..5, 4..6)
    @test setdiff(X, Y) == [ IntervalBox(3..4, 3..4),
                              IntervalBox(2..3, 3..5) ]

    @test setdiff(X.v, Y) == [ IntervalBox(3..4, 3..4),
                              IntervalBox(2..3, 3..5) ]

    @test setdiff(X, Y.v) == [ IntervalBox(3..4, 3..4),
                              IntervalBox(2..3, 3..5) ]

    X = IntervalBox(2..5, 3..6)
    Y = IntervalBox(-10..10, 4..5)
    @test setdiff(X, Y) == [ IntervalBox(2..5, 3..4),
                              IntervalBox(2..5, 5..6) ]

    X = IntervalBox(2..5, 3..6)
    Y = IntervalBox(4..6, 4..5)
    @test setdiff(X, Y) == [ IntervalBox(4..5, 3..4),
                              IntervalBox(4..5, 5..6),
                              IntervalBox(2..4, 3..6) ]

    X = IntervalBox(2..5, 3..6)
    Y = IntervalBox(3..4, 4..5)
    @test setdiff(X, Y) == [ IntervalBox(3..4, 3..4),
                              IntervalBox(3..4, 5..6),
                              IntervalBox(2..3, 3..6),
                              IntervalBox(4..5, 3..6) ]

    X = IntervalBox(2..5, 3..6)
    Y = IntervalBox(2..4, 10..20)
    @test setdiff(X, Y) == typeof(X)[X]

    X = IntervalBox(2..5, 3..6)
    Y = IntervalBox(-10..10, -10..10)
    @test setdiff(X, Y) == typeof(X)[]

    X = IntervalBox(1..4, 3..6, 7..10)
    Y = IntervalBox(2..3, 4..5, 8..9)
    @test setdiff(X, Y) == [ IntervalBox(2..3, 4..5, 7..8),
                              IntervalBox(2..3, 4..5, 9..10),
                              IntervalBox(2..3, 3..4, 7..10),
                              IntervalBox(2..3, 5..6, 7..10),
                              IntervalBox(1..2, 3..6, 7..10),
                              IntervalBox(3..4, 3..6, 7..10) ]

end

@testset "mid, diam, × for IntervalBox" begin
    X = (0..2) × (3..5)
    @test length(X) == 2
    @test X == IntervalBox(Interval(0, 2), Interval(3, 5))

    @test diam(X) == 2
    @test mid(X) == [1, 4]

    Y = X × (4..8)
    @test isa(Y, IntervalBox)
    @test length(Y) == 3
    @test Y == IntervalBox(Interval(0, 2), Interval(3, 5), Interval(4, 8))
    @test diam(Y) == 4

    Z = X × Y
    @test isa(Z, IntervalBox)
    @test length(Z) == 5
    @test Z == IntervalBox(Interval(0, 2), Interval(3, 5), Interval(0, 2), Interval(3, 5), Interval(4, 8))
    @test diam(Z) == 4

    Z = X × Y.v
    @test isa(Z, IntervalBox)
    @test length(Z) == 5
    @test Z == IntervalBox(Interval(0, 2), Interval(3, 5), Interval(0, 2), Interval(3, 5), Interval(4, 8))
    @test diam(Z) == 4

    Z = X.v × Y
    @test isa(Z, IntervalBox)
    @test length(Z) == 5
    @test Z == IntervalBox(Interval(0, 2), Interval(3, 5), Interval(0, 2), Interval(3, 5), Interval(4, 8))
    @test diam(Z) == 4

    @test mid(IntervalBox(0..1, 3), 0.75) == [0.75, 0.75, 0.75]
end

@testset "Constructing multidimensional IntervalBoxes" begin
    @test IntervalBox(1..2, Val(1)) == IntervalBox(1..2)
    @test IntervalBox(1..2, Val(2)) == (1..2) × (1..2)
    @test IntervalBox(1..2, Val(5)) == (1..2) × (1..2) × (1..2) × (1..2) × (1..2)

    @test IntervalBox(1..2, 3) == IntervalBox(1..2, Val(3))
    @test IntervalBox((1..2, 2..3)) == IntervalBox(1..2, 2..3)
    @test IntervalBox((1, 2)) == IntervalBox(1..1, 2..2)
    @test IntervalBox( (1, 2, 3) ) == IntervalBox(1..1, 2..2, 3..3)
    @test IntervalBox( (1, 2, 3.1) ) == IntervalBox(1..1, 2..2, interval(3.1))
    @test IntervalBox( SVector(1, 2, 3.1) ) == IntervalBox(1..1, 2..2, interval(3.1))
    @test IntervalBox( interval.(SVector(1, 2, 3.1)) ) == IntervalBox(1..1, 2..2, interval(3.1))
    @test IntervalBox(3) == IntervalBox(3..3)
    @test IntervalBox(1:5) == IntervalBox(1..1, 2..2, 3..3, 4..4, 5..5)
    @test IntervalBox([1:5...]) == IntervalBox(1..1, 2..2, 3..3, 4..4, 5..5)
    @test IntervalBox((1..2) × (2..3), 2) == (1..2) × (2..3) × (1..2) × (2..3)

    # construct from corners:
    @test IntervalBox(SVector(1, 2), SVector(3, 4)) == (1..3) × (2..4)
    @test IntervalBox(SVector(3, 4), SVector(1, 2)) == (1..3) × (2..4)

end

@testset "getindex and setindex" begin
    X = IntervalBox(3..4, 5..6)
    @test X[1] == 3..4
    @test X[2] == 5..6
    @test_throws BoundsError X[3]

    @test setindex(X, 5..5, 2) == IntervalBox(3..4, 5..5)
    @test_throws BoundsError setindex(X, 5..5, 3)
end

@testset "Iteration" begin
    X = IntervalBox(3..4, 5..6)
    Y = collect(X)
    @test Y == [3..4, 5..6]
    @test eltype(Y) == Interval{Float64}
end

@testset "Broadcasting" begin
    X = IntervalBox(3..4, 5..6)

    @test sin.(X) == IntervalBox(sin(X[1]), sin(X[2]))
    @test mid.(X) == SVector(mid(X[1]), mid(X[2]))
    @test diam.(X) == SVector(diam(X[1]), diam(X[2]))
end

@testset "∈" begin
    X = IntervalBox(3..4, 5..6)
    @test mid(X) ∈ X
    @test mid(X, 0.75) ∈ X

    @test (zero(mid(X)) ∈ X) == false
    @test zero(mid(X)) ∉ X

    @test_throws ArgumentError (3..4) ∈ X

    @test_throws ArgumentError [3..4, 5..6] ∈ X

end

@testset "Multiplication by a matrix" begin
    A = [1 2; 3 4]
    X = IntervalBox(1..2, 3..4)

    @test A * X == IntervalBox(7..10, 15..22)

    B = SMatrix{2,2}(A)

    @test B * X == IntervalBox(7..10, 15..22)

end

end


text:
module Search

include("beam_search.jl")
export beam_search, greedy_search

end # module


text:
using Test
import ClimaCore: Domains, Meshes, Topologies
import ClimaCore.Geometry: Geometry
using StaticArrays
using IntervalSets

function rectangular_grid(
    n1,
    n2,
    x1periodic,
    x2periodic;
    x1min = 0.0,
    x1max = 1.0,
    x2min = 0.0,
    x2max = 1.0,
)
    domain = Domains.RectangleDomain(
        Geometry.XPoint(x1min) .. Geometry.XPoint(x1max),
        Geometry.YPoint(x2min) .. Geometry.YPoint(x2max),
        x1periodic = x1periodic,
        x2periodic = x2periodic,
        x1boundary = x1periodic ? nothing : (:west, :east),
        x2boundary = x2periodic ? nothing : (:south, :north),
    )
    mesh = Meshes.RectilinearMesh(domain, n1, n2)
    return Topologies.Topology2D(mesh)
end

@testset "opposing face tests" begin
    @testset "assert correct element numbering" begin
        topology = rectangular_grid(1, 1, true, true)
        @test_throws BoundsError Topologies.opposing_face(topology, 0, 1)
        @test_throws BoundsError Topologies.opposing_face(topology, 2, 1)
    end

    @testset "assert correct face numbering" begin
        topology = rectangular_grid(1, 1, true, true)
        @test_throws AssertionError Topologies.opposing_face(topology, 1, 5)
        @test_throws AssertionError Topologies.opposing_face(topology, 1, 0)
    end

    @testset "1×1 element quad mesh with all periodic boundries" begin
        topology = rectangular_grid(1, 1, true, true)

        @test Topologies.opposing_face(topology, 1, 1) == (1, 3, true)
        @test Topologies.opposing_face(topology, 1, 2) == (1, 4, true)
        @test Topologies.opposing_face(topology, 1, 3) == (1, 1, true)
        @test Topologies.opposing_face(topology, 1, 4) == (1, 2, true)
    end

    @testset "1×1 element quad mesh with 1 periodic boundary" begin
        topology = rectangular_grid(1, 1, true, false)
        @test Topologies.opposing_face(topology, 1, 1) == (0, 1, false)
        @test Topologies.opposing_face(topology, 1, 2) == (1, 4, true)
        @test Topologies.opposing_face(topology, 1, 3) == (0, 2, false)
        @test Topologies.opposing_face(topology, 1, 4) == (1, 2, true)
    end

    @testset "1×1 element quad mesh with non-periodic boundaries" begin
        topology = rectangular_grid(1, 1, false, false)
        @test Topologies.opposing_face(topology, 1, 1) == (0, 3, false)
        @test Topologies.opposing_face(topology, 1, 2) == (0, 2, false)
        @test Topologies.opposing_face(topology, 1, 3) == (0, 4, false)
        @test Topologies.opposing_face(topology, 1, 4) == (0, 1, false)
    end

    @testset "2×2 element quad mesh with non-periodic boundaries" begin
        topology = rectangular_grid(2, 2, false, false)
        @test Topologies.opposing_face(topology, 1, 1) == (0, 3, false)
        @test Topologies.opposing_face(topology, 1, 2) == (2, 4, true)
        @test Topologies.opposing_face(topology, 1, 3) == (3, 1, true)
        @test Topologies.opposing_face(topology, 1, 4) == (0, 1, false)
        @test Topologies.opposing_face(topology, 2, 1) == (0, 3, false)
        @test Topologies.opposing_face(topology, 2, 2) == (0, 2, false)
        @test Topologies.opposing_face(topology, 2, 3) == (4, 1, true)
        @test Topologies.opposing_face(topology, 2, 4) == (1, 2, true)
    end
end

@testset "simple rectangular mesh interior faces iterator" begin

    @testset "1×1 element quad mesh with all periodic boundries" begin
        topology = rectangular_grid(1, 1, true, true)
        @test length(Topologies.interior_faces(topology)) == 2
        faces = collect(Topologies.interior_faces(topology))
        @test sort(faces) == sort([(1, 4, 1, 2, true), (1, 3, 1, 1, true)])
    end

    @testset "1×1 element quad mesh with 1 periodic boundary" begin
        topology = rectangular_grid(1, 1, true, false)
        @test length(Topologies.interior_faces(topology)) == 1
        faces = collect(Topologies.interior_faces(topology))
        @test faces == [(1, 4, 1, 2, true)]
    end

    @testset "1×1 element quad mesh with non-periodic boundaries" begin
        topology = rectangular_grid(1, 1, false, false)
        @test length(Topologies.interior_faces(topology)) == 0
        faces = collect(Topologies.interior_faces(topology))
        @test isempty(faces)
    end

    @testset "2×2 element quad mesh with non-periodic boundaries" begin
        topology = rectangular_grid(2, 2, false, false)
        @test length(Topologies.interior_faces(topology)) == 4
        faces = collect(Topologies.interior_faces(topology))
        @test sort(faces) == sort([
            (2, 4, 1, 2, true),
            (3, 1, 1, 3, true),
            (4, 4, 3, 2, true),
            (4, 1, 2, 3, true),
        ])
    end
end

@testset "simple rectangular mesh boundry faces iterator" begin
    @testset "1×1 element quad mesh with all periodic boundries" begin
        topology = rectangular_grid(1, 1, true, true)
        @test isempty(Topologies.boundary_tags(topology))
    end

    @testset "1×1 element quad mesh with 1 periodic boundary" begin
        topology = rectangular_grid(1, 1, true, false)
        @test length(Topologies.boundary_faces(topology, 1)) == 1
        @test length(Topologies.boundary_faces(topology, 2)) == 1
        @test collect(Topologies.boundary_faces(topology, 1)) == [(1, 1)]
        @test collect(Topologies.boundary_faces(topology, 2)) == [(1, 3)]
    end

    @testset "1×1 element quad mesh with non-periodic boundaries" begin
        topology = rectangular_grid(1, 1, false, false)
        @test length(Topologies.boundary_faces(topology, 1)) == 1
        @test length(Topologies.boundary_faces(topology, 2)) == 1
        @test length(Topologies.boundary_faces(topology, 3)) == 1
        @test length(Topologies.boundary_faces(topology, 4)) == 1
        @test collect(Topologies.boundary_faces(topology, 1)) == [(1, 4)]
        @test collect(Topologies.boundary_faces(topology, 2)) == [(1, 2)]
        @test collect(Topologies.boundary_faces(topology, 3)) == [(1, 1)]
        @test collect(Topologies.boundary_faces(topology, 4)) == [(1, 3)]
        @test Topologies.boundary_tag(topology, :west) == 1
        @test Topologies.boundary_tag(topology, :east) == 2
        @test Topologies.boundary_tag(topology, :south) == 3
        @test Topologies.boundary_tag(topology, :north) == 4
    end

    @testset "2×3 element quad mesh with non-periodic boundaries" begin
        topology = rectangular_grid(2, 3, false, false)
        @test length(Topologies.boundary_faces(topology, 1)) == 3
        @test length(Topologies.boundary_faces(topology, 2)) == 3
        @test length(Topologies.boundary_faces(topology, 3)) == 2
        @test length(Topologies.boundary_faces(topology, 4)) == 2

        @test collect(Topologies.boundary_faces(topology, 1)) ==
              [(1, 4), (3, 4), (5, 4)]
        @test collect(Topologies.boundary_faces(topology, 2)) ==
              [(2, 2), (4, 2), (6, 2)]
        @test collect(Topologies.boundary_faces(topology, 3)) ==
              [(1, 1), (2, 1)]
        @test collect(Topologies.boundary_faces(topology, 4)) ==
              [(5, 3), (6, 3)]
    end
end

@testset "simple rectangular mesh vertex iterator" begin

    @testset "1×1 element quad mesh with all periodic boundries" begin
        topology = rectangular_grid(1, 1, true, true)
        # this has 1 global vertex
        @test length(Topologies.vertices(topology)) == 1
        V = collect(Topologies.vertices(topology))
        @test V[1] isa Topologies.Vertex
        @test sort(collect(V[1])) == [(1, 1), (1, 2), (1, 3), (1, 4)]
    end

    @testset "1×1 element quad mesh with 1 periodic boundary" begin
        topology = rectangular_grid(1, 1, true, false)
        @test length(Topologies.vertices(topology)) == 2
        V = collect(Topologies.vertices(topology))
        @test V[1] isa Topologies.Vertex
        @test sort(collect(V[1])) == [(1, 1), (1, 2)]
        @test sort(collect(V[2])) == [(1, 3), (1, 4)]
    end

    @testset "1×1 element quad mesh with non-periodic boundaries" begin
        topology = rectangular_grid(1, 1, false, false)
        @test length(Topologies.vertices(topology)) == 4
        V = collect(Topologies.vertices(topology))
        @test V[1] isa Topologies.Vertex
        @test collect(V[1]) == [(1, 1)]
        @test collect(V[2]) == [(1, 2)]
        @test collect(V[3]) == [(1, 3)]
        @test collect(V[4]) == [(1, 4)]
    end

    @testset "2×3 element quad mesh with non-periodic boundaries" begin
        topology = rectangular_grid(2, 3, false, false)
        @test length(Topologies.vertices(topology)) == 3 * 4
        V = collect(Topologies.vertices(topology))
        @test length(V[1]) == 1
        @test collect(V[1]) == [(1, 1)]
    end

    @testset "2×3 element quad mesh with periodic boundaries" begin
        topology = rectangular_grid(2, 3, true, true)
        @test length(Topologies.vertices(topology)) == 2 * 3
        V = collect(Topologies.vertices(topology))
        @test length(V) == 6
        @test sort([sort(collect(vert)) for vert in V]) == sort([
            sort([(1, 1), (2, 2), (5, 4), (6, 3)]),
            sort([(2, 1), (1, 2), (6, 4), (5, 3)]),
            sort([(3, 1), (4, 2), (1, 4), (2, 3)]),
            sort([(4, 1), (3, 2), (2, 4), (1, 3)]),
            sort([(5, 1), (6, 2), (3, 4), (4, 3)]),
            sort([(6, 1), (5, 2), (4, 4), (3, 3)]),
        ])
    end
end

@testset "simple rectangular mesh coordinates" begin

    @testset "1×1 element quad mesh with all periodic boundries" begin
        topology = rectangular_grid(1, 1, true, true)
        c1, c2, c3, c4 = Topologies.vertex_coordinates(topology, 1)
        @test c1 == Geometry.XYPoint(0.0, 0.0)
        @test c2 == Geometry.XYPoint(1.0, 0.0)
        @test c3 == Geometry.XYPoint(1.0, 1.0)
        @test c4 == Geometry.XYPoint(0.0, 1.0)
    end

    @testset "1×1 element quad mesh with non-periodic boundries" begin
        topology =
            rectangular_grid(1, 1, false, false; x1min = -1.0, x2min = -1.0)
        c1, c2, c3, c4 = Topologies.vertex_coordinates(topology, 1)
        @test c1 == Geometry.XYPoint(-1.0, -1.0)
        @test c2 == Geometry.XYPoint(1.0, -1.0)
        @test c3 == Geometry.XYPoint(1.0, 1.0)
        @test c4 == Geometry.XYPoint(-1.0, 1.0)
    end

    @testset "2×4 element quad mesh with non-periodic boundaries" begin
        topology = rectangular_grid(2, 4, false, false)
        c1, c2, c3, c4 = Topologies.vertex_coordinates(topology, 1)
        @test c1 == Geometry.XYPoint(0.0, 0.0)
        @test c2 == Geometry.XYPoint(0.5, 0.0)
        @test c3 == Geometry.XYPoint(0.5, 0.25)
        @test c4 == Geometry.XYPoint(0.0, 0.25)

        c1, c2, c3, c4 = Topologies.vertex_coordinates(topology, 8)
        @test c1 == Geometry.XYPoint(0.5, 0.75)
        @test c2 == Geometry.XYPoint(1.0, 0.75)
        @test c3 == Geometry.XYPoint(1.0, 1.0)
        @test c4 == Geometry.XYPoint(0.5, 1.0)
    end

    @testset "check coordinate type accuracy" begin
        topology = rectangular_grid(
            3,
            1,
            false,
            false,
            x1min = big(0.0),
            x1max = big(1.0),
            x2min = big(0.0),
            x2max = big(1.0),
        )
        c1, c2, c3, c4 = Topologies.vertex_coordinates(topology, 1)
        @test eltype(c2) == BigFloat
        @test getfield(c2, 1) ≈ big(1.0) / big(3.0) rtol = eps(BigFloat)
        @test getfield(c2, 2) == 0.0
    end
end


text:
# ---
# title: 1680. Concatenation of Consecutive Binary Numbers
# id: problem1680
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Math
# link: <https://leetcode.com/problems/concatenation-of-consecutive-binary-numbers/description/>
# hidden: true
# ---
#
# Given an integer `n`, return _the **decimal value** of the binary string
# formed by concatenating the binary representations of _`1` _to_`n` _in order,
# **modulo**_`109 + 7`.
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 1
#     Output: 1
#     Explanation: "1" in binary corresponds to the decimal value 1.
#
#
# **Example 2:**
#
#
#
#     Input: n = 3
#     Output: 27
#     Explanation: In binary, 1, 2, and 3 corresponds to "1", "10", and "11".
#     After concatenating them, we have "11011", which corresponds to the decimal value 27.
#
#
# **Example 3:**
#
#
#
#     Input: n = 12
#     Output: 505379714
#     **Explanation** : The concatenation results in "1101110010111011110001001101010111100".
#     The decimal value of that is 118505380540.
#     After modulo 109 + 7, the result is 505379714.
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 105`
#
#
## @lc code=start
using LeetCode

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


text:
using BernsteinExpansions, Compat, BenchmarkTools

SUITE = BenchmarkGroup()  # parent BenchmarkGroup to contain our suite

SUITE["Univariate"] = SUni = BenchmarkGroup()
SUni["Univariate monomial, floating point input, degree 3"] = @benchmarkable univariate(3, 3, 1.0, 2.0)
SUni["Univariate monomial, floating point input, degree 5"] = @benchmarkable univariate(1, 5, 0.0, 1.0)

SUITE["Multivariate"] = SMulti = BenchmarkGroup()
k = [1,2,3,4,5,2,3,4]
l = [5,5,5,5,5,5,5,5]

low = [1//1,1//1,1//1,1//1,1//1,1//1,1//1,1//1]
high = [2//1,2//1,2//1,2//1,2//1,2//1,2//1,2//1]
SMulti["Rational multivariate expansion"] = @benchmarkable multivariate($k, $l, $low, $high)

low = convert(Vector{Float64}, low)
high = convert(Vector{Float64}, high)
SMulti["Float64 multivariate expansion"] = @benchmarkable multivariate($k, $l, $low, $high)


text:
module ForwardBackwardOptim

using Reexport
@reexport using EmpiricalRisks

import Base.LinAlg: axpy!

# source files

include("callback.jl")
include("fbs.jl")
include("fista.jl")
include("fasta.jl")

export  fbs!, fbs,
        fista!, fista,
        fasta!, fasta

end # module


text:
using BinDeps
using Libdl

@BinDeps.setup

version = "1.20"
url = "https://github.com/google/leveldb/archive/v$(version).tar.gz"

libleveldb = library_dependency("libleveldb")

provides(Sources, URI(url), libleveldb, unpacked_dir="leveldb-$(version)")

leveldbbuilddir = BinDeps.builddir(libleveldb)
leveldbsrcdir = joinpath(BinDeps.depsdir(libleveldb),"src", "leveldb-$version")
leveldblibdir = BinDeps.libdir(libleveldb)
leveldblibfile = joinpath(leveldblibdir,libleveldb.name*".$(Libdl.dlext)")

provides(BuildProcess,
    (@build_steps begin
        GetSources(libleveldb)
        CreateDirectory(leveldblibdir)
        @build_steps begin
            ChangeDirectory(leveldbsrcdir)
            FileRule(leveldblibfile, @build_steps begin
                     `make`
                     `cp out-shared/libleveldb.$(Libdl.dlext).$(version) $(leveldblibfile)`
            end)
        end
    end), libleveldb, os = :Unix)

@BinDeps.install Dict(:libleveldb => :libleveldb)


text:
# Periodic Function

@doc """
# Description
Constructor for the Periodic kernel (covariance)

k(x,x') = σ²exp(-2sin²(π|x-x'|/p)/ℓ²)
# Arguments:
* `ll::Vector{Float64}`: Log of length scale ℓ
* `lσ::Float64`        : Log of the signal standard deviation σ
* `lp::Float64`        : Log of the period
""" ->
type Periodic <: Isotropic
    ℓ2::Float64
    σ2::Float64
    p::Float64      # Log of period
    Periodic(ll::Float64, lσ::Float64, lp::Float64) = new(exp(2*ll), exp(2*lσ), exp(lp))
end

get_params(pe::Periodic) = Float64[log(pe.ℓ2)/2.0, log(pe.σ2)/2.0, log(pe.p)]
get_param_names(pe::Periodic) = [:ll, :lσ, :lp]
num_params(pe::Periodic) = 3
metric(pe::Periodic) = Euclidean()

function set_params!(pe::Periodic, hyp::Vector{Float64})
    length(hyp) == 3 || throw(ArgumentError("Periodic function has only three parameters"))
    pe.ℓ2, pe.σ2 = exp(2.0*hyp[1:2])
    pe.p = exp(hyp[3])
end

cov(pe::Periodic, r::Float64) = pe.σ2*exp(-2.0/pe.ℓ2*sin(π*r/pe.p)^2)

function grad_kern(pe::Periodic, x::Vector{Float64}, y::Vector{Float64})
    r = distance(pe, x, y)

    dK_dlℓ   = 4.0*pe.σ2*(sin(π*r/pe.p)^2/pe.ℓ2)*exp(-2/pe.ℓ2*sin(π*r/pe.p)^2)
    dK_dlσ = 2.0*pe.σ2*exp(-2/pe.ℓ2*sin(π*r/pe.p)^2)
    dK_dlp     = 4.0/pe.ℓ2*pe.σ2*(π*r/pe.p)*sin(π*r/pe.p)*cos(π*r/pe.p)*exp(-2/pe.ℓ2*sin(π*r/pe.p)^2)
    dK_theta = [dK_dlℓ, dK_dlσ, dK_dlp]
    return dK_theta
end

function grad_stack!(stack::AbstractArray,  pe::Periodic, X::Matrix{Float64}, data::IsotropicData)
    d, nobsv = size(X)
    R = distance(pe, X, data)

    for i in 1:nobsv, j in 1:i
        @inbounds stack[i,j,1] = 4.0*pe.σ2*(sin(π*R[i,j]/pe.p)^2/pe.ℓ2)*exp(-2/pe.ℓ2*sin(π*R[i,j]/pe.p)^2)  # dK_dlogℓ
        @inbounds stack[j,i,1] = stack[i,j,1]

        @inbounds stack[i,j,2] = 2.0*pe.σ2*exp(-2/pe.ℓ2*sin(π*R[i,j]/pe.p)^2)        # dK_dlogσ
        @inbounds stack[j,i,2] = stack[i,j,2]

        @inbounds stack[i,j,3] = 4.0/pe.ℓ2*pe.σ2*(π*R[i,j]/pe.p)*sin(π*R[i,j]/pe.p)*cos(π*R[i,j]/pe.p)*exp(-2/pe.ℓ2*sin(π*R[i,j]/pe.p)^2)    # dK_dlogp
        @inbounds stack[j,i,3] = stack[i,j,3]
    end
    return stack
end


text:
include("test_filewriter.jl")
include("test_kind0.jl")
include("test_operators.jl")
include("test_show.jl")
include("test_state.jl")
include("test_dicke.jl")
include("test_collective.jl")