text:
"""
    PersistentDict{K,V} <: AbstractDict{K,V}
PersistentDict is the abstract subtype for all implementations of persistent maps such as
PersistentHashDict, PersistentAutoDict, and PersistentArrayDict.

To construct a PersistentDict use the PersistentDict constructor.
"""
abstract type PersistentDict{K, V} <: AbstractDict{K, V} end

struct NotFound end

struct PersistentArrayDict{K, V} <: PersistentDict{K, V}
    kvs::Vector{Pair{K, V}}
end
PersistentArrayDict{K, V}() where {K, V} =
    PersistentArrayDict{K, V}(Pair{K, V}[])
PersistentArrayDict(kvs::(Union{Tuple{K, V}, Pair{K, V}})...) where {K, V} =
    PersistentArrayDict{K, V}(Pair{K, V}[Pair(k, v) for (k, v) in kvs])
PersistentArrayDict(; kwargs...) = PersistentArrayDict(kwargs...)

Base.isequal(m1::PersistentArrayDict, m2::PersistentArrayDict) =
    isequal(Set(m1.kvs), Set(m2.kvs))
==(m1::PersistentArrayDict, m2::PersistentArrayDict) =
    Set(m1.kvs) == Set(m2.kvs)

Base.length(m::PersistentArrayDict)  = length(m.kvs)
Base.isempty(m::PersistentArrayDict) = length(m) == 0

findkeyidx(m::PersistentArrayDict, k) = findfirst(kv -> kv[1] == k, m.kvs)

function _get(m::PersistentArrayDict, k, default, hasdefault::Bool)
    for kv in m.kvs
        kv[1] == k && return kv[2]
    end
    hasdefault ? default : default()
end

Base.get(m::PersistentArrayDict, k) =
    _get(m, k, ()->error("key not found: $k"), false)
Base.get(m::PersistentArrayDict, k, default) =
    _get(m, k, default, true)
Base.getindex(m::PersistentArrayDict, k) = get(m, k)

Base.haskey(m::PersistentArrayDict, k) = get(m, k, NotFound()) != NotFound()

function assoc(m::PersistentArrayDict{K, V}, k, v) where {K, V}
    idx = findkeyidx(m, k)
    idx === nothing && return PersistentArrayDict{K, V}(push!(m.kvs[1:end], Pair{K,V}(k, v)))

    kvs = m.kvs[1:end]
    kvs[idx] = Pair{K,V}(k, v)
    PersistentArrayDict{K, V}(kvs)
end

function dissoc(m::PersistentArrayDict{K, V}, k) where {K, V}
    idx = findkeyidx(m, k)
    idx === nothing && return m

    kvs = m.kvs[1:end]
    splice!(kvs, idx)
    PersistentArrayDict{K, V}(kvs)
end

function Base.iterate(m::PersistentArrayDict, i = 1)
    if i > length(m)
        return nothing
    else
        return (m.kvs[i], i + 1)
    end
end

Base.map(f::( Union{DataType, Function}), m::PersistentArrayDict) =
    PersistentArrayDict([f(kv) for kv in m]...)

Base.show(io::IO, ::MIME"text/plain", m::PersistentArrayDict{K, V}) where {K, V} =
    print(io, "Persistent{$K, $V}$(m.kvs)")

# Persistent Hash Maps
# ====================

struct PersistentHashDict{K, V} <: PersistentDict{K, V}
    trie::SparseBitmappedTrie{PersistentArrayDict{K, V}}
    length::Int
end
PersistentHashDict{K, V}() where {K, V} =
    PersistentHashDict{K, V}(SparseNode(PersistentArrayDict{K, V}), 0)

function PersistentHashDict(itr)
    if length(itr) == 0
        return PersistentHashDict()
    end
    K, V = typejoin(map(typeof, itr)...).types
    m = PersistentHashDict{K, V}()
    for (k, v) in itr
        m = assoc(m, k, v)
    end
    m
end

function PersistentHashDict(kvs::(Tuple{Any, Any})...)
    PersistentHashDict([kvs...])
end

function PersistentHashDict(kvs::(Pair)...)
    PersistentHashDict([kvs...])
end

function PersistentHashDict(; kwargs...)
    isempty(kwargs) ?
    PersistentHashDict{Any, Any}() :
    PersistentHashDict(kwargs...)
end

Base.length(m::PersistentHashDict) = m.length
Base.isempty(m::PersistentHashDict) = length(m) == 0

zipd(x,y) = map(p -> p[1] => p[2], zip(x,y))
Base.isequal(m1::PersistentHashDict, m2::PersistentHashDict) =
    length(m1) == length(m2) && all(x -> isequal(x...), zipd(m1, m2))

tup_eq(x) = x[1] == x[2]
==(m1::PersistentHashDict, m2::PersistentHashDict) =
    length(m1) == length(m2) && all(x -> x[1] == x[2], zipd(m1, m2))

function _update(f::Function, m::PersistentHashDict{K, V}, key) where {K, V}
    keyhash = reinterpret(Int, hash(key))
    arraymap = get(m.trie, keyhash, PersistentArrayDict{K, V}())
    newmap = f(arraymap)
    newtrie, _ = update(m.trie, keyhash, newmap)
    PersistentHashDict{K, V}(newtrie,
                            m.length + (length(newmap) < length(arraymap) ? -1 :
                                        length(newmap) > length(arraymap) ? 1 :
                                        0))
end

function assoc(m::PersistentHashDict{K, V}, key, value) where {K, V}
    _update(m, key) do arraymap
        assoc(arraymap, key, value)
    end
end

function dissoc(m::PersistentHashDict, key)
    _update(m, key) do arraymap
        dissoc(arraymap, key)
    end
end

function Base.getindex(m::PersistentHashDict, key)
    val = get(m.trie, reinterpret(Int, hash(key)), NotFound())
    (val === NotFound()) && error("key not found")
    val[key]
end

Base.get(m::PersistentHashDict, key) = m[key]
function Base.get(m::PersistentHashDict, key, default)
    val = get(m.trie, reinterpret(Int, hash(key)), NotFound())
    (val === NotFound()) && return default
    val[key]
end

function Base.haskey(m::PersistentHashDict, key)
    val = get(m.trie, reinterpret(Int, hash(key)), NotFound())
    (val != NotFound()) && haskey(val, key)
end

function Base.iterate(m::PersistentHashDict)
    trie_iter_result = iterate(m.trie)
    if trie_iter_result === nothing
        return nothing
    else
        arrmap, triestate = trie_iter_result
        kvs = arrmap.kvs
        return iterate(m, (kvs, triestate))
    end
end

function Base.iterate(m::PersistentHashDict, (kvs, triestate))
    if isempty(kvs) && isempty(triestate)
        return nothing
    else
        if isempty(kvs)
            arrmap, triestate = iterate(m.trie, triestate)
            return iterate(m, (arrmap.kvs, triestate))
        else
            return (kvs[1], (kvs[2:end], triestate))
        end
    end
end

function Base.map(f::( Union{Function, DataType}), m::PersistentHashDict)
    PersistentHashDict([f(kv) for kv in m]...)
end

function Base.filter(f::Function, m::PersistentHashDict{K, V}) where {K, V}
    arr = Array{Pair{K, V},1}()
    for el in m
        f(el) && push!(arr, el)
    end
    isempty(arr) ? PersistentHashDict{K, V}() : PersistentHashDict(arr...)
end

# Suppress ambiguity warning while allowing merging with array
function _merge(d::PersistentHashDict, others...)
    acc = d
    for other in others
        for (k, v) in other
            acc = assoc(acc, k, v)
        end
    end
    acc
end

# This definition suppresses ambiguity warning
Base.merge(d::PersistentHashDict, others::AbstractDict...) =
    _merge(d, others...)
Base.merge(d::PersistentHashDict, others...) =
    _merge(d, others...)

 function Base.show(io::IO, ::MIME"text/plain", m::PersistentHashDict{K, V}) where {K, V}
    print(io, "PHash{$K, $V}(")
    print(io, join(["$k => $v" for (k, v) in m], ", "))
    print(io, ")")
end

# Persistent Default Maps
# =======================
# The point of these "default" maps is this: it is generally fastest to use array maps for
# small/short maps and best to use hash-maps for large maps. Usually, we don't want to
# worry too much about tracking when things are large or small, though. Enter the default
# maps. Default maps are just a struct that wraps and imitates a PersistentArrayDict with
# one exception: when the size of an array map would grow to 64, a hash-map is instead
# returned, so maps get promoted automatically.

const _max_array_map_size = 64

struct PersistentAutoDict{K,V} <: PersistentDict{K,V}
    impl::PersistentArrayDict{K,V}
end

PersistentAutoDict{K, V}() where {K, V} =
    PersistentAutoDict{K, V}(PersistentArrayDict(Pair{K, V}[]))
PersistentAutoDict(arr::Array{Pair{K,V}}) where {K,V} =
    PersistentAutoDict{K,V}(PersistentArrayDict{K,V}(arr))
PersistentAutoDict(arr::Array{Tuple{K,V}}) where {K,V} =
    PersistentAutoDict{K,V}(PersistentArrayDict{K,V}([k=>v for (k,v) in arr]))
PersistentAutoDict(kvs::(Union{Tuple{K, V}, Pair{K, V}})...) where {K, V} =
    PersistentAutoDict{K,V}(PersistentArrayDict{K, V}(Pair{K, V}[Pair(k, v) for (k, v) in kvs]))
Base.isequal(m1::PersistentAutoDict, m2::PersistentAutoDict) = isequal(m1.impl, m2.impl)
==(m1::PersistentAutoDict, m2::PersistentAutoDict) = (m1.impl == m2.impl)

Base.length(m::PersistentAutoDict)  = length(m.impl)
Base.isempty(m::PersistentAutoDict) = isempty(m.impl)

findkeyidx(m::PersistentAutoDict, k) = findkeyidx(m.impl)

Base.get(m::PersistentAutoDict, k) = get(m.impl, k)
Base.get(m::PersistentAutoDict, k, default) = get(m.impl, k, default)
Base.getindex(m::PersistentAutoDict, k) = getindex(m.impl, k)

Base.haskey(m::PersistentAutoDict, k) = haskey(m.impl, k)

function FunctionalCollections.assoc(m::PersistentAutoDict{K, V}, k, v) where {K, V}
    if length(m.impl) + 1 >= _max_array_map_size && !haskey(m.impl, k)
        assoc(PersistentHashDict(collect(m)), k, v)
    else
        q = assoc(m.impl, k, v)
        (q === m.impl) ? m : PersistentAutoDict{K,V}(q)
    end
end

function FunctionalCollections.dissoc(m::PersistentAutoDict{K, V}, k) where {K, V}
    q = dissoc(m.impl, k)
    (q === m.impl) ? m : PersistentAutoDict{K,V}(q)
end

function Base.iterate(m::PersistentAutoDict, i = 1)
    iterate(m.impl, i)
end

Base.map(f::( Union{DataType, Function}), m::PersistentAutoDict) = PersistentAutoDict(map(f, m.impl))

Base.show(io::IO, q::MIME"text/plain", m::PersistentAutoDict{K, V}) where {K, V} =
    print(io, "PDict{$K, $V}(")
    print(io, join(["$k => $v" for (k, v) in m.impl], ", "))
    print(io, ")")

max_issues_repo_name: noahbenson/FunctionalUtilities

max_forks_repo_name: noahbenson/FunctionalUtilities


text:
using Makie

res = 500
game = Scene(resolution=(res, res))
campixel!(game)
cell = 10
cells = res ÷ 10
snakelen = 5
middle = cells ÷ 2
snakestart = middle - (snakelen ÷ 2)
segments = Point2f0.(middle, range(snakestart, step=cell, length=snakelen))
snake = scatter!(game, segments, markersize=cell, raw=true, color=:blue)[end]
speed = 1 / 5
global last_dir = (0, 1)
dir = lift(game.events.keyboardbuttons) do but
    global last_dir
    ispressed(but, Keyboard.left) && return (last_dir = (-1, 0))
    ispressed(but, Keyboard.up) && return (last_dir = (0, 1))
    ispressed(but, Keyboard.right) && return (last_dir = (1, 0))
    ispressed(but, Keyboard.down) && return (last_dir = (0, -1))
    last_dir
end
display(game)
newsnake = copy(segments)
while isopen(game)
    curr_snake = snake[1][]
    circshift!(newsnake, curr_snake, 1)
    newsnake[1] = newsnake[2] .+ Point2f0(dir[] .* cell)
    snake[1] = newsnake
    newsnake = curr_snake
    sleep(speed)
end

ext: jl


text:
module JLBoost

#using DataFrames
using SortingLab
using StatsBase: sample
using Base.Iterators: drop
using LossFunctions: LogitProbLoss, deriv, deriv2, SupervisedLoss
#using Zygote: gradient, hessian
#using ForwardDiff: gradient, hessian
#using Flux: logitcrossentropy, logitbinarycrossentropy
#using RCall

export jlboost, find_best_split, _find_best_split, predict, fit_tree, logloss, jlboost!
export JLBoostTree, show, *, print, println
export LogitLogloss, value, deriv, deriv2, trees
export JLBoostTreeModel, JLBoostTree, WeightedJLBoostTree, features, feature_importance, vcat
export getproperty, AbstractJLBoostTree, predict
export max_depth_stopping_criterion, max_leaves_stopping_criterion

include("JLBoostTrees/JLBoostTrees.jl");
using ..JLBoostTrees: JLBoostTree, AbstractJLBoostTree, WeightedJLBoostTree,
    JLBoostTreeModel, trees, vcat, getproperty, get_leaf_nodes, treedepth, keeprow_vec,
    is_left_child, is_right_child

include("tree-growth.jl")

include("column_sampling_strategies.jl")
include("row_sampling_strategies.jl")

include("diagnostics.jl")
include("g_h.jl")
include("find_best_split.jl")

include("fit_tree_stopping_criterion.jl")
using ..TreeFitStoppingCriterion: max_depth_stopping_criterion, max_leaves_stopping_criterion

include("fit_tree.jl")

include("predict.jl")
include("save.jl")
include("jlboost-fit.jl")
include("get-features.jl")
include("feature-importance.jl")
include("tables.jl")

end # module


text:
using LinearAlgebra: BlasInt, BlasFloat, LU, UnitLowerTriangular, ldiv!, BLAS, checknonsingular

function lu(A::AbstractMatrix, pivot::Union{Val{false}, Val{true}} = Val(true);
            check::Bool = true, blocksize::Integer = 16)
    lu!(copy(A), pivot; check = check, blocksize = blocksize)
end

function lu!(A, pivot::Union{Val{false}, Val{true}} = Val(true);
             check::Bool = true, blocksize::Integer = 16)
    lu!(A, Vector{BlasInt}(undef, min(size(A)...)), pivot;
        check = check, blocksize = blocksize)
end

function lu!(A::AbstractMatrix{T}, ipiv::AbstractVector{<:Integer},
             pivot::Union{Val{false}, Val{true}} = Val(true);
             check::Bool=true, blocksize::Integer=16) where T
    info = Ref(zero(BlasInt))
    m, n = size(A)
    mnmin = min(m, n)
    if T <: BlasFloat && A isa StridedArray
        reckernel!(A, pivot, m, mnmin, ipiv, info, blocksize)
        if m < n # fat matrix
            # [AL AR]
            AL = @view A[:, 1:m]
            AR = @view A[:, m+1:n]
            apply_permutation!(ipiv, AR)
            ldiv!(UnitLowerTriangular(AL), AR)
        end
      else # generic fallback
        _generic_lufact!(A, pivot, ipiv, info)
    end
    check && checknonsingular(info[])
    LU{T, typeof(A)}(A, ipiv, info[])
end

function nsplit(::Type{T}, n) where T
    k = 128 ÷ sizeof(T)
    k_2 = k ÷ 2
    return n >= k ? ((n + k_2) ÷ k) * k_2 : n ÷ 2
end

Base.@propagate_inbounds function apply_permutation!(P, A)
    for i in axes(P, 1)
        i′ = P[i]
        i′ == i && continue
        @simd for j in axes(A, 2)
            A[i, j], A[i′, j] = A[i′, j], A[i, j]
        end
    end
    nothing
end

function reckernel!(A::AbstractMatrix{T}, pivot::Val{Pivot}, m, n, ipiv, info, blocksize)::Nothing where {T,Pivot}
    @inbounds begin
        if n <= max(blocksize, 1)
            _generic_lufact!(A, pivot, ipiv, info)
            return nothing
        end
        n1 = nsplit(T, n)
        n2 = n - n1
        m2 = m - n1

        # ======================================== #
        # Now, our LU process looks like this
        # [ P1 ] [ A11 A21 ]   [ L11 0 ] [ U11 U12  ]
        # [    ] [         ] = [       ] [          ]
        # [ P2 ] [ A21 A22 ]   [ L21 I ] [ 0   A′22 ]
        # ======================================== #

        # ======================================== #
        # Partition the matrix A
        # [AL AR]
        AL = @view A[:, 1:n1]
        AR = @view A[:, n1+1:n]
        #  AL  AR
        # [A11 A12]
        # [A21 A22]
        A11 = @view A[1:n1, 1:n1]
        A12 = @view A[1:n1, n1+1:n]
        A21 = @view A[n1+1:m, 1:n1]
        A22 = @view A[n1+1:m, n1+1:n]
        # [P1]
        # [P2]
        P1 = @view ipiv[1:n1]
        P2 = @view ipiv[n1+1:n]
        # ========================================

        #   [ A11 ]   [ L11 ]
        # P [     ] = [     ] U11
        #   [ A21 ]   [ L21 ]
        reckernel!(AL, pivot, m, n1, P1, info, blocksize)
        # [ A12 ]    [ P1 ] [ A12 ]
        # [     ] <- [    ] [     ]
        # [ A22 ]    [ 0  ] [ A22 ]
        Pivot && apply_permutation!(P1, AR)
        # A12 = L11 U12  =>  U12 = L11 \ A12
        ldiv!(UnitLowerTriangular(A11), A12)
        # Schur complement:
        # We have A22 = L21 U12 + A′22, hence
        # A′22 = A22 - L21 U12
        BLAS.gemm!('N', 'N', -one(T), A21, A12, one(T), A22)
        # record info
        previnfo = info[]
        # P2 A22 = L22 U22
        reckernel!(A22, pivot, m2, n2, P2, info, blocksize)
        # A21 <- P2 A21
        Pivot && apply_permutation!(P2, A21)

        info[] != previnfo && (info[] += n1)
        @simd for i in 1:n2
            P2[i] += n1
        end
        return nothing
    end # inbounds
end

#=
    Modified from https://github.com/JuliaLang/julia/blob/b56a9f07948255dfbe804eef25bdbada06ec2a57/stdlib/LinearAlgebra/src/lu.jl
    License is MIT: https://julialang.org/license
=#
function _generic_lufact!(A, ::Val{Pivot}, ipiv, info) where Pivot
    m, n = size(A)
    minmn = length(ipiv)
    @inbounds begin
        for k = 1:minmn
            # find index max
            kp = k
            if Pivot
              amax = abs(zero(eltype(A)))
              for i = k:m
                  absi = abs(A[i,k])
                  if absi > amax
                      kp = i
                      amax = absi
                  end
              end
            end
            ipiv[k] = kp
            if !iszero(A[kp,k])
                if k != kp
                    # Interchange
                    @simd for i = 1:n
                        tmp = A[k,i]
                        A[k,i] = A[kp,i]
                        A[kp,i] = tmp
                    end
                end
                # Scale first column
                Akkinv = inv(A[k,k])
                @simd for i = k+1:m
                    A[i,k] *= Akkinv
                end
            elseif info[] == 0
                info[] = k
            end
            # Update the rest
            for j = k+1:n
                @simd for i = k+1:m
                    A[i,j] -= A[i,k]*A[k,j]
                end
            end
        end
    end
    return nothing
end


text:
module TestJSON
import JSON

mutable struct S
	a
	b
end

mutable struct T
	a
	b
end

s = S(1,T(10, 25))

JSON.print(s)

n = (a=1, b=(a=10, b=25))

JSON.print(n)

@show JSON.print(s) == JSON.print(n)

end


text:
# Experimental threads example. WIP!

using StanSample
using Statistics
using Distributions

ProjDir = @__DIR__
cd(ProjDir) #do

bernoullimodel = "
data {
  int<lower=1> N;
  int<lower=0,upper=1> y[N];
}
parameters {
  real<lower=0,upper=1> theta;
}
model {
  theta ~ beta(1,1);
  y ~ bernoulli(theta);
}
";

#isdir(ProjDir * "/tmp") && rm(ProjDir * "/tmp", recursive=true)
tmpdir = ProjDir * "/tmp"

p1 = 12 # p1 is the number of models to fit
p = range(0, stop=1, length=p1)

n = 100;
observeddata = [Dict("N" => n, "y" => rand(Bernoulli(p[i]), n)) for i in 1:p1]

sm = [SampleModel("bernoulli_m$i", bernoullimodel; tmpdir=tmpdir) for i in 1:p1];

println("\nThreads loop\n")

estimates = Vector(undef, p1)
Threads.@threads for i in 1:p1
    rc = stan_sample(sm[i]; data=observeddata[i]);
    if success(rc)
      samples = read_samples(sm[i]; output_format=:namedtuple)

      estimates[i] = [mean(reshape(samples.theta, 4000)), std(reshape(samples.theta, 4000))]
    end
end

estimates |> display

#end


text:
#Define the absolute Gini Index, see [4] for a complete definition
function GI(x::Array)

    @assert eltype(x)<:Number

    #Reshape the incoming matrix into a sorted 1xn list of positive numbers
    x=abs(x[:])
    sort!(x)

    n=length(x)

    #s is the sum of all terms, g is the gini index
    g=0
    s=0
    for i = 1:n
        s += x[i]
        g += x[i]*i
    end

    g *= 2/n/s
    g -= (n+1)/n

    #return G
    g
end;

max_stars_repo_name: dahlend/CompressedSensing

max_issues_repo_name: dahlend/CompressedSensing

max_forks_repo_name: dahlend/CompressedSensing


text:
using LabVIEW0

try
    include(get_script_path(raw"C:\_LabView_projects\ZMQ\LabVIEW0.jl\src\LabVIEW\jl-scripts\simple_examples.jl"))
    LabVIEW0.setglobals(isOK=true)
catch excep
    LabVIEW0.setglobals(isOK=false, excpn=(;excep, stack_trace=backtrace()))
end

server_0mq4lv(fns)


text:
### Special constructors for HSS matrices
#
# Written by Boris Bonev, Feb. 2020

function lowrank2hss(U::Matrix{T}, V::Matrix{T}, rcl::ClusterTree, ccl::ClusterTree) where T
  (k = size(U,2)) == size(V,2) || throw(ArgumentError("second dimension of U and V must agree"))
  return _lowrank2hss(U, V, rcl, ccl, k;  rootnode=true)
end

function _lowrank2hss(U::Matrix{T}, V::Matrix{T}, rcl::ClusterTree, ccl::ClusterTree, k::Int; rootnode=false) where T
  if isleaf(rcl) && isleaf(ccl)
    HssMatrix(U[rcl.data,:]*V[ccl.data,:]', U[rcl.data,:], V[ccl.data,:], rootnode)
  elseif isbranch(rcl) && isbranch(ccl)
    A11 = _lowrank2hss(U, V, rcl.left, ccl.left, k)
    A22 = _lowrank2hss(U, V, rcl.right, ccl.right, k)
    HssMatrix(A11, A22, Matrix{T}(I,k,k), Matrix{T}(I,k,k), Matrix{T}(I,k,k), Matrix{T}(I,k,k), Matrix{T}(I,k,k), Matrix{T}(I,k,k), rootnode)
  else
    throw(ArgumentError("row and column clusters are not compatible"))
  end
end


text:
_themes[:cooper] = PlotTheme(linewidth = 4,
                             markersize = 10,
                             fontfamily = "sans-serif",
                             titlefontsize = 24,
                             guidefontsize = 18,
                             tickfontsize = 18,
                             colorbar_tickfontsize = 18,
                             legend_font_pointsize = 18,
                             legend=:outertopright,
#                             legend_font = Plots.Font("sans-serif",
#                                                      18, :hcenter, :vcenter,
#                                                      0.0, RGB{N0f8}(0.0,0.0,0.0)),
                             grid = false
                            )


text:
module Wigxjpf

using Wigxjpf_jll

export wig_table_init, wig_temp_init, wig_thread_temp_init, wig3jj, wig6jj, wig9jj, wig_table_free, wig_temp_free

@doc raw"""
```
wig_table_init(max_two_j::Integer, wigner_type::Integer)
```

Initialize the calculation table. Must be called before evaluating any Wigner symbols.

`max_two_j` is twice the highest absolute value of all `j` values to be evaluated.

`wigner_type` should be 3, 6, or 9. When multiple types are to be used, you should use the highest value for this parameter.

> When used in a multi-threaded environment, this function should be called **globally**.
"""
function wig_table_init(max_two_j::Integer, wigner_type::Integer)
    @assert max_two_j >= 0
    @assert wigner_type ∈ [3, 6, 9]

    ccall((:wig_table_init, wigxjpf), Cvoid, (Int32, Int32), convert(Int32, max_two_j), convert(Int32, wigner_type))
end

@doc raw"""
```
wig_temp_init(max_two_j::Integer)
```

Initialize the temporary array. Must be called before evaluating any Wigner symbols.

> When used in a multi-threaded environment, call `wig_thread_temp_init` instead.
"""
function wig_temp_init(max_two_j::Integer)
    @assert max_two_j >= 0

    ccall((:wig_temp_init, wigxjpf), Cvoid, (Int32,), convert(Int32, max_two_j))
end

@doc raw"""
```
wig_thread_temp_init(max_two_j::Integer)
```

Initialize the temporary array. Must be called before evaluating any Wigner symbols.

> When used in a multi-threaded environment, this function should be called in each thread.
"""
function wig_thread_temp_init(max_two_j::Integer)
    @assert max_two_j >= 0

    ccall((:wig_thread_temp_init, wigxjpf), Cvoid, (Int32,), convert(Int32, max_two_j))
end

@doc raw"""
```
wig3jj(two_j1::Integer, two_j2::Integer, two_j3::Integer, two_m1::Integer, two_m2::Integer, two_m3::Integer)
```

Evaluate Wigner 3-j symbols. Note that all parameter are twice the actual j or m parameters to allow representation of half-integers.
"""
function wig3jj(two_j1::Integer, two_j2::Integer, two_j3::Integer, two_m1::Integer, two_m2::Integer, two_m3::Integer)
    return ccall((:wig3jj, wigxjpf), Cdouble, (Int32, Int32, Int32, Int32, Int32, Int32), convert(Int32, two_j1), convert(Int32, two_j2), convert(Int32, two_j3), convert(Int32, two_m1), convert(Int32, two_m2), convert(Int32, two_m3))
end

@doc raw"""
```
wig6jj(two_j1::Integer, two_j2::Integer, two_j3::Integer, two_j4::Integer, two_j5::Integer, two_j6::Integer)
```

Evaluate Wigner 6-j symbols. Note that all parameter are twice the actual j parameters to allow representation of half-integers.
"""
function wig6jj(two_j1::Integer, two_j2::Integer, two_j3::Integer, two_j4::Integer, two_j5::Integer, two_j6::Integer)
    return ccall((:wig6jj, wigxjpf), Cdouble, (Int32, Int32, Int32, Int32, Int32, Int32), convert(Int32, two_j1), convert(Int32, two_j2), convert(Int32, two_j3), convert(Int32, two_j4), convert(Int32, two_j5), convert(Int32, two_j6))
end

@doc raw"""
```
wig9jj(two_j1::Integer, two_j2::Integer, two_j3::Integer, two_j4::Integer, two_j5::Integer, two_j6::Integer, two_j7::Integer, two_j8::Integer, two_j9::Integer)
```

Evaluate Wigner 9-j symbols. Note that all parameter are twice the actual j parameters to allow representation of half-integers.
"""
function wig9jj(two_j1::Integer, two_j2::Integer, two_j3::Integer, two_j4::Integer, two_j5::Integer, two_j6::Integer, two_j7::Integer, two_j8::Integer, two_j9::Integer)
    return ccall((:wig9jj, wigxjpf), Cdouble, (Int32, Int32, Int32, Int32, Int32, Int32, Int32, Int32, Int32,), convert(Int32, two_j1), convert(Int32, two_j2), convert(Int32, two_j3), convert(Int32, two_j4), convert(Int32, two_j5), convert(Int32, two_j6), convert(Int32, two_j7), convert(Int32, two_j8), convert(Int32, two_j9))
end

@doc raw"""
```
wig_temp_free()
```

Free the temporary array.

> When used in a multi-threaded environment, this function should be called in each thread.
"""
function wig_temp_free()
    ccall((:wig_temp_free, wigxjpf), Cvoid, ())
end

@doc raw"""
```
wig_table_free()
```

Free the calculation table.

> When used in a multi-threaded environment, this function should be called **globally**.
"""
function wig_table_free()
    ccall((:wig_table_free, wigxjpf), Cvoid, ())
end

end


text:
function mul_revDR(a::Interval,b::Interval,c::Interval) # DONE
    if (in(0,c) && ~in(0,b))
        c = c ∩ (a / b)
    elseif (~in(0,c) && in(0,b))
        b = b ∩ (a / c)
    elseif (~in(0,c) && ~in(0,b))
        b = b ∩ (a / c)
        c = c ∩ (a / b)
    end
    return a, b, c
end
mul_revDR(a,b,c) = div_revDR(promote(a,b,c)...)

function exp_revDR(y::Interval, x::Interval) # DONE
    y_new = y ∩ (0..∞)
    x_new = x ∩ log(y_new)

    return y_new, x_new
end
exp_revDR(a,b) = exp_revDR(promote(a,b)...)

function acos_rev(y::Interval, x::Interval) # DONE
    y_new = y ∩ (0.0..2*pi)
    x_new = x ∩ cos(y_new)

    return y_new, x_new
end
acos_rev(a,b) = acos_rev(promote(a,b)...)

function atan_rev(y::Interval, x::Interval) # DONE
    y_new = y ∩ (-pi/2.0..pi/2.0)
    x_new = x ∩ tan(y_new)

    return y_new, x_new
end
atan_rev(a,b) = atan_rev(promote(a,b)...)

function sinh_rev(y::Interval,x::Interval) # DONE
    x_new = x ∩ asinh(y_new)

    return y, x_new
end
sinh_rev(a,b) = sinh_rev(promote(a,b)...)

function cosh_rev(y::Interval,x::Interval) # DONE
    y_new = y ∩ (1.0..∞)
    x_new = x ∩ acosh(y_new)

    return y_new, x_new
end
cosh_rev(a,b) = cosh_rev(promote(a,b)...)

function tanh_rev(y::Interval,x::Interval) # DONE
    y_new = y ∩ ((-1.0)..(1.0))
    x_new = x ∩ atanh(y_new)

    return y_new, x_new
end
tanh_rev(a,b) = tanh_rev(promote(a,b)...)

function asinh_rev(y::Interval,x::Interval) # DONE
    x_new = x ∩ sinh(y_new)

    return y, x_new
end
asinh_rev(a,b) = asinh_rev(promote(a,b)...)

function acosh_rev(y::Interval,x::Interval) # DONE
    y_new = y ∩ (0.0..∞)
    x_new = x ∩ cosh(y_new)

    return y_new, x_new
end
acosh_rev(a,b) = acosh_rev(promote(a,b)...)

function atanh_rev(y::Interval,x::Interval) # DONE
    x_new = x ∩ tanh(y_new)

    return y, x_new
end
atanh_rev(a,b) = atanh_rev(promote(a,b)...)

function step_rev(y::Interval,x::Interval) # DONE
    y_new = y ∩ (0..1)

    return y_new, x
end
step_rev(a,b) = step_rev(promote(a,b)...)

function sign_rev(y::Interval,x::Interval) # DONE
    y_new = y ∩ ((-1.0)..(1.0))

    return y_new, x
end
sign_rev(a,b) = sign_rev(promote(a,b)...)

function exp2_rev(y::Interval,x::Interval) # DONE
    y_new = y ∩ (0..∞)
    x_new = x ∩ log10(y_new)

    return y_new, x_new
end
exp2_rev(a,b) = exp2_rev(promote(a,b)...)

function exp10_rev(y::Interval,x::Interval) # DONE
    y_new = y ∩ (0..∞)
    x_new = x ∩ log10(y_new)

    return y_new, x_new
end
exp10_rev(a,b) = exp10_rev(promote(a,b)...)

function log2_rev(y::Interval,x::Interval) # DONE
    y_new = y ∩ (0..∞)
    x_new = x ∩ exp2(y_new)

    return y_new, x_new
end
log2_rev(a,b) = log2_rev(promote(a,b)...)

function log10_rev(y::Interval,x::Interval) # DONE
    y_new = y ∩ (0..∞)
    x_new = x ∩ exp10(y_new)

    return y_new, x_new
end
log10_rev(a,b) = log10_rev(promote(a,b)...)

function one_rev(y::Interval,x::Interval) # DONE
    y_new = y ∩ ((1.0)..(1.0))

    return y_new, x
end
one_rev(a,b) = one_rev(promote(a,b)...)

function zero_rev(y::Interval,x::Interval) # DONE
    y_new = y ∩ ((0.0)..(0.0))

    return y_new, x
end
zero_rev(a,b) = zero_rev(promote(a,b)...)

function min_rev(z::Interval,x::Interval,y::Interval)
    y_new = y ∩ (0..∞)
    x_new = x ∩ log(y_new)

    return y_new, x_new
end
function max_rev(z::Interval,x::Interval,y::Interval)
    y_new = y ∩ (0..∞)
    x_new = x ∩ log(y_new)

    return y_new, x_new
end


text:
# A hopefully correct implementation of "non-blocking unordered single-reader channel."
#
# Important requirements:
# * Lock-free (in fact wait-free) and yield-free on the writer side (so that it can be used
#   inside the finalizer).
# * Reader can wait cooperatively (so that this can be used in a "background" task).
# * Implementable without atomics on GC-manged Julia object fields (for Julia < 1.7).
#
# The implementation can be simplified because of:
# * The insertion order is irrelevant.
# * Single-reader.

baremodule WorkerStates
using Base: UInt8
const Kind = UInt8
const ACTIVE = Kind(0)
const WAITING = Kind(1)
const NOTIFYING = Kind(2)
end  # baremodule WorkerStates

"""
    WriterState{T}

Let `w = state.writable[]` where `state::WriterState`.  The lowest bit ("index bit") of `w`
indicates 0-origin index of the writable buffer; i.e.,

    i = (w & 0b1) + 1
    state.buffers[i]

is writable (if locked; see blow). The other buffer, i.e.,

    j = ~(w & 0b1) + 1
    state.buffers[j]

is readable.  If the second bit ("lock bit") of `w` is set, `state.buffers[i]` is considered
locked.

The writer obtains the exclusive access by setting the lock bit while the reader obtains the
exclusive access by flipping the index bit.

This is loosely inspired by the left-right technique.
"""
struct WriterState{T}
    writable::Threads.Atomic{UInt32}
    buffers::NTuple{2,Vector{T}}
end

WriterState{T}() where {T} = WriterState{T}(Threads.Atomic{UInt32}(), (T[], T[]))

"""
    SingleReaderDualBag{T}()

"Non-blocking unordered single-reader channel."
"""
mutable struct SingleReaderDualBag{T}
    # TODO: Use a proper thread-local storage for `states`
    @const states::Vector{WriterState{T}}
    @const readerstate::Threads.Atomic{WorkerStates.Kind}
    waiter::Union{Task,Nothing}
end

function SingleReaderDualBag{T}() where {T}
    states = [WriterState{T}() for _ in 1:Threads.nthreads()]
    readerstate = Threads.Atomic{WorkerStates.Kind}(WorkerStates.ACTIVE)
    return SingleReaderDualBag{T}(states, readerstate, nothing)
end

Base.eltype(::Type{SingleReaderDualBag{T}}) where {T} = T

const LOCKBITMASK = UInt32(0b10)

"""
    Base.put!(bag::SingleReaderDualBag, x)

Insert `x` to `bag`.  Notify the reader if required.  This method is wait-free.  It also has
no yield points provided that `convert(eltype(bag), x)` does not yield.
"""
function Base.put!(bag::SingleReaderDualBag{T}, x) where {T}
    @record(:put_begin, bag, x, {yield = false})
    x = convert(T, x)

    @yield_unsafe let state, writable  # [^crucially_yield_unsafe]
        state = bag.states[Threads.threadid()]

        # Acquire a writable buffer
        writable = atomic_or!(state.writable, LOCKBITMASK)
        @assert iszero(writable & LOCKBITMASK)
        @assert writable in (0, 1)

        push!(state.buffers[writable+1], x)

        # Release the buffer
        state.writable[] = writable
    end

    atomic_fence()  # [^store_buffering_1]
    if bag.readerstate[] == WorkerStates.WAITING
        old = atomic_cas!(bag.readerstate, WorkerStates.WAITING, WorkerStates.NOTIFYING)
        if old == WorkerStates.WAITING
            @record(:put_schedule, bag, x, {yield = false})
            schedule(bag.waiter::Task)
        end
    end

    @record(:put_end, bag, x, {yield = false})
    return bag
end
# [^crucially_yield_unsafe]: According to the docstring, everything in this function except
# the `convert` call is yield-unsafe.  However, the explicit `@yield_unsafe` block is "much
# more" yield-unsafe since the yield must not exist inside this block for correctness even
# without the specification in the docstring.

function trytakemanyto!(output, bag::SingleReaderDualBag)
    took = false
    locked = false
    for state in bag.states
        # Acquire a readable buffer
        writable = state.writable[]
        local readable::Int = -1
        while true
            if iszero(writable & LOCKBITMASK)
                new = (writable + one(writable)) & one(writable)
                old = atomic_cas!(state.writable, writable, new)
                if old == writable
                    readable = writable
                    @goto acquired
                end
            else
                locked = true
                break
            end
        end
        continue
        @label acquired
        @assert readable in (0, 1)

        buffer = state.buffers[readable+1]
        isempty(buffer) && continue
        append!(output, buffer)
        empty!(buffer)
        took = true
    end
    return (; took = took, locked = locked)
end

"""
    takemanyto!(output, bag::SingleReaderDualBag)

Take one or many items from `bag` and `append!` them into `output`.

The items that has been `put!` into `bag` are taken out from the `bag` if it is non-empty.
Otherwise, it waits until at least one item is available.

This method is lock-free if `bag` is not empty.
"""
function takemanyto!(output, bag::SingleReaderDualBag)
    @record(:takemanyto_begin, bag)
    chaos_takemanyto()
    result = trytakemanyto!(output, bag)
    if result.took
        @record(:takemanyto_end, bag, waited = false)
        return result
    end

    bag.waiter = current_task()
    while true
        @record(:takemanyto_wait_trying, bag)
        @yield_unsafe begin
            bag.readerstate[] = WorkerStates.WAITING
            atomic_fence()  # [^store_buffering_1]
            result = trytakemanyto!(output, bag)
            if result.took
                old =
                    atomic_cas!(bag.readerstate, WorkerStates.WAITING, WorkerStates.ACTIVE)
                if old == WorkerStates.WAITING
                    # It is yield-safe once the state transition is successfully cancelled.
                    @record(:takemanyto_wait_cancelled, bag, success = true)
                    break
                end
                @record(:takemanyto_wait_cancelled, bag, success = false, {yield = false})
                # Otherwise, it means that the reader (this task) has failed to cancel the
                # state transition.  It has to "receive" the `schedule` by calling `wait()`
                # and then return immediately. [^failed_to_cancel].
            end
            @record(:takemanyto_wait_begin, bag, {yield = false})
        end
        wait()
        @record(:takemanyto_wait_end, bag)
        @assert bag.readerstate[] == WorkerStates.NOTIFYING
        bag.readerstate[] = WorkerStates.ACTIVE
        result.took && break  # [^failed_to_cancel]: if it had failed to cancel

        result = trytakemanyto!(output, bag)
        result.took && break
    end
    bag.waiter = nothing
    @record(:takemanyto_end, bag, waited = true)
    return result
end
# [^store_buffering_1]: These fences are used to avoid the cycle 2b -> 1a -> 1b -> 2a -> 2b
# where
# * Reader:
#   * 1a: `bag.readerstate[] = WorkerStates.WAITING`
#   * 1b: `trytakemanyto!(output, bag).took` is `false` (acquiring of state locks)
# * Writer:
#   * 2a: state lock release in `put!(bag, x)`
#   * 2b: `bag.readerstate[]` in `put!` is not `WAITING`
# i.e., the reader misses the insertion and writer misses the sleep state transition.
#
# The edges of the cycle are:
# * 2b -> 1a: reads-before, due to the value read in 2b
# * 1a -> 1b: sequenced-before
# * 1b -> 2a: reads-before, due to the value read in 1b
# * 2a -> 2b: sequenced-before
#
# Ref:
# https://github.com/JuliaLang/julia/pull/43418#discussion_r790206434

function takemany!(bag::SingleReaderDualBag{T}) where {T}
    output = T[]
    takemanyto!(output, bag)
    return output
end


text:
# contains the information required for running a task
struct GWTask
    tid::String
    task_dir::String

    # Sys dat
    dat::Dict{Any, Any}

    GWTask(tid, task_dir) = new(tid, task_dir, Dict{Any, Any}())
    GWTask(;tid, task_dir) = GWTask(tid, task_dir)

end

# overwrite base
import Base.show

function Base.show(io::IO, gwt::GWTask)
    println(io, "GWTask(;")
    println(io, "   task_id = \"", task_id(gwt), "\",")
    println(io, "   task_dir = \"", task_dir(gwt), "\"")
    println(io, ")")
end

import Base.get!
get!(f::Function, gwt::GWTask, key) = get!(f, gwt.dat, key)
get!(gwt::GWTask, key, val) = get!(gwt.dat, key, val)

import Base.get
get(f::Function, gwt::GWTask, key) = get(f, gwt.dat, key)
get(gwt::GWTask, key, val) = get(gwt.dat, key, val)

set!(gwt::GWTask, key, val) = (gwt.dat[key] = val)
set!(f::Function, gwt::GWTask, key) = (gwt.dat[key] = f())

task_id(gwt::GWTask) = gwt.tid
task_dir(gwt::GWTask) = gwt.task_dir

const _GIT_WORKER_KEY = :_gitworker
function gitworker(gwt::GWTask)
    get(gwt, _GIT_WORKER_KEY) do
        gwfile = _find_worker_file(task_dir(gwt))
        gw = _gw_from_toml(gwfile)
        isnothing(gw) && error("GitWorker not found!")
        return gw
    end
end

remote_url(gwt::GWTask) = remote_url(gitworker(gwt))
sys_root(gwt::GWTask) = sys_root(gitworker(gwt))

import Base.lock
lock(f::Function, gwt::GWTask; kwargs...) = lock(f, gitworker(gw); kwargs...)

max_stars_repo_name: josePereiro/GitWorkers

max_issues_repo_name: josePereiro/GitWorkers

max_forks_repo_name: josePereiro/GitWorkers


text:
using SimpleSolvers
using Test

n = 1
T = Float64
x = [T(π),]
j = reshape(2x, 1, 1)

function F!(f::Vector, x::Vector)
    f .= x.^2
end

function J!(g::Matrix, x::Vector)
    g .= 0
    for i in eachindex(x)
        g[i,i] = 2x[i]
    end
end

JPAD = JacobianParameters{T}(F!, n; mode = :autodiff, diff_type = :forward)
JPFD = JacobianParameters{T}(F!, n; mode = :autodiff, diff_type = :finite)
JPUS = JacobianParameters{T}(J!, n; mode = :user)

@test typeof(JPAD) <: JacobianParametersAD
@test typeof(JPFD) <: JacobianParametersFD
@test typeof(JPUS) <: JacobianParametersUser

function test_jac(j1, j2, atol)
    for i in eachindex(j1,j2)
        @test j1[i] ≈ j2[i] atol=atol
    end
end

jad = zero(j)
jfd = zero(j)
jus = zero(j)

compute_jacobian!(jad, x, JPAD)
compute_jacobian!(jfd, x, JPFD)
compute_jacobian!(jus, x, JPUS)

test_jac(jad, j, eps())
test_jac(jfd, j, 1E-7)
test_jac(jus, j, 0)

jad1 = zero(j)
jfd1 = zero(j)
jus1 = zero(j)

compute_jacobian!(jad1, x, F!; mode = :autodiff, diff_type = :forward)
compute_jacobian!(jfd1, x, F!; mode = :autodiff, diff_type = :finite)
compute_jacobian!(jus1, x, J!; mode = :user)

test_jac(jad, jad1, 0)
test_jac(jfd, jfd1, 0)
test_jac(jus, jus1, 0)

jad2 = zero(j)
jfd2 = zero(j)

compute_jacobian_ad!(jad2, x, F!)
compute_jacobian_fd!(jfd2, x, F!)

test_jac(jad, jad2, 0)
test_jac(jfd, jfd2, 0)


text:
using DataStructures: OrderedDict

abstract type Tabulation{V} end

value(t::Tabulation) = getfield(t, :value)
Base.adjoint(t::Tabulation) = value(t)
Base.getproperty(t::Tabulation, k::Symbol) = value(t)[k]
Base.getindex(t::Tabulation, k::Symbol) = getproperty(t, k)
Base.iterate(t::Tabulation) = iterate(value(t))
Base.iterate(t::Tabulation, i) = iterate(value(t), i)

struct TabulationCols{V} <: Tabulation{V}
    value::OrderedDict{Symbol,V}
    TabulationCols{V}(x...) where V = new{V}(OrderedDict(x...))
end

struct TabulationRows{V} <: Tabulation{V}
    value::OrderedDict{Symbol,TabulationCols{V}}
    TabulationRows{V}(x...) where V = new{V}(OrderedDict(x...))
end

import DataFrames: DataFrames, DataFrame
DataFrames.DataFrame(t::TabulationCols) = DataFrame(value(t))
DataFrames.DataFrame(t::TabulationRows; index=true) = begin
    v = value(t)
    I = DataFrame("" => collect(keys(v)))
    C = DataFrame(value.(values(v)))
    index ? [I C] : C
end
Base.Matrix(t::TabulationRows) = Matrix(DataFrame(t; index=false))

Base.show(io::IO, t::TabulationCols) = show(io, DataFrame(t); summary=false, eltypes=false, show_row_number=false, vlines=:none)
Base.show(io::IO, t::TabulationRows) = show(io, DataFrame(t); summary=false, eltypes=false, show_row_number=false, vlines=[1])

tabulation(m, R, C, V) = TabulationRows{V}(R .=> [TabulationCols{V}(zip(C, V.(m[i,:]))) for i in 1:size(m, 1)])

struct Tabulate{V} <: State{V}
    value::TabulationRows{V}
    rows::Tuple{Vararg{Symbol}}
    columns::Tuple{Vararg{Symbol}}
end

Base.getproperty(t::Tabulate, k::Symbol) = t[k]
Base.show(io::IO, t::Tabulate) = show(io, Matrix(value(t)))
Base.show(io::IO, ::MIME"text/plain", t::Tabulate) = show(io, value(t))

Tabulate(; unit, rows, columns=(), _value, _type, _...) = begin
    U = value(unit)
    V = valuetype(_type, U)
    columns = isempty(columns) ? rows : columns
    v = tabulation(_value, rows, columns, V)
    Tabulate{V}(v, rows, columns)
end

constructortags(::Val{:Tabulate}) = (:unit, :rows, :columns)

genvartype(v::VarInfo, ::Val{:Tabulate}; V, _...) = @q Tabulate{$V}

gendefault(v::VarInfo, ::Val{:Tabulate}) = gendefaultvalue(v, parameter=true)

genupdate(v::VarInfo, ::Val{:Tabulate}, ::MainStep; kw...) = nothing


text:
include("../src/classification.jl")
using Main.Classification
using Random
using PyCall
using Statistics
#using Gadfly
import PyPlot; const plt = PyPlot

Random.seed!(4)

sklsets = pyimport("sklearn.datasets")
observables, target = sklsets.make_classification(n_samples=2000,
                                                  n_features=25,
                                                  n_informative=24,
                                                  n_redundant=1,
                                                  n_repeated=0,
                                                  n_classes=2,
                                                  random_state=42,
                                                  flip_y=0.05,
                                                  class_sep=1.0)

(Xtrain, Ytrain), (Xtest, Ytest) = traintestsplit(observables, target, 0.8, 2000)

Xtrain, Ytrain = (Xtrain, Ytrain) .|> collect ∘ transpose
Xtest, Ytest = (Xtest, Ytest) .|> collect ∘ transpose

Y = onehot(Ytrain)

gd_optim = GradientDescent(GDContext(learningrate=0.1, maxiterations=500))
optim = StochasticGradientDescent(SGDContext(learningrate=0.1, maxiterations=500))
optim2 = StochasticGradientDescent(SGDContext(learningrate=0.1, maxiterations=500))
addvalidationset!(optim2, Xtest, vec(Ytest))
addvalidationset!(gd_optim, Xtest, vec(Ytest))
addvalidationset!(optim, Xtest, vec(Ytest))
net = NeuralNet([(5, :relu), (2, :sigmoid)])
#addSoftmaxLayer!(net)
@time fit!(net, gd_optim, copy(Xtrain), copy(Y), copy(Ytrain))

net = NeuralNet([(5, :leakyrelu), (2, :sigmoid)])
@time fit!(net, optim2, copy(Xtrain), copy(Y), copy(Ytrain))

net = NeuralNet([(5, :relu), (2, :sigmoid)])
@time fit!(net, optim, copy(Xtrain), copy(Y), copy(Ytrain))

@show mean(predict(net, Xtrain) .== Ytrain)
@show mean(predict(net, Xtest) .== Ytest)

# plt = plot(x = 1:length(optim.loss), y = optim.loss,
     # Geom.line, Guide.xlabel("iterations"),
     # Guide.ylabel("loss"))
     #
fig, ax = plt.subplots()
line = ax.plot(gd_optim.loss, label="Gradient Descent")
ax.plot(gd_optim.validationloss, "--", c=line[1].get_color())
line = ax.plot(optim.loss, label="Stochastic")
ax.plot(optim.validationloss, "--", c=line[1].get_color())
line = ax.plot(optim2.loss, label="Stochastic Leaky")
ax.plot(optim2.validationloss, "--", c=line[1].get_color())
ax.set_ylabel("Accuracy")
ax.set_xlabel("Epochs")
ax.legend()
plt.show()


text:
""" Spectral Newton algorithm using a line search supporting stopping
"""
function Newton_G_StopLS(nlp :: AbstractNLPModel;
                       x :: Vector{T}=copy(nlp.meta.x0),
                       stp :: NLPStopping = NLPStopping(nlp,
                                                        NLPAtX(nlp.meta.x0)),
                       LS_algo   :: Function = bracket,
                       LS_logger :: AbstractLogger = Logging.NullLogger(),
                       NwtDirection :: Function = NwtdirectionSpectral,
                       hessian_rep :: Function = hessian_dense,
                       kwargs...      # eventually options for the line search
                       ) where T

    n = nlp.meta.nvar

    @info log_header([:iter, :f, :dual, :step, :slope], [Int, T, T, T, T],
                     hdr_override=Dict(:f=>"f(x)", :dual=>"‖∇f‖", :slope=>"∇fᵀd"))
    f ::T  = obj(nlp,x)
    ∇f :: Vector{T} = grad(nlp, x)

    OK = update_and_start!(stp, x = x, fx = f, gx = ∇f)

    τ₀ = 0.0005
    τ₁ = 0.9999
    d = similar(x)
    ϕ, ϕstp = prepare_LS(stp, x, d, τ₀, f, ∇f)

    @info log_row(Any[0, f, norm(∇f)])

    while !OK
        #Hx = hess(nlp, x)
        #H = Matrix(Symmetric(Hx,:L))

        H = hessian_rep(nlp, x)

        #Stopping.update!(stp.current_state; Hx = H)

        #Δ, O = eigen(H)
        # Boost negative values of Δ to 1e-8
        # devise an adaptative value for γ
        #γ = 1e-6
        #D = abs.(Δ) + max.((γ .- abs.(Δ)), 0.0) .*ones(n)
        #d = - O*diagm(1.0 ./ D)*O'*∇f

        d = NwtDirection(H, ∇f; kwargs...)

        #hp0 = ∇f'*d
        # Simple line search call

        t, x, f, ∇f = linesearch(ϕ, ϕstp, x, d, f, ∇f, τ₀, τ₁, logger = LS_logger, algo = LS_algo)#; kwargs...)
        fail_sub_pb = ~ϕstp.meta.optimal

        if fail_sub_pb
            OK = true
            stp.meta.fail_sub_pb = true
        else
            OK = update_and_stop!(stp, x = x, gx = ∇f, fx=f)
        end
        @info log_row(Any[stp.meta.nb_of_stop, f, norm(∇f), t, hp0])
    end

    if !stp.meta.optimal
        #@warn "Optimalité non atteinte"
        @warn status(stp, list=true)
    end

    return stp
end

max_stars_repo_name: vepiteski/LSDescent

max_issues_repo_name: vepiteski/LSDescent

max_forks_repo_name: vepiteski/LSDescent


text:
# Static Parameters
f_c = 20
x_max = 1.0
num_x = 10*f_c
filter_x = ones(2*f_c+1)

# Dynamic Parameters
K  = 2
tau = 1.0/(K*1.0)
v_max = 0.5
num_v = 10

# Variability of weights
minWeights = 0.9
maxWeights = 1.1

# Variability of number of particles
min_number_part = 20
max_number_part = 30

# Number of generated examples
num_trials = 1000

# test case
test_case = () -> TestCases.cloud_1d_full(x_max, v_max, minWeights, maxWeights, K, tau, rand(min_number_part:max_number_part))
# WARNING: changing anything in the test_case will generate troubles if the rejection_sampling algorithm
# in Utils.js is not either updated for it, or disabled.

# Noise parameters
noises_data = linspace(0,0.1,5)
noises_data = noises_data[2:end]
noises_position = linspace(0,0.06,5)
noises_position = noises_position[2:end]

### Location of data folder
dataFolder = "data/1Dsimulations"
# loaded data for rejection sampling of simulations
bins = np.load(dataFolder*"/separationDistribBins_20particles.npy")
density = np.load(dataFolder*"/separationDistribVal_20particles.npy")

println(" WARNING: The rejection sampling algorithm is being used, with interval of interest: [",bins[1]," ",bins[end],"]. To change it check the separationDistribution.jl file.")

### Cases to be tested

noiseless_dynamic = true
noiseless_static = true
noise_dynamic = false
noise_static = false
curvature_static = false

cases = (noiseless_dynamic, noiseless_static, noise_dynamic, noise_static, curvature_static)


text:
using TightBindingLattice
using JLD2

POINT_SYMMETRY_DATABASE = [PointSymmetryDatabase.get(gn) for gn in 1:32]
@save "PointGroup3D.jld2" POINT_SYMMETRY_DATABASE


text:
using Random

"""
    function information_set(rng::AbstractRNG, j::Node, n_I::Int)

Generates random information sets for chance and decision nodes.
"""
function information_set(rng::AbstractRNG, j::Node, n_I::Int)
    m = min(rand(rng, 0:n_I), j-1)
    I_j = shuffle(rng, 1:(j-1))[1:m]
    return sort(I_j)
end

"""
    function information_set(rng::AbstractRNG, leaf_nodes::Vector{Node}, n::Int)

Generates random information sets for value nodes.
"""
function information_set(rng::AbstractRNG, leaf_nodes::Vector{Node}, n::Int)
    @assert n ≥ 1
    non_leaf_nodes = shuffle(rng, setdiff(1:n, leaf_nodes))
    l = length(non_leaf_nodes)
    if isempty(leaf_nodes)
        m = rand(rng, 1:l)
    else
        m = rand(rng, 0:l)
    end
    I_v = [leaf_nodes; non_leaf_nodes[1:m]]
    return sort(I_v)
end

"""
    random_diagram!(rng::AbstractRNG, diagram::InfluenceDiagram, n_C::Int, n_D::Int, n_V::Int, m_C::Int, m_D::Int, states::Vector{Int})

Generate random decision diagram with `n_C` chance nodes, `n_D` decision nodes, and `n_V` value nodes.
Parameter `m_C` and `m_D` are the upper bounds for the size of the information set.

# Arguments
- `rng::AbstractRNG`: Random number generator.
- `diagram::InfluenceDiagram`: The (empty) influence diagram structure that is filled by this function
- `n_C::Int`: Number of chance nodes.
- `n_D::Int`: Number of decision nodes.
- `n_V::Int`: Number of value nodes.
- `m_C::Int`: Upper bound for size of information set for chance nodes.
- `m_D::Int`: Upper bound for size of information set for decision nodes.
- `states::Vector{State}`: The number of states for each chance and decision node
    is randomly chosen from this set of numbers.

# Examples
```julia
rng = MersenneTwister(3)
diagram = InfluenceDiagram()
random_diagram!(rng, diagram, 5, 2, 3, 2, 2, [2,3])
```
"""
function random_diagram!(rng::AbstractRNG, diagram::InfluenceDiagram, n_C::Int, n_D::Int, n_V::Int, m_C::Int, m_D::Int, states::Vector{Int})
    n = n_C + n_D
    n_C ≥ 0 || throw(DomainError("There should be `n_C ≥ 0` chance nodes."))
    n_D ≥ 0 || throw(DomainError("There should be `n_D ≥ 0` decision nodes"))
    n ≥ 1 || throw(DomainError("There should be at least one chance or decision node `n_C+n_D≥1`."))
    n_V ≥ 1 || throw(DomainError("There should be `n_V ≥ 1` value nodes."))
    m_C ≥ 1 || throw(DomainError("Maximum size of information set should be `m_C ≥ 1`."))
    m_D ≥ 1 || throw(DomainError("Maximum size of information set should be `m_D ≥ 1`."))
    all(s > 1 for s in states) || throw(DomainError("Minimum number of states possible should be 2."))

    # Create node indices
    U = shuffle(rng, 1:n)
    diagram.C = [Node(c) for c in sort(U[1:n_C])]
    diagram.D = [Node(d) for d in sort(U[(n_C+1):n])]
    diagram.V = [Node(v) for v in collect((n+1):(n+n_V))]

    diagram.I_j = Vector{Vector{Node}}(undef, n+n_V)
    # Create chance and decision nodes
    for c in diagram.C
        diagram.I_j[c] = information_set(rng, c, m_C)
    end
    for d in diagram.D
        diagram.I_j[d] = information_set(rng, d, m_D)
    end

    # Assign each leaf node to a random value node
    leaf_nodes = setdiff(1:n, (diagram.I_j[c] for c in diagram.C)..., (diagram.I_j[d] for d in diagram.D)...)
    leaf_nodes_v = Dict(v=>Node[] for v in diagram.V)
    for j in leaf_nodes
        v = rand(rng, diagram.V)
        push!(leaf_nodes_v[v], j)
    end

    # Create values nodes
    for v in diagram.V
        diagram.I_j[v] = information_set(rng, leaf_nodes_v[v], n)
    end

    diagram.S = States(State[rand(rng, states, n)...])
    diagram.X = Vector{Probabilities}(undef, n_C)
    diagram.Y = Vector{Utilities}(undef, n_V)

    diagram.Names = ["$(i)" for i in 1:(n+n_V)]
    statelist = []
    for i in 1:n
        push!(statelist, ["$(j)" for j in 1:diagram.S[i]])
    end
    diagram.States = statelist

    for c in diagram.C
        random_probabilities!(rng, diagram, c)
    end

    for v in diagram.V
        random_utilities!(rng, diagram, v)
    end

    diagram.P = DefaultPathProbability(diagram.C, diagram.I_j[diagram.C], diagram.X)
    diagram.U = DefaultPathUtility(diagram.I_j[diagram.V], diagram.Y)

    return diagram
end

"""
    function random_probabilities!(rng::AbstractRNG, diagram::InfluenceDiagram, c::Node; n_inactive::Int=0)

Generate random probabilities for chance node `c`.

# Examples
```julia
rng = MersenneTwister(3)
diagram = InfluenceDiagram()
random_diagram!(rng, diagram, 5, 2, 3, 2, 2, [2,3])
c = diagram.C[1]
random_probabilities!(rng, diagram, c)
```
"""
function random_probabilities!(rng::AbstractRNG, diagram::InfluenceDiagram, c::Node; n_inactive::Int=0)
    if !(c in diagram.C)
        throw(DomainError("Probabilities can only be added for chance nodes."))
    end
    I_c = diagram.I_j[c]
    states = diagram.S[I_c]
    state = diagram.S[c]
    if !(0 ≤ n_inactive ≤ prod([states...; (state - 1)]))
        throw(DomainError("Number of inactive states must be < prod([S[I_j]...;, S[j]-1])"))
    end

    # Create the random probabilities
    data = rand(rng, states..., state)

    # Create inactive chance states
    if n_inactive > 0
        # Count of inactive states per chance stage.
        a = zeros(Int, states...)
        # There can be maximum of (state - 1) inactive chance states per stage.
        b = repeat(vec(CartesianIndices(a)), state - 1)
        # Uniform random sample of n_inactive states.
        r = shuffle(rng, b)[1:n_inactive]
        for s in r
            a[s] += 1
        end
        for s in CartesianIndices(a)
            indices = CartesianIndices(data[s, :])
            i = shuffle(rng, indices)[1:a[s]]
            data[s, i] .= 0.0
        end
    end

    # Normalize the probabilities
    for s in CartesianIndices((states...,))
        data[s, :] /= sum(data[s, :])
    end

    index_c = findfirst(j -> j==c, diagram.C)
    diagram.X[index_c] = Probabilities(c, data)
end

scale(x::Utility, low::Utility, high::Utility) = x * (high - low) + low

"""
    function random_utilities!(rng::AbstractRNG, diagram::InfluenceDiagram, v::Node; low::Float64=-1.0, high::Float64=1.0)

Generate random utilities between `low` and `high` for value node `v`.

# Examples
```julia
rng = MersenneTwister(3)
diagram = InfluenceDiagram()
random_diagram!(rng, diagram, 5, 2, 3, 2, 2, [2,3])
v = diagram.V[1]
random_utilities!(rng, diagram, v)
```
"""
function random_utilities!(rng::AbstractRNG, diagram::InfluenceDiagram, v::Node; low::Float64=-1.0, high::Float64=1.0)
    if !(v in diagram.V)
        throw(DomainError("Utilities can only be added for value nodes."))
    end
    if !(high > low)
        throw(DomainError("high should be greater than low"))
    end
    I_v = diagram.I_j[v]
    data = rand(rng, Utility, diagram.S[I_v]...)
    data = scale.(data, Utility(low), Utility(high))

    index_v = findfirst(j -> j==v, diagram.V)
    diagram.Y[index_v] = Utilities(v, data)
end

"""
    function LocalDecisionStrategy(rng::AbstractRNG, diagram::InfluenceDiagram, d::Node)

Generate random decision strategy for decision node `d`.

# Examples
```julia
rng = MersenneTwister(3)
diagram = InfluenceDiagram()
random_diagram!(rng, diagram, 5, 2, 3, 2, 2, rand(rng, [2,3], 5))
LocalDecisionStrategy(rng, diagram, diagram.D[1])
```
"""
function LocalDecisionStrategy(rng::AbstractRNG, diagram::InfluenceDiagram, d::Node)
    I_d = diagram.I_j[d]
    states = diagram.S[I_d]
    state = diagram.S[d]
    data = zeros(Int, states..., state)
    for s in CartesianIndices((states...,))
        s_j = rand(rng, 1:state)
        data[s, s_j] = 1
    end
    LocalDecisionStrategy(d, data)
end


text:
using FastGaussQuadrature
using BenchmarkTools
import DoubleFloats: Double64
using HypergeometricFunctions: drummond2F0

function L₀truncated(x::T, μ) where T
    N = 30
	Σ = zero(T)
	for n in 1:N
        s = one(T)
        for i in 0:n-1
            s *= (1/2 + μ + i)/(i+1)
            s *= (1/2 - μ + i) / x
        end
		if iseven(n)
			Σ += s
		else
			Σ -= s
		end
	end
	return one(T) + Σ
end

# U(a, b, z::T) = z^(-a) * pFq([a, one(T) + a - b], T[], -one(T)/z)
function L₀(x::T, μ) where T
    return drummond2F0(1/2+μ, 1/2-μ, -1/x; kmax=10_000Int(ceil(μ)))
end

W₀(μ, z) = exp(-z/2) * L₀(z, μ)

function QS_H1(f, a, b, ω::T, ν, n) where T
    i_ω = im / ω  # this is used a lot
    x_gl, w_gl = gausslaguerre(n);
    Σa = complex(zero(T))
    for k in 1:n
        s = w_gl[k] * f(a + i_ω * x_gl[k]) *
            L₀(-2im * ω * (a + i_ω * x_gl[k]), ν) / √(a + i_ω * x_gl[k])
        Σa += s
    end
    Σb = complex(zero(T))
    for k in 1:n
        s = w_gl[k] * f(b + i_ω * x_gl[k]) *
            L₀(-2im * ω * (b + i_ω * x_gl[k]), ν) / √(b + i_ω * x_gl[k])
        Σb += s
    end
    diff = i_ω * exp(im * ω * a) * Σa - i_ω * exp(im * ω * b) * Σb
    res = √(2 / π / ω) * exp(-im * π/2 * (ν+1/2)) * diff
    return real(res)
end

function QS_H2(f, a, b, ω::T, ν, n) where T
    i_ω = im / ω  # this is used a lot
    x_gl, w_gl = gausslaguerre(n);
    Σa = complex(zero(T))
    for k in 1:n
        s = w_gl[k] * f(a - i_ω * x_gl[k]) *
            L₀(2im * ω * (a - i_ω * x_gl[k]), ν) / √(a - i_ω * x_gl[k])
        Σa += s
    end
    Σb = complex(zero(T))
    for k in 1:n
        s = w_gl[k] * f(b - i_ω * x_gl[k]) *
            L₀(2im * ω * (b - i_ω * x_gl[k]), ν) / √(b - i_ω * x_gl[k])
        Σb += s
    end
    diff = i_ω * exp(-im * ω * b) * Σb - i_ω * exp(-im * ω * a) * Σa
    res = √(2 / π / ω) * exp(im * π/2 * (ν+1/2)) * diff
    return real(res)
end

##
Tn = Double64
f(z) = z^(-2)
         #  a    b             ω        ν     n
QS_H1(f, Tn(1.0), Tn(2.0),  Tn(100.0), Tn(400.0), 1000)

##
# z = 0.5
# μ = 3.0
# κ = 4.0
# W(k, m, z) = exp(-z/2) * z^(m + 1/2) * U(1/2 + m - k, 1 + 2m, z)


text:
include("smear_FD.jl")

using PyPlot

const plt = PyPlot

function test_plot(kT)

    mu = -2.0

    @printf("kT = %f\n", kT)

    xmin = -5.0
    xmax =  5.0

    NptsPlot = 500
    x = Array{Float64}(NptsPlot)
    f = Array{Float64}(NptsPlot)
    for i = 1:NptsPlot
        x[i] = xmin + (i-1)*(xmax-xmin)/(NptsPlot-1)
        f[i] = smear_FD(x[i], mu, kT, is_spinpol=true)
    end

    plt.plot(x, f, marker="o", label="kT="*string(kT))

end

function test_main()
    plt.clf()
    test_plot(0.001)
    test_plot(0.01)
    test_plot(0.03)
    test_plot(0.1)
    test_plot(0.2)
    test_plot(0.3)
    test_plot(0.4)
    plt.legend()
    plt.grid()
    plt.savefig("plot_smear_FD.png", dpi=300)
end

test_main()


text:
# Helper data for CIE observer functions
include("cie_data.jl")

# Linear interpolation in [a, b] where x is in [0,1],
# or coerced to be if not.
function lerp(x, a, b)
    a + (b - a) * max(min(x, one(x)), zero(x))
end

"""
    hex(c)

Print a color as a RGB hex triple, or a transparent paint as an ARGB
hex quadruplet.
"""
function hex(c::RGB)
    @sprintf("%02X%02X%02X",
             round(Int, lerp(c.r, 0.0, 255.0)),
             round(Int, lerp(c.g, 0.0, 255.0)),
             round(Int, lerp(c.b, 0.0, 255.0)))
end
function hex(c::ARGB)
    @sprintf("%02X%02X%02X%02X",
             round(Int, lerp(alpha(c), 0.0, 255.0)),
             round(Int, lerp(red(c), 0.0, 255.0)),
             round(Int, lerp(green(c), 0.0, 255.0)),
             round(Int, lerp(blue(c), 0.0, 255.0)))
end

hex(c::Color) = hex(convert(RGB, c))
hex(c::Colorant) = hex(convert(ARGB, c))

"""
    weighted_color_mean(w1, c1, c2)

Returns the color `w1*c1 + (1-w1)*c2` that is the weighted mean of `c1` and
`c2`, where `c1` has a weight 0 ≤ `w1` ≤ 1.
"""
function weighted_color_mean(w1::Real, c1::Colorant, c2::Colorant)
    weight1 = convert(promote_type(eltype(c1), eltype(c2)),w1)
    weight2 = weight1 >= 0 && weight1 <= 1 ? oftype(weight1,1-weight1) : throw(DomainError())
    mapc((x,y)->weight1*x+weight2*y, c1, c2)
end
function weighted_color_mean(w1::Real, c1::Gray{Bool}, c2::Gray{Bool})
    # weighting of two Gray{Bool} would return different color type and therefore omitted
    throw(DomainError())
end

"""
    linspace(c1::Color, c2::Color, n=100)

Generates `n`>2 colors in a linearly interpolated ramp from `c1` to`c2`,
inclusive, returning an `Array` of colors.
"""
function linspace{T<:Colorant}(c1::T, c2::T, n::Integer=100)
    return T[weighted_color_mean(w1, c1, c2) for w1 in linspace(1.0,0.0,n)]
end

#Double quadratic Bezier curve
function Bezier{T<:Real}(t::T, p0::T, p2::T, q0::T, q1::T, q2::T)
    B(t,a,b,c)=a*(1.0-t)^2.0+2.0b*(1.0-t)*t+c*t^2.0
    if t <= 0.5
        return B(2.0t, p0, q0, q1)
    else #t > 0.5
        return B(2.0(t-0.5), q1, q2, p2)
    end
end

#Inverse double quadratic Bezier curve
function invBezier{T<:Real}(t::T, p0::T, p2::T, q0::T, q1::T, q2::T)
    invB(t,a,b,c)=(a-b+sqrt(b^2.0-a*c+(a-2.0b+c)*t))/(a-2.0b+c)
    if t < q1
        return 0.5*invB(t,p0,q0,q1)
    else #t >= q1
        return 0.5*invB(t,q1,q2,p2)+0.5
    end
end


text:
using Mocha

############################################################
# Based on the paper "Auto-Encoding Variational Bayes"
# by Diederik P Kingma, Max Welling
# http://arxiv.org/abs/1312.6114
#
# We perform unsupervised learning of a generative model of
# binary MNIST data, using the "Variational Auto-Encoder" architecture.
#
# 50-dim Gaussian latent variable z, with a standard normal prior
#
# Encoder MLP:
# - one hidden layer with 500 hidden units and tanh()
# - two linear readout layers for mu(z) and logvar(z)
#
# (We'll just use elementwise layers and the RandomNormal layer to
# implement the "reparameterization trick".)
#
# Decoder MLP matches Encoder:
# - one hidden layer with 500 hidden units and tanh()
# - one linear readout layer for p(x)
#
############################################################

function make_vae(backend, N_LATENT=50, N_HIDDEN_DEC=500, N_HIDDEN_ENC=500, N_OUT=784)

  data_layer      = HDF5DataLayer(name="train-data", source="data/train.txt", batch_size=100)

  enc1_layer      = InnerProductLayer(name="enc1", output_dim=N_HIDDEN_ENC, neuron=Neurons.Tanh(),
                                      weight_init = GaussianInitializer(std=0.01),
                                      bottoms=[:data], tops=[:enc1])
  enc_split_layer = SplitLayer(name="enc-split", bottoms=[:enc1], tops=[:enc1_mu_in, :enc1_sigma_in])

  enc_mu_layer    = InnerProductLayer(name="enc1-mu", output_dim=N_LATENT, neuron=Neurons.Identity(),
                                      weight_init = GaussianInitializer(std=0.01),
                                      bottoms=[:enc1_mu_in], tops=[:z_mu])
  enc_sigma_layer = InnerProductLayer(name="enc1-sigma", output_dim=N_LATENT,
                                      # This layer outputs log(sigma^2) in the original paper
                                      neuron=Neurons.Exponential(), # still has large jumps in KL on CPU
                                      weight_init = GaussianInitializer(std=0.01),
                                      bottoms=[:enc1_sigma_in], tops=[:z_sigma])
  zm_split_layer = SplitLayer(name="zm-split", bottoms=[:z_mu],
                              tops=[:z_mu_1, :z_mu_2])
  zs_split_layer = SplitLayer(name="zs-split", bottoms=[:z_sigma],
                              tops=[:z_sigma_1, :z_sigma_2])

  eps_layer       = RandomNormalLayer(name="eps",
                                      batch_sizes=[data_layer.batch_size], output_dims=[N_LATENT],
                                      tops=[:eps])
  z_noise_layer   = ElementWiseLayer(name="z-noise", bottoms=[:z_sigma_1, :eps], tops=[:z_noise],
                                     operation=ElementWiseFunctors.Multiply())
  z_layer         = ElementWiseLayer(name="z", bottoms=[:z_noise, :z_mu_1], tops=[:z],
                                     operation=ElementWiseFunctors.Add())

  dec1_layer      = InnerProductLayer(name="dec1", output_dim=N_HIDDEN_DEC, neuron=Neurons.Tanh(),
                                      weight_init = GaussianInitializer(std=0.01),
                                      bottoms=[:z], tops=[:dec1])

  dec_out_layer   = InnerProductLayer(name="dec-out", output_dim=N_OUT, neuron=Neurons.Sigmoid(),
                                      weight_init = GaussianInitializer(std=0.01),
                                      bottoms=[:dec1], tops=[:out])

  bce_loss_layer  = BinaryCrossEntropyLossLayer(name="bce-loss", bottoms=[:out, :data])

  kl_loss_layer   = GaussianKLLossLayer(name="kl-loss", bottoms=[:z_mu_2, :z_sigma_2])

  enc_layers = [enc1_layer, enc_split_layer, enc_mu_layer, enc_sigma_layer, zm_split_layer, zs_split_layer]
  z_layers = [eps_layer, z_noise_layer, z_layer]
  dec_layers = [dec1_layer, dec_out_layer]
  common_layers = [enc_layers, z_layers, dec_layers]
  loss_layers = [bce_loss_layer, kl_loss_layer]
  # put training net together, note that the correct ordering will automatically be established by the constructor
  net = Net("MNIST-VAE-train", backend, [data_layer, common_layers..., loss_layers...])
  return net
end

# This function reproduces the feed-forward part of the model above, layers 12-13,
# as a simple function

asarray{T}(blob::CuTensorBlob{T}) = begin
  arr = zeros(T, blob.shape)
  copy!(arr, blob)
  return arr
end

asarray(blob::CPUBlob) = blob.data

function make_latent_to_output(net::Net, dec1_index, dec_out_index)
  sigmoid(x) = 1 ./ (1 + exp(-x))

  function latent_to_output(z)
    # pad z with zeros to match the expected latent size
    N_Z = size(net.states[dec1_index].parameters[1].blob)[1]
    if length(z) < N_Z
      z_ = zeros(eltype(z), N_Z)
      z_[1:length(z)] = z
      z = z_
    end

    w1 = asarray(net.states[dec1_index].parameters[1].blob)
    b1 = asarray(net.states[dec1_index].parameters[2].blob)

    w2 = asarray(net.states[dec_out_index].parameters[1].blob)
    b2 = asarray(net.states[dec_out_index].parameters[2].blob)

    dec1 = tanh(w1' * z + b1)
    dec_out = sigmoid(w2' * dec1 + b2)
    return reshape(dec_out, 28, 28)'
  end
  return latent_to_output
end


text:
test_string ="abc

a
b
c

ab
ac

a
a
a
a

b"

test_list = split(test_string,'\n')

list = readlines("day6_input.txt")

function get_union(list)
    groups = ""
    total = 0
    for n in 1:length(list)
        if list[n] == ""
            groups *= ","
        else
            groups *= list[n]
        end
    end
    return groups
end

function get_union_total_from_groups(list)
    return sum(length.(Set.(split(get_union(list), ','))))
end

function get_union_and_total(list)
    groups = Set()
    group = ""
    total = 0
    for n in 1:length(list)
        if list[n] == ""
            push!(groups, group)
            total += length(group)
            group = ""
        else
            group = union(group, list[n])
        end
    end
    push!(groups, group)
    total += length(group)
    return groups, total
end

function get_union_total(list)
    group = ""
    total = 0
    for n in 1:length(list)
        if list[n] == ""
            total += length(group)
            group = ""
        else
            group = union(group, list[n])
        end
    end
    total += length(group)
    return total
end

get_union_total_oneliner(list) = sum(length.(Set.(split(join(map(
        x -> x == "" ? "," : x, list)), ','
    ))))

get_intersection_total_oneliner(list) = sum(length.(Set.(split(join(map(
        x -> x == "" ? "," : x, list)), ','
    ))))

# No gain in making a map to broadcast directly to the composition
# of length and Set. I thought that would avoid one more loop through
# the list, but it either doesn't avoid or it does, but makes no difference
# somehow

get_union_total_oneliner2(list) = sum((length∘Set).(split(join(map(
    x -> x == "" ? "," : x, list)),',')))

@show get_union_total_from_groups(test_list)

@show get_union_total_from_groups(list)

@show get_union_total(list)

#=
```julia
julia> @btime get_union_total_from_groups(list)
  3.438 ms (8115 allocations: 18.38 MiB)
6809

julia> @btime get_union_total(list)
  1.860 ms (18572 allocations: 1.64 MiB)
6809

julia> @btime get_union_total_oneliner(list)
  660.330 μs (4136 allocations: 431.29 KiB)
6809

julia> @btime get_union_total_oneliner2(list)
  757.013 μs (5137 allocations: 454.79 KiB)
6809
```
 =#
function get_intersection_and_total(list)
    nmax = length(list)
    groups = Set()
    abc = "abcdefghijklmnopqrstuvwxyz"
    group = abc
    total = 0
    for n in 1:nmax
        if list[n] == ""
            push!(groups, group)
            total += length(group)
            group = abc
        else
            group = intersect(group, list[n])
        end
    end
    push!(groups, group)
    total += length(group)
    return groups, total
end

function get_intersection_total(list)
    nmax = length(list)
    abc = "abcdefghijklmnopqrstuvwxyz"
    group = abc
    total = 0
    for n in 1:nmax
        if list[n] == ""
            total += length(group)
            group = abc
        else
            group = intersect(group, list[n])
        end
    end
    total += length(group)
    return total
end

@show get_intersection_total(list)

#=
```julia
julia> @btime get_intersection_total(list)
  4.145 ms (29922 allocations: 2.80 MiB)
3394
```
 =#

function get_totals(list)
    nmax = length(list)
    abc = "abcdefghijklmnopqrstuvwxyz"
    group_union = ""
    total_union = 0
    group_intersection = abc
    total_intersection = 0
    for n in 1:nmax
        if list[n] == ""
            total_union += length(group_union)
            total_intersection += length(group_intersection)
            group_union = ""
            group_intersection = abc
        else
            group_union = union(group_union, list[n])
            group_intersection = intersect(group_intersection, list[n])
        end
    end
    total_union += length(group_union)
    total_intersection += length(group_intersection)
    return total_union, total_intersection
end

# No gain in time doing both at once:
#=
```julia
julia> @btime get_union_total(list)
  1.860 ms (18572 allocations: 1.64 MiB)
6809

julia> @btime get_intersection_total(list)
  4.145 ms (29922 allocations: 2.80 MiB)
3394

julia> @btime get_totals(list)
  6.185 ms (48493 allocations: 4.44 MiB)
(6809, 3394)
``` =#


text:
# -------------- Install packages --------------
using Pkg
Pkg.add(
    ["Parameters", "RecursiveArrayTools", "Statistics", "NLsolve",
    "DifferentialEquations", "PyPlot"]
    )

max_stars_repo_name: McCannLab/NDID

max_issues_repo_name: McCannLab/NDID

max_forks_repo_name: McCannLab/NDID


text:
#= ######################################################
inside functions
###################################################### =#

const sigma = 5.670374419E-8

function create_identity_matrix(n::T) where T<:Integer
    id_mat = zeros(T,n,n)
    for i = 1:n
        id_mat[i,i] = 1
    end
    return id_mat
end

#= ######################################################
final functions
###################################################### =#

function tempsolver(m, vfmat, temp, epsilon)
    # old algorithm used in matlab (slightly optimized)
    # prepare matrix for LGS
    # J --> outgoing radiation of element
    # G --> incoming radiation of element
    mat = ((epsilon .- 1) .* vfmat)
    # id_mat = Matrix{Int64}(I, m.nelem, m.nelem) # LinearAlgebra Pkg
    id_mat = create_identity_matrix(m.nelem)
    mat .= mat .+ id_mat
    # solve LGS
    B = (temp[:,1].^4) .* sigma .* epsilon
    J = mat \ B
    # Bestimmung der einfallenden Strahlung G
    G = zeros(Float64,m.nelem,1)
    for i = 1:m.nelem
        G[i,1] = sum(J[:,1] .* vfmat[i, :])
    end
    # Qp Radiation per element
    area = [m.elem[i].area for i = m.elem2par[1].first:m.elem2par[end].last]
    Qp = (J[:,1] - G[:,1]) .* area[:,1]
    Ga = G[:,1] .* area[:,1]
    return Qp, Ga
end

function set_bc_part!(m, bcv, parts, bcp)
    # give all face elements one temperature
    for i = parts
        e1 = m.elem2par[i].first
        e2 = m.elem2par[i].last
        bcv[e1:e2,1] .= bcp
    end
end


text:
# ---------------------------------------------------------------
# Parts of the Julia code in this package is ported from the Lisp
# and C++ code in 'Calendrical Calculations' by Nachum Dershowitz
# and Edward M. Reingold, Software---Practice & Experience,
# vol. 20, no. 9 (September, 1990), pp. 899--928.
# This code is in the public domain, but any use of it should
# publicly acknowledge the above sources, which is done herewith.
# Copyright of this port Peter Luschny, 2022.
# Licensed under the MIT license.
# ---------------------------------------------------------------

module Calendars

export EC, CE, JD, AM, AH, ID, EN, JN, DN, XX
export Calendar, Year, Month, Day, Date, CDate, YMD
export CalendarDates, CalendarName, CDateStr, DateStr, DateTable
export DayNumberFromDate, DateFromDayNumber, ConvertDate
export isValidDate, Duration, DayOfYear, ConvertOrdinalDate
export PrintDateLine, PrintDateTable, PrintEuropeanMonth
export SaveEuropeanMonth, WeekDay, WeekDays, PrintIsoWeek, IDate
export ProfileYearAsEuropean

# Export only functions from CalendarUtils and from this file,
# which basically is a wrapper around a collection of calendars.
include("CalendarUtils.jl")

# Consider the following included files as 'intern', and
# do not export their functions directly. Their algorithms
# follow Dershowitz/Reingold and use RataDie-numbers
# whereas here we use EuroDay-numbers.
include("DayNumbers.jl")
include("GregorianCalendar.jl")
include("JulianCalendar.jl")
include("EuropeanCalendar.jl")
include("HebrewCalendar.jl")
include("IslamicCalendar.jl")
include("IsoCalendar.jl")
include("InteractiveCalendar.jl")

"""

```julia
DayNumberFromDate(date::Cdate, string::Bool, show::Bool)
```

Return the day number correspondig to the calendar date.
The parts of the date can also be given separately.

If the optional parameter `string` is `true` the function returns
the day number in string format, otherwise as an integer.
`string` is 'false' by default.

If the optional parameter `show` is 'true', date
and number are printed. `show` is 'false' by default.

For example:
```julia
julia> DayNumberFromDate("Gregorian", 1756, 1, 27)
```

Using the acronym CE for 'Common Epoch' which is synonymous with
'Gregorian', this can also be written as

```julia
julia> DayNumberFromDate(CE, 1756, 1, 27)
```

This returns the European day number 641027 as an integer by default.
If `string` is `true` the string "EN#641029" is returned. If `show`
is `true` the line "CE-1756-01-27 -> EN#641029" is printed.

If an error occurs 0 (representing the invalid day
number) is returned.
"""
function DayNumberFromDate(d::CDate, string=false, show=false)
    # Don't process invalid dates
    ! isValidDate(d) && return 0

    calendar = d[1]
    if calendar == CE
        dn = DayNumberGregorian(d)
    elseif calendar == JD
        dn = DayNumberJulian(d)
    elseif calendar == EC
        dn = DayNumberEuropean(d)
    elseif calendar == AM
        dn = DayNumberHebrew(d)
    elseif calendar == AH
        dn = DayNumberIslamic(d)
    elseif calendar == ID
        dn = DayNumberIso(d)
    elseif calendar == DN
        dn = d[4]
    elseif calendar == JN
        dn = JulianNumToFixNum(d[4])
    elseif calendar == EN
        dn = EuroNumToFixNum(d[4])
    else
        @warn("Unknown calendar: $calendar")
        return InvalidDayNumber
    end

    enum = FixNumToEuroNum(dn)

    if show
        println(CDateStr(d), " -> ", CDateStr(EN, enum))
    end

    return string ? CDateStr(EN, enum) : enum
end

DayNumberFromDate(cal::String, year::DPart, month::DPart, day::DPart, string=false,
show=false) = DayNumberFromDate((StringToSymbol(cal), year, month, day), string, show)

DayNumberFromDate(cal::String, d::YMD, string=false, show=false) =
DayNumberFromDate((StringToSymbol(cal), d[1], d[2], d[3]), string, show)

DayNumberFromDate(cal::DPart, year::DPart, month::DPart, day::DPart, string=false,
show=false) = DayNumberFromDate((cal, year, month, day), string, show)

"""

```julia
DateFromDayNumber(calendar::DPart, enum::DPart, string::Bool, show::Bool)
```

Return the date from a European day number expressed in the calendar.

The day number `enum` must be an integer >= 1.

If the optional parameter `string` is `true` the function returns the date
as a string, otherwise as an integer tuple. `string` is `false` by default.

If the optional parameter `show` is set to `true`, date and number
are printed. `show` is `false` by default.

For example:
```julia
julia> DateFromDayNumber("Gregorian", 641029)
```

Alternatively you can also write
```julia
julia> DateFromDayNumber(CE, 641029)
```

By default this returns the calendar representation (2, 1756, 1, 27).
If `string` is true the string "CE-1756-01-27" is returned.
If `show` is 'true' the line "EN#0641029 -> CE-1756-01-27"
is printed.

If an error occurs (0, 0, 0, 0) (representing the invalid date) is returned.
"""
function DateFromDayNumber(calendar::DPart, enum::DPart, string=false, show=false)
    # Use a symbol for the calendar name.
    dn = EuroNumToFixNum(enum)

    if calendar == CE
        date = DateGregorian(dn)
    elseif calendar == JD
        date = DateJulian(dn)
    elseif calendar == EC
        date = DateEuropean(dn)
    elseif calendar == AM
        date = DateHebrew(dn)
    elseif calendar == AH
        date = DateIslamic(dn)
    elseif calendar == ID
        date = DateIso(dn)
    elseif calendar == DN
        date = (DN, 0, 0, EuroNumToFixNum(enum))
    elseif calendar == EN
        date = (EN, 0, 0, enum)
    elseif calendar == JN
        date = (JN, 0, 0, EuroNumToJulianNum(enum))
    else
        @warn("Unknown calendar: $calendar")
        return InvalidDate
    end

    if show
        println(CDateStr(EN, enum), " -> ", CDateStr(date))
    end

    return string ? CDateStr(date) : date
end

DateFromDayNumber(calendar::String, enum::DPart, string=false, show=false) =
DateFromDayNumber(StringToSymbol(calendar, true), enum, string, show)

"""

```julia
ConvertDate(date::CDate, calendar::DPart, string::Bool, show::Bool)
```

Convert the given date to a date in the calendar `calendar`.
Admissible names for the calendar are listed in the CDate docstring.

If the optional parameter `string` is `true` the function returns the date
as a string, otherwise as an integer tuple. `string` is `false` by default.

If the optional parameter `show` is set to 'true', both dates
are printed. `show` is 'false' by default.

For example:
```julia
julia> ConvertDate(CE, 1756, 1, 27, AM)
```

Alternatively this can also be written as

```julia
julia> ConvertDate("Gregorian", 1756, 1, 27, "Hebrew")
```

This converts the Gregorian date (1756, 1, 27) to the Hebrew date (5516, 11, 25).
By default the calendar representation (4, 5516, 11, 25) is returned.

If `string` is true the string "AM-5516-11-25" is returned.
If `show` is 'true' the following line is printed:

```julia
    "CE-1756-01-27 -> AM-5516-11-25"
```

If an error occurs (0, 0, 0, 0) (representing the invalid date) is returned.
"""
function ConvertDate(date::CDate, calto::DPart, string=false, show=false)
    ! isValidDate(date) && return InvalidDate

    # TODO There should be a general validation method in Utils.
    if StringToSymbol(SymbolToString[calto]) != calto
        @warn("Unknown calendar: $calto")
        return InvalidDate
    end

    enum = DayNumberFromDate(date)
    rdate = DateFromDayNumber(calto, enum)
    show && println(CDateStr(date), " -> ", CDateStr(rdate))

    return string ? CDateStr(rdate) : rdate
end

ConvertDate(calendar::DPart, year::DPart, month::DPart, day::DPart, to::DPart,
string=false, show=false) = ConvertDate((calendar, year, month, day), to, string, show)

ConvertDate(calendar::String, year::DPart, month::DPart, day::DPart, to::String,
string=false, show=false) = ConvertDate((StringToSymbol(calendar), year, month, day),
StringToSymbol(to), string, show)

ConvertDate(calendar::String, d::YMD, to::String, string=false, show=false) =
ConvertDate((StringToSymbol(calendar), d[1], d[2], d[3]), StringToSymbol(to), string, show)

ConvertDate(date::CDate, calto::String, string=false, show=false) =
ConvertDate(date, StringToSymbol(calto), string, show)

"""

```julia
CalendarDates(date::CDate, show::Bool)
```

Return a table of the dates of all supported calendars
corresponding to `date`. If the optional parameter `show`
is set to 'true', the date table is printed. `show` is
'false' by default.

```julia
julia> CalendarDates("Gregorian", 1756, 1, 27, true)
```

Alternatively you can also write

```julia
julia> CalendarDates(CE, 1756, 1, 27, true)
```

This computes a `DateTable`, which is a tuple of six calendar
dates plus the Euro day number and the Julian day number.
If `show` is 'true' the table below will be printed.
```julia
    European  EC-1756-01-27
    Common    CE-1756-01-27
    Julian    JD-1756-01-16
    Hebrew    AM-5516-11-25
    Islamic   AH-1169-04-24
    IsoDate   ID-1756-05-02
    EuroNum   EN#0641029
    JulianNum JN#2362452
```
"""
function CalendarDates(date::CDate, show=false)::DateTable

    dn = DayNumberFromDate(date)
    Table = (
        DateFromDayNumber(EC, dn),
        DateFromDayNumber(CE, dn),
        DateFromDayNumber(JD, dn),
        DateFromDayNumber(AM, dn),
        DateFromDayNumber(AH, dn),
        DateFromDayNumber(ID, dn),
        DateFromDayNumber(EN, dn),
        DateFromDayNumber(JN, dn)
    )
    show && PrintDateTable(Table)
    return Table
end

CalendarDates(calendar::String, year::DPart, month::DPart, day::DPart,
show=false) = CalendarDates((StringToSymbol(calendar), year, month, day), show)

CalendarDates(calendar::String, d::YMD, show=false) =
CalendarDates((StringToSymbol(calendar), d[1], d[2], d[3]), show)

CalendarDates(cal::DPart, year::DPart, month::DPart, day::DPart,
show=false) = CalendarDates((cal, year, month, day), show)

"""

```julia
isValidDate(date::CDate, warn=true)::Bool
```

Query if the given `date`` is a valid calendar date.

For example:
```julia
julia> isValidDate("Gregorian", 1756, 1, 27)
```

Alternatively you can also write

```julia
julia> isValidDate(CE, 1756, 1, 27)
```

This query affirms that 1756-01-27 is a valid Gregorian date.
But the next two queries return 'false' and 'true', respectively.

```julia
julia> isValidDate(CE, 1900, 2, 29)

julia> isValidDate(JD, 1900, 2, 29)
```
"""
function isValidDate(d::CDate, warn=true)
    calendar = d[1]

    if calendar == CE
        val = isValidDateGregorian(d, warn)
    elseif calendar == JD
        val = isValidDateJulian(d, warn)
    elseif calendar == EC
        val = isValidDateEuropean(d, warn)
    elseif calendar == AM
        val = isValidDateHebrew(d, warn)
    elseif calendar == AH
        val = isValidDateIslamic(d, warn)
    elseif calendar == ID
        val = isValidDateIso(d, warn)
    elseif calendar == DN
        val = d[4] > 0
    elseif calendar == EN
        val = d[4] > 0
    elseif calendar == JN
        val = d[4] > 0
    else
        @warn("Unknown calendar: $calendar")
        return false
    end
    return val
end

isValidDate(calendar::String, year::DPart, month::DPart, day::DPart) =
isValidDate((StringToSymbol(calendar), year, month, day))

isValidDate(calendar::DPart, year::DPart, month::DPart, day::DPart) =
isValidDate((calendar, year, month, day))

isValidDate(calendar::String, d::YMD) =
isValidDate((StringToSymbol(calendar), d[1], d[2], d[3]))

"""

```julia
DayOfYear(date::CDate)
```

Return the ordinal of the day in the given calendar. Also known as the 'ordinal date' relative to the epoch of the calendar.

```julia
julia> DayOfYear(EC, 1756, 1, 27)
```

Alternatively you can write

```julia
julia> DayOfYear("European", 1756,  1, 27)
```

From the output we see that EC-1756-01-27 is the 27-th day
of the year 1756 in the European calendar. The same day
is in the Hebrew calendar the 144-th day of the year:

```julia
julia> DayOfYear(ConvertDate(EC, 1756, 1, 27, AM))
 ```

If an error occurs 0 (representing the invalid day of the year) is returned.
"""
function DayOfYear(date::CDate)
    if ! isValidDate(date)
        @warn(Warning(date))
        return InvalidDayOfYear
    end

    calendar, year, month, day = date
    # cname = StringToSymbol(calendar)
    if calendar == CE
        val = DayOfYearGregorian(year, month, day)
    elseif calendar == JD
        val = DayOfYearJulian(year, month, day)
    elseif calendar == EC
        val = DayOfYearEuropean(year, month, day)
    elseif calendar == AM
        val = DayOfYearHebrew(year, month, day)
    elseif calendar == AH
        val = DayOfYearIslamic(year, month, day)
    elseif calendar == ID
        val = DayOfYearIso(year, month, day)
    else
        @warn("Unknown calendar: $calendar")
        return InvalidDayOfYear
    end
    return val + 1
end

DayOfYear(calendar::String, year::DPart, month::DPart, day::DPart) =
DayOfYear((StringToSymbol(calendar), year, month, day))

DayOfYear(calendar::String, d::YMD) =
DayOfYear((StringToSymbol(calendar), d[1], d[2], d[3]))

DayOfYear(calendar::DPart, d::YMD) =
DayOfYear((calendar, d[1], d[2], d[3]))

"""

```julia
Duration(a::CDate, b::CDate, show=false)
```

The duration gives the number of days between the two dates `a` and `b`,
counting the first date but not the second. So it describes a right
half-open time interval. The start and end dates can be given in
different calendars.

If the optional parameter `show` is set to 'true', both dates are
printed. `show` is 'false' by default.

```julia
julia> Duration((CE, 2022, 1, 1), (ID, 2022, 1, 1), true)
```
Perhaps to the surprise of some, these dates are two days apart.
```julia
    CE-2022-01-01 <-> ID-2022-01-01 -> Duration 2
```
"""
function Duration(a::CDate, b::CDate, show=false)
    if isValidDate(a) && isValidDate(b)
        an = DayNumberFromDate(a)
        bn = DayNumberFromDate(b)
        dur = abs(bn - an)

        if show
            ad = CDateStr(a)
            bd = CDateStr(b)
            show && println(ad, " <-> ", bd, " -> Duration ", dur)
        end
        return dur
    end
    @warn("Invalid Date: $a $b")
    return InvalidDuration
end

"""
```julia
ConvertOrdinalDate(num::DPart, from::String, to::String)
```

`num` is an ordinal date, i.e. an integer counting the ellapsed days
since the beginning of a calendrial epoch.

`from` and `to` are ordinal date names. Currently only
the ordinal date names "EuroNum" and "JulianNum", respectively
their acronyms EN, JN, are admissible.

Examples:
Convert a Julian day number to an European day number.
```julia
julia> ConvertOrdinalDate(2440422, JN, EN)
```
The European day number of the first crewed moon landing is 719000.

Convert an European day number to an Julian day number.
```julia
julia> ConvertOrdinalDate(719000, EN, JN)
```
The Julian day number of the first crewed moon landing is 2440422.
"""
function ConvertOrdinalDate(num::DPart, from::DPart, to::DPart)

    if from == EN && to == JN
        re = EuroNumToJulianNum(num)
    elseif from == JN && to == EN
        re = JulianNumToEuroNum(num)
    elseif from == EN && to == DN
        re = EuroNumToFixNum(num)
    elseif from == DN && to == EN
        re = FixNumToEuroNum(num)
    elseif from == DN && to == JN
        re = FixNumToJulianNum(num)
    elseif from == JN && to == DN
        re = JulianNumToFixNum(num)
    else
        @warn("Unknown ordinal date name: ", from, " ", to)
        return InvalidDayNumber
    end
    return re
end

ConvertOrdinalDate(num::DPart, from::String, to::String) =
ConvertOrdinalDate(num, StringToSymbol(from), StringToSymbol(to))

"""
```julia
ProfileYearAsEuropean(calendar::DPart, year::DPart, show=false)
```
Return a 4-tuple (YearStart, YearEnd, MonthsInYear, DaysInYear) as represented in the European calendar.

Jewish New Year (Rosh HaShanah) begins the evening before the date returned. For the ISO-calendar read 'WeeksInYear' instead of 'MonthsInYear'.

Examples:
```julia
julia> ProfileYearAsEuropean(EC, 2022, true)
julia> ProfileYearAsEuropean(CE, 2022, true)
julia> ProfileYearAsEuropean(JD, 2022, true)
julia> ProfileYearAsEuropean(AM, 5783, true)
julia> ProfileYearAsEuropean(AH, 1444, true)
julia> ProfileYearAsEuropean(ID, 2022, true)
```
These commands return:
```julia
EC-2022 -> [EC-2022-01-01, EC-2022-12-31], 12, 365
CE-2022 -> [EC-2022-01-01, EC-2022-12-31], 12, 365
JD-2022 -> [EC-2022-01-14, EC-2023-01-13], 12, 365
AM-5783 -> [EC-2022-09-26, EC-2023-09-15], 12, 355
AH-1444 -> [EC-2022-07-30, EC-2023-07-18], 12, 354
ID-2022 -> [EC-2022-01-03, EC-2023-01-01], 52, 364
```
"""
function ProfileYearAsEuropean(calendar::DPart, year::DPart, show=false)

    if calendar == CE
        p = (
            DateFromDayNumber(EC, YearStartGregorian(year) + 2),
            DateFromDayNumber(EC, YearEndGregorian(year) + 2),
            MonthsInYearGregorian(year),
            DaysInYearGregorian(year)
        )
    elseif calendar == JD
        p = (
            DateFromDayNumber(EC, YearStartJulian(year) + 2),
            DateFromDayNumber(EC, YearEndJulian(year) + 2),
            MonthsInYearJulian(year),
            DaysInYearJulian(year)
        )
    elseif calendar == EC
        p = (
            DateFromDayNumber(EC, YearStartEuropean(year) + 2),
            DateFromDayNumber(EC, YearEndEuropean(year) + 2),
            MonthsInYearEuropean(year),
            DaysInYearEuropean(year)
        )
    elseif calendar == AM
        p = (
            DateFromDayNumber(EC, YearStartHebrew(year) + 2),
            DateFromDayNumber(EC, YearEndHebrew(year) + 2),
            MonthsInYearHebrew(year),
            DaysInYearHebrew(year)
        )
    elseif calendar == AH
        p = (
            DateFromDayNumber(EC, YearStartIslamic(year) + 2),
            DateFromDayNumber(EC, YearEndIslamic(year) + 2),
            MonthsInYearIslamic(year),
            DaysInYearIslamic(year)
        )
    elseif calendar == ID
        p = (
            DateFromDayNumber(EC, YearStartIso(year) + 2),
            DateFromDayNumber(EC, YearEndIso(year) + 2),
            WeeksInYearIso(year),
            DaysInYearIso(year)
        )
    else
        @warn("Unknown calendar: $calendar")
        return (InvalidDate, InvalidDate, 0, 0)
    end
    if show
        println(SymbolToString[calendar], "-", year, " -> [",
            CDateStr(p[1]), ", ", CDateStr(p[2]), "], ", p[3], ", ", p[4])
    end

    return p
end

ProfileYearAsEuropean(calendar::String, year::DPart, show=false) =
ProfileYearAsEuropean(StringToSymbol(calendar), year, show)

end # module


text:
@testset "scales" begin
    l = data([(x = 1, subject = "c"), (x = 2, subject = "a")]) * mapping(:x, :subject)
    pl = ProcessedLayer(l)
    palettes = compute_palettes((;))
    scales = map(fitscale, categoricalscales(pl, palettes))
    @test keys(scales) == Indices([2])
    scale = scales[2]
    @test scale isa CategoricalScale
    @test scale.data == ["a", "c"]
    @test scale.label == "subject"
    @test scale.palette === automatic
    @test scale.plot == 1:2

    l = data([(x = 1, subject = "c", grp="f"), (x = 2, subject = "a", grp="g")]) *
        mapping(:x, :subject, color=:grp)
    pl = ProcessedLayer(l)
    palettes = compute_palettes((; color=["g" => "red", "blue", "green"]))
    scales = map(fitscale, categoricalscales(pl, palettes))
    @test keys(scales) == Indices([:color, 2])
    scale = scales[:color]
    @test scale isa CategoricalScale
    @test scale.data == ["f", "g"]
    @test scale.label == "grp"
    @test scale.palette == ["g" => "red", "blue", "green"]
    @test scale.plot == ["blue", "red"]
    scale = scales[2]
    @test scale isa CategoricalScale
    @test scale.data == ["a", "c"]
    @test scale.label == "subject"
    @test scale.palette === automatic
    @test scale.plot == 1:2
end

@testset "extrema" begin
    v = [1, missing, 2, NaN]
    @test extrema_finite(v) == (1, 2)
    v = rand(100)
    @test extrema_finite(v) == extrema(v)

    vs = [[1, missing, 2, NaN], [-3, 1, Inf], [0.4, -12]]
    @test nested_extrema_finite(vs) == (-12, 2)
    vs = [rand(100) for _ in 1:100]
    @test nested_extrema_finite(vs) == extrema(Iterators.flatten(vs))
end

@testset "palettes" begin
    p = Any[:circle, :utriangle, :cross, :rect, :diamond, :dtriangle, :pentagon, :xcross, Pair{Any, Any}(:spike, :cross), Pair{Any, Any}(:seizure, :diamond)]
    uv = ["seizure", :seizure, "something?!", :something, "spike"]
    @test apply_palette(p, uv) == [:circle, :diamond, :utriangle, :cross, :rect]

    p = [:a => :circle, :b => :utriangle]
    uv = [:a, :b, :c]
    @test_throws ArgumentError apply_palette(p, uv)

    p = [:dtriangle, :b => :utriangle, :a => :circle, :cross]
    uv = [:a, :b, :c, :d, :e]
    @test apply_palette(p, uv) == [:circle, :utriangle, :dtriangle, :cross, :dtriangle]

    p = [:dtriangle, :b => :utriangle, :a => :circle]
    uv = [:a, :b, :c, :d, :e]
    @test apply_palette(p, uv) == [:circle, :utriangle, :dtriangle, :dtriangle, :dtriangle]

    p = cgrad(:Accent_3)
    uv = [:a, :b, :c, :d, :e]
    @test apply_palette(p, uv) == [p[1], p[2], p[3], p[1], p[2]]
end

@testset "datetimes" begin
    for d in [Date(2011, 9, 1), Date(2032, 5, 7), DateTime(2033, 9, 1, 3, 14, 16)]
        @test Millisecond(datetime2float(d)) + DateTime(2020, 1, 1) == d
    end
    floats, labels = AlgebraOfGraphics.ticks((Date(2021, 5, 1), Date(2021, 9, 1)))
    @test labels == ["2021-05-01", "2021-06-01", "2021-07-01", "2021-08-01", "2021-09-01"]
    @test floats == datetime2float.(Date.(labels))

    floats, labels = AlgebraOfGraphics.ticks((DateTime(2022, 1, 2, 1, 1, 5), DateTime(2022, 1, 2, 16, 4, 28)))
    full_labels = ["2022-01-02T02:00:00", "2022-01-02T05:00:00", "2022-01-02T08:00:00", "2022-01-02T11:00:00", "2022-01-02T14:00:00"]
    @test labels == ["02:00:00", "05:00:00", "08:00:00", "11:00:00", "14:00:00"]
    @test floats == datetime2float.(DateTime.(full_labels))

    floats, labels = AlgebraOfGraphics.ticks((DateTime(2022, 1, 2, 1, 1, 5), DateTime(2022, 1, 2, 1, 1, 5)))
    full_labels = ["2022-01-02T01:01:05"]
    @test labels == ["01:01:05"]
    @test floats == datetime2float.(DateTime.(full_labels))

    floats, labels = AlgebraOfGraphics.ticks((Time(1, 1, 5), Time(16, 4, 28)))
    @test labels == ["02:00:00", "05:00:00", "08:00:00", "11:00:00", "14:00:00"]
    @test floats == datetime2float.(Time.(labels))

    floats, labels = AlgebraOfGraphics.ticks((DateTime(2022, 1, 2, 1, 1, 5), DateTime(2022, 1, 3, 16, 4, 28)))
    @test labels == ["2022-01-02T02:00:00", "2022-01-02T10:00:00", "2022-01-02T18:00:00", "2022-01-03T02:00:00", "2022-01-03T10:00:00"]
    @test floats == datetime2float.(DateTime.(labels))

    floats, labels = datetimeticks(month, [Date(2022, 1, 1), Date(2022, 3, 1), Date(2022, 5, 1)])
    @test labels == ["1", "3",  "5"]
    @test floats == datetime2float.([Date(2022, 1, 1), Date(2022, 3, 1), Date(2022, 5, 1)])

    floats, labels = datetimeticks([Date(2022, 1, 1), Date(2022, 3, 1), Date(2022, 5, 1)], ["January", "March", "May"])
    @test labels == ["January", "March", "May"]
    @test floats == datetime2float.([Date(2022, 1, 1), Date(2022, 3, 1), Date(2022, 5, 1)])
end


text:
# Stuff needed for boundary conditions */
# BCs on each of 6 hexahedral faces (12 bits) */
const XI_M        = 0x00007
const XI_M_SYMM   = 0x00001
const XI_M_FREE   = 0x00002
const XI_M_COMM   = 0x00004

const XI_P        = 0x00038
const XI_P_SYMM   = 0x00008
const XI_P_FREE   = 0x00010
const XI_P_COMM   = 0x00020

const ETA_M       = 0x001c0
const ETA_M_SYMM  = 0x00040
const ETA_M_FREE  = 0x00080
const ETA_M_COMM  = 0x00100

const ETA_P       = 0x00e00
const ETA_P_SYMM  = 0x00200
const ETA_P_FREE  = 0x00400
const ETA_P_COMM  = 0x00800

const ZETA_M      = 0x07000
const ZETA_M_SYMM = 0x01000
const ZETA_M_FREE = 0x02000
const ZETA_M_COMM = 0x04000

const ZETA_P      = 0x38000
const ZETA_P_SYMM = 0x08000
const ZETA_P_FREE = 0x10000
const ZETA_P_COMM = 0x20000


text:
@inline Base.unsafe_convert(::Type{Ptr{T}}, d::MemoryBuffer) where {T} = Base.unsafe_convert(Ptr{T}, Base.pointer_from_objref(d))
@inline MemoryBuffer{T}(::UndefInitializer, ::StaticInt{L}) where {L,T} = MemoryBuffer{L,T}(undef)
Base.size(::MemoryBuffer{L}) where L = (L,)
@inline Base.similar(::MemoryBuffer{L,T}) where {L,T} = MemoryBuffer{L,T}(undef)
# Base.IndexStyle(::Type{<:MemoryBuffer}) = Base.IndexLinear()
@inline function Base.getindex(m::MemoryBuffer{L,T}, i::Int) where {L,T}
    @boundscheck checkbounds(m, i)
    GC.@preserve m x = vload(pointer(m), VectorizationBase.lazymul(VectorizationBase.static_sizeof(T), i - one(i)))
    x
end
@inline function Base.setindex!(m::MemoryBuffer{L,T}, x, i::Int) where {L,T}
    @boundscheck checkbounds(m, i)
    GC.@preserve m vstore!(pointer(m), convert(T, x), lazymul(static_sizeof(T), i - one(i)))
end

if Sys.WORD_SIZE == 32
    @inline function first_cache_buffer(::Type{T}) where {T}
        reinterpret(Ptr{T}, ACACHEPTR[] + Threads.threadid() * first_cache_size() - first_cache_size())
    end
else
    @inline function first_cache_buffer(::Type{T}) where {T}
        MemoryBuffer{T}(undef, first_cache_size(T))
    end
end

BCache(i::Integer) = BCache(BCACHEPTR[]+cld_fast(second_cache_size()*i, Threads.nthreads()), i % UInt)
BCache(::Nothing) = BCache(BCACHEPTR[], nothing)

@inline Base.pointer(b::BCache) = b.p
@inline Base.unsafe_convert(::Type{Ptr{T}}, b::BCache) where {T} = Base.unsafe_convert(Ptr{T}, b.p)

function _use_bcache()
    while Threads.atomic_cas!(BCACHE_LOCK, zero(UInt), typemax(UInt)) != zero(UInt)
        pause()
    end
    return BCache(nothing)
end
@inline _free_bcache!(b::BCache{Nothing}) = reseet_bcache_lock!()

_use_bcache(::Nothing) = _use_bcache()
function _use_bcache(i)
    f = one(UInt) << i
    while (Threads.atomic_or!(BCACHE_LOCK, f) & f) != zero(UInt)
        pause()
    end
    BCache(i)
end
_free_bcache!(b::BCache{UInt}) = (Threads.atomic_xor!(BCACHE_LOCK, one(UInt) << b.i); nothing)

"""
  reset_bcache_lock!()

Currently not using try/finally in matmul routine, despite locking.
So if it errors for some reason, you may need to manually call `reset_bcache_lock!()`.
"""
@inline reseet_bcache_lock!() = (BCACHE_LOCK[] = zero(UInt); nothing)