text:
export hessian_backback

@i function backback(f, args...; index::Int, iloss::Int, kwargs...)
    # forward
    Grad(f)(args...; kwargs..., iloss=iloss)

    for i = 1:length(args)
        GVar(args |> tget(i) |> grad)
        GVar(args |> tget(i) |> value)
    end
    (args |> tget(index) |> grad |> grad) += 1
    # backward#2
    (~Grad(f))(args...; kwargs..., iloss=iloss)
end

"""
    hessian_backback(f, args; iloss::Int, kwargs...)

Obtain the Hessian matrix of `f(args..., kwargs...)` by back propagating adjoint program.
"""
function hessian_backback(f, args; iloss::Int, kwargs...)
    N = length(args)
    hmat = zeros(N, N)
    for i=1:N
        if !(args[i] isa Integer || args[i] isa AbstractVector)
            res = backback(f, args...; kwargs..., index=i, iloss=iloss)
            hmat[:,i] .= map(x->grad(value(x)), res[2:end])
        end
    end
    hmat
end

function hessian_numeric(f, args; iloss::Int, η=1e-5, kwargs...)
    narg = length(args)
    res = zeros(narg, narg)
    largs = [args...]
    for i = 1:narg
        if nparams(args[i]) == 1
            @instr (largs[i] |> value) += η/2
            gpos = gradient(f, (largs...,); iloss=iloss, kwargs...)
            @instr (largs[i] |> value) -= η
            gneg = gradient(f, (largs...,); iloss=iloss, kwargs...)
            @instr (largs[i] |> value) += η/2
            res[:,i] .= (gpos .- gneg)./η
        end
    end
    return res
end

function local_hessian_numeric(f, args; kwargs...)
    nargs = length(args)
    hes = zeros(nargs,nargs,nargs)
    for j=1:nargs
        if nparams(args[j]) == 1
            hes[:,:,j] .= hessian_numeric(f, args; kwargs..., iloss=j)
        end
    end
    mask = BitArray(nparams.(args) .== 1)
    hes[mask, mask, mask]
end


text:
abstract type Healing end
abstract type NoHealing <: Healing end
abstract type AnnealingOnly <: Healing end
abstract type IsotropicDomainOnly <: Healing end
abstract type AnnealingIsotropicDomain <: Healing end
abstract type FabricDestructionOnly <: Healing end
abstract type AnnealingFabricDestruction <: Healing end

# Define isotropic domain
struct IsotropicDomain{T}
    r::T  # depth above which Ω is isotropic
end

# Rate of annealing
struct Annealing{T}
    rate::T # rate of annealing
end

# Destruction of fabric
struct FabricDestruction{T}
    ϵ::T # ratio a1/a2 at which aggregate goes back to its anisotropic form
end

struct FiniteStrain{T} <: Healing
    fse::Matrix{FiniteStrainEllipsoid{Float64}} # matrix containing the FSE at every ip
    annealing::Annealing # rate of annealing
    destruction::FabricDestruction # ratio a1/a2 at which aggregate goes back to its anisotropic form
    isotropic_domain::IsotropicDomain # depth above which Ω is isotropic

    function FiniteStrain(nel; nip = 6, ϵ = 1e3, annealing_rate = 0, r_iso = 0)
        # instantiate isotropic finite strain ellipsoid
        fse = [FiniteStrainEllipsoid(1.0, 0.0, 0.0, 1.0, 1.0, 1.0) for _ in 1:nel, _ in 1:nip]
        # define active healing mechanism(s)
        if (annealing_rate != 0) && (ϵ != 0)
            type = AnnealingFabricDestruction

        elseif (annealing_rate != 0) && (ϵ == 0)
            type = AnnealingOnly

        elseif (annealing_rate == 0) && (ϵ != 0)
            type = FabricDestructionOnly

        else
            type = NoHealing

        end
        # create object
        new{type}(
            fse,
            Annealing(annealing_rate),
            IsotropicDomain(r_iso)
        )

    end
end

*(a::Number, annealing::Annealing) = a*annealing.rate
*(annealing::Annealing, a::Number) = a*annealing.rate

#=
Unstrecth the FSE semi-axes by a factor of s.
The area of the ellipse must remain the same
after unstretching
=#
function unstretch_axes(a1, a2, s)
    # a2, a1 are the old semiaxes
    # s is the unstretching factor
    r0 = a2/a1 # old ratio
    r = @muladd r0 + (1-r0)*s # new aspect ratio
    A = π*a2*a1 # Area of the FSE
    # find new unstrecthed axes such that the
    # aspect ratio is r and A remains untouched
    a1u = √(A/(π*r))
    a2u = r*a1u

    return a1u, a2u
end

#=
Recover the new F tensor after unstretching the FSE semi-axes.
We know that F = LR, where L is the left-stretch tensor and R⁻¹=Rᵀ
is the orthogonal rotation tensor. First we recover the new L using
the old eigenvalues (we do not want to rotate the FSE, only unstretch
it) and new eigenvalues:
    L = PDP⁻¹
where P = [λ₁ λ₃ λ₃] and D = I⋅[a₁ a₃ a₃]ᵀ. R is obtained using the
old F and L0 tensors:
    R = F\L0
and finally we obtain the new deformation gradient tensor:
    Fn = LR
=#
function recover_F(F, L0, λ11, λ12, λ21, λ22, a1, a2)
    P = @SMatrix [
        λ11 λ21
        λ12 λ22
    ]
    D = @SMatrix [
        a1*a1  0.0
        0.0    a2*a2
    ]
    # recover left-stretch tensor
    L = P*D\P
    # rotation tensor R = F\(F*Fᵀ) = F\l0
    R = F\L0
    # unstretched tensor of deformation
    Fu = L*R
end

ext: jl

max_stars_repo_name: albert-de-montserrat/Persephone

max_issues_repo_name: albert-de-montserrat/Persephone

max_forks_repo_name: albert-de-montserrat/Persephone


text:
push!(LOAD_PATH,"../")

using Documenter
using Downhill

makedocs(
    ;
    sitename="Downhill.jl documentation",
    modules=[Downhill],
    pages = [
        "index.md",
        "Optimization Methods" => "core_types.md",
        "Basic Functions" => "functions.md",
        "Customization" => "wrappers.md"
    ],
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    )
)

deploydocs(
    repo = "github.com/vvpisarev/Downhill.jl.git",
)


text:
using Pkg
packages = keys(Pkg.installed())
if !in("Vec", packages)
    Pkg.add(PackageSpec(url="https://github.com/sisl/Vec.jl.git"))
end
if !in("Records", packages)
    Pkg.add(PackageSpec(url="https://github.com/sisl/Records.jl.git"))
end
if !in("AutomotiveDrivingModels", packages)
    Pkg.add(PackageSpec(url="https://github.com/sisl/AutomotiveDrivingModels.jl.git"))
end
if !in("AutoViz", packages)
    Pkg.add(PackageSpec(url="https://github.com/sisl/AutoViz.jl.git"))
end


text:
abstract type PDMPSampler
end

"""
    ContinuousDynamics

Abstract type for the deterministic dynamics of PDMPs
"""
abstract type ContinuousDynamics end

"""
    struct ZigZag(Γ, μ) <: ContinuousDynamics

Local ZigZag sampler which exploits
any independence structure of the target measure,
in form the argument `Γ`, a sparse precision matrix approximating
target precision. `μ` is the approximate target mean.
"""
struct ZigZag{T,S,S2,R} <: ContinuousDynamics
    Γ::T
    μ::S
    σ::S2
    λref::R
    ρ::R
    ρ̄::R
end
ZigZag(Γ, μ, σ=(Vector(diag(Γ))).^(-0.5); λref=0.0, ρ=0.0) = ZigZag(Γ, μ, σ, λref, ρ, sqrt(1-ρ^2))

"""
    BouncyParticle(λ) <: ContinuousDynamics
Input: argument `Γ`, a sparse precision matrix approximating target precision.
Bouncy particle sampler,  `λ` is the refreshment rate, which has to be
strictly positive.
"""
struct BouncyParticle{T, S, R, V, LT} <: ContinuousDynamics
    Γ::T
    μ::S
    λref::R
    ρ::R
    U::V
    L::LT
end
BouncyParticle(Γ, μ, λ; ρ=0.0) = BouncyParticle(Γ, μ, λ, ρ, nothing, (cholesky(Symmetric(Γ)).L))
# simple constructor for first experiments
BouncyParticle(λ, d) = BouncyParticle(1.0I(d), zeros(d), λ, 0.0, nothing)
"""
    Boomerang(μ, λ) <: ContinuousDynamics

Dynamics preserving the `N(μ, Σ)` measure (Boomerang)
with refreshment time `λ`
"""
struct Boomerang{U, T, S, LT} <: ContinuousDynamics
    Γ::U
    μ::T
    λref::S
    ρ::S
    L::LT
end
Boomerang(Γ, μ, λ; ρ=0.0) = Boomerang(Γ, μ, λ, ρ, (cholesky(Symmetric(Γ)).L))
"""
    FactBoomerang(Γ, μ, λ) <: ContinuousDynamics

Factorized Boomerang dynamics preserving the `N(μ, inv(Diagonal(Γ)))` measure
with refreshment time `λ`.

Exploits the conditional independence structure of the target measure,
in form the argument `Γ`, a sparse precision matrix approximating
target precision. `μ` is the approximate target mean.
"""
struct FactBoomerang{R, T, S} <: ContinuousDynamics
    Γ::R
    μ::T
    σ::T
    λref::S
    ρ::S
    ρ̄::S
end
FactBoomerang(Γ, μ, λ, σ=(Vector(diag(Γ))).^(-0.5); ρ=0.0) = FactBoomerang(Γ, μ, σ, λ, ρ, sqrt(1-ρ^2))

"""
    ZigZag1d <: ContinuousDynamics

1-d toy ZigZag sampler, dynamics preserving the Lebesgue measure.
"""
struct ZigZag1d <: ContinuousDynamics  end

"""
    Boomerang1d(Σ, μ, λ) <: ContinuousDynamics

1-d toy boomerang samper. Dynamics preserving the `N(μ, Σ)` measure
with refreshment time `λ`.
"""
struct Boomerang1d{S,T} <: ContinuousDynamics
    Σ::S
    μ::T
    λref::T
end
Boomerang1d(λ) = Boomerang1d(1.0, 0.0, λ)
Boomerang1d(μstar, λ) = Boomerang1d(1.0, μstar, λ)

"""
    ExtendedForm()

Indicates as `args[1]` that `∇ϕ`
depends on the extended arguments

    ∇ϕ(t, x, θ, i, t′, Z, args...)

instead of

    ∇ϕ(x, i, args...)

Can be used to implement `∇ϕ` depending on random coefficients.
"""
struct ExtendedForm
end

abstract type Bound end
struct LocalBound{T} <: Bound
    c::T
end
struct GlobalBound{T} <: Bound
    c::T
end

Base.:*(C::T, a) where {T <: Bound} = T(C.c*a)


text:
struct SymmetricSinkhornSolver{A<:Sinkhorn,M,CT,E<:Real,T<:Real,R<:Real,C1,C2}
    source::M
    C::CT
    eps::E
    alg::A
    atol::T
    rtol::R
    maxiter::Int
    check_convergence::Int
    cache::C1
    convergence_cache::C2
end

struct SymmetricSinkhornGibbs <: Sinkhorn end

Base.show(io::IO, ::SymmetricSinkhornGibbs) = print(io, "Symmetric Sinkhorn algorithm")

struct SymmetricSinkhornGibbsCache{U,KT}
    u::U
    K::KT
    Kv::U
end

function build_cache(
    ::Type{T},
    ::SymmetricSinkhornGibbs,
    size2::Tuple,
    μ::AbstractVecOrMat,
    C::AbstractMatrix,
    ε::Real,
) where {T}
    K = similar(C, T)
    @. K = exp(-C / ε)

    u = similar(μ, T, size(μ, 1), size2...)
    fill!(u, one(T))

    Kv = similar(u)

    return SymmetricSinkhornGibbsCache(u, K, Kv)
end

function check_convergence(solver::SymmetricSinkhornSolver)
    A_batched_mul_B!(solver.cache.Kv, solver.cache.K, solver.cache.u)
    return OptimalTransport.check_convergence(
        solver.source,
        solver.cache.u,
        solver.cache.Kv,
        solver.convergence_cache,
        solver.atol,
        solver.rtol,
    )
end

function build_solver(
    μ::AbstractVecOrMat,
    C::AbstractMatrix,
    ε::Real,
    alg::SymmetricSinkhornGibbs;
    atol=nothing,
    rtol=nothing,
    check_convergence=1,
    maxiter::Int=25,
)
    size2 = size(μ)[2:end]

    # compute type
    T = float(Base.promote_eltype(μ, one(eltype(C)) / ε))

    # build caches using SinkhornGibbsCache struct (since there is no dependence on ν)
    cache = build_cache(T, alg, size2, μ, C, ε)
    convergence_cache = build_convergence_cache(T, size2, μ)

    # set tolerances
    _atol = atol === nothing ? 0 : atol
    _rtol = rtol === nothing ? (_atol > zero(_atol) ? zero(T) : sqrt(eps(T))) : rtol

    # create solver
    solver = SymmetricSinkhornSolver(
        μ, C, ε, alg, _atol, _rtol, maxiter, check_convergence, cache, convergence_cache
    )
    return solver
end

function init_step!(solver::SymmetricSinkhornSolver{SymmetricSinkhornGibbs})
    source = solver.source
    cache = solver.cache
    u = cache.u
    K = cache.K
    Kv = cache.Kv
    return A_batched_mul_B!(Kv, K, u)
end

prestep!(::SymmetricSinkhornSolver{SymmetricSinkhornGibbs}, ::Int) = nothing

function step!(solver::SymmetricSinkhornSolver{SymmetricSinkhornGibbs}, iter::Int)
    source = solver.source
    cache = solver.cache
    u = cache.u
    K = cache.K
    Kv = cache.Kv
    @. u = sqrt(source * u / Kv)
    return A_batched_mul_B!(Kv, K, u)
end

function sinkhorn_plan(solver::SymmetricSinkhornSolver{SymmetricSinkhornGibbs})
    cache = solver.cache
    return sinkhorn_plan(cache.u, cache.u, cache.K)
end

function sinkhorn(μ, C, ε, alg::SymmetricSinkhornGibbs; kwargs...)
    # build solver
    solver = build_solver(μ, C, ε, alg; kwargs...)

    # perform Sinkhorn algorithm
    solve!(solver)

    # compute optimal transport plan
    γ = sinkhorn_plan(solver)

    return γ
end

function sinkhorn2(
    μ, C, ε, alg::SymmetricSinkhornGibbs; regularization=false, plan=nothing, kwargs...
)
    cost = if regularization && plan === nothing
        # special case where we can take advantage of dual objective formula
        # build solver
        solver = build_solver(μ, C, ε, alg; kwargs...)
        # perform Sinkhorn algorithm
        solve!(solver)
        # return loss
        cache = solver.cache
        sinkhorn_dual_objective(cache.u, cache.u, cache.Kv, cache.K, solver.eps)
    else
        γ = if plan === nothing
            sinkhorn(μ, C, ε, alg; kwargs...)
        else
            size(plan) == size(C) || error(
                "optimal transport plan `plan` and cost matrix `C` must be of the same size",
            )
            plan
        end
        sinkhorn_cost_from_plan(γ, C, ε; regularization=regularization)
    end
    return cost
end


text:
module JuliennedArrays

import Base: axes, getindex, setindex!, size
using Base: promote_op, @pure, @propagate_inbounds, tail

map_unrolled(call, variables::Tuple{}) = ()
map_unrolled(call, variables) =
    call(first(variables)), map_unrolled(call, tail(variables))...

map_unrolled(call, variables1::Tuple{}, variables2::Tuple{}) = ()
map_unrolled(call, variables1, variables2) =
    call(first(variables1), first(variables2)),
    map_unrolled(call, tail(variables1), tail(variables2))...

is_in(needle::Needle, straw1::Needle, straws...) where {Needle} = True()
is_in(needle, straw1, straws...) = is_in(needle, straws...)
is_in(needle) = False()

in_unrolled(straws, needle1, needles...) =
    is_in(needle1, straws...),
    in_unrolled(straws, needles...)...
in_unrolled(straws) = ()

@pure as_vals(them::Int...) = map(Val, them)

abstract type TypedBool end
"""
    struct True

Typed `true`
"""
struct True <: TypedBool end
"""
    struct False

Typed `false`
"""
struct False <: TypedBool end

@inline untyped(::True) = true
@inline untyped(::False) = false

not(::False) = True()
not(::True) = False()

export True
export False

getindex_unrolled(into::Tuple{}, switches::Tuple{}) = ()
function getindex_unrolled(into, switches)
    next = getindex_unrolled(tail(into), tail(switches))
    if untyped(first(switches))
        (first(into), next...)
    else
        next
    end
end

setindex_unrolled(old::Tuple{}, something, ::Tuple{}) = ()
setindex_unrolled(old, new, switches) =
    if untyped(first(switches))
        first(new),
        setindex_unrolled(tail(old), tail(new), tail(switches))...
    else
        first(old),
        setindex_unrolled(tail(old), new, tail(switches))...
    end

struct Slices{Item, Dimensions, Whole, Alongs} <: AbstractArray{Item, Dimensions}
    whole::Whole
    alongs::Alongs
end
Slices{Item, Dimensions}(whole::Whole, alongs::Alongs) where {Item, Dimensions, Whole, Alongs} =
    Slices{Item, Dimensions, Whole, Alongs}(whole, alongs)

axes(slices::Slices) =
    getindex_unrolled(axes(slices.whole), map_unrolled(not, slices.alongs))
size(slices::Slices) = map_unrolled(length, axes(slices))

slice_index(slices, indices) = setindex_unrolled(
    axes(slices.whole),
    indices,
    map_unrolled(not, slices.alongs)
)
@propagate_inbounds getindex(slices::Slices, indices::Int...) =
    view(slices.whole, slice_index(slices, indices)...)
@propagate_inbounds setindex!(slices::Slices, value, indices::Int...) =
    slices.whole[slice_index(slices, indices)...] = value

axis_or_1(switch, axis) =
    if untyped(switch)
        axis
    else
        1
    end
"""
    Slices(whole, alongs::TypedBool...)

Slice `whole` into `view`s.

`alongs`, made of [`True`](@ref) and [`False`](@ref) objects, shows which dimensions will be replaced with `:` when slicing.

```jldoctest
julia> using JuliennedArrays

julia> whole = [1 2; 3 4];

julia> slices = Slices(whole, False(), True())
2-element Slices{SubArray{Int64,1,Array{Int64,2},Tuple{Int64,Base.OneTo{Int64}},true},1,Array{Int64,2},Tuple{False,True}}:
 [1, 2]
 [3, 4]

julia> slices[1] == whole[1, :]
true

julia> slices[1] = [2, 1];

julia> whole
2×2 Array{Int64,2}:
 2  1
 3  4

julia> larger = rand(5, 5, 5);

julia> larger_slices = Slices(larger, True(), False(), False());

julia> size(first(larger_slices))
(5,)
```
"""
Slices(whole::AbstractArray, alongs::TypedBool...) =
    Slices{
        typeof(@inbounds view(
            whole,
            map_unrolled(axis_or_1, alongs, axes(whole))...
        )),
        length(getindex_unrolled(alongs, map_unrolled(not, alongs)))
    }(whole, alongs)

"""
    Slices(whole, alongs::Int...)

Alternative syntax: `alongs` is which dimensions will be replaced with `:` when slicing.

```jldoctest
julia> using JuliennedArrays

julia> input = reshape(1:8, 2, 2, 2)
2×2×2 reshape(::UnitRange{Int64}, 2, 2, 2) with eltype Int64:
[:, :, 1] =
 1  3
 2  4

[:, :, 2] =
 5  7
 6  8

julia> Slices(input, 1, 3)
2-element Slices{SubArray{Int64,2,Base.ReshapedArray{Int64,3,UnitRange{Int64},Tuple{}},Tuple{Base.OneTo{Int64},Int64,Base.OneTo{Int64}},false},1,Base.ReshapedArray{Int64,3,UnitRange{Int64},Tuple{}},Tuple{True,False,True}}:
 [1 5; 2 6]
 [3 7; 4 8]
```
"""
Slices(whole::AbstractArray{Item, NumberOfDimensions}, alongs::Int...) where {Item, NumberOfDimensions} =
    Slices(whole, in_unrolled(
        as_vals(alongs...),
        ntuple(Val, NumberOfDimensions)...
    )...)
export Slices

struct Align{Item, Dimensions, Sliced, Alongs} <: AbstractArray{Item, Dimensions}
    slices::Sliced
    alongs::Alongs
end
Align{Item, Dimensions}(slices::Sliced, alongs::Alongs) where {Item, Dimensions, Sliced, Alongs} =
    Align{Item, Dimensions, Sliced, Alongs}(slices, alongs)

axes(aligned::Align) = setindex_unrolled(
    setindex_unrolled(
        aligned.alongs,
        axes(aligned.slices),
        map_unrolled(not, aligned.alongs)
    ),
    axes(first(aligned.slices)),
    aligned.alongs
)
size(aligned::Align) = map_unrolled(length, axes(aligned))

split_indices(aligned, indices) =
    getindex_unrolled(indices, map_unrolled(not, aligned.alongs)),
    getindex_unrolled(indices, aligned.alongs)
@propagate_inbounds function getindex(aligned::Align, indices::Int...)
    outer, inner = split_indices(aligned, indices)
    aligned.slices[outer...][inner...]
end
@propagate_inbounds function setindex!(aligned::Align, value, indices::Int...)
    outer, inner = split_indices(aligned, indices)
    aligned.slices[outer...][inner...] = value
end

"""
    Align(slices, alongs::TypedBool...)

`Align` an array of arrays, all with the same size.

`alongs`, made of [`True`](@ref) and [`False`](@ref) objects, shows which dimensions will be taken up by the inner arrays. Inverse of [`Slices`](@ref).

```jldoctest
julia> using JuliennedArrays

julia> slices = [[1, 2], [3, 4]];

julia> aligned = Align(slices, False(), True())
2×2 Align{Int64,2,Array{Array{Int64,1},1},Tuple{False,True}}:
 1  2
 3  4

julia> aligned[1, :] == slices[1]
true

julia> aligned[1, 1] = 0;

julia> slices
2-element Array{Array{Int64,1},1}:
 [0, 2]
 [3, 4]
```
"""
Align(slices::AbstractArray{<:AbstractArray{Item, InnerDimensions}, OuterDimensions}, alongs::TypedBool...) where {Item, InnerDimensions, OuterDimensions} =
    Align{Item, OuterDimensions + InnerDimensions}(slices, alongs)
export Align

"""
    Along(slices, alongs::Int...)

Alternative syntax: `alongs` is which dimensions will be taken up by the inner arrays.

```jldoctest
julia> using JuliennedArrays

julia> input = reshape(1:8, 2, 2, 2)
2×2×2 reshape(::UnitRange{Int64}, 2, 2, 2) with eltype Int64:
[:, :, 1] =
 1  3
 2  4

[:, :, 2] =
 5  7
 6  8

julia> slices = collect(Slices(input, 1, 3))
2-element Array{SubArray{Int64,2,Base.ReshapedArray{Int64,3,UnitRange{Int64},Tuple{}},Tuple{Base.OneTo{Int64},Int64,Base.OneTo{Int64}},false},1}:
 [1 5; 2 6]
 [3 7; 4 8]

julia> Align(slices, 1, 3)
2×2×2 Align{Int64,3,Array{SubArray{Int64,2,Base.ReshapedArray{Int64,3,UnitRange{Int64},Tuple{}},Tuple{Base.OneTo{Int64},Int64,Base.OneTo{Int64}},false},1},Tuple{True,False,True}}:
[:, :, 1] =
 1  3
 2  4

[:, :, 2] =
 5  7
 6  8
```
"""
Align(slices::AbstractArray{<:AbstractArray{Item, InnerDimensions}, OuterDimensions}, alongs::Int...) where {Item, InnerDimensions, OuterDimensions} =
    Align(slices, in_unrolled(
        as_vals(alongs...),
        ntuple(Val, InnerDimensions + OuterDimensions)...
    )...)

end


text:
using Gtk, Graphics, Colors

function drawline(ctx, p1, p2, color, width)
    move_to(ctx, p1.x, p1.y)
    set_source(ctx, color)
    line_to(ctx, p2.x, p2.y)
    set_line_width(ctx, width)
    stroke(ctx)
end

const can = @GtkCanvas()
const win = GtkWindow(can, "Colour pinstripe/Display", 400, 400)
const colors = [colorant"black", colorant"red", colorant"green", colorant"blue",
          colorant"magenta", colorant"cyan", colorant"yellow", colorant"white"]
const numcolors = length(colors)

@guarded draw(can) do widget
    ctx = getgc(can)
    h = height(can)
    w = width(can)
    deltaw = 1.0
    for (i, x) in enumerate(0:deltaw:w)
        drawline(ctx, Point(x, 0.25*h), Point(x, 0), colors[i % numcolors + 1], deltaw)
    end
    for (i, x) in enumerate(0:deltaw*2:w)
        drawline(ctx, Point(x, 0.5*h), Point(x, 0.25*h), colors[i % numcolors + 1], deltaw*2)
    end
    for (i, x) in enumerate(0:deltaw*3:w)
        drawline(ctx, Point(x, 0.75*h), Point(x, 0.5*h), colors[i % numcolors + 1], deltaw*3)
    end
    for (i, x) in enumerate(0:deltaw*4:w)
        drawline(ctx, Point(x, h), Point(x, 0.75*h), colors[i % numcolors + 1], deltaw*4)
    end
end

show(can)
const cond = Condition()
endit(w) = notify(cond)
signal_connect(endit, win, :destroy)
wait(cond)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

using ..OpticSim, .Geometry
using LinearAlgebra
using Distributions
using StaticArrays

import Makie

using .Emitters
using .Emitters.Spectrum
using .Emitters.Directions
using .Emitters.Origins
using .Emitters.AngularPower
using .Emitters.Sources

const ARRROW_LENGTH = 0.5
const ARRROW_SIZE = 0.01
const MARKER_SIZE = 1

#-------------------------------------
# draw debug information - local axes and positions
#-------------------------------------
function maybe_draw_debug_info(scene::Makie.LScene, o::Origins.AbstractOriginDistribution; transform::Geometry.Transform = Transform(), debug::Bool=false, kwargs...) where {T<:Real}

    dir = forward(transform)
    uv = SVector{3}(right(transform))
    vv = SVector{3}(up(transform))
    pos = origin(transform)

    if (debug)
        # this is a stupid hack to force makie to render in 3d - for some scenes, makie decide with no apperent reason to show in 2d instead of 3d
        Makie.scatter!(scene, [pos[1], pos[1]+0.1], [pos[2], pos[2]+0.1], [pos[3], pos[3]+0.1], color=:red, markersize=0)

        # draw the origin and normal of the surface
        Makie.scatter!(scene, pos, color=:blue, markersize = MARKER_SIZE * visual_size(o))

        # normal
        arrow_size = ARRROW_SIZE * visual_size(o)
        arrow_start = pos
        arrow_end = dir * ARRROW_LENGTH * visual_size(o)
        Makie.arrows!(scene.scene, [Makie.Point3f0(arrow_start)], [Makie.Point3f0(arrow_end)], arrowsize=arrow_size, linewidth=arrow_size * 0.5, linecolor=:blue, arrowcolor=:blue)
        arrow_end = uv * 0.5 * ARRROW_LENGTH * visual_size(o)
        Makie.arrows!(scene.scene, [Makie.Point3f0(arrow_start)], [Makie.Point3f0(arrow_end)], arrowsize= 0.5 * arrow_size, linewidth=arrow_size * 0.5, linecolor=:red, arrowcolor=:red)
        arrow_end = vv * 0.5 * ARRROW_LENGTH * visual_size(o)
        Makie.arrows!(scene.scene, [Makie.Point3f0(arrow_start)], [Makie.Point3f0(arrow_end)], arrowsize= 0.5 * arrow_size, linewidth=arrow_size * 0.5, linecolor=:green, arrowcolor=:green)

        # draw all the samples origins
        positions = map(x -> transform*x, collect(o))
        positions = collect(Makie.Point3f0, positions)
        Makie.scatter!(scene, positions, color=:green, markersize = MARKER_SIZE * visual_size(o))

        # positions = collect(Makie.Point3f0, o)
        # Makie.scatter!(scene, positions, color=:green, markersize = MARKER_SIZE * visual_size(o))
    end

end

#-------------------------------------
# draw point origin
#-------------------------------------
function OpticSim.Vis.draw!(scene::Makie.LScene, o::Origins.Point; transform::Geometry.Transform = Transform(), kwargs...) where {T<:Real}
        maybe_draw_debug_info(scene, o; transform=transform, kwargs...)
end

#-------------------------------------
# draw RectGrid and RectUniform origins
#-------------------------------------
function OpticSim.Vis.draw!(scene::Makie.LScene, o::Union{Origins.RectGrid, Origins.RectUniform}; transform::Geometry.Transform = Transform(), kwargs...) where {T<:Real}
    dir = forward(transform)
    uv = SVector{3}(right(transform))
    vv = SVector{3}(up(transform))
    pos = origin(transform)

    # @info "RECT: transform $(pos)"

    plane = OpticSim.Plane(dir, pos)
    rect = OpticSim.Rectangle(plane, o.width / 2, o.height / 2, uv, vv)

    OpticSim.Vis.draw!(scene, rect;  kwargs...)

    maybe_draw_debug_info(scene, o; transform=transform, kwargs...)
end

#-------------------------------------
# draw hexapolar origin
#-------------------------------------
function OpticSim.Vis.draw!(scene::Makie.LScene, o::Origins.Hexapolar; transform::Geometry.Transform = Transform(), kwargs...) where {T<:Real}
    dir = forward(transform)
    uv = SVector{3}(right(transform))
    vv = SVector{3}(up(transform))
    pos = origin(transform)

    plane = OpticSim.Plane(dir, pos)
    ellipse = OpticSim.Ellipse(plane, o.halfsizeu, o.halfsizev, uv, vv)

    OpticSim.Vis.draw!(scene, ellipse;  kwargs...)

    maybe_draw_debug_info(scene, o; transform=transform, kwargs...)
end

#-------------------------------------
# draw source
#-------------------------------------
function OpticSim.Vis.draw!(scene::Makie.LScene, s::Sources.Source{T}; parent_transform::Geometry.Transform = Transform(), debug::Bool=false, kwargs...) where {T<:Real}

    OpticSim.Vis.draw!(scene, s.origins;  transform=parent_transform * s.transform, debug=debug, kwargs...)

    if (debug)
        m = zeros(T, length(s), 7)
        for (index, optical_ray) in enumerate(s)
            ray = OpticSim.ray(optical_ray)
            ray = parent_transform * ray
            m[index, 1:7] = [ray.origin... ray.direction... OpticSim.power(optical_ray)]
        end

        m[:, 4:6] .*= m[:, 7] * ARRROW_LENGTH * visual_size(s.origins)

        # Makie.arrows!(scene, [Makie.Point3f0(origin(ray))], [Makie.Point3f0(rayscale * direction(ray))]; kwargs..., arrowsize = min(0.05, rayscale * 0.05), arrowcolor = color, linecolor = color, linewidth = 2)
        color = :yellow
        arrow_size = ARRROW_SIZE * visual_size(s.origins)
        Makie.arrows!(scene, m[:,1], m[:,2], m[:,3], m[:,4], m[:,5], m[:,6]; kwargs...,  arrowcolor=color, linecolor=color, arrowsize=arrow_size, linewidth=arrow_size*0.5)
    end

    # for ray in o
    #     OpticSim.Vis.draw!(scene, ray)
    # end
end

#-------------------------------------
# draw optical rays
#-------------------------------------
function OpticSim.Vis.draw!(scene::Makie.LScene, rays::AbstractVector{OpticSim.OpticalRay{T, 3}}; kwargs...) where {T<:Real}
    m = zeros(T, length(rays)*2, 3)
    for (index, optical_ray) in enumerate(rays)
        ray = OpticSim.ray(optical_ray)
        m[(index-1)*2+1, 1:3] = [origin(optical_ray)...]
        m[(index-1)*2+2, 1:3] = [(OpticSim.origin(optical_ray) + OpticSim.direction(optical_ray) * 1 * OpticSim.power(optical_ray))... ]
    end

    color = :green
    Makie.linesegments!(scene, m[:,1], m[:,2], m[:,3]; kwargs...,  color = color, linewidth = 2, )
end

#-------------------------------------
# draw composite source
#-------------------------------------
function OpticSim.Vis.draw!(scene::Makie.LScene, s::Sources.CompositeSource{T}; parent_transform::Geometry.Transform = Transform(), kwargs...) where {T<:Real}
    for source in s.sources
        OpticSim.Vis.draw!(scene, source; parent_transform=parent_transform*s.transform, kwargs...)
    end
end


text:
#=
The number 3797 has an interesting property. Being prime itself, it is possible to continuously remove digits from left to right, and remain prime at each stage: 3797, 797, 97, and 7. Similarly we can work from right to left: 3797, 379, 37, and 3.

Find the sum of the only eleven primes that are both truncatable from left to right and right to left.

NOTE: 2, 3, 5, and 7 are not considered to be truncatable primes.
=#
include("projecteulerutils.jl")
using Primes

function istruncatableprime(a)
  da = reverse(digits(a))
  f1,f2 = copy(da), copy(da)
  for i in 1:length(da)-1
    deleteat!(f1,1)
    pop!(f2)
    if !isprime(nd(f1)) || !isprime(nd(f2))
      return false
    end
  end
  true
end

function calc()
  i,count,total = 10,0,0
  while count < 11
    if isprime(i) && istruncatableprime(i)
      count += 1
      total += i
    end
    i += 1
  end
  total
end
@time println(calc())

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
@doc raw"""
    init_aux_variables([rng::AbstractRNG], ::Likelihood, n::Int) -> TupleVector

Initialize collections of `n` auxiliary variables in  a `TupleVector` to be used
in the context of sampling.
`n` should be the number of data inputs.

See also [`init_aux_posterior`](@ref) for variational inference.
"""
init_aux_variables

@doc raw"""
    init_aux_posterior([T::DataType], lik::Likelihood, n::Int) -> AbstractProductMeasure

Initialize collections of `n` (independent) posteriors for the auxiliary
variables in the context of variational inference.
`n` should be the size of the data used every iteration.
The real variational parameters will be given type `T` (`Float64` by default)

See also [`init_aux_variables`](@ref) for sampling.
"""
init_aux_posterior

@doc raw"""
    aux_full_conditional(lik::AbstractLikelihood, yᵢ, fᵢ) -> AbstractNTDist

Return the full conditional distribution of `Ωᵢ`.
"""
aux_full_conditional

@doc raw"""
    aux_sample!([rng::AbstractRNG], Ω, lik::Likelihood, y, f) -> TupleVector

Sample the auxiliary variables `Ω` in-place based on the full-conditional [`aux_full_conditional`](@ref)
associated with the augmented likelihood:
```math
    p(\Omega|y,f) \propto p(\Omega,y|f).
```

See also [`aux_sample`](@ref) for an allocating version and [`aux_posterior!`](@ref)
for variational inference.
"""
aux_sample!

@doc raw"""
    aux_sample([rng::AbstractRNG], lik::Likelihood, y, f) -> TupleVector

Sample and allocate the auxiliary variables `Ω` in a `TupleVector` based
on the full-conditional associated with the likelihood.

See [`aux_sample!`](@ref) for an in-place version.
"""
aux_sample

@doc raw"""
    aux_full_conditional(lik::Likelihood, y, f::Real) -> Distribution

Given the observation `y` and latent `f`, returns the full conditional on the
auxiliary variables `Ω`.
"""

@doc raw"""
    aux_posterior!(qΩ, lik::Likelihood, y, qf) -> AbstractProductMeasure

Compute the optimal posterior of the auxiliary variables ``q^*(\Omega)`` given
the marginal distributions `qf` by updating the variational parameters
in-place using the formula
```math
    q^*(\Omega) \propto \exp\left(E_{q(f)}\left[p(\Omega|f,y)\right]\right)
```

See also [`aux_posterior`](@ref) and [`aux_sample!`](@ref)
"""
aux_posterior!

@doc raw"""
    aux_posterior(lik::Likelihood, y, qf) -> AbstractProductMeasure

Compute the optimal posterior of the auxiliary variables in a new
`AbstractProductMeasure` (`For` by default).

See [`aux_posterior!`](@ref) for more details
"""
aux_posterior

@doc raw"""
    auglik_potential(lik::Likelihood, Ω, y) -> Tuple

Given the augmented likelihood ``l(\Omega,y,f) \propto \exp(\beta(\Omega,y) f + \gamma(\Omega,y)f^2)``,
return the potential, ``\beta(\Omega,y)``.
Note that this equivalent to the shift of the first natural parameter ``\eta_1 = \Sigma^{-1}\mu``.
The `Tuple` contains a `Vector` for each latent.

See also [`expected_auglik_potential`](@ref) for variational inference.
"""
auglik_potential

@doc raw"""
    auglik_precision(lik::Likelihood, Ω, y) -> Tuple

Given the augmented likelihood ``l(\Omega,y,f) \propto \exp(\beta(\Omega,y) f + \frac{\gamma(\Omega,y)}{2}f^2)``,
return the precision, $$\gamma(\Omega,y)$$, note that this equivalent to the
shift of the precision ``\Lambda = \Sigma^{-1}``.
The `Tuple` contains a `Vector` for each latent.

See also [`expected_auglik_precision`](@ref) for variational inference.
"""
auglik_precision

@doc raw"""
    auglik_potential_and_precision(lik::Likelihood, Ω, y) -> Tuple{Tuple, Tuple}

Returns both [`auglik_potential`](@ref) and [`auglik_precision`](@ref) when some
computation can be saved by doing both at the same time.
"""
auglik_potential_and_precision

@doc raw"""
    expected_auglik_potential(lik::Likelihood, qΩ, y) -> Tuple

Given the augmented likelihood ``l(\Omega,y,f) \propto \exp(\beta(\Omega,y) f + \frac{\gamma(\Omega,y)}{2}f^2)``,
return the expected potential, ``E_{q(\Omega)}[\beta(\Omega,y)]``, note that this equivalent to the
shift of the first variational natural parameter ``\eta_1 = \Sigma^{-1}\mu``.
The `Tuple` contains a `Vector` for each latent.

See also [`auglik_potential`](@ref) for sampling.
"""
expected_auglik_potential

@doc raw"""
    expected_auglik_precision(lik::Likelihood, qΩ, y) -> Tuple

Given the augmented likelihood $$l(\Omega,y,f) \propto \exp(\beta(\Omega,y) f + \frac{\gamma(\Omega,y)}{2}f^2)$$,
return the expected precision, $$E_{q(\Omega)}[\gamma(\Omega,y)]$$, note that this equivalent to the
shift of the variational precision $\Lambda = \Sigma^{-1}$.
The `Tuple` contains a `Vector` for each latent.

See also [`auglik_precision`](@ref) for sampling.
"""
expected_auglik_precision

@doc raw"""
    expected_auglik_potential_and_precision(lik::Likelihood, Ω, y) -> Tuple{Tuple, Tuple}

Returns both [`expected_auglik_potential`](@ref) and [`expected_auglik_precision`](@ref) when some
computation can be saved by doing both at the same time.
"""
expected_auglik_potential_and_precision

@doc raw"""
    aug_loglik(lik::Likelihood, Ω, y, f) -> Real

Return the augmented log-likelihood with the given parameters.
The augmented log-likelihood is of the form
```math
    \log p(y,\Omega|f) = \log l(y,\Omega,f) + \log p(\Omega|y).
```
To only obtain the $$p(\Omega|y)$$ part see [`aux_prior`](@ref) and see [`logtilt`](@ref)
for $$\log l(y, \Omega, f)$$.

A generic fallback exists based on [`logtilt`](@ref) and [`aux_prior`](@ref) but
specialized implementations are encouraged.
"""
aug_loglik

@doc raw"""
    aux_kldivergence(lik::Likelihood, qΩ::For, pΩ::For) -> Real
    aux_kldivergence(lik::Likelihood, qΩ::For, y) -> Real

Compute the analytical KL divergence between the auxiliary variables posterior
``q(\Omega)``, obtained with [`aux_posterior`](@ref) and prior
``p(\Omega)``, obtained with [`aux_prior`](@ref).
"""
aux_kldivergence

@doc raw"""
    aux_prior(lik::Likelihood, y) -> AbstractProductMeasure

Returns a `NamedTuple` of distributions with the same structure as [`aux_posterior`](@ref),
[`init_aux_posterior`](@ref) and [`init_aux_variables`](@ref).
"""
aux_prior

@doc raw"""
    logtilt(lik::Likelihood, Ω, y, f) -> Real

Compute the quadratic part on ``f`` of the augmented likelihood:
```math
    \log C(\Omega, y) + \alpha(\Omega,y) + \beta(\Omega,y)f + \frac{\gamma(\Omega,y)}{2}f^2.
```

See also [`expected_logtilt`](@ref) for variational inference.
"""
logtilt

@doc raw"""
    expected_logtilt(lik::Likelihood, qΩ, y, qf) -> Real

Compute the expectation of the quadratic part on $$f$$ of the augmented likelihood.
```math
    E_{q(\Omega,f)}\left[\log C(\Omega, y) + \alpha(\Omega,y) + \beta(\Omega,y)f + \frac{\gamma(\Omega,y)}{2}f^2\right].
```

See also [`logtilt`](@ref).
"""
function expected_logtilt(lik::AbstractLikelihood, qΩ, y, qf::AbstractVector{<:Normal})
    return mapreduce(+, y, qf, @ignore_derivatives marginals(qΩ)) do yᵢ, qfᵢ, qΩᵢ
        expected_logtilt(lik, qΩᵢ, yᵢ, qfᵢ)
    end
end


text:
@testset "Test DelayPECE method" begin
    function ϕ(x)
        if x == 0
            return 19.00001
        else
            return 19.0
        end
    end

    function f(t, y, ϕ)
        return 3.5*y*(1-ϕ/19)
    end

    h = 0.5
    α = 0.97
    τ = 0.8
    T = 1
    fddeprob = FDDEProblem(f, ϕ, α, τ, T)
    V, y = solve(fddeprob, h, DelayPECE())

    @test V≈[19.0, 19.0, 1.0]
    @test y≈[19.00001, 19.00001, 37.274176448220274]
end

@testset "Test DelayPI method" begin
    function ϕ(x)
        if x == 0
            return 19.00001
        else
            return 19.0
        end
    end
    function f(t, y, ϕ)
        return 3.5*y*(1-ϕ/19)
    end
    prob = FDDEProblem(f, ϕ, 0.97, 0.8, (0, 2))
    result = solve(prob, 0.5, DelayPI())
    @test result≈[19.00001, 19.00001, 19.00001, 18.99999190949352, 18.99997456359874]
end

@testset "Test DelayABM method" begin
    h=0.5; T=5; α=0.97; τ=2; q=0.5
    delayabmfun(t, ϕ, y) = 2*ϕ/(1+ϕ^9.65)-y
    prob = FDDEProblem(delayabmfun, q, α, τ, (0, T))
    x, y=solve(prob, h, DelayABM())

    @test isapprox(x, [1.078559863692747, 1.175963999045738, 1.1661317460354588, 1.128481756921719, 1.0016061526083417, 0.7724564325042358, 0.5974978685646778]; atol=1e-3)
    @test isapprox(y, [0.8889787467894421, 0.9404487875504524, 0.9667449499617093, 0.9803311436135411, 1.078559863692747, 1.175963999045738, 1.1661317460354588]; atol=1e-3)
end

@testset "Test DelayABM for FDDESystem" begin
    function EnzymeKinetics!(dy, y, ϕ, t)
        dy[1] = 10.5-y[1]/(1+0.0005*ϕ[4]^3)
        dy[2] = y[1]/(1+0.0005*ϕ[4]^3)-y[2]
        dy[3] = y[2]-y[3]
        dy[4] = y[3]-0.5*y[4]
    end
    q = [60, 10, 10, 20]; α = [0.95, 0.95, 0.95, 0.95]
    prob = FDDESystem(EnzymeKinetics!, q, α, 4, 6)
    sold, sol = solve(prob, 0.1, DelayABM())

    @test isapprox(sold, [ 56.3     11.3272   11.4284  22.4025
    56.3229  11.2293   11.4106  22.4194
    56.3468  11.141    11.3849  22.4332
    56.3731  11.0592   11.3533  22.4434
    56.4026  10.9824   11.3171  22.4495
    56.4358  10.9096   11.2773  22.4516
    56.4729  10.8398   11.2346  22.4493
    56.5144  10.7725   11.1896  22.4429
    56.5604  10.707    11.1429  22.4322
    56.6111  10.6431   11.0946  22.4175
    56.6668  10.5803   11.0452  22.3987
    56.7275  10.5186   10.9947  22.376
    56.7932  10.4577   10.9435  22.3495
    56.864   10.3975   10.8916  22.3194
    56.9398  10.3382   10.8392  22.2857
    57.0205  10.2795   10.7863  22.2488
    57.1059  10.2217   10.7332  22.2086
    57.1961  10.1647   10.6799  22.1653
    57.2908  10.1087   10.6265  22.1191
    57.3897  10.0537   10.5731  22.0702
    57.4928   9.99987  10.5198  22.0187]; atol=1e-3)

    @test isapprox(sol, [60.5329  10.8225  10.6355  20.6245
    60.3612  10.9555  10.6638  20.6616
    60.1978  11.0679  10.7015  20.6988
    60.0407  11.1632  10.7456  20.7379
    59.8889  11.2442  10.7935  20.7792
    59.7419  11.3129  10.8437  20.8231
    59.5991  11.3711  10.8946  20.8697
    59.4603  11.4202  10.9452  20.9187
    59.3251  11.4613  10.9947  20.9701
    59.1935  11.4954  11.0425  21.0235
    59.0652  11.5234  11.0883  21.0786
    58.9401  11.5461  11.1315  21.1351
    58.818   11.5642  11.1722  21.1927
    58.6988  11.5781  11.21    21.2511
    58.5824  11.5883  11.245   21.3098
    58.4686  11.5954  11.2772  21.3687
    58.3575  11.5997  11.3065  21.4274
    58.2489  11.6016  11.3331  21.4857
    58.1427  11.6013  11.357   21.5433
    58.0389  11.5992  11.3783  21.6001
    57.9373  11.5953  11.3971  21.6557]; atol=1e-3)
end

@testset "Test FDDE Matrix Form method" begin
    τ=3.1416; h=0.5; α=0.4; tspan = (0, 1)
    fx0(t) = [sin(t)*cos(t); sin(t)*cos(t); cos(t)^2-sin(t)^2; cos(t)^2-sin(t)^2]
    A=[0 0 1 0; 0 0 0 1; 0 -2 0 0; -2 0 0 0]
    B=[0 0 0 0; 0 0 0 0; -2 0 0 0; 0 -2 0 0]
    f=[0; 0; 0; 0]
    prob = FDDEMatrixProblem(α, τ, A, B, f, fx0, tspan)
    result=solve(prob, h, MatrixForm())

    @test isapprox(result, [0.139708   0.139708   0.96017    0.96017
    0.479462   0.479462   0.283662   0.283662
    0.378401   0.378401  -0.653644  -0.653644
   -0.07056   -0.07056   -0.989992  -0.989992
   -0.454649  -0.454649  -0.416147  -0.416147
   -0.420735  -0.420735   0.540302   0.540302
    0.0        0.0        1.0        1.0
    0.242142   0.242142   0.932342   0.932342
    0.401838   0.401838   0.60134    0.60134]; atol=1e-4)
end


text:
module FittingPWALikelihood

using LinearAlgebra
using PWAHelper

export createLLHandGRAD, minimize

using NLopt
using ForwardDiff

function createLLHandGRAD(PsiDT, form, block_masks)

    Nd, Nwaves = size(PsiDT);
    Tmap = get_parameter_map(block_masks,Nwaves)
    Np = size(Tmap,2)::Int

    pblocks = make_pblock_inds(block_masks)
    vblocks = [let v = fill(0.0,Np)
            v[bl] .= 1.0
            v
        end for bl in pblocks]

    BM = (real.(contract_to_intensity(form,block_masks)))::Array{Float64,2};

    COHTS!(X) = cohts!(X,pblocks)
    COHSQ(X) = cohsq(X,pblocks)
    EXTND(Ψ) = extnd(Ψ,Tmap)
    EXTND!(X,Ψ) = extnd!(X,Ψ,Tmap)

    transposedExtendedPsiDT = Array{Complex{Float64}}(undef,Np,Nd)
    for e in 1:Nd
        EXTND!(@view(transposedExtendedPsiDT[:,e]),@view(PsiDT[e,:]))
    end

    # function LLH(pars)
    #     v = EXTND(@view(transposedPsiDT[:,1]))
    #     res = sum(log,
    #         let
    #             EXTND!(v,@view(transposedPsiDT[:,e]))
    #             v .*= pars
    #             COHSQ(v)
    #         end for e in 1:Nd)
    #     bilin = Nd * (pars' * BM * pars)
    #     return -res + bilin
    # end
    # function LLH(pars)
    #     res = sum(log,
    #         let
    #             v = transposedExtendedPsiDT[:,e]
    #             v .*= pars
    #             COHSQ(v)
    #         end for e in 1:Nd)
    #     bilin = Nd * (pars' * BM * pars)
    #     return -res + bilin
    # end
    function LLH(pars)
        res = sum(log,
                @views sum(abs2,
                    transposedExtendedPsiDT[bl,e] ⋅ pars[bl]
                for bl in pblocks)
                    for e in 1:Nd)
        bilin = Nd * (pars' * BM * pars)
        return -res + bilin
    end

    # function GETDV!(X, psi, pars)
    #     ExtΨ = Array{Complex{Float64}}(undef, length(pars))
    #     EXTND!(ExtΨ, psi)
    #     Y = ExtΨ.*pars
    #     COHTS!(Y)
    #     Y .*= conj.(ExtΨ)
    #     X .= real.(Y)
    #     return
    # end
    # function GETDV!(X, ExtΨ, pars)
    #     Y = ExtΨ.*pars
    #     COHTS!(Y)
    #     Y .*= conj.(ExtΨ)
    #     X .= real.(Y)
    # end
    function GETDV!(X, ExtΨ, Ypars)
        @inbounds Ypars .*= ExtΨ
        @views for bl in pblocks
            @inbounds Ypars[bl] .= sum(Ypars[bl])
        end
        Ypars .*= conj.(ExtΨ)
        X .= real.(Ypars)
    end
    # function GETDV!(X, ExtΨ, pars)
    #     @views for bl in pblocks
    #         X[bl] .= real.(
    #             conj.(ExtΨ[bl]) .*
    #             (pars[bl]' * ExtΨ[bl])
    #         )
    #     end
    # end

    function GETHS(ExtΨ, pars)
        return sum(
            let v = ExtΨ.*bl
                real.(2*v*v')
            end for bl in vblocks)
    end

    function HESSIAN(pars)
        Np = size(pars,1)
        hes = fill(0.0, Np, Np)
        v = Array{Float64}(undef, Np)
        Yp = Array{Complex{Float64}}(undef, Np)
        for e in 1:Nd
            Yp .= pars
            GETDV!(v, @view(transposedExtendedPsiDT[:,e]),Yp)
            deriv = 2*v
            vale = pars' * v
            num = GETHS(@view(transposedExtendedPsiDT[:,e]),pars)*vale - deriv * deriv'
            hes .-= num ./ vale^2
        end
        hes .+= BM* (2Nd);
        return hes;
    end

    function LLH_and_GRAD!(pars, grad)
        val = 0.0; grad .= 0.0
        temp = Array{Float64}(undef,Np)
        Yp = Array{Complex{Float64}}(undef, Np)
        for e in 1:Nd
            Yp .= pars
            GETDV!(temp, @view(transposedExtendedPsiDT[:,e]), Yp)
            vale = (pars' * temp)::Float64
            temp .*= (1.0/vale)
            grad .-= temp
            val -= log(vale);
        end
        grad .*= 2.0
        BB = BM*pars;
        val += Nd * (pars' * BB);
        grad .+= BB .* (2Nd);
        return val;
    end

    function GRAD(pars)
        grad = fill(0.0,Np)
        temp = Array{Float64}(undef,Np)
        Yp = Array{Complex{Float64}}(undef, Np)
        for e in 1:Nd
            Yp .= pars
            GETDV!(temp, @view(transposedExtendedPsiDT[:,e]),Yp)
            vale = (pars' * temp)::Float64
            temp .*= 1/vale
            grad .-= temp
        end
        grad .*= 2.0
        BB = BM*pars;
        grad .+= BB* (2Nd);
        return grad;
    end

    return LLH, GRAD, LLH_and_GRAD!, HESSIAN
end

function minimize(minusLogLikelihood, andDerive!;
                  algorithm::Symbol = :LD_LBFGS,
                  verbose::Int = 0,
                  maxeval::Int = 50000,
                  # parsprecision::Float64 = 1e-4,
                  llhtolerance::Float64 = 1e-4,
                  starting_pars::Vector{Float64} = error("I need the starting parameters! Say starting_pars = [...]"))
    i_iter = 1;
    function to_minimize(x::Vector, grad::Vector)
        if length(grad) > 0
            v = andDerive!(x,grad)
        else
            v = minusLogLikelihood(x)
        end
        i_iter += 1
        verbose==1 && println(i_iter,": ",v)
        verbose==2 && @show v,x
        verbose==3 && @show v,grad
        return v;
    end
    opt = Opt(algorithm, length(starting_pars)) # try LD_LBFGS || LD_MMA || LD_SLSQP
    # xtol_rel!(opt,parsprecision)
    maxeval!(opt,maxeval)
    ftol_abs!(opt,llhtolerance);

    min_objective!(opt, to_minimize)

    (minf,pars,ret) = optimize(opt, starting_pars)#rand(size(TT,2)))
    println("got $minf at $pars after some iterations (returned $ret)")

    pars
end

end


text:
@with_kw struct FirstOrderGeodesicPoint{T,P} <: AbstractGeodesicPoint{T}
    retcode::Symbol
    t::T
    u::AbstractVector{T}
    v::AbstractVector{T}
    p::P
end

function geodesic_point_type(m::CarterMethodBL{T}) where {T}
    p_type = typeof(make_parameters(T(0), T(0), 1, T))
    FirstOrderGeodesicPoint{T,p_type}
end

make_parameters(L, Q, sign_θ, T) =
    (L = L, Q = Q, r = -1, θ = convert(Int, sign_θ), changes = T[0.0, 0.0])

function metric_callback(m::CarterMethodBL{T}) where {T}
    (
        DiscreteCallback(ensure_domain(m), terminate!),
        DiscreteCallback(radial_negative_check(m), flip_radial_sign!),
        DiscreteCallback(angular_negative_check(m), flip_angular_sign!),
    )
end

convert_velocity_type(u::StaticVector{S,T}, v) where {S,T} = convert(SVector{S,T}, v)
convert_velocity_type(u::AbstractVector{T}, v) where {T} = convert(typeof(u), collect(v))

function get_endpoint(
    m::CarterMethodBL{T},
    sol::SciMLBase.AbstractODESolution{T,N,S},
) where {T,N,S}
    us, ts, p = unpack_solution(sol)
    u = us[end]
    v = carter_velocity(u, m.E, m.M, m.a, p)
    t = ts[end]
    FirstOrderGeodesicPoint(sol.retcode, t, u, convert_velocity_type(u, v), p)
end

Vr(m::AbstractFirstOrderMetricParams{T}, u, p) where {T} =
    error("Not implmented for $(typeof(m)).")
Vθ(m::AbstractFirstOrderMetricParams{T}, u, p) where {T} =
    error("Not implmented for $(typeof(m)).")

export FirstOrderGeodesicPoint


text:
### A Pluto.jl notebook ###
# v0.16.0

using Markdown
using InteractiveUtils

# ╔═╡ c17cdcb5-6b4c-4528-a471-304383bd736a
begin
    import Pkg
    # activate a clean environment
    Pkg.activate(mktempdir())
	pkgs =
    Pkg.add(
		["Cascadia","Gumbo","HTTP","DataFrames"] .|>
		name -> Pkg.PackageSpec(name=name))
	using Cascadia, Gumbo, HTTP, DataFrames
end

# ╔═╡ 215ff1a2-d65a-49c3-9a6e-aa394f420342
@eval(PlutoRunner, table_row_display_limit = 270)

# ╔═╡ ac9490eb-1447-4662-a955-b76e22a79f47
md"# Completed"

# ╔═╡ db306e32-94c5-4360-bdcf-19b6346f38d3
md"# All"

# ╔═╡ 3a72039c-97c4-4b64-87f5-eba76879cd24
read_url(url) = url |> HTTP.get |> x -> x.body |> String |> parsehtml

# ╔═╡ e4bf36f4-8382-437d-8c86-14a6808227dd
const daily_schedule_url = "https://www.webtoons.com/en/dailySchedule"

# ╔═╡ 0b7c8808-c6af-4f04-a82b-297cc4e69654
function single_selector(page, selector)
	matches = eachmatch(Selector(selector), page.root)
	@assert(length(matches) == 1, "$(length(matches))")
	matches[1].children |>
	x->filter(x->typeof(x) != HTMLElement{:br}, x) .|>
	child -> child.text
end

# ╔═╡ e453acd5-804a-47b3-b595-d6835582db33
function get_name(page)
	l = single_selector(page, ".detail_header>.info>.subj")
	join(l, " ")
end

# ╔═╡ 0eb299fa-d46d-4166-904d-5f212f39bf91
function get_genre(page)
	l = single_selector(page, ".detail_header>.info>.genre")
	@assert length(l) == 1
	l[1]
end

# ╔═╡ fa6fec9b-3c29-4e6b-97d7-2a7961bb3806
function get_score(page)
	l = single_selector(page, "#_starScoreAverage")
	@assert length(l) == 1
	parse(Float64, l[1])
end

# ╔═╡ 52bf2979-08d8-41ec-81a3-08880961a314
get_is_completed(page) = !all(isempty.(
							[".txt_ico_completed", ".txt_ico_completed2"] .|>
							selector -> eachmatch(Selector(selector), page.root)
						))

# ╔═╡ cd5f059d-44d5-466e-804f-63d8c516283f
function get_last_chapter(page)
	l = single_selector(page, ".detail_lst>#_listUl>:first-child>a>.tx")
	@assert length(l) == 1
	@assert(l[1][1] == '#', l[1])
	parse(Int, l[1][2:end])
end

# ╔═╡ a534f08a-cca2-4456-b77e-9e78c4d3b44f
values_getters = [
	("Score", get_score),
	("Genre", get_genre),
	("Name", get_name),
	("Completed", get_is_completed),
	("Last Chapter", get_last_chapter),
]

# ╔═╡ 3a615405-7861-467f-b29a-df0ae80a6789
function get_data_frame(series_urls)
	series_pages = asyncmap(read_url, series_urls)
	df = DataFrame()
	for (name, getter) ∈ values_getters
		df[:, name] = getter.(series_pages)
	end
	df[:, "Url"] = series_urls .|> url -> HTML("Link")
	sort!(df, rev=true)
end

# ╔═╡ 3486e367-1e67-4e78-9693-8e53e6907e1f
function get_all_series_urls(;completed::Bool=false)
	h = read_url(daily_schedule_url)
	series_selector = ".daily_section>.daily_card>li>a"
	completed && (series_selector = ".comp>" * series_selector)
	as = eachmatch(Selector(series_selector), h.root)
	map(q->q.attributes["href"], as) |> unique!
end

# ╔═╡ 2aaaa8e5-3c82-4c46-b1d4-564885531d1f
get_all_series_urls(completed=true) |> get_data_frame

# ╔═╡ 1111481e-7705-402b-92a9-14e39933ca51
get_all_series_urls() |> get_data_frame

# ╔═╡ Cell order:
# ╠═215ff1a2-d65a-49c3-9a6e-aa394f420342
# ╟─ac9490eb-1447-4662-a955-b76e22a79f47
# ╠═2aaaa8e5-3c82-4c46-b1d4-564885531d1f
# ╠═db306e32-94c5-4360-bdcf-19b6346f38d3
# ╠═1111481e-7705-402b-92a9-14e39933ca51
# ╠═3a72039c-97c4-4b64-87f5-eba76879cd24
# ╠═e4bf36f4-8382-437d-8c86-14a6808227dd
# ╠═c17cdcb5-6b4c-4528-a471-304383bd736a
# ╠═0b7c8808-c6af-4f04-a82b-297cc4e69654
# ╠═e453acd5-804a-47b3-b595-d6835582db33
# ╠═0eb299fa-d46d-4166-904d-5f212f39bf91
# ╠═fa6fec9b-3c29-4e6b-97d7-2a7961bb3806
# ╠═52bf2979-08d8-41ec-81a3-08880961a314
# ╠═cd5f059d-44d5-466e-804f-63d8c516283f
# ╠═a534f08a-cca2-4456-b77e-9e78c4d3b44f
# ╠═3a615405-7861-467f-b29a-df0ae80a6789
# ╠═3486e367-1e67-4e78-9693-8e53e6907e1f


text:
module TestMINRES

using IterativeSolvers
using Test
using Random
using SparseArrays
using LinearAlgebra
using LinearMaps

@testset "MINRES" begin

function hermitian_problem(T, n)
    B = rand(T, n, n) + n * I
    A = B + B'
    x = ones(T, n)
    b = B * x
    A, x, b
end

function skew_hermitian_problem(T, n)
    B = rand(T, n, n) + n * I
    A = B - B'
    x = ones(T, n)
    b = A * x
    A, x, b
end

Random.seed!(123)
n = 15

@testset "Hermitian Matrix{$T}" for T in (Float32, Float64, ComplexF32, ComplexF64)
    A, x, b = hermitian_problem(T, n)
    reltol = sqrt(eps(real(T)))
    x0 = rand(T, n)

    x1, hist1 = minres(A, b, maxiter = 10n, reltol = reltol, log = true)
    x2, hist2 = minres!(x0, A, b, maxiter = 10n, reltol = reltol, log = true)

    @test isa(hist1, ConvergenceHistory)
    @test norm(b - A * x1) / norm(b) ≤ reltol
    @test hist1.isconverged
    @test norm(b - A * x2) / norm(b) ≤ reltol
    @test x2 == x0
end

@testset "Skew-Hermitian Matrix{$T}" for T in (Float32, Float64, ComplexF32, ComplexF64)
    A, x, b = skew_hermitian_problem(T, n)
    reltol = sqrt(eps(real(T)))
    x_approx, hist = minres(A, b, skew_hermitian = true, maxiter = 10n, reltol = reltol, log = true)

    @test norm(b - A * x_approx) / norm(b) ≤ reltol
    @test hist.isconverged
end

@testset "SparseMatrixCSC{$T, $Ti}" for T in (Float32, Float64, ComplexF32, ComplexF64), Ti in (Int64, Int32)
    A = let
        B = sprand(n, n, 2 / n)
        B + B' + I
    end

    x = ones(T, n)
    b = A * x
    reltol = sqrt(eps(real(T)))

    x_approx, hist = minres(A, b, maxiter = 10n, reltol = reltol, log = true)

    @test norm(b - A * x_approx) / norm(b) ≤ reltol
    @test hist.isconverged
end

@testset "Termination criterion" begin
    for T in (Float32, Float64, ComplexF32, ComplexF64)
        A = T[ 2 -1  0
              -1  2 -1
               0 -1  2]
        n = size(A, 2)
        b = ones(T, n)
        x0 = A \ b
        perturbation = 10 * sqrt(eps(real(T))) * T[(-1)^i for i in 1:n]

        # If the initial residual is small and a small relative tolerance is used,
        # many iterations are necessary
        x = x0 + perturbation
        initial_residual = norm(A * x - b)
        x, ch = minres!(x, A, b, log=true)
        @test 2 ≤ niters(ch) ≤ n

        # If the initial residual is small and a large absolute tolerance is used,
        # no iterations are necessary
        x = x0 + perturbation
        initial_residual = norm(A * x - b)
        x, ch = minres!(x, A, b, abstol=2*initial_residual, reltol=zero(real(T)), log=true)
        @test niters(ch) == 0
    end
end

end

end # module


text:
using MakieThemes
using Test
using Makie
using CSV
using DataFrames
using StatsMakie

@testset "simple" begin
  for i in ggthemr_colorthemes()
        scene = show_ggthemr(i)
  end
end

@testset "original" begin
  for dataset ∈ (:www, :drivers, :mtcars, :diamonds)
    path = joinpath(@__DIR__, "..", "data", string(dataset, ".tsv"))
    @eval $(dataset) = CSV.read($path, delim = '\t', allowmissing = :none)
  end

  AbstractPlotting.set_theme!(ggthemr(:fresh))

  p1 = scatterlines(Data(www), :Minute, :Users,
    Group(color = :Measure, marker = :Measure),
    markersize = 6, marker = [:rect, :circle]);

  p2 = plot(density, Data(mtcars),
    :mpg, Group(color = :cyl));

  p3 = plot(Position.stack, histogram, Data(diamonds),
    :price, Group(color = :cut));

  p4 = boxplot(Data(drivers), :Year, :Deaths);

  vbox(hbox(p1, p2), hbox(p3, p4))
end


text:
function [Points,Faces,Cell]=Grid2DCFP(nx,ny,GridType,lenx,leny)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Creation of the grid, with point coordinates, point numeration, face
% numeration, face definition, cell numeration and cell definition.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
option=GridType;
switch option
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% rectangular grid
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  case 'quad'
    EpsRel=0.0;
    Points((nx+1)*(ny+1)).P.x=0;
    iP=0;
    dx=1/nx;
    dy=1/ny;
    for iy=1:ny+1
      for ix=1:nx+1
        iP=iP+1;
        Points(iP).iP=iP;
        Points(iP).x=(ix-1)/nx;
        Points(iP).y=(iy-1)/ny;
        Points(iP).ix=ix;
        Points(iP).iy=iy;
        if ix==1 || ix==nx+1 || iy==1 || iy==ny+1
          Points(iP).Boundary=1;
        else
          Points(iP).Boundary=0;
        end
      end
    end
    for iP=1:size(Points,2)
      if Points(iP).Boundary==0
        Pert=rand()-0.5;
        Points(iP).x=Points(iP).x+Pert*dx*EpsRel;
        Pert=rand()-0.5;
        Points(iP).y=Points(iP).y+Pert*dy*EpsRel;
      end
    end

    Faces(nx*(ny+1)+(nx+1)*ny).P(2)=0;
    iF=0;
    for iy=1:ny
      for ix=1:nx+1
        iF=iF+1;
        Faces(iF).iF=iF;
        Faces(iF).Type='x';
        Faces(iF).P(1)=ix+(nx+1)*(iy-1);
        Faces(iF).P(2)=ix+(nx+1)*iy;
        if ix==1
          Faces(iF).Boundary=1;
        elseif ix==nx+1
          Faces(iF).Boundary=1;
        else
          Faces(iF).Boundary=0;
        end
        Faces(iF).uE=zeros(2,2);
      end
    end
    for iy=1:ny+1
      for ix=1:nx
        iF=iF+1;
        Faces(iF).iF=iF;
        Faces(iF).Type='y';
        Faces(iF).P(1)=ix+(nx+1)*(iy-1);
        Faces(iF).P(2)=ix+1+(nx+1)*(iy-1);
        if iy==1
          Faces(iF).Boundary=1;
        elseif iy==ny+1
          Faces(iF).Boundary=1;
        else
          Faces(iF).Boundary=0;
        end
        Faces(iF).uE=zeros(2,2);
      end
    end

    Cell(nx*ny).P(4)=0;
    iC=0;
    for iy=1:ny
      for ix=1:nx
        iC=iC+1;
        Cell(iC).iC=iC;
        Cell(iC).nP=4;
        Cell(iC).P(1)=ix+(nx+1)*(iy-1);
        Cell(iC).P(2)=ix+1+(nx+1)*(iy-1);
        Cell(iC).P(3)=ix+1+(nx+1)*iy;
        Cell(iC).P(4)=ix+(nx+1)*iy;
        Cell(iC).P(5)=ix+(nx+1)*(iy-1);
        Cell(iC).F(1)=(nx+1)*ny+ix+nx*(iy-1);
        Cell(iC).F(2)=ix+1+(nx+1)*(iy-1);
        Cell(iC).F(3)=(nx+1)*ny+ix+nx*iy;
        Cell(iC).F(4)=ix+(nx+1)*(iy-1);
      end
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% triangular grid on base of the rectangular grid with one diagonal put in
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  case 'quadtri'
    EpsRel=0.0;
    Points((nx+1)*(ny+1)).P.x=0;
    iP=0;
    dx=1/nx;
    dy=1/ny;
    for iy=1:ny+1
      for ix=1:nx+1
        iP=iP+1;
        Points(iP).iP=iP;
        Points(iP).x=(ix-1)/nx;
        Points(iP).y=(iy-1)/ny;
        Points(iP).ix=ix;
        Points(iP).iy=iy;
        if ix==1 || ix==nx+1 || iy==1 || iy==ny+1
          Points(iP).Boundary=1;
        else
          Points(iP).Boundary=0;
        end
      end
    end
    for iP=1:size(Points,2)
      if Points(iP).Boundary==0
        Pert=rand()-0.5;
        Points(iP).x=Points(iP).x+Pert*dx*EpsRel;
        Pert=rand()-0.5;
        Points(iP).y=Points(iP).y+Pert*dy*EpsRel;
      end
    end

    Faces(nx*(ny+1)+(nx+1)*ny).P(2)=0;
    iF=0;
    for iy=1:ny
      for ix=1:nx+1
        iF=iF+1;
        Faces(iF).iF=iF;
        Faces(iF).Type='x';
        Faces(iF).P(1)=ix+(nx+1)*(iy-1);
        Faces(iF).P(2)=ix+(nx+1)*iy;
        if ix==1
          Faces(iF).Boundary=1;
        elseif ix==nx+1
          Faces(iF).Boundary=1;
        else
          Faces(iF).Boundary=0;
        end
        Faces(iF).uE=zeros(2,2);
      end
    end
    for iy=1:ny+1
      for ix=1:nx
        iF=iF+1;
        Faces(iF).iF=iF;
        Faces(iF).Type='y';
        Faces(iF).P(1)=ix+(nx+1)*(iy-1);
        Faces(iF).P(2)=ix+1+(nx+1)*(iy-1);
        if iy==1
          Faces(iF).Boundary=1;
        elseif iy==ny+1
          Faces(iF).Boundary=1;
        else
          Faces(iF).Boundary=0;
        end
        Faces(iF).uE=zeros(2,2);
      end
    end
    for iy=1:ny
      for ix=1:nx
        iF=iF+1;
        Faces(iF).iF=iF;
        Faces(iF).Type='diag';
        Faces(iF).P(1)=ix+(nx+1)*(iy-1);
        Faces(iF).P(2)=ix+1+(nx+1)*(iy);
        Faces(iF).Boundary=0;
        Faces(iF).uE=zeros(2,2);
      end
    end

    Cell(nx*ny).P(4)=0;
    iC=0;
    for iy=1:ny
      for ix=1:nx
        iC=iC+1;
        Cell(iC).iC=iC;
        Cell(iC).nP=3;
        Cell(iC).P(1)=ix+(nx+1)*(iy-1);
        Cell(iC).P(2)=ix+1+(nx+1)*(iy-1);
        Cell(iC).P(3)=ix+1+(nx+1)*iy;
        Cell(iC).P(4)=ix+(nx+1)*(iy-1);
        Cell(iC).F(1)=(nx+1)*ny+ix+nx*(iy-1);
        Cell(iC).F(2)=ix+1+(nx+1)*(iy-1);
        Cell(iC).F(3)=(nx+1)*ny+nx*(ny+1)+ceil(iC/2);
        iC=iC+1;
        Cell(iC).iC=iC;
        Cell(iC).nP=3;
        Cell(iC).P(1)=ix+(nx+1)*(iy-1);
        Cell(iC).P(2)=ix+1+(nx+1)*iy;
        Cell(iC).P(3)=ix+(nx+1)*iy;
        Cell(iC).P(4)=ix+(nx+1)*(iy-1);
        Cell(iC).F(1)=(nx+1)*ny+nx*(ny+1)+ceil(iC/2);
        Cell(iC).F(2)=ix+(nx+1)*(iy-1);
        Cell(iC).F(3)=(nx+1)*ny+ix+nx*iy;
      end
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% triangular grid
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  case 'tri'
    EpsRel=0.0;
    Points(nx*ny+nx+ny+2).P.x=0;
    iP=0;
    dx=1/nx;
    dy=1/ny;
    for iy=1:ny+1
      for ix=1:nx+1
        iP=iP+1;
        Points(iP).iP=iP;
        if mod((iy-1),2)==0
          Points(iP).x=(ix-1)/nx;
          Points(iP).y=(iy-1)/ny;
          Points(iP).ix=ix;
          Points(iP).iy=iy;
        else
          if ix==1
            Points(iP).x=(ix-1)/nx;
            Points(iP).y=(iy-1)/ny;
            Points(iP).ix=ix;
            Points(iP).iy=iy;
            iP=iP+1;
            Points(iP).iP=iP;
            Points(iP).x=(2*ix-1)/(2*nx);
            Points(iP).y=(iy-1)/ny;
            Points(iP).ix=ix;
            Points(iP).iy=iy;
          elseif ix==nx+1
            Points(iP).x=(ix-1)/nx;
            Points(iP).y=(iy-1)/ny;
            Points(iP).ix=ix;
            Points(iP).iy=iy;
          else
            Points(iP).x=(2*ix-1)/(2*nx);
            Points(iP).y=(iy-1)/ny;
            Points(iP).ix=ix;
            Points(iP).iy=iy;
          end
        end
        if ix==1 || ix==nx+1 || iy==1 || iy==ny+1
          Points(iP).Boundary=1;
        else
          Points(iP).Boundary=0;
        end
      end
    end
    for iP=1:size(Points,2)
      if Points(iP).Boundary==0
        Pert=rand()-0.5;
        Points(iP).x=Points(iP).x+Pert*dx*EpsRel;
        Pert=rand()-0.5;
        Points(iP).y=Points(iP).y+Pert*dy*EpsRel;
      end
    end

    Faces(2*nx*ny+2*ny+ceil((ny+1)/2)*nx+floor((ny+1)/2)*(nx+1)).P(2)=0;
    iF=0;
    for iP=1:size(Points,2)-nx
      if and(Points(iP).x==0,Points(iP).y<1)
        iF=iF+1;
        Faces(iF).iF=iF;
        Faces(iF).Type='x';
        if mod(Points(iP).iy-1,2)==0
          Faces(iF).P(1)=iP;
          Faces(iF).P(2)=iP+nx+1;
        else
          Faces(iF).P(1)=iP;
          Faces(iF).P(2)=iP+nx+2;
        end
        Faces(iF).Boundary=1;
        Faces(iF).uE=zeros(2,2);
      elseif and(Points(iP).x==1,Points(iP).y<1)
        iF=iF+1;
        Faces(iF).iF=iF;
        Faces(iF).Type='x';
        if mod(Points(iP).iy-1,2)==0
          Faces(iF).P(1)=iP;
          Faces(iF).P(2)=iP+nx+2;
        else
          Faces(iF).P(1)=iP;
          Faces(iF).P(2)=iP+nx+1;
        end
        Faces(iF).Boundary=1;
        Faces(iF).uE=zeros(2,2);
      end
    end
    iF=iF+1;
    for iP=1:size(Points,2)-1
      Faces(iF).iF=iF;
      Faces(iF).Type='y';
      if Points(iP).x~=1
        Faces(iF).P(1)=iP;
        Faces(iF).P(2)=iP+1;
        if or(Points(iP).y==0,Points(iP).y==1)
          Faces(iF).Boundary=1;
        else
          Faces(iF).Boundary=0;
        end
        Faces(iF).uE=zeros(2,2);
        iF=iF+1;
      end
    end
    iF=iF-1;
    for iP=1:size(Points,2)
      if and(and(mod(Points(iP).iy-1,2)==1,Points(iP).x>0),...
          Points(iP).x<1)
        iF=iF+1;
        Faces(iF).iF=iF;
        Faces(iF).Type='diag';
        Faces(iF).P(1)=iP-nx-2;
        Faces(iF).P(2)=iP;
        Faces(iF).Boundary=0;
        Faces(iF).uE=zeros(2,2);
        iF=iF+1;
        Faces(iF).iF=iF;
        Faces(iF).Type='diag';
        Faces(iF).P(1)=iP-nx-1;
        Faces(iF).P(2)=iP;
        Faces(iF).Boundary=0;
        Faces(iF).uE=zeros(2,2);
      end
    end
    for iP=1:size(Points,2)
      if and(and(and(mod(Points(iP).iy-1,2)==1,Points(iP).x>0),...
          Points(iP).x<1),Points(iP).y<1)
        iF=iF+1;
        Faces(iF).iF=iF;
        Faces(iF).Type='diag';
        Faces(iF).P(1)=iP;
        Faces(iF).P(2)=iP+nx+1;
        Faces(iF).Boundary=0;
        Faces(iF).uE=zeros(2,2);
        iF=iF+1;
        Faces(iF).iF=iF;
        Faces(iF).Type='diag';
        Faces(iF).P(1)=iP;
        Faces(iF).P(2)=iP+nx+2;
        Faces(iF).Boundary=0;
        Faces(iF).uE=zeros(2,2);
      end
    end

    Cell((2*nx+1)*ny).P(4)=0;
    iC=0;
    for iP=1:size(Points,2)
      if and(mod(Points(iP).iy-1,2)==0,Points(iP).y<1)
        if Points(iP).x==1
          iC=iC+1;
          Cell(iC).iC=iC;
          Cell(iC).nP=3;
          Cell(iC).P(1)=iP;
          Cell(iC).P(2)=iP+nx+2;
          Cell(iC).P(3)=iP+nx+1;
          Cell(iC).P(4)=iP;
        else
          iC=iC+1;
          Cell(iC).iC=iC;
          Cell(iC).nP=3;
          Cell(iC).P(1)=iP;
          Cell(iC).P(2)=iP+nx+2;
          Cell(iC).P(3)=iP+nx+1;
          Cell(iC).P(4)=iP;
          iC=iC+1;
          Cell(iC).iC=iC;
          Cell(iC).nP=3;
          Cell(iC).P(1)=iP;
          Cell(iC).P(2)=iP+1;
          Cell(iC).P(3)=iP+nx+2;
          Cell(iC).P(4)=iP;
        end
      elseif and(mod(Points(iP).iy-1,2)==1,Points(iP).y<1)
        if Points(iP).x==0
          iC=iC+1;
          Cell(iC).iC=iC;
          Cell(iC).nP=3;
          Cell(iC).P(1)=iP;
          Cell(iC).P(2)=iP+1;
          Cell(iC).P(3)=iP+nx+2;
          Cell(iC).P(4)=iP;
        elseif and(Points(iP).x<1,Points(iP).x>0)
          iC=iC+1;
          Cell(iC).iC=iC;
          Cell(iC).nP=3;
          Cell(iC).P(1)=iP;
          Cell(iC).P(2)=iP+nx+2;
          Cell(iC).P(3)=iP+nx+1;
          Cell(iC).P(4)=iP;
          iC=iC+1;
          Cell(iC).iC=iC;
          Cell(iC).nP=3;
          Cell(iC).P(1)=iP;
          Cell(iC).P(2)=iP+1;
          Cell(iC).P(3)=iP+nx+2;
          Cell(iC).P(4)=iP;
        end
      end
    end

    for iC=1:size(Cell,2)
      for iF=1:size(Faces,2)
        for i=1:Cell(iC).nP
          if and(Cell(iC).P(i)==Faces(iF).P(1),Cell(iC).P(i+1)==Faces(iF).P(2))
            Cell(iC).F(1,i)=Faces(iF).iF;
          elseif and(Cell(iC).P(i)==Faces(iF).P(2),Cell(iC).P(i+1)==Faces(iF).P(1))
            Cell(iC).F(1,i)=Faces(iF).iF;
          end
        end
      end
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Hexagonal Grid in Progress
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  case 'hex'
    ny=2*ny;
    EpsRel=0.0;
    if mod(nx,2)==1
      Points((ny/2+1)*(nx+1)+(ny/2-1)*nx+nx+2).P.x=0;
    else
      Points((ny/2)*(nx+1)+(ny/2-1)*nx+2*(nx+2)).P.x=0;
    end
    iP=0;
    dx=lenx/nx;
    dy=leny/ny;
    for iy=1:ny+1
      for ix=1:nx+1
        iP=iP+1;
        Points(iP).iP=iP;
        if or(iy==1,iy==ny+1)
          if or(mod((iy-1),4)==1,mod(iy-1,4)==2)
            Points(iP).x=(ix-1)/nx;
            Points(iP).y=(iy-1)/(ny);
            Points(iP).ix=ix;
            Points(iP).iy=iy;
          else
            if ix==1
              Points(iP).x=(ix-1)/nx;
              Points(iP).y=(iy-1)/ny;
              Points(iP).ix=ix;
              Points(iP).iy=iy;
              iP=iP+1;
              Points(iP).iP=iP;
              Points(iP).x=(2*ix-1)/(2*nx);
              Points(iP).y=(iy-1)/ny;
              Points(iP).ix=ix;
              Points(iP).iy=iy;
            elseif ix==nx+1
              Points(iP).x=(ix-1)/nx;
              Points(iP).y=(iy-1)/ny;
              Points(iP).ix=ix;
              Points(iP).iy=iy;
            else
              Points(iP).x=(2*ix-1)/(2*nx);
              Points(iP).y=(iy-1)/ny;
              Points(iP).ix=ix;
              Points(iP).iy=iy;
            end
          end
        else
          if or(mod((iy-1),4)==1,mod(iy-1,4)==2)
            Points(iP).x=(ix-1)/nx;
            Points(iP).y=(iy-1)/(ny);
            Points(iP).ix=ix;
            Points(iP).iy=iy;
          else
            if (2*ix-1)/(2*nx)>1
              iP=iP-1;
              continue
            else
              Points(iP).iP=iP;
              Points(iP).x=(2*ix-1)/(2*nx);
              Points(iP).y=(iy-1)/ny;
              Points(iP).ix=ix;
              Points(iP).iy=iy;
            end
          end
        end
        if ix==1 || ix==nx+1 || iy==1 || iy==ny+1
          Points(iP).Boundary=1;
        else
          Points(iP).Boundary=0;
        end
      end
    end
    for iP=1:size(Points,2)
      if Points(iP).Boundary==0
        Pert=rand()-0.5;
        Points(iP).x=Points(iP).x+Pert*dx*EpsRel;
        Pert=rand()-0.5;
        Points(iP).y=Points(iP).y+Pert*dy*EpsRel;
      end
    end

    Faces(nx*ny).P(2)=0;
    iF=1;
    for iP=1:size(Points,2)-1
      for iPP=iP+1:size(Points,2)
        if Points(iP).x==Points(iPP).x
          if and(and(Points(iPP).y-Points(iP).y>=2/ny,Points(iP).x>0)...
              ,Points(iP).x<1)
            break
          else
            Faces(iF).iF=iF;
            Faces(iF).Type='x';
            Faces(iF).P(1)=iP;
            Faces(iF).P(2)=iPP;
            if and(Points(iP).x==0,Points(iPP).x==0)
              Faces(iF).Boundary=1;
            elseif and(Points(iP).x==1,Points(iPP).x==1)
              Faces(iF).Boundary=1;
            else
              Faces(iF).Boundary=0;
            end
            iF=iF+1;
            break;
          end
        end
      end
    end
    for iP=1:size(Points,2)-1
      if Points(iP).y==Points(iP+1).y
        if or(Points(iP).y==1,Points(iP).y==0)
          Faces(iF).iF=iF;
          Faces(iF).Type='y';
          Faces(iF).P(1)=iP;
          Faces(iF).P(2)=iP+1;
          Faces(iF).Boundary=1;
          iF=iF+1;
        end
      end
    end
    iF=iF-1;
    for iP=1:size(Points,2)
      if mod(Points(iP).iy-1,4)==1
        if Points(iP).iy==2
          if Points(iP).x==0
            iF=iF+1;
            Faces(iF).iF=iF;
            Faces(iF).Type='diag';
            Faces(iF).P(1)=iP-nx-1;
            Faces(iF).P(2)=iP;
            Faces(iF).Boundary=0;
            Faces(iF).uE=zeros(2,2);
          elseif Points(iP).x==1
            iF=iF+1;
            Faces(iF).iF=iF;
            Faces(iF).Type='diag';
            Faces(iF).P(1)=iP-nx-2;
            Faces(iF).P(2)=iP;
            Faces(iF).Boundary=0;
            Faces(iF).uE=zeros(2,2);
          else
            iF=iF+1;
            Faces(iF).iF=iF;
            Faces(iF).Type='diag';
            Faces(iF).P(1)=iP-nx-2;
            Faces(iF).P(2)=iP;
            Faces(iF).Boundary=0;
            Faces(iF).uE=zeros(2,2);
            iF=iF+1;
            Faces(iF).iF=iF;
            Faces(iF).Type='diag';
            Faces(iF).P(1)=iP-nx-1;
            Faces(iF).P(2)=iP;
            Faces(iF).Boundary=0;
            Faces(iF).uE=zeros(2,2);
          end
        else
          if Points(iP).x==0
            iF=iF+1;
            Faces(iF).iF=iF;
            Faces(iF).Type='diag';
            Faces(iF).P(1)=iP-nx;
            Faces(iF).P(2)=iP;
            Faces(iF).Boundary=0;
            Faces(iF).uE=zeros(2,2);
          elseif Points(iP).x==1
            iF=iF+1;
            Faces(iF).iF=iF;
            Faces(iF).Type='diag';
            Faces(iF).P(1)=iP-nx-1;
            Faces(iF).P(2)=iP;
            Faces(iF).Boundary=0;
            Faces(iF).uE=zeros(2,2);
          else
            iF=iF+1;
            Faces(iF).iF=iF;
            Faces(iF).Type='diag';
            Faces(iF).P(1)=iP-nx-1;
            Faces(iF).P(2)=iP;
            Faces(iF).Boundary=0;
            Faces(iF).uE=zeros(2,2);
            iF=iF+1;
            Faces(iF).iF=iF;
            Faces(iF).Type='diag';
            Faces(iF).P(1)=iP-nx;
            Faces(iF).P(2)=iP;
            Faces(iF).Boundary=0;
            Faces(iF).uE=zeros(2,2);
          end
        end
      end
    end
    for iP=1:size(Points,2)
      if mod(Points(iP).iy-1,4)==3
        iF=iF+1;
        Faces(iF).iF=iF;
        Faces(iF).Type='diag';
        Faces(iF).P(1)=iP;
        Faces(iF).P(2)=iP-nx-1;
        Faces(iF).Boundary=0;
        Faces(iF).uE=zeros(2,2);
        iF=iF+1;
        Faces(iF).iF=iF;
        Faces(iF).Type='diag';
        Faces(iF).P(1)=iP;
        Faces(iF).P(2)=iP-nx;
        Faces(iF).Boundary=0;
        Faces(iF).uE=zeros(2,2);
      end
    end

    Cell(nx*ny/2).P(4)=0;
    iC=0;
    for iP=1:size(Points,2)-1
      if and(Points(iP).y==0,Points(iP+1).y==0)
        iC=iC+1;
        Cell(iC).iC=iC;
        Cell(iC).nP=3;
        Cell(iC).P(1)=iP;
        Cell(iC).P(2)=iP+1;
        Cell(iC).P(3)=iP+nx+2;
        Cell(iC).P(4)=iP;
      elseif and(Points(iP).y==1,Points(iP+1).y==1)
        if and(Points(iP).x==0,mod(ny,4)==0)
          iC=iC+1;
          Cell(iC).iC=iC;
          Cell(iC).nP=4;
          Cell(iC).P(1)=iP;
          Cell(iC).P(2)=iP-2*nx-1;
          Cell(iC).P(3)=iP-nx;
          Cell(iC).P(4)=iP+1;
          Cell(iC).P(5)=iP;
        elseif and(Points(iP+1).x==1,mod(ny,4)==0)
          iC=iC+1;
          Cell(iC).iC=iC;
          Cell(iC).nP=4;
          Cell(iC).P(1)=iP;
          Cell(iC).P(2)=iP-nx-1;
          Cell(iC).P(3)=iP-2*nx-1;
          Cell(iC).P(4)=iP+1;
          Cell(iC).P(5)=iP;
        else
          iC=iC+1;
          Cell(iC).iC=iC;
          Cell(iC).nP=5;
          Cell(iC).P(1)=iP;
          Cell(iC).P(2)=iP-nx-1;
          Cell(iC).P(3)=iP-2*nx-1;
          Cell(iC).P(4)=iP-nx;
          Cell(iC).P(5)=iP+1;
          Cell(iC).P(6)=iP;
        end
      end
      if iP<size(Points,2)-nx
        if and(and(Points(iP).x==Points(iP+nx).x,Points(iP).y~=1),...
            Points(iP).x<1-1/(2*nx))
          iC=iC+1;
          Cell(iC).iC=iC;
          Cell(iC).nP=6;
          Cell(iC).P(1)=iP;
          Cell(iC).P(2)=iP-nx;
          Cell(iC).P(3)=iP+1;
          Cell(iC).P(4)=iP+nx+1;
          Cell(iC).P(5)=iP+2*nx+1;
          Cell(iC).P(6)=iP+nx;
          Cell(iC).P(7)=iP;
        end
      end
      if iP<size(Points,2)-nx-1
        if and(and(and(and(Points(iP).x==Points(iP+nx+1).x,...
            Points(iP).y~=1),Points(iP).y~=0),Points(iP).x~=1),...
            Points(iP).iy<ny-1)
          if Points(iP).iy==2
            iC=iC+1;
            Cell(iC).iC=iC;
            Cell(iC).nP=6;
            Cell(iC).P(1)=iP;
            Cell(iC).P(2)=iP-nx-1;
            Cell(iC).P(3)=iP+1;
            Cell(iC).P(4)=iP+nx+2;
            Cell(iC).P(5)=iP+2*nx+2;
            Cell(iC).P(6)=iP+nx+1;
            Cell(iC).P(7)=iP;
          else
            iC=iC+1;
            Cell(iC).iC=iC;
            Cell(iC).nP=6;
            Cell(iC).P(1)=iP;
            Cell(iC).P(2)=iP-nx;
            Cell(iC).P(3)=iP+1;
            Cell(iC).P(4)=iP+nx+2;
            Cell(iC).P(5)=iP+2*nx+2;
            Cell(iC).P(6)=iP+nx+1;
            Cell(iC).P(7)=iP;
          end
        end
      end
      if and(iP<size(Points,2)-3*nx-1,ny~=4)
        if and(and(Points(iP).x==0,Points(iP+3*nx+1).x==0),Points(iP).iy<ny-1)
          iC=iC+1;
          Cell(iC).iC=iC;
          Cell(iC).nP=4;
          Cell(iC).P(1)=iP;
          Cell(iC).P(2)=iP+nx+1;
          Cell(iC).P(3)=iP+2*nx+1;
          Cell(iC).P(4)=iP+3*nx+1;
          Cell(iC).P(5)=iP;
        elseif and(Points(iP).x==1,Points(iP+3*nx+1).x==1)
          iC=iC+1;
          Cell(iC).iC=iC;
          Cell(iC).nP=4;
          Cell(iC).P(1)=iP;
          Cell(iC).P(2)=iP+3*nx+1;
          Cell(iC).P(3)=iP+2*nx;
          Cell(iC).P(4)=iP+nx;
          Cell(iC).P(5)=iP;
        end
      end
    end

    for iC=1:size(Cell,2)
      for iF=1:size(Faces,2)
        for i=1:Cell(iC).nP
          if and(Cell(iC).P(i)==Faces(iF).P(1),Cell(iC).P(i+1)==Faces(iF).P(2))
            Cell(iC).F(1,i)=Faces(iF).iF;
          elseif and(Cell(iC).P(i)==Faces(iF).P(2),Cell(iC).P(i+1)==Faces(iF).P(1))
            Cell(iC).F(1,i)=Faces(iF).iF;
          end
        end
      end
    end
end
end


text:
export Sequential,learn!

mutable struct Sequential <: Models
    common
    window_size
    Sequential(gpu_flg=false)=new(ModelCommon(gpu_flg),0)
end

function learn!(model::Models;data,t_data,window_size,max_epoch,verification_params=nothing,verification=nothing,notebook=false)
    D = size(data,3) #データ次元数
    T = window_size #内部ユニット数
    N = size(data,1) #バッチ数
    max_ite = size(t_data,2) #イテレーション数
    continue_flg = true
    p = Progress(max_epoch*max_ite,1,"Progress : ")
    model.common.epoch = max_epoch
    model.common.loss = zeros(max_epoch)
    model.window_size = window_size

    @inbounds for epoch in 1:max_epoch
        ite_total_loss = 0 #損失合計
        avg_loss = 0 #1エポックの平均損失
        model.common.now_epoch = epoch
        @inbounds for ite in 1:max_ite
            #ミニバッチ作成
            st = Int(1+(ite-1)*T) #data切り取り位置
            ed = Int(T*ite)
            xs = view(data,:,st:ed,:)
            t = view(t_data,:,ite,:)
            #順伝播
            model.common.layers[end].t = t #教師データ挿入
            loss = predict!(model,xs,true)
            #ite毎の平均損失を加算
            ite_total_loss += sum(loss)/length(loss)
            #逆伝播
            backward!(model)
            #更新
            update!(model.common.optimizer,model)

            next!(p)
        end

        reset!(model)

        avg_loss = ite_total_loss/max_ite
        model.common.loss[epoch] = avg_loss

        if notebook
            print(model.common.learn_io,"\n","\e[0F","\e[2K","epoch: ",epoch," | loss: ",avg_loss," | ")
        else
            print(model.common.learn_io,"\n","epoch:",epoch," | ","loss:",avg_loss," | ")
        end

        #検証
        if verification!=nothing
            continue_flg = verification(model,verification_params)
            reset!(model)
        end

        if notebook
            IJulia.clear_output(true)
            plot(model.common.learn_plot) |> display
        else
            print(model.common.learn_io,"\e[2F") #上(プログレスバーの先頭へ)
        end

        take!(model.common.learn_io) |> String |> println

        # min_avg_loss，その時のep,今のep,avg_loss
        if model.common.min_loss > avg_loss || epoch==1
            model.common.min_loss = avg_loss
            model.common.min_epoch = epoch
        end

    end

    return 0
end

max_stars_repo_name: catmanium/Spike

max_issues_repo_name: catmanium/Spike

max_forks_repo_name: catmanium/Spike


text:
module BVHFiles

using LinearAlgebra, Rotations, LightGraphs, Flux

include("BVHGraphs.jl")
include("utils.jl")
include("file_io.jl")
include("global_positions.jl")
include("transformation.jl")
include("optimization.jl")

end


text:
using Transformers.Basic
using Transformers.BidirectionalEncoder:WordPiece
using PreSumm.Model

include("loss.jl")

function preprocess(text::String, wordpiece::WordPiece, tokenizer::Function)::AbstractVector{String}
    tokens = text |> tokenizer |> wordpiece
    max_length = min(4096, length(tokens)) # Truncate to 4096 tokens
    tokens = tokens[1:max_length]
    return ["[CLS]"; tokens; "[SEP]"]
end

label_smoothing_α = 0.0 # Label smoothing doesn't work yet.
function loss(
    inputs::AbstractVector{String},
    outputs::AbstractVector{String},
    ground_truth::AbstractMatrix{<:Number},
    model::Translator,
    vocabulary::Vocabulary,
)::AbstractFloat
    prediction = model.transformers(vocabulary, inputs, outputs)
    loss = logtranslationloss(prediction, ground_truth, α=label_smoothing_α)
    return loss
end


text:
module TestVector

import ..repl_output
using AMGX, Defer, Test, JSON, CUDA
using AMGX: Config, Resources, AMGXVector, dDDI, dFFI

@scope @testset "Vector" begin
    c = @! Config("")
    r = @! Resources(c)

    @scope @testset "upload/download" begin
        v = @! AMGXVector(r, dDDI)
        @test occursin("dDDI", repl_output(v))
        v_h = [1.0, 2.0, 3.0]
        AMGX.upload!(v, v_h)
        @test occursin("of length 3", repl_output(v))
        @test length(v) == 3
        @test AMGX.download(v) == v_h
        @test Vector(v) == v_h
        @test Array(v) == v_h

        # Copy in place
        v_h_buffer = zeros(Float64, 3)
        copy!(v_h_buffer, v)
        @test v_h_buffer == v_h

        # Copy in place wrong size
        v_h_buffer_wrong_length = zeros(Float64, 4)
        @test_throws ArgumentError copy!(v_h_buffer_wrong_length, v)

        # Wrong element type
        v_h = Float32[1.0, 2.0, 3.0]
        @test_throws ArgumentError AMGX.upload!(v, v_h)

        # Float32
        v = @! AMGXVector(r, dFFI)
        v_h = Float32[1.0, 2.0, 3.0]
        AMGX.upload!(v, v_h)
        @test length(v) == 3
        @test AMGX.download(v) == v_h

        # Set zero
        v = @! AMGXVector(r, dDDI)
        AMGX.set_zero!(v, 5)
        @test length(v) == 5
        @test Vector(v) == zeros(Float64, 5)

        # Blocked
        v = @! AMGXVector(r, dDDI)
        v_h = [1.0, 2.0, 3.0, 4.0]
        AMGX.upload!(v, v_h; block_dim = 2)
        @test occursin("of length 2⋅2", repl_output(v))
        @test length(v_h) == 4
        @test AMGX.vector_get_size(v) == (2, 2)
        v_h = [1.0, 2.0, 3.0]
        @test_throws ArgumentError AMGX.upload!(v, v_h; block_dim = 2)
    end

    @scope @testset "upload/download cuarray" begin
        v = @! AMGXVector(r, dDDI)
        v_h = [1.0, 2.0, 3.0]
        v_d = CuArray(v_h)
        AMGX.upload!(v, v_d)
        @test length(v) == 3
        @test AMGX.download(v) == v_h

        # Copy in place
        v_d_buffer = CUDA.zeros(Float64, 3)
        copy!(v_d_buffer, v)
        @test Vector(v_d_buffer) == v_h
        v_h_buffer = zeros(Float64, 3)
        copy!(v_h_buffer, v)
        @test v_h_buffer == v_h
        v_d_buffer = CUDA.zeros(Float32, 3)
        @test_throws ArgumentError copy!(v_d_buffer, v)

        # Float32
        v = @! AMGXVector(r, dFFI)
        v_h = Float32[1.0, 2.0, 3.0]
        v_d = CuArray(v_h)
        AMGX.upload!(v, v_d)
        @test Vector(v) == v_h
    end
end
end # module


text:
using Parameters

abstract type AbstractViscosity end
export AbstractViscosity

struct EmptyViscosity <: AbstractViscosity end
export EmptyViscosity

@with_kw struct ConstantViscosity{T<:AbstractFloat} <: AbstractViscosity
    mu::T
end

(cv::ConstantViscosity)(T::Q) where {Q<:AbstractFloat} = cv.mu
export ConstantViscosity

@with_kw struct RiedelViscosity{N<:AbstractFloat} <: AbstractViscosity
    A::N
    B::N
    C::N
    D::N
    E::N
end

(rv::RiedelViscosity)(T::Q) where {Q<:Number} = @fastmath exp(rv.A+rv.B/T+rv.C*log(T)+rv.D*T^rv.E)
export RiedelViscosity


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -33.59531116254298 <= q <= 31.312008728451346)

                     add_NL_constraint(m, :(sigmoid(0.7464435003747418 + -0.7676428720221384*sigmoid(-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4])) + -0.5937285886485864*sigmoid(0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4])) + 0.70936536775501*sigmoid(-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4])) + 0.9647291374973435*sigmoid(0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4])) + 0.527810271934432*sigmoid(-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))) + sigmoid(-0.35636357579153266 + 0.8665784119884141*sigmoid(-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4])) + 0.03175667478129274*sigmoid(0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4])) + -0.7799384758999897*sigmoid(-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4])) + 0.4104070349373172*sigmoid(0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4])) + -0.7416668360909755*sigmoid(-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))) + sigmoid(0.5038932713674513 + 0.7223377417739489*sigmoid(-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4])) + -0.17110213406275054*sigmoid(0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4])) + 0.7726394930743403*sigmoid(-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4])) + 0.6118790711423934*sigmoid(0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4])) + 0.16250007841990932*sigmoid(-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))) + sigmoid(-0.7190092552996492 + 0.6952008118484776*sigmoid(-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4])) + 0.975536624348428*sigmoid(0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4])) + 0.03367235210941821*sigmoid(-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4])) + 0.8059867473566515*sigmoid(0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4])) + 0.5874603423029647*sigmoid(-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))) + sigmoid(-0.5715134646750415 + 0.6325466436773857*sigmoid(-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4])) + -0.2633742432221431*sigmoid(0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4])) + -0.3786600783719858*sigmoid(-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4])) + 0.5923927477083897*sigmoid(0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4])) + 0.9114166052727355*sigmoid(-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
mutable struct BinarySearchTree{T,K} <: AbstractDictionary{T,K}
    root::BinarySearchNode{T}
    size::Int
    function BinarySearchTree{T,K}() where {T,K}
        root = BinarySearchNode{T}()
        size = 0
        new(root, size)
    end
end

#FIXME You need a fucking EmptyNode Type

### Helper methods

function _insert(r::BinarySearchNode, key, item)

end

function _remove(r::BinarySearchNode, key)

end

function _find(r::BinarySearchNode, key)
    if (key < r.key) && (r.leftchild !== r)
        _find(r.leftchild, key)
    elseif (key > r.key) && (r.rightchild !== r)
        _find(r.rightchild, key)
    elseif key == r.key
        return r.data
    else
        return nothing
    end
end

# Non-Null approach avoids reassigning every child to every
# node along the path from the root to the parent of the new
# node. However, at less simple insertion process (see book).
function _insert!(r::BinarySearchNode, key, item)
    if (key < r.key)
        if r.leftchild === r
            setleftchild!(r, BinarySearchNode(key, item))
        else
            _insert!(r.leftchild, key, item)
        end
    else
        if r.rightchild === r
            setrightchild!(r, BinarySearchNode(key, item))
        else
            _insert!(r.rightchild, key, item)
        end
    end
end

### Interface

isempty(d::BinarySearchTree) = d.size == 0
length(d::BinarySearchTree) = d.size

function find(t::BinarySearchTree, key)
    (t.size == 0) && throw(ArgumentError("Empty tree"))
    return _find(t.root, key)
end

function insert!(t::BinarySearchTree, key, item)
    if t.size == 0
        setkey!(t.root, key)
        setvalue!(t.root, item)
    else
        _insert!(t.root, key, item)
    end
    t.size += 1
end


text:
const MOD_REPORT_CHANNEL_NAME = "mod-report"
const MOD_REPORT_CHANNEL = Ref{DiscordChannel}()
const BAD_WORDS = Vector{String}()

function handler(c::Client, e::MessageCreate, ::Val{:mod})
    # @info "mod handler: $(e.message.content)"

    mod_report_channel = mod_get_report_channel(c, e.message.guild_id)
    if mod_report_channel === nothing
        return # config error; nothing we can do about...
    end
    if e.message.channel_id == mod_report_channel.id
        return # Do not moderate the mod channel itself
    end

    mod_init() # lazy initialization
    result, word_involed = mod_check(e.message.content)
    if result != :good
        username = e.message.author.username
        report = mod_report(username, e.message.content, result, word_involed)
        if mod_report_channel !== nothing
            create(c, Message, mod_report_channel; content = report)
        end
    end
    return nothing
end

"Get the mod-report channel."
function mod_get_report_channel(c::Client, guild_id::UInt64)
    if !isassigned(MOD_REPORT_CHANNEL)
        channels = @discord get_guild_channels(c, guild_id)
        channel = filter(c -> c.name == MOD_REPORT_CHANNEL_NAME, channels)
        if length(channel) == 1
            MOD_REPORT_CHANNEL[] = channel[1]
            return MOD_REPORT_CHANNEL[]
        else
            @error "Unable to find mod report channel: $MOD_REPORT_CHANNEL_NAME"
            return nothing
        end
    end
    return MOD_REPORT_CHANNEL[]
end

"Initialize mod function e.g. read bad words list."
function mod_init()
    if isempty(BAD_WORDS)
        bad_words_file = joinpath(@__DIR__, "..", "..", "config", "mod", "bad-words.txt")
        bad_words = readlines(bad_words_file)
        append!(BAD_WORDS, bad_words)
    end
end

"Returns true if the word string is questionably bad."
mod_is_questionable(w::AbstractString) = startswith(w, "?")

"Return true if the string contains any non-alpha/numeric character."
mod_contains_symbols(w::AbstractString) = occursin(r"[^0-9a-zA-Z]", w)

"Return true if `content` contains `word` considering word boundaries."
function mod_contains(content::AbstractString, word::AbstractString)
    content = lowercase(content)
    try
        regex = Regex("\\b" * word * "\\b")
        return match(regex, content) !== nothing
    catch ex
        @error "cannot compile: $word"
        Base.showerror(stdout, ex)
        return false
    end
end

"Return true if `word` appears anywhere in `content` without considering word boundaries."
function mod_contains_exactly(content::AbstractString, word::AbstractString)
    content = lowercase(content)
    return occursin(word, content)
end

"""
    mod_check(content::AbstractString)

Check the content against a profanity set of words. Returns status as
one of the symbols - `:good`, `:bad`, or `:questionable`.
"""
function mod_check(content::AbstractString)
    for bad_word in BAD_WORDS
        questionable = mod_is_questionable(bad_word)
        word = questionable ? bad_word[2:end] : bad_word
        if (mod_contains_symbols(word) && mod_contains_exactly(content, word)) ||
                mod_contains(content, word)
            @info "Bad word `$word` detected: $content"
            return questionable ? (:questionable, word) : (:bad, word)
        end
    end
    return (:good, "")
end

function mod_report(
    username::AbstractString,
    content::AbstractString,
    result::Symbol,
    word_involved::AbstractString,
)
    return if result == :bad
        "Bad message from $username with `$word_involved`: $(content)"
    elseif result == :questionable
        "Questionable message from $username with `$word_involved`: $(content)"
    else # should never happen
        "Unclear how bad is the message from $username with `$word_involved`: $(content)"
    end
end