text:
@doc raw"""
    AffTransformation2

Represents two-dimensional affine transformations.

The general form of an affine transformation is based on a homogeneous
representation of points. Thereby all transformations can be realized by matrix
multiplications.

Multiplying the transformation matrix by a scalar does not change the
represented transformation. Therefore, any transformation represented by a
matrix with rational entries can be represented by a transformation matrix with
integer entries as well. (Multiply the matrix with the common denominator of
the rational entries.) Hence, it is sufficient to use the number type
[`RT`](@ref) to represent the entries of the transformation matrix.

CGAL offers several specialized affine transformations. Different constructors
are provided to create them. They are paremeterized with a symbolic name to
denote the transformation type, followed by additional parameters. The symbolic
name tags solve ambiguities in the function overloading and they make the code
more readable, i.e., what type of transformation is created.

Since two-dimensional points have three homogeneous coordinates, we have a
``3 × 3`` matrix ``(m_{ij})_{i,j=0…2}``.

If the homogeneous representations are normalized (the homogenizing coordinate
is 1), then the upper left ``2 × 2`` matrix realizes linear transformations. In
the matrix form of a translation, the translation vector ``(v_0,\,v_1,\,1)``
appears in the last column of the matrix. The entries ``m_{20}`` and ``m_{21}``
are always zero and therefore do not appear in the constructors.

# Examples
```julia-repl
julia> for (A, T) ∈ ((:Transformation, AffTransformation2),
                     (:Point, Point2),
                     (:Vector, Vector2),
                     (:Direction, Direction2))
           @eval const $A = $T
       end

julia> rotate = Transformation(ROTATION, sin(π), cos(π))
Aff_TransformationC2(0, -1)

julia> rational_rotate = Transformation(ROTATION, Direction(1,1), 1, 100)
Aff_TransformationC2(0.710059, 0.704142)

julia> translate = Transformation(TRANSLATION, Vector(-2, 0))
Aff_TransformationC2(VectorC2(-2, 0))

julia> scale = Transformation(SCALING, 3)
Aff_TransformationC2(3)

julia> q = Point(0, 1)
PointC2(0, 1)

julia> q = rational_rotate(q)
PointC2(-0.710059, 0.704142)

julia> p = Point(1, 1)
PointC2(1, 1)

julia> p = rotate(p)
PointC2(-1, -1)

julia> p = translate(p)
PointC2(-3, -1)

julia> p = scale(p)
PointC2(-9, -3)
```

The same would have been achieved with

```julia-repl
julia> transform = scale * (translate * rotate)
Aff_TransformationC2(-3 0 -6
                     0 -2 0)

julia> p = transform(Point(1.0, 1.0))
PointC2(-9, -3)
```

See also: [`IdentityTransformation`](@ref), [`Rotation`](@ref),
[`Scaling`](@ref), [`Translation`](@ref), [`Reflection`](@ref),
[`rational_rotation_approximation`](@ref)
"""
AffTransformation2

"""
    AffTransformation2(::IdentityTransformation)

Introduces an identity transformation.
"""
AffTransformation2(::IdentityTransformation)

"""
    AffTransformation2(::Translation, v::Vector2)

Introduces a translation by vector `v`.
"""
AffTransformation2(::Translation, v::Vector2)

"""
    AffTransformation2(::Rotation, d::Direction2, num::Real, den::Real = 1)

Approximates the rotation over the angle indicated by direction `d`, such that
the differences between the sines and cosines of the rotation given by `d` and
the approximating rotation are at most `num/den` each.

!!! info "Precondition"

    ``num/den>0`` and ``d ≠ 0``.
"""
AffTransformation2(::Rotation, d::Direction2, num::Real, den::Real = 1)

@cxxdereference AffTransformation2(R::Union{Rotation,Ref{Rotation}}, d::Direction2, num::Real, den::Real = 1) =
    AffTransformation2(R, d, convert(FT, num), convert(FT, den))

@doc raw"""
    AffTransformation2(::Rotation, sineᵨ::Real, cosineᵨ::Real, hw::Real = 1)

Introduces a rotation by the angle `ρ`.

!!! info "Precondition"

    ``sine_ρ^2 + cosine_ρ^2 = hw^2``.
"""
AffTransformation2(::Rotation, sinᵨ::Real, cosᵨ::Real, hw::Real = 1)

AffTransformation2(R::Union{Rotation,Ref{Rotation}}, sinᵨ::Real, cosᵨ::Real, hw::Real = 1) =
    AffTransformation2(R, convert(FT, sinᵨ), convert(FT, cosᵨ), convert(FT, hw))

"""
    AffTransformation2(::Scaling, s::Real, hw::Real = 1)

Introduces a scaling by a scale factor ``s/hw``.
"""
AffTransformation2(::Scaling, s::Real, hw::Real = 1)

AffTransformation2(S::Union{Scaling,Ref{Scaling}}, s::Real, hw::Real = 1) =
    AffTransformation2(S, convert(FT, s), convert(FT, hw))

"""
    AffTransformation2(::Reflection, l::Line2)

Introduces a reflection by a line ``l``.
"""
AffTransformation2(::Reflection, l::Line2)

@doc raw"""
    AffTransformation2(m₀₀::Real, m₀₁::Real, m₀₂::Real,
                       m₁₀::Real, m₁₁::Real, m₁₂::Real,
                                                        hw::Real = 1)

Introduces a general affine transformation in the ``3 × 3`` matrix form ``
\small\left(\begin{array}{ccc}
    m_{00} & m_{01} & m_{02} \\ m_{10} & m_{11} & m_{12} \\ 0 & 0 & hw
\end{array}\right)
``.

The sub-matrix ``
\frac{1}{hw}\small\left(\begin{array}{cc}
    m_{00} & m_{01} \\ m_{10} & m_{11}
\end{array}\right)
``
contains the scaling and rotation information, the vector ``
\small\left(\begin{array}{c}
    m_{02}\\ m_{12}
\end{array}\right)
``
contains the translational part of the transformation.
"""
AffTransformation2(m₀₀::Real, m₀₁::Real, m₀₂::Real,
                   m₁₀::Real, m₁₁::Real, m₁₂::Real,
                                                    hw::Real = 1) =
    AffTransformation2(convert.(FT, (m₀₀, m₀₁, m₀₂,
                                     m₁₀, m₁₁, m₁₂,
                                                    hw))...)

@doc raw"""
    AffTransformation2(m₀₀::Real, m₀₁::Real,
                       m₁₀::Real, m₁₁::Real,
                                             hw::Real = 1)

Introduces a general linear transformation ``
\small\left(\begin{array}{ccc}
    m_{00} & m_{01} & 0 \\ m_{10} & m_{11} & 0 \\ 0 & 0 & hw
\end{array}\right)
``,
i.e. there is no translational part.
"""
AffTransformation2(m₀₀::Real, m₀₁::Real,
                   m₁₀::Real, m₁₁::Real,
                                         hw::Real = 1) =
    AffTransformation2(convert.(FT, (m₀₀, m₀₁, m₁₀, m₁₁, hw))...)

"""
    transform(t::AffTransformation2, x::Union{Point2,Vector2,Direction2,Line2})
    (::AffTransformation2)(x::Union{Point2,Vector2,Direction2,Line2})

The [`transform`](@ref) function has methods for points, vectors, directions,
and lines. The same functionality is available through the invocation of
[`AffTransformation2`](@ref) instances.
"""
transform(t::AffTransformation2, p::Union{Point2,Vector2,Direction2,Line2}),
(::AffTransformation2)(p::Union{Point2,Vector2,Direction2,Line2})

"""
    *(t::AffTransformation2, s::AffTransformation2)

Composes two affine transformations.
"""
*(t::AffTransformation2, s::AffTransformation2)

"""
    inverse(t::AffTransformation2)

Gives the inverse transformation.
"""
inverse(t::AffTransformation2)

"""
    ==(t::AffTransformation2, s::AffTransformation2)

Compares two affine transformations.
"""
==(t::AffTransformation2, s::AffTransformation2)

"""
    is_even(t::AffTransformation2)

Returns `true`, if the transformation is not reflecting, i.e., the determinant
of the involved linear transformation is non-negative.
"""
is_even(t::AffTransformation2)

"""
    is_odd(t::AffTransformation2)

Returns `true`, if the transformation is reflecting.
"""
is_odd(t::AffTransformation2)

"""
    is_scaling(t::AffTransformation2)

Returns `true`, if the object was constructed using the tag
[`CGAL.Scaling`](@ref), or is the result of the composition of only such scaling
transformation objects.
"""
is_scaling(t::AffTransformation2)

"""
    is_translation(t::AffTransformation2)

Returns `true`, if the object was constructed using the tag
[`CGAL.Translation`](@ref), or is the result of the composition of only such
translation transformation objects.
"""
is_translation(t::AffTransformation2)

"""
    is_rotation(t::AffTransformation2)

Returns `true`, if the object was constructed using the tag
[`CGAL.Rotation`](@ref), or is the result of the composition of only such
rotation transformation objects.
"""
is_rotation(t::AffTransformation2)

"""
    is_reflection(t::AffTransformation2)

Returns `true`, if the object was constructed using the tag
[`CGAL.Reflection`](@ref), or is the result of the composition of only such
reflection transformation objects.
"""
is_reflection(t::AffTransformation2)

"""
    m(t::AffTransformation2, i::Integer, j::Integer)

Returns entry ``m_{ij}`` in a matrix representation in which ``m_{22}`` is 1.
"""
m(t::AffTransformation2, i::Integer, j::Integer)

"""
    hm(t::AffTransformation2, i::Integer, j::Integer)

Returns entry ``m_{ij}`` in some fixed matrix representation.
"""
hm(t::AffTransformation2, i::Integer, j::Integer)


text:
#!/usr/bin/env julia

using SailRoute, Test
@testset "Discretization error functions" begin
    value_1 = 0.970500
    value_2 = 0.968540
    value_3 = 0.961780
    h1 = 1.0
    h2 = 2.0
    h3 = 4.0
    ratio_21 = h2/h1
    ratio_32 = h3/h2
    ooc = SailRoute.order_of_convergence(value_1, value_2, value_3, ratio_21, ratio_32)
    @test ooc ≈ 1.786170 atol=4
    f_exact_21 = SailRoute.richardson_extrapolate(value_1, value_2, ratio_21, ooc)
    f_exact_32 = SailRoute.richardson_extrapolate(value_2, value_3, ratio_32, ooc)
    e21a, e21ext = SailRoute.error_estimates(value_1, value_2, f_exact_21)
    @test e21a ≈ 0.002020 atol = 5
    @test e21ext ≈ 0.000824 atol = 5
    e32a, e32ext = SailRoute.error_estimates(value_2, value_3, f_exact_32)
    gci_fine_21, gci_coarse_21 = SailRoute.gci(ratio_21, e21a, ooc)
    gci_fine_32, gci_coarse_32 = SailRoute.gci(ratio_32, e32a, ooc)
    ratio = SailRoute.asymptotic_ratio(gci_fine_21, gci_fine_32, ratio_21, ooc)
    @test ratio ≈ 0.997980 atol = 5
    # @test SailRoute.perf_interp(test_performance, 60.0, 10.0) ≈ 5.0
    # @test SailRoute.perf_interp(test_performance, 60.0, 15.0) ≈ 7.5
    # @test SailRoute.perf_interp(test_performance, 30.0, 10.0) ≈ 0.0
end

ext: jl


text:
"""Struct that holds the quantised array as UInts with an additional
field for the min, max of the original range."""
struct LogQuantArray{T,N} <: AbstractArray{Unsigned,N}
    A::Array{T,N}       # array of UInts
    min::Float64        # min of value range
    max::Float64        # max of value range
end

Base.size(QA::LogQuantArray) = size(QA.A)
Base.getindex(QA::LogQuantArray,i...) = getindex(QA.A,i...)
Base.eltype(Q::LogQuantArray{T,N}) where {T,N} = T

"""Return minimum, ignoring negatives and zeroes."""
function minpos(A::AbstractArray{T}) where T
    o = zero(T)
    mi = foldl((x,y) -> y > o ? min(x,y) : x, A; init=typemax(T))
    mi == typemax(T) && return o
    return mi
end

"""Quantize elements of an array logarithmically into UInts with
either round to nearest in linear or Logarithmic space."""
function LogQuantization(   ::Type{T},
                            A::AbstractArray,
                            round_nearest_in::Symbol=:linspace) where {T<:Unsigned}

    (any(A .< zero(eltype(A))) || ~all(isfinite.(A))) &&
        throw(DomainError("Logarithmic  quantization only for positive&zero entries."))

    # min/max of non-zero entries
    mi = Float64(minpos(A))
    logmin = log(Float64(mi))
    logmax = log(Float64(maximum(A)))

    # throw error in case the range is zero.
    if logmin == logmax
        Δ = 0.0
        c = 0.0
    else
        # inverse log spacing
        # map min to 1 and max to ff..., reserve 0 for 0.
        Δ = (2^(sizeof(T)*8)-2)/(logmax-logmin)

        # shift to round-to-nearest in lin or log-space
        if round_nearest_in == :linspace
            c = 1/2 - Δ*log(mi*(exp(1/Δ)+1)/2)
        elseif round_nearest_in == :logspace
            c = -logmin*Δ
        else
            throw(ArgumentError("Round-to-nearest either :linspace or :logspace"))
        end
    end

    # preallocate output
    Q = similar(A,T)

    @inbounds for i in eachindex(A)
        # store 0 as 0x00...
        # store positive numbers via convert to logpacking as 0x1-0xff..
        Q[i] = iszero(A[i]) ? zero(T) : convert(T,round(c + Δ*log(Float64(A[i]))))+one(T)
    end

    return LogQuantArray{T,ndims(Q)}(Q,Float64(logmin),Float64(logmax))
end

# define for 8, 16, 24 and 32 bit uints
LogQuant8Array(A::AbstractArray{T,N},rn::Symbol=:linspace) where {T,N} = LogQuantization(UInt8,A,rn)
LogQuant16Array(A::AbstractArray{T,N},rn::Symbol=:linspace) where {T,N} = LogQuantization(UInt16,A,rn)
LogQuant24Array(A::AbstractArray{T,N},rn::Symbol=:linspace) where {T,N} = LogQuantization(UInt24,A,rn)
LogQuant32Array(A::AbstractArray{T,N},rn::Symbol=:linspace) where {T,N} = LogQuantization(UInt32,A,rn)

"""De-quantise a LogQuantArray into floats."""
function Base.Array{T}(n::Integer,Q::LogQuantArray) where {T<:AbstractFloat}
    Qlogmin = Q.min                 # log(min::Float64)
    Qlogmax = Q.max                 # log(max::Float64)

    # spacing in logspace ::Float64
    Δ = (Qlogmax-Qlogmin)/(2^n-2)   # -2 as 0x00.. is reserved for 0

    A = similar(Q,T)                # preallocate

    @inbounds for i in eachindex(A)
        # 0x0 is unpack as 0
        # exp in Float64 then convert to T at assignment =
        A[i] = iszero(Q[i]) ? zero(T) : A[i] = exp(Qlogmin + (Q[i]-1)*Δ)
    end

    return A
end

Base.Array{T}(Q::LogQuantArray{UInt8,N}) where {T,N} = Array{T}(8,Q)
Base.Array{T}(Q::LogQuantArray{UInt16,N}) where {T,N} = Array{T}(16,Q)
Base.Array{T}(Q::LogQuantArray{UInt24,N}) where {T,N} = Array{T}(24,Q)
Base.Array{T}(Q::LogQuantArray{UInt32,N}) where {T,N} = Array{T}(32,Q)

Base.Array(Q::LogQuantArray{UInt8,N}) where N = Array{Float32}(8,Q)
Base.Array(Q::LogQuantArray{UInt16,N}) where N = Array{Float32}(16,Q)
Base.Array(Q::LogQuantArray{UInt24,N}) where N = Array{Float32}(24,Q)
Base.Array(Q::LogQuantArray{UInt32,N}) where N = Array{Float64}(32,Q)

# one quantization per layer
"""Logarithmic quantization independently for every element along dimension
dim in array A. Returns a Vector{LogQuantArray}."""
function LogQuantArray(::Type{TUInt},A::AbstractArray{T,N},dim::Int) where {TUInt,T,N}
    @assert dim <= N   "Can't quantize a $N-dimensional array in dim=$dim"
    n = size(A)[dim]
    L = Vector{LogQuantArray}(undef,n)
    t = [if j == dim 1 else Colon() end for j in 1:N]
    for i in 1:n
        t[dim] = i
        L[i] = LogQuantization(TUInt,A[t...])
    end
    return L
end

# for 8,16,24 and 32 bit
LogQuant8Array(A::AbstractArray{T,N},dim::Int) where {T,N} = LogQuantArray(UInt8,A,dim)
LogQuant16Array(A::AbstractArray{T,N},dim::Int) where {T,N} = LogQuantArray(UInt16,A,dim)
LogQuant24Array(A::AbstractArray{T,N},dim::Int) where {T,N} = LogQuantArray(UInt24,A,dim)
LogQuant32Array(A::AbstractArray{T,N},dim::Int) where {T,N} = LogQuantArray(UInt32,A,dim)

"""Undo the logarithmic quantisation independently along one dimension, and returns
an array whereby the dimension always comes last. Hence, might be permuted compared
to the uncompressed array."""
function Base.Array{T}(L::Vector{LogQuantArray}) where T
    N = ndims(L[1])
    n = length(L)
    s = size(L[1])
    t = axes(L[1])
    A = Array{T,N+1}(undef,s...,length(L))
    for i in 1:n
        A[t...,i] = Array{T}(L[i])
    end
    return A
end


text:
export SampleList, SampleListMeta

import Base: show, ndims, length, size, precision, getindex, broadcasted, map
import Distributions: mean, var, cov, std
import StatsBase: Weights

using StaticArrays
using LoopVectorization

abstract type AbstractSampleListSamplingMethod end

struct BootstrapImportanceSampling <: AbstractSampleListSamplingMethod end

mutable struct SampleListCache{M, C}
    mean :: M
    cov  :: C
    is_mean_cached :: Bool
    is_cov_cached  :: Bool
end

SampleListCache(::Type{T}, dims::Tuple{})         where T = SampleListCache(zero(T), zero(T), false, false)
SampleListCache(::Type{T}, dims::Tuple{Int})      where T = SampleListCache(zeros(T, first(dims)), zeros(T, first(dims), first(dims)), false, false)
SampleListCache(::Type{T}, dims::Tuple{Int, Int}) where T = SampleListCache(zeros(T, dims), zeros(T, prod(dims), prod(dims)), false, false)

is_mean_cached(cache::SampleListCache) = cache.is_mean_cached
is_cov_cached(cache::SampleListCache)  = cache.is_cov_cached

get_mean_storage(cache::SampleListCache) = cache.mean
get_cov_storage(cache::SampleListCache)  = cache.cov

cache_mean!(cache::SampleListCache, mean) = begin cache.mean = mean; cache.is_mean_cached = true; mean end
cache_cov!(cache::SampleListCache, cov)   = begin cache.cov  = cov; cache.is_cov_cached = true; cov end

struct SampleListMeta{W, E, LP, LI}
    unnormalisedweights :: W
    entropy             :: E
    logproposal         :: LP
    logintegrand        :: LI
end

get_unnormalised_weights(meta::SampleListMeta) = meta.unnormalisedweights
get_entropy(meta::SampleListMeta)              = meta.entropy
get_logproposal(meta::SampleListMeta)          = meta.logproposal
get_logintegrand(meta::SampleListMeta)         = meta.logintegrand

call_logproposal(logproposal::Function, x) = logproposal(x)
call_logproposal(logproposal::Any, x)      = logpdf(logproposal, x)

call_logintegrand(logintegrand::Function, x) = logintegrand(x)
call_logintegrand(logintegrand::Any, x)      = logpdf(logintegrand, x)

"""
    SampleList

Generic distribution represented as a list of weighted samples.

# Arguments
- `samples::S`
- `weights::W`: optional, equivalent to `fill(1 / N, N)` by default, where `N` is the length of `samples` container
"""
struct SampleList{D, S, W, C, M}
    samples :: S
    weights :: W
    cache   :: C
    meta    :: M

    function SampleList(::Val{D}, samples::S, weights::W, meta::M = nothing) where { D, S, W, M }
        @assert div(length(samples), prod(D)) === length(weights) "Invalid sample list samples and weights lengths. `samples` has length $(length(samples)), `weights` has length $(length(weights))"
        @assert eltype(samples) <: Number "Invalid eltype of samples container. Should be a subtype of `Number`, but $(eltype(samples)) has been found. Samples should be stored in a linear one dimensional vector even for multivariate and matrixvariate cases."
        @assert eltype(weights) <: Number "Invalid eltype of weights container. Should be a subtype of `Number`, but $(eltype(weights)) has been found."
        cache = SampleListCache(promote_type(eltype(samples), eltype(weights)), D)
        return new{D, S, W, typeof(cache), M}(samples, weights, cache, meta)
    end
end

Base.show(io::IO, sl::SampleList)         = sample_list_show(io, variate_form(sl), sl)
Base.similar(sl::SampleList{ D }) where D = SampleList(Val(D), similar(sl.samples), similar(sl.weights))

sample_list_show(io::IO, ::Type{Univariate}, sl::SampleList) = print(io, "SampleList(Univariate, ", length(sl), ")")
sample_list_show(io::IO, ::Type{Multivariate}, sl::SampleList) = print(io, "SampleList(Multivariate(", ndims(sl),"), ", length(sl), ")")
sample_list_show(io::IO, ::Type{Matrixvariate}, sl::SampleList) = print(io, "SampleList(Matrixvariate", ndims(sl), ", ", length(sl), ")")

function SampleList(samples::S) where { S <: AbstractVector }
    N = length(samples)
    return SampleList(samples, fill(one(deep_eltype(S)) / N, N))
end

function SampleList(samples::S, weights::W, meta::M = nothing) where { S, W, M }
    nsamples = length(samples)
    @assert nsamples !== 0 "Empty samples list"
    @assert sum(weights) ≈ one(eltype(weights)) "Weights must sum up to one. sum(weights) = $(sum(weights))"
    D = size(first(samples))
    return SampleList(Val(D), sample_list_linearize(samples, nsamples, prod(D)), weights, meta)
end

const DEFAULT_SAMPLE_LIST_N_SAMPLES = 5000

## Utility functions

Base.eltype(::Type{ <: SampleList{D, S, W} }) where { D, S, W } = Tuple{ sample_list_eltype(SampleList, D, S), eltype(W) }

sample_list_eltype(::Type{ SampleList }, ndims::Tuple{}, ::Type{S}) where S         = eltype(S)
sample_list_eltype(::Type{ SampleList }, ndims::Tuple{Int}, ::Type{S}) where S      = SVector{ ndims[1], eltype(S) }
sample_list_eltype(::Type{ SampleList }, ndims::Tuple{Int, Int}, ::Type{S}) where S = SMatrix{ ndims[1], ndims[2], eltype(S), ndims[1] * ndims[2] }

deep_eltype(::Type{ <: SampleList{ D, S } }) where { D, S } = eltype(S)

## Variate forms

variate_form(::SampleList{ D }) where { D } = sample_list_variate_form(D)

sample_list_variate_form(::Tuple{})         = Univariate
sample_list_variate_form(::Tuple{Int})      = Multivariate
sample_list_variate_form(::Tuple{Int, Int}) = Matrixvariate

## Getters

get_weights(sl::SampleList) = get_linear_weights(sl)
get_samples(sl::SampleList) = SamplesOnlyIterator(sl)

get_linear_weights(sl::SampleList) = sl.weights
get_linear_samples(sl::SampleList) = sl.samples
get_cache(sl::SampleList)          = sl.cache
get_meta(sl::SampleList)           = sample_list_check_meta(sl.meta)
is_meta_present(sl::SampleList)    = sl.meta !== nothing

get_data(sl::SampleList) = (length(sl), get_linear_samples(sl), get_linear_weights(sl))

sample_list_check_meta(meta::Any)     = meta
sample_list_check_meta(meta::Nothing) = error("SampleList object has not associated meta information with it.")

get_unnormalised_weights(sl::SampleList) = get_unnormalised_weights(get_meta(sl))
get_entropy(sl::SampleList)              = get_entropy(get_meta(sl))
get_logproposal(sl::SampleList)          = get_logproposal(get_meta(sl))
get_logintegrand(sl::SampleList)         = get_logintegrand(get_meta(sl))

call_logproposal(sl::SampleList, x)  = call_logproposal(get_logproposal(sl), x)
call_logintegrand(sl::SampleList, x) = call_logintegrand(get_logintegrand(sl), x)

Base.length(sl::SampleList) = div(length(get_linear_samples(sl)), prod(ndims(sl)))
Base.ndims(sl::SampleList)  = sample_list_ndims(variate_form(sl), sl)
Base.size(sl::SampleList)   = (length(sl), )

sample_list_ndims(::Type{ Univariate }, sl::SampleList{ D }) where { D }    = 1
sample_list_ndims(::Type{ Multivariate }, sl::SampleList{ D }) where { D }  = first(D)
sample_list_ndims(::Type{ Matrixvariate }, sl::SampleList{ D }) where { D } = D

## Statistics

# Returns a zeroed container for mean
function sample_list_zero_element(sl::SampleList)
    T = promote_type(eltype(get_linear_weights(sl)), eltype(get_linear_samples(sl)))
    return sample_list_zero_element(variate_form(sl), T, sl)
end

sample_list_zero_element(::Type{ Univariate }, ::Type{T}, sl::SampleList) where T    = zero(T)
sample_list_zero_element(::Type{ Multivariate }, ::Type{T}, sl::SampleList) where T  = zeros(T, ndims(sl))
sample_list_zero_element(::Type{ Matrixvariate }, ::Type{T}, sl::SampleList) where T = zeros(T, ndims(sl))

# Generic mean_cov

mean_cov(sl::SampleList) = sample_list_mean_cov(sl, Val(true))
mean_var(sl::SampleList) = sample_list_mean_var(variate_form(sl), sl)

##

Distributions.mean(sl::SampleList)      = sample_list_mean(sl, Val(true))
Distributions.var(sl::SampleList)       = last(mean_var(sl))
Distributions.cov(sl::SampleList)       = last(mean_cov(sl))
Distributions.invcov(sl::SampleList)    = cholinv(cov(sl))
Distributions.std(sl::SampleList)       = cholsqrt(cov(sl))
Distributions.logdetcov(sl::SampleList) = logdet(cov(sl))

Base.precision(sl::SampleList) = invcov(sl)

function mean_precision(sl::SampleList)
    μ, Σ = mean_cov(sl)
    return μ, cholinv(Σ)
end

function weightedmean_precision(sl::SampleList)
    μ, Λ = mean_precision(sl)
    return Λ * μ, Λ
end

weightedmean(sl::SampleList)    = first(weightedmean_precision(sl))
logmean(sl::SampleList)         = sample_list_logmean(variate_form(sl), sl)
meanlogmean(sl::SampleList)     = sample_list_meanlogmean(variate_form(sl), sl)
mirroredlogmean(sl::SampleList) = sample_list_mirroredlogmean(variate_form(sl), sl)

##

# Differential entropy for SampleList
# Entropy is pre-computed during computation of the marginal in `approximate_prod_with_sample_list` function
Distributions.entropy(sl::SampleList) = get_entropy(get_meta(sl))

##

vague(::Type{ SampleList }; nsamples::Int = DEFAULT_SAMPLE_LIST_N_SAMPLES)                        = sample_list_vague(Univariate, nsamples)
vague(::Type{ SampleList }, dims::Int; nsamples::Int = DEFAULT_SAMPLE_LIST_N_SAMPLES)             = sample_list_vague(Multivariate, dims, nsamples)
vague(::Type{ SampleList }, dims::Tuple{Int, Int}; nsamples::Int = DEFAULT_SAMPLE_LIST_N_SAMPLES) = sample_list_vague(Matrixvariate, dims, nsamples)
vague(::Type{ SampleList }, dim1::Int, dim2::Int; nsamples::Int = DEFAULT_SAMPLE_LIST_N_SAMPLES)  = sample_list_vague(Matrixvariate, (dim1, dim2), nsamples)

##

sample_list_default_prod_strategy() = BootstrapImportanceSampling()

## prod related stuff
function approximate_prod_with_sample_list(x, y, nsamples::Int = DEFAULT_SAMPLE_LIST_N_SAMPLES)
    return approximate_prod_with_sample_list(Random.GLOBAL_RNG, x, y, nsamples)
end

function approximate_prod_with_sample_list(rng::AbstractRNG, x, y, nsamples::Int = DEFAULT_SAMPLE_LIST_N_SAMPLES)
    return approximate_prod_with_sample_list(rng, sample_list_default_prod_strategy(), x, y, nsamples)
end

# `x` is proposal distribution
# `y` is integrand distribution
function approximate_prod_with_sample_list(rng::AbstractRNG, ::BootstrapImportanceSampling, x::Any, y::Any, nsamples::Int = DEFAULT_SAMPLE_LIST_N_SAMPLES)
    @assert nsamples >= 1 "Number of samples should be non-positive"

    xlogpdf, xsample = logpdf_sample_friendly(x)
    ylogpdf, ysample = logpdf_sample_friendly(y)

    T            = promote_type(eltype(x), eltype(y))
    U            = variate_form(x)
    xsize        = size(x)
    preallocated = preallocate_samples(T, xsize, nsamples)
    samples      = rand!(rng, xsample, reshape(preallocated, (xsize..., nsamples)))

    raw_weights  = Vector{T}(undef, nsamples) # un-normalised
    norm_weights = Vector{T}(undef, nsamples) # normalised

    H_x         = zero(T)
    weights_sum = zero(T)

    for i in 1:nsamples
        # Static indexing from reshaped array
        sample_i = static_getindex(U, xsize, samples, i)
        # Apply log-pdf functions to the samples
        log_sample_x = logpdf(xlogpdf, sample_i)
        log_sample_y = logpdf(ylogpdf, sample_i)

        raw_weight = exp(log_sample_y)

        raw_weights[i]  = raw_weight
        norm_weights[i] = raw_weight # will be renormalised later

        weights_sum += raw_weight
        H_x         += raw_weight * (log_sample_x + log_sample_y)
    end

    # Normalise weights
    @turbo for i in 1:nsamples
        norm_weights[i] /= weights_sum
    end

    # Renormalise H_x
    H_x /= weights_sum

    # Compute the separate contributions to the entropy
    H_y = log(weights_sum) - log(nsamples)
    H_x = -H_x

    entropy = H_x + H_y

    # Inform next step about the proposal and integrand to be used in entropy calculation in smoothing
    logproposal  = xlogpdf
    logintegrand = ylogpdf

    meta = SampleListMeta(raw_weights, entropy, logproposal, logintegrand)

    return SampleList(Val(xsize), preallocated, norm_weights, meta)
end

function approximate_prod_with_sample_list(rng::AbstractRNG, ::AbstractSampleListSamplingMethod, x::SampleList, y::SampleList, nsamples::Int = DEFAULT_SAMPLE_LIST_N_SAMPLES)
    error("Unsupported SampleList × SampleList prod operation.")
end

function approximate_prod_with_sample_list(rng::AbstractRNG, method::AbstractSampleListSamplingMethod, x::SampleList, y::Any, nsamples::Int = DEFAULT_SAMPLE_LIST_N_SAMPLES)
    return approximate_prod_with_sample_list(rng, method, y, x, nsamples)
end

# prod of a pdf (or distribution) message and a SampleList message
# this function is capable to calculate entropy with SampleList messages in VMP setting
function approximate_prod_with_sample_list(rng::AbstractRNG, ::BootstrapImportanceSampling, x::Any, y::SampleList{ D }, nsamples::Int = DEFAULT_SAMPLE_LIST_N_SAMPLES) where { D }

    # TODO: In principle it is possible to implement different prod approximation for different nsamples
    # TODO: This feature would be probably super rare in use so lets postpone it and mark as todo
    @assert length(y) === nsamples "Unsupported SampleList prod approximation: nsamples should match"

    # TODO: Is it possible to support different variate forms here?
    @assert variate_form(x) === variate_form(y) "Unsupported SampleList prod approximation: variate forms should match"

    # Suppose in the previous time step m1(pdf) and m2(pdf) messages collided.
    # The resulting collision m3 (sampleList) = m1*m2 is supposed to carry
    # the proposal (m1) and integrand (m2) distributions. m1 is the message from which
    # the samples are drawn. m2 is the message on which the samples are evaluated and
    # weights are calculated. In case Particle Filtering (BP), entropy will not be calculated
    # and in the first step there won't be any integrand information.

    xlogpdf, xsample = logpdf_sample_friendly(x)

    log_integrand = if is_meta_present(y) && get_logintegrand(y) !== nothing
        # recall that we are calculating m3*m4. If m3 consists of integrand information
        # update it: new_integrand = m2*m3. This allows us to collide arbitrary number of beliefs
        # to approximate posterior and yet estimate the entropy.
        let y_logintegrand = get_logintegrand(y)
            (sample) -> call_logintegrand(y_logintegrand, sample) + logpdf(xlogpdf, sample)
        end
    else
        # If there is no integrand information before, set it to m4
        xlogpdf
    end

    samples = get_samples(y) # samples come from proposal (m1)
    weights = get_weights(y)

    # Resulting samples and weights will go here
    rcontainer = similar(y)
    # vec here just to convert unusual containers into an array just in case
    # does nothing if `get_weights` returns an array
    rsamples, rweights = get_samples(rcontainer), vec(get_weights(rcontainer))

    rweights_raw      = similar(rweights)
    rweights_prod_sum = zero(eltype(get_weights(y)))

    H_x = zero(eltype(rweights_raw))

    # Compute sample weights
    @inbounds for i in 1:nsamples
        log_sample_x    = logpdf(xlogpdf, samples[i]) # evaluate samples in logm4, i.e. logm4(s)
        raw_weight      = exp(log_sample_x)       # m4(s)
        raw_weight_prod = raw_weight * weights[i] # update the weights of posterior w_unnormalized = m4(s)*w_prev

        rweights_raw[i] = raw_weight
        rweights[i]     = raw_weight_prod

        rweights_prod_sum += raw_weight_prod

        H_x += raw_weight_prod * log_sample_x
    end

    H_x /= rweights_prod_sum

    # Normalize prod weights
    @turbo for i in 1:nsamples
        rweights[i] /= rweights_prod_sum
    end

    # Effective number of particles (Base.Generator is allocation free version of map)
    # neff = sum(rweights .^ 2)
    neff = 1 / mapreduce(abs2, +, rweights)

    # Resample and readjust entropy approximation if required
    if neff < nsamples / 10
        sample!(rng, samples, Weights(weights), rsamples)
        fillv = one(eltype(rweights)) / nsamples
        fill!(rweights, fillv)
        # Readjust H_x after resampling
        H_x = zero(eltype(rweights_raw))
        for i in 1:nsamples
            log_sample_x = logpdf(xlogpdf, rsamples[i])
            H_x += (exp(log_sample_x) * fillv) * log_sample_x
        end
        # H_x /= 1 rweights are normalized at this point, no need to track of its sum
    else
        # Just copy the existing samples instead
        copyto!(get_linear_samples(rsamples), get_linear_samples(samples))
    end

    meta = if is_meta_present(y) && get_logproposal(y) !== nothing && get_unnormalised_weights(y) !== nothing

        y_unnormalised_weights     = get_unnormalised_weights(y)
        r_unnormalised_weights     = similar(y_unnormalised_weights)
        r_unnormalised_weights_sum = zero(eltype(r_unnormalised_weights))
        @inbounds for i in 1:nsamples
            r_unnormalised_weights_prod = y_unnormalised_weights[i] * rweights_raw[i]
            r_unnormalised_weights[i]   = r_unnormalised_weights_prod
            r_unnormalised_weights_sum += r_unnormalised_weights_prod

            H_x += rweights[i] * (call_logproposal(y, rsamples[i]) + log(y_unnormalised_weights[i]))
        end

        H_y = log(r_unnormalised_weights_sum) - log(nsamples)
        H_x = -H_x

        entropy = H_x + H_y

        SampleListMeta(r_unnormalised_weights, entropy, get_logproposal(y), log_integrand)
    else
        SampleListMeta(nothing, nothing, nothing, log_integrand)
    end

    return SampleList(Val(D), get_linear_samples(rsamples), rweights, meta)
end

############################################################################################
## Everything below this comment is a low-level implementation of SampleList
## It consists of various routines to compute statistics from internal list implementation
## Here we also implement iteration utilities
## ##########################################################################################

## Lowlevel implementation below...

@inline static_getindex(::Type{ Univariate }, ndims::Tuple{}, samples, i)            = samples[i]
@inline static_getindex(::Type{ Multivariate }, ndims::Tuple{Int}, samples, i)       = view(samples, :, i)
@inline static_getindex(::Type{ Matrixvariate }, ndims::Tuple{Int, Int}, samples, i) = view(samples, :, :, i)

## Preallocation utilities

preallocate_samples(::Type{T}, dims::Tuple, length::Int) where T = Vector{T}(undef, length * prod(dims))

## Linearization functions

# Here we cast an array of arrays into a single flat array of floats for better performance
# There is a package for this called ArraysOfArrays.jl, but the performance of handwritten version is way better
# We provide custom optimized mean/cov function for our implementation with LoopVectorization.jl package
function sample_list_linearize end

function sample_list_linearize(samples::AbstractVector{T}, nsamples, size) where { T <: Number }
    return samples
end

function sample_list_linearize(samples::AbstractVector, nsamples, size)
    T = deep_eltype(samples)
    alloc = Vector{T}(undef, nsamples * size)
    for i in 1:nsamples
        copyto!(alloc, (i - 1) * size + 1, samples[i], 1, size)
    end
    return alloc
end

## Cache utilities

sample_list_mean(sl::SampleList, cached)     = sample_list_mean(variate_form(sl), sl, cached)
sample_list_mean_cov(sl::SampleList, cached) = sample_list_mean_cov(variate_form(sl), sl, cached)

# Cache ignoring versions
function sample_list_mean(::Type{ U }, sl::SampleList, ::Val{ false }) where U
    return sample_list_mean!(fill!(similar(get_mean_storage(get_cache(sl))), zero(deep_eltype(sl))), U, sl)
end

function sample_list_mean_cov(::Type{ U }, sl::SampleList, ::Val{ false }) where U
    mean = sample_list_mean(U, sl, Val(false))
    cov = fill!(similar(get_cov_storage(get_cache(sl))), zero(deep_eltype(sl)))
    sample_list_covm!(cov, mean, U, sl)
    return mean, cov
end

# By default we try to save mean in an internal cache
function sample_list_mean(::Type{ U }, sl::SampleList, ::Val{ true }) where U
    cache = get_cache(sl)
    mean  = get_mean_storage(cache)
    # If no cache present, compute and save
    if !is_mean_cached(cache)
        mean = cache_mean!(cache, sample_list_mean!(mean, U, sl))
    end
    return mean
end

# By default we try to save cov in an internal cache
function sample_list_mean_cov(::Type{ U }, sl::SampleList, ::Val{ true }) where U
    cache = get_cache(sl)
    mean  = sample_list_mean(U, sl, Val(true))
    cov   = get_cov_storage(cache)
    # If no cache present, compute and save
    if !is_cov_cached(cache)
        cov = cache_cov!(cache, sample_list_covm!(cov, mean, U, sl))
    end
    return (mean, cov)
end

## Specific implementations

# Compute mean in a preallocated container and return it
function sample_list_mean! end
# Compute covariance with known mean in a preallocated container and return it
function sample_list_covm! end
# Compute mean and variance
function sample_list_mean_var end
# Compute E[log(x)]
function sample_list_logmean end
# Compute E[xlog(x)]
function sample_list_meanlogmean end
# Compute E[log(1 - x)]
function sample_list_mirroredlogmean end
# Return vague weak-informative sample list
function sample_list_vague end

## Univariate

function sample_list_mean!(μ, ::Type{ Univariate }, sl::SampleList)
    n, samples, weights = get_data(sl)
    @turbo for i in 1:n
        μ += weights[i] * samples[i]
    end
    return μ
end

function sample_list_covm!(σ², μ, ::Type{ Univariate }, sl::SampleList)
    n, samples, weights = get_data(sl)
    @turbo for i in 1:n
        σ² += weights[i] * abs2(samples[i] - μ)
    end
    σ² = (n / (n - 1)) * σ²
    return σ²
end

function sample_list_mean_var(::Type{ Univariate }, sl::SampleList)
    return sample_list_mean_cov(Univariate, sl, Val(true))
end

function sample_list_logmean(::Type{ Univariate }, sl::SampleList)
    n, samples, weights = get_data(sl)
    logμ = sample_list_zero_element(sl)
    @turbo for i in 1:n
        logμ += weights[i] * log(samples[i])
    end
    return logμ
end

function sample_list_meanlogmean(::Type{ Univariate }, sl::SampleList)
    n, samples, weights = get_data(sl)
    μlogμ = sample_list_zero_element(sl)
    @turbo for i in 1:n
        μlogμ += weights[i] * samples[i] * log(samples[i])
    end
    return μlogμ
end

function sample_list_mirroredlogmean(::Type{ Univariate }, sl::SampleList)
    @assert all(0 .<= sl .< 1) "mirroredlogmean does not apply to variables outside of the range [0, 1]"
    n, samples, weights = get_data(sl)
    mirμ = sample_list_zero_element(sl)
    @turbo for i in 1:n
        mirμ += weights[i] * log(1 - samples[i])
    end
    return mirμ
end

function sample_list_vague(::Type{ Univariate }, nsamples::Int)
    targetdist   = vague(Uniform)
    preallocated = preallocate_samples(Float64, (), nsamples)
    rand!(targetdist, preallocated)
    return SampleList(Val(()), preallocated, fill(one(Float64) / nsamples, nsamples), nothing)
end

## Multivariate

function sample_list_mean!(μ, ::Type{ Multivariate }, sl::SampleList)
    n, samples, weights = get_data(sl)
    k = length(μ)
    @turbo for i in 1:n, j in 1:k
        μ[j] += (weights[i] * samples[(i - 1) * k + j])
    end
    return μ
end

function sample_list_covm!(Σ, μ, ::Type{ Multivariate }, sl::SampleList)
    n, samples, weights = get_data(sl)
    tmp = similar(μ)
    k   = length(tmp)

    @inbounds for i in 1:n
        for j in 1:k
            tmp[j] = samples[(i - 1) * k + j] - μ[j]
        end
        # Fast equivalent of Σ += w .* (tmp * tmp')
        for h in 1:k, l in 1:k
            Σ[(h - 1) * k + l] += weights[i] * tmp[h] * tmp[l]
        end
    end
    s = n / (n - 1)
    @turbo for i in 1:length(Σ)
        Σ[i] *= s
    end
    return Σ
end

function sample_list_mean_var(::Type{ Multivariate }, sl::SampleList)
    μ, Σ = sample_list_mean_cov(Multivariate, sl, Val(true))
    return μ, diag(Σ)
end

function sample_list_logmean(::Type{ Multivariate }, sl::SampleList)
    n, samples, weights = get_data(sl)
    logμ = sample_list_zero_element(sl)
    k = length(logμ)
    @turbo for i in 1:n, j in 1:k
        logμ[j] += (weights[i] * log(samples[(i - 1) * k + j]))
    end
    return logμ
end

function sample_list_meanlogmean(::Type{ Multivariate }, sl::SampleList)
    n, samples, weights = get_data(sl)
    μlogμ = sample_list_zero_element(sl)
    k = length(μlogμ)
    @turbo for i in 1:n, j in 1:k
        cs = samples[(i - 1) * k + j]
        μlogμ[j] += (weights[i] * cs * log(cs))
    end
    return μlogμ
end

function sample_list_vague(::Type{ Multivariate }, dims::Int, nsamples::Int)
    targetdist   = vague(Uniform)
    preallocated = preallocate_samples(Float64, (dims, ), nsamples)
    rand!(targetdist, preallocated)
    return SampleList(Val((dims, )), preallocated, fill(one(Float64) / nsamples, nsamples), nothing)
end

## Matrixvariate

function sample_list_mean!(μ, ::Type{ Matrixvariate }, sl::SampleList)
    n, samples, weights = get_data(sl)
    k = length(μ)
    @turbo for i in 1:n, j in 1:k
        μ[j] += (weights[i] * samples[(i - 1) * k + j])
    end
    return μ
end

function sample_list_covm!(Σ, μ, ::Type{ Matrixvariate }, sl::SampleList)
    n, samples, weights = get_data(sl)
    k   = length(μ)
    rμ  = reshape(μ, k)
    tmp = similar(rμ)
    @inbounds for i in 1:n
        for j in 1:k
            tmp[j] = samples[(i - 1) * k + j] - μ[j]
        end
        # Fast equivalent of Σ += w .* (tmp * tmp')
        for h in 1:k, l in 1:k
            Σ[(h - 1) * k + l] += weights[i] * tmp[h] * tmp[l]
        end
    end
    s = n / (n - 1)
    @turbo for i in 1:length(Σ)
        Σ[i] *= s
    end
    return Σ
end

function sample_list_mean_var(::Type{ Matrixvariate }, sl::SampleList)
    μ, Σ = sample_list_mean_cov(Matrixvariate, sl, Val(true))
    return μ, reshape(diag(Σ), size(μ))
end

function sample_list_logmean(::Type{ Matrixvariate }, sl::SampleList)
    n, samples, weights = get_data(sl)
    logμ = sample_list_zero_element(sl)
    k = length(logμ)
    @turbo for i in 1:n, j in 1:k
        logμ[j] += (weights[i] * log(samples[(i - 1) * k + j]))
    end
    return logμ
end

function sample_list_meanlogmean(::Type{ Matrixvariate }, sl::SampleList)
    n, samples, weights = get_data(sl)
    μlogμ = sample_list_zero_element(sl)
    k = length(μlogμ)
    @turbo for i in 1:n, j in 1:k
        cs = samples[(i - 1) * k + j]
        μlogμ[j] += (weights[i] * cs * log(cs))
    end
    return μlogμ
end

function sample_list_vague(::Type{ Matrixvariate }, dims::Tuple{Int, Int}, nsamples::Int)
    targetdist   = vague(Uniform)
    preallocated = preallocate_samples(Float64, dims, nsamples)
    rand!(targetdist, preallocated)
    return SampleList(Val(dims), preallocated, fill(one(Float64) / nsamples, nsamples), nothing)
end

## Array operations, broadcasting and mapping

struct SamplesOnlyIterator{T, L} <: AbstractVector{T}
    samplelist :: L

    function SamplesOnlyIterator(samplelist::L) where { L <: SampleList }
        return new{ samples_type(eltype(samplelist)), L }(samplelist)
    end
end

samples_type(::Type{ T }) where { L, R, T <: Tuple{L, R} } = L

@inline Base.size(iter::SamplesOnlyIterator) = (length(iter.samplelist), )
@inline Base.getindex(iter::SamplesOnlyIterator, i::Int)    = first(getindex(iter.samplelist, i))

@inline function Base.setindex!(iter::SamplesOnlyIterator, v, i::Int)
    samples    = get_linear_samples(iter.samplelist)
    sample_len = prod(ndims(iter.samplelist))
    left  = (i - 1) * sample_len + 1
    right = left + sample_len - 1
    copyto!(view(samples, left:right), v)
    v
end

get_linear_samples(iter::SamplesOnlyIterator) = get_linear_samples(iter.samplelist)

##

Base.iterate(sl::SampleList)             = (sl[1], 2)
Base.iterate(sl::SampleList, state::Int) = state <= length(sl) ? (sl[state], state + 1) : nothing

@inline Base.getindex(sl::SampleList, i::Int) = sample_list_get_index(variate_form(sl), ndims(sl), sl, i)

@inline function sample_list_get_index(::Type{ Univariate }, ndims, sl, i)
    return (get_linear_samples(sl)[i], get_linear_weights(sl)[i])
end

@inline function sample_list_get_index(::Type{ Multivariate }, ndims, sl, i)
    samples = get_linear_samples(sl)
    left  = (i - 1) * ndims + 1
    right = left + ndims - 1
    # ndims is compile-time here
    return (SVector{ndims}(view(samples, left:right)), get_linear_weights(sl)[i])
end

@inline function sample_list_get_index(::Type{ Matrixvariate }, ndims, sl, i)
    p = prod(ndims)
    samples = get_linear_samples(sl)
    left  = (i - 1) * p + 1
    right = left + p - 1
    # ndims are compile-time here
    return (SMatrix{ ndims[1], ndims[2] }(reshape(view(samples, left:right), ndims)), get_linear_weights(sl)[i])
end

## Transformation routines

transform_samples(f::Function, sl::SampleList) = sample_list_transform_samples(variate_form(sl), f, sl)

@inline input_for_transform(::Type{ Univariate }, samples, size, left, right)    = samples[left]
@inline input_for_transform(::Type{ Multivariate }, samples, size, left, right)  = SVector{size}(view(samples, left:right))
@inline input_for_transform(::Type{ Matrixvariate }, samples, size, left, right) = SMatrix{size[1],size[2]}(reshape(view(samples, left:right), size))

function sample_list_transform_samples(::Type{ U }, f::Function, sl::SampleList) where U
    n, samples, weights = get_data(sl)
    input_size = ndims(sl)
    input_len  = prod(input_size)

    # Here we simulate an original implementation of map function from Julia Base
    # Trick here is to compute the first value so compiler may infer the actual output type
    # Later on output_size and output_len are compile-time constants (given that `f` is type-stable)
    first_item  = f(input_for_transform(U, samples, input_size, 1, input_len))

    # After computing first value compiler knows the output type and size of this type
    output_size = size(first_item)
    output_len  = prod(output_size)

    preallocated = preallocate_samples(promote_type(eltype(first_item), eltype(samples)), output_size, n)
    copyto!(view(preallocated, 1:output_len), first_item)

    # We then compute all values from 2 to n into a preallocated buffer
    @views for i in 2:n
        input_left  = (i - 1) * input_len + 1
        input_right = input_left + input_len - 1
        output_left = (i - 1) * output_len + 1
        output_right = output_left + output_len - 1
        # We use static matrix size to ensure that we do not allocate extra memory on a heap
        # Instead we try to do all computations on stack as much as possible
        # If `f` function is "bad" and still allocates matrices this optimisation doesn't really work well
        copyto!(preallocated[output_left:output_right], f(input_for_transform(U, samples, input_size, input_left, input_right)))
    end

    # if `f` is type stable and returns StaticArray `output_size` is a compile-time constant
    return SampleList(Val(output_size), preallocated, weights)
end

function transform_weights!(f::Function, sl::SampleList)
    n, _, weights = get_data(sl)
    map!(f, weights, weights)
    norm = sum(weights)
    @turbo for i in 1:n
        weights[i] /= norm
    end
    return sl
end


text:
module Lingua

export stem

function stem(word, lang)
    word = lowercase(word)
    chars = length(word)

    if chars < 4
        return word
    end

    if chars > 5
        local stem = replace(word, r"ища$", "")
        if stem != word
            return stem
        end
    end

    word = remove_article(word)
    word = remove_plural(word)
    chars = length(word)

    if chars > 3
        word = replace(word, r"я $"x, "")  # masculine

        # normalization (e.g., -а could be a definite article or plural form)
        word = replace(word, r"[аео] $"x, "")

        chars = length(word)
    end

    if chars > 4
        word = replace(word, r"е (?= н $)"x, "")  # -ен → -н

        chars = length(word)
    end

    if chars > 5
        word = replace(word, r"ъ (?= \p{Cyrillic} $)"x, "")  # -ъ� → -�
    end

    return word
end

function remove_article(word)
    chars = length(word)

    if chars > 6
        # definite article with adjectives and masculine
        stem = replace(word, r"ият $"x, "")
        if stem != word
            return stem
        end
    end

    if chars > 5
        stem = replace(word, r"(?:
              ия  # definite articles for nouns:
            | ът  # ∙ masculine
            | та  # ∙ feminine
            | то  # ∙ neutral
            | те  # ∙ plural
        ) $"x, "")
        if stem != word
            return stem
        end
    end

    if chars > 4
        # article for masculine
        stem = replace(word, r"ят $"x, "")
        if stem != word
            return stem
        end
    end

    return word
end

function remove_plural(word)
    chars = length(word)

    if chars > 6
        # specific plural rules for some words (masculine)
        stem = replace(word, r"ове $"x, "")
        if stem != word
            return stem
        end
        stem = replace(word, r"еве $"x, "й")
        if stem != word
            return stem
        end
        stem = replace(word, r"овци $"x, "о")
        if stem != word
            return stem
        end
    end

    if chars > 5
        stem = replace(word, r"зи $"x, "г")
        if stem != word
            return stem
        end
        # -е�и → -я�
        m = match(r"е ( \p{Cyrillic} ) и $"x, word)
        if m != nothing
            return replace(word, Regex("$(m.match)\$"), "я$(m.captures[1])")
        end
        stem = replace(word, r"ци $"x, "к")
        if stem != word
            return stem
        end
        stem = replace(word, r"(?: та | ища ) $"x, "")
        if stem != word
            return stem
        end
    end

    if chars > 4
        stem = replace(word, r"си $"x, "х")
        if stem != word
            return stem
        end
        # plural for various nouns and adjectives
        stem = replace(word, r"и $"x, "")
        if stem != word
            return stem
        end
    end

    return word
end

end

max_stars_repo_name: patch-orphan/lingua-jl

max_issues_repo_name: patch-orphan/lingua-jl

max_forks_repo_name: patch-orphan/lingua-jl


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

function _gcirc{T<:AbstractFloat}(units::Integer, ra1::T, dec1::T, ra2::T, dec2::T)
    # Convert all quantities to radians.
    if units == 0
        # All radians
        λ_1 = ra1
        λ_2 = ra2
        φ_1 = dec1
        φ_2 = dec2
    elseif units == 1
        # Right ascensions are in hours, declinations in degrees.
        λ_1 = ra1*pi/12.0
        λ_2 = ra2*pi/12.0
        φ_1 = deg2rad(dec1)
        φ_2 = deg2rad(dec2)
    elseif units == 2
        # Right ascensions and declinations are in degrees.
        λ_1 = deg2rad(ra1)
        λ_2 = deg2rad(ra2)
        φ_1 = deg2rad(dec1)
        φ_2 = deg2rad(dec2)
    else
        # In any other case throw an error.
        error("units must be 0 (radians), 1 (hours, degrees) or 2 (degrees)")
    end
    Δφ_2 = (φ_2 - φ_1) * 0.5
    Δλ_2 = (λ_2 - λ_1) * 0.5
    Δσ = 2asin(sqrt(abs2(sin(Δφ_2)) + cos(φ_1) * cos(φ_2) * abs2(sin(Δλ_2))))
    if units == 0
        return Δσ
    else
        return rad2sec(Δσ)
    end
end

"""
    gcirc(units, ra1, dec1, ra2, dec2) -> angular_distance

### Purpose ###

Computes rigorous great circle arc distances.

### Explanation ###

Input position can be either radians, sexagesimal right ascension and
declination, or degrees.

### Arguments ###

* `units`: integer, can be either 0, or 1, or 2.  Describes units of inputs and
 output:
    * 0: everything (input right ascensions and declinations, and output
      distance) is radians
    * 1: right ascensions are in decimal hours, declinations in decimal degrees,
      output distance in arc seconds
    * 2: right ascensions and declinations are in degrees, output distance in arc
      seconds
* `ra1`:  right ascension or longitude of point 1
* `dec1`: declination or latitude of point 1
* `ra2`: right ascension or longitude of point 2
* `dec2`: declination or latitude of point 2

Both `ra1` and `dec1`, and `ra2` and `dec2` can be given as 2-tuples `(ra1,
dec1)` and `(ra2, dec2)`.

### Output ###

Angular distance on the sky between points 1 and 2, as a `AbstractFloat`.  See
`units` argument above for the units.

### Method ###

"Haversine formula" see http://en.wikipedia.org/wiki/Great-circle_distance.

### Example ###

``` julia
gcirc(0, 120, -43, 175, +22)
# => 1.590442261600714
```

### Notes ###

* If `ra1`, `dec1` are scalars, and `ra2`, `dec2` are vectors, then the output
 is a vector giving the distance of each element of `ra2`, `dec2` to `ra1`,
 `dec1`.  Similarly, if `ra1`,`de1` are vectors, and `ra2`,` dec2` are scalars,
 then the output is a vector giving the distance of each element of `ra1`,
 `dec1` to `ra2`, `dec2`.  If both `ra1`, `dec1` and `ra2`, `dec2` are vectors
 then the output is a vector giving the distance of each element of `ra1`,
 `dec1` to the corresponding element of `ra2`, `dec2`.
* The function `sphdist` provides an alternate method of computing a spherical
 distance.
* The Haversine formula can give rounding errors for antipodal points.

Code of this function is based on IDL Astronomy User's Library.
"""
gcirc(units::Integer, ra1::Real, dec1::Real, ra2::Real, dec2::Real) =
    _gcirc(units, promote(float(ra1), float(dec1), float(ra2), float(dec2))...)

function gcirc{R1<:Real, D1<:Real}(units::Integer,
                                   ra1::AbstractArray{R1},
                                   dec1::AbstractArray{D1},
                                   ra2::Real,
                                   dec2::Real)
    @assert length(ra1) == length(dec1)
    dist = similar(ra1, typeof(float(one(R1))))
    for i in eachindex(ra1)
        dist[i] = gcirc(units, ra1[i], dec1[i], ra2, dec2)
    end
    return dist
end

function gcirc{R2<:Real, D2<:Real}(units::Integer,
                                   ra1::Real,
                                   dec1::Real,
                                   ra2::AbstractArray{R2},
                                   dec2::AbstractArray{D2})
    @assert length(ra2) == length(dec2)
    dist = similar(ra2, typeof(float(one(R2))))
    for i in eachindex(ra2)
        dist[i] = gcirc(units, ra1, dec1, ra2[i], dec2[i])
    end
    return dist
end

function gcirc{R1<:Real, D1<:Real, R2<:Real, D2<:Real}(units::Integer,
                                                       ra1::AbstractArray{R1},
                                                       dec1::AbstractArray{D1},
                                                       ra2::AbstractArray{R2},
                                                       dec2::AbstractArray{D2})
    @assert length(ra1) == length(dec1) == length(ra2) == length(dec2)
    dist = similar(ra1, typeof(float(one(R1))))
    for i in eachindex(ra1)
        dist[i] = gcirc(units, ra1[i], dec1[i], ra2[i], dec2[i])
    end
    return dist
end

### Tuples input
gcirc(units::Integer, radec1::Tuple{Real, Real}, ra2::Real, dec2::Real) =
    gcirc(units, radec1..., ra2, dec2)

gcirc(units::Integer, ra1::Real, dec1::Real, radec2::Tuple{Real, Real}) =
    gcirc(units, ra1, dec1, radec2...)

gcirc(units::Integer, radec1::Tuple{Real, Real}, radec2::Tuple{Real, Real}) =
    gcirc(units, radec1..., radec2...)


text:
# 1x1 convolution operator using Householder matrices.
# Adapted from Putzky and Welling (2019): https://arxiv.org/abs/1911.10914
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020
#

export Conv1x1

"""
    C = Conv1x1(k; logdet=false)

 or

    C = Conv1x1(v1, v2, v3; logdet=false)

 Create network layer for 1x1 convolutions using Householder reflections.

 *Input*:

 - `k`: number of channels

 - `v1`, `v2`, `v3`: Vectors from which to construct matrix.

 - `logdet`: if true, returns logdet in forward pass (which is always zero)

 *Output*:

 - `C`: Network layer for 1x1 convolutions with Householder reflections.

 *Usage:*

 - Forward mode: `Y, logdet = C.forward(X)`

 - Backward mode: `ΔX, X = C.backward((ΔY, Y))`

 *Trainable parameters:*

 - Householder vectors `C.v1`, `C.v2`, `C.v3`

 See also: [`get_params`](@ref), [`clear_grad!`](@ref)
"""
struct Conv1x1 <: NeuralNetLayer
    k::Integer
    v1::Parameter
    v2::Parameter
    v3::Parameter
    logdet::Bool
end

@Flux.functor Conv1x1

# Constructor with random initializations
function Conv1x1(k; logdet=false)
    v1 = Parameter(glorot_uniform(k))
    v2 = Parameter(glorot_uniform(k))
    v3 = Parameter(glorot_uniform(k))
    return Conv1x1(k, v1, v2, v3, logdet)
end

function Conv1x1(v1, v2, v3; logdet=false)
    k = length(v1)
    v1 = Parameter(v1)
    v2 = Parameter(v2)
    v3 = Parameter(v3)
    return Conv1x1(k, v1, v2, v3, logdet)
end

function partial_derivative_outer(v::AbstractArray{T, 1}) where T
    k = length(v)
    out1 = v * v'
    n = v' * v
    outer = cuzeros(v, k, k, k)
    for i=1:k
        copyto!(view(outer, i, :, :), out1)
    end
    broadcast!(*, outer, v, outer)
    broadcast!(*, outer, -2/n, outer)
    for j=1:k
        v1 = view(outer,j, :, j)
        broadcast!(+, v1, v1, v)
        v1 = view(outer,j, j, :)
        broadcast!(+, v1, v1, v)
    end
    broadcast!(*, outer, 1/n, outer)
    return outer
end

function partial_derivative_outer(v::CuArray{T, 1}) where T
    k = length(v)
    out1 = v * v'
    n = v' * v
    outer = cuzeros(v, k, k, k)
    for i=1:k
        copyto!(view(outer, i, :, :), out1)
    end
    broadcast!(*, outer, v, outer)
    broadcast!(*, outer, -2/n, outer)
    for j=1:k
        v1 = view(outer,j, :, j)
        broadcast!(+, v1, v1, v)
        v1 = view(outer,j, j, :)
        broadcast!(+, v1, v1, v)
    end
    broadcast!(*, outer, 1/n, outer)
    return outer
end

function mat_tens_i(out::AbstractArray{T, 3}, Mat::AbstractArray{T, 2},
                    Tens::AbstractArray{T, 3}, Mat2::AbstractArray{T, 2}) where T
    tmp = cuzeros(out, size(out, 2), size(out, 3))
    for i=1:size(out, 1)
        mul!(tmp, Mat, Tens[i, :, :])
        broadcast!(*, tmp, tmp, Mat2)
        @views copyto!(out[i, :, :], tmp)
    end
    return out
end

function custom_sum(a::AbstractArray{T, 3}, dims::Tuple{Integer, Integer}) where T
    summed = sum(a, dims=dims)
    return dropdims(summed, dims = (findall(size(summed) .== 1)...,))
end

function conv1x1_grad_v(X::AbstractArray{T, N}, ΔY::AbstractArray{T, N},
                        C::Conv1x1; adjoint=false) where {T, N}

    # Reshape input
    n_in, batchsize = size(X)[N-1:N]
    v1 = C.v1.data
    v2 = C.v2.data
    v3 = C.v3.data
    k = length(v1)

    dv1 = cuzeros(X, k)
    dv2 = cuzeros(X, k)
    dv3 = cuzeros(X, k)

    V1 = v1*v1'/(v1'*v1)
    V2 = v2*v2'/(v2'*v2)
    V3 = v3*v3'/(v3'*v3)

    dV1 = partial_derivative_outer(v1)
    dV2 = partial_derivative_outer(v2)
    dV3 = partial_derivative_outer(v3)

    ∂V1 = deepcopy(dV1)
    ∂V2 = deepcopy(dV2)
    ∂V3 = deepcopy(dV3)

    M1 = (I - 2 * (V2 + V3) + 4*V2*V3)
    M3 = (I - 2 * (V1 + V2) + 4*V1*V2)
    tmp = cuzeros(X, k, k)
    for i=1:k
        # ∂V1
        mul!(tmp, ∂V1[i, :, :], M1)
        @views adjoint ? adjoint!(∂V1[i, :, :], tmp) : copyto!(∂V1[i, :, :], tmp)
        # ∂V2
        v2 = ∂V2[i, :, :]
        broadcast!(+, tmp, v2, 4 * V1 * v2 * V3 - 2 * (V1 * v2 + v2 * V3))
        @views adjoint ? adjoint!(∂V2[i, :, :], tmp) : copyto!(∂V2[i, :, :], tmp)
        # ∂V3
        mul!(tmp, M3, ∂V3[i, :, :])
        @views adjoint ? adjoint!(∂V3[i, :, :], tmp) : copyto!(∂V3[i, :, :], tmp)
    end

    prod_res = cuzeros(X, size(∂V1, 1), prod(size(X)[1:N-2]), n_in)
    inds = [i<N ? (:) : 1 for i=1:N]
    for i=1:batchsize
        inds[end] = i
        Xi = -2f0*reshape(view(X, inds...), :, n_in)
        ΔYi = reshape(view(ΔY, inds...), :, n_in)
        broadcast!(+, dv1, dv1, custom_sum(mat_tens_i(prod_res, Xi, ∂V1, ΔYi), (3, 2)))
        broadcast!(+, dv2, dv2, custom_sum(mat_tens_i(prod_res, Xi, ∂V2, ΔYi), (3, 2)))
        broadcast!(+, dv3, dv3, custom_sum(mat_tens_i(prod_res, Xi, ∂V3, ΔYi), (3, 2)))
    end
    return dv1, dv2, dv3
end

# Forward pass
function forward(X::AbstractArray{T, N}, C::Conv1x1; logdet=nothing) where {T, N}
    isnothing(logdet) ? logdet = C.logdet : logdet = logdet
    Y = cuzeros(X, size(X)...)
    n_in = size(X, N-1)

    v1 = C.v1.data
    v2 = C.v2.data
    v3 = C.v3.data

    inds = [i<N ? (:) : 1 for i=1:N]
    for i=1:size(X, N)
        inds[end] = i
        Xi = reshape(view(X, inds...), :, n_in)
        Yi = chain_lr(Xi, v1, v2, v3)
        view(Y, inds...) .= reshape(Yi, size(view(Y, inds...))...)
    end
    logdet == true ? (return Y, 0) : (return Y)   # logdet always 0
end

# Forward pass and update weights
function forward(X_tuple::Tuple, C::Conv1x1; set_grad::Bool=true)
    ΔX = X_tuple[1]
    X = X_tuple[2]
    ΔY = forward(ΔX, C; logdet=false)    # forward propagate residual
    Y = forward(X, C; logdet=false)  # recompute forward state
    Δv1, Δv2, Δv3 = conv1x1_grad_v(Y, ΔX, C; adjoint=true)  # gradient w.r.t. weights
    if set_grad
        isnothing(C.v1.grad) ? (C.v1.grad = Δv1) : (C.v1.grad += Δv1)
        isnothing(C.v2.grad) ? (C.v2.grad = Δv2) : (C.v2.grad += Δv2)
        isnothing(C.v3.grad) ? (C.v3.grad = Δv3) : (C.v3.grad += Δv3)
    else
        Δθ = [Parameter(Δv1), Parameter(Δv2), Parameter(Δv3)]
    end
    set_grad ? (return ΔY, Y) : (return ΔY, Δθ, Y)
end

# Inverse pass
function inverse(Y::AbstractArray{T, N}, C::Conv1x1; logdet=nothing) where {T, N}
    isnothing(logdet) ? logdet = C.logdet : logdet = logdet
    X = cuzeros(Y, size(Y)...)
    n_in = size(X, N-1)

    v1 = C.v1.data
    v2 = C.v2.data
    v3 = C.v3.data

    inds = [i<N ? (:) : 1 for i=1:N]
    for i=1:size(Y, N)
        inds[end] = i
        Yi = reshape(view(Y, inds...), :, n_in)
        Xi = chain_lr(Yi, v3, v2, v1)
        view(X, inds...) .= reshape(Xi, size(view(X, inds...))...)
    end
    logdet == true ? (return X, 0) : (return X)   # logdet always 0
end

# Inverse pass and update weights
function inverse(Y_tuple::Tuple, C::Conv1x1; set_grad::Bool=true)
    ΔY = Y_tuple[1]
    Y = Y_tuple[2]
    ΔX = inverse(ΔY, C; logdet=false)    # derivative w.r.t. input
    X = inverse(Y, C; logdet=false)  # recompute forward state
    Δv1, Δv2, Δv3 =  conv1x1_grad_v(X, ΔY, C)  # gradient w.r.t. weights
    if set_grad
        isnothing(C.v1.grad) ? (C.v1.grad = Δv1) : (C.v1.grad += Δv1)
        isnothing(C.v2.grad) ? (C.v2.grad = Δv2) : (C.v2.grad += Δv2)
        isnothing(C.v3.grad) ? (C.v3.grad = Δv3) : (C.v3.grad += Δv3)
    else
        Δθ = [Parameter(Δv1), Parameter(Δv2), Parameter(Δv3)]
    end
    set_grad ? (return ΔX, X) : (return ΔX, Δθ, X)
end

## Jacobian-related functions

function jacobian(ΔX::AbstractArray{T, N}, Δθ::Array{Parameter, 1}, X::AbstractArray{T, N}, C::Conv1x1) where {T, N}
    Y = cuzeros(X, size(X)...)
    ΔY = cuzeros(ΔX, size(ΔX)...)
    n_in = size(X, N-1)

    v1 = C.v1.data
    v2 = C.v2.data
    v3 = C.v3.data
    dv1 = Δθ[1].data
    dv2 = Δθ[2].data
    dv3 = Δθ[3].data

    inds = [i<N ? (:) : 1 for i=1:N]
    for i=1:size(X, N)
        inds[end] = i
        Xi = reshape(view(X, inds...), :, n_in)
        Yi = chain_lr(Xi, v1, v2, v3)
        view(Y, inds...) .= reshape(Yi, size(view(Y, inds...))...)

        ΔXi = reshape(view(ΔX, inds...), :, n_in)
        ΔYi = chain_lr(Xi, v1, v2, v3)
        # this is a lot of outer products of 1D vecotrs, need to be cleaned up that's overkill computationnaly
        n1 = norm(v1); n2 = norm(v2); n3 = norm(v3);
        c1 = I - 2f0*v1*v1'/n1^2f0; c2 = I - 2f0*v2*v2'/n2^2f0; c3 = I - 2f0*v3*v3'/n3^2f0;
        ΔYi = chain_lr(ΔXi, v1, v2, v3)
        ΔYi += -2f0*Xi*((dv1*v1'+v1*dv1'-2f0*dot(v1,dv1)*v1*v1'/n1^2f0)/n1^2f0*c2*c3+
                       c1*(dv2*v2'+v2*dv2'-2f0*dot(v2,dv2)*v2*v2'/n2^2f0)/n2^2f0*c3+
                       c1*c2*(dv3*v3'+v3*dv3'-2f0*dot(v3,dv3)*v3*v3'/n3^2f0)/n3^2f0)
        view(ΔY, inds...) .= reshape(ΔYi, size(view(ΔY, inds...))...)
    end

    return ΔY, Y
end

function adjointJacobian(ΔY::AbstractArray{T, N}, Y::AbstractArray{T, N}, C::Conv1x1) where {T, N}
    return inverse((ΔY, Y), C; set_grad=false)
end

function jacobianInverse(ΔY::AbstractArray{T, N}, Δθ::Array{Parameter, 1}, Y::AbstractArray{T, N}, C::Conv1x1) where {T, N}
    return inverse(C).jacobian(ΔY, Δθ[end:-1:1], Y)
end

function adjointJacobianInverse(ΔX::AbstractArray{T, N}, X::AbstractArray{T, N}, C::Conv1x1) where {T, N}
    ΔX, Δθinv, X = inverse(C).adjointJacobian(ΔX, X)
    return ΔX, Δθinv[end:-1:1], X
end

function inverse(C::Conv1x1)
    return Conv1x1(C.k, C.v3, C.v2, C.v1, C.logdet)
end

## Other utils

# Clear gradients
function clear_grad!(C::Conv1x1)
    C.v1.grad = nothing
    C.v2.grad = nothing
    C.v3.grad = nothing
end

# Get parameters
get_params(C::Conv1x1) = [C.v1, C.v2, C.v3]


text:
using ConfigEnv
using CSV, DataFrames
using AlphaVantage

dotenv()

df1 = time_series_daily_adjusted("SPY", outputsize = "full", datatype = "csv", parser = x -> CSV.File(IOBuffer(x.body))) |> DataFrame
sort!(df1, :timestamp)
rename!(df1, :timestamp => :Date)
rename!(df1, :adjusted_close => "Adj Close")
rename!(df1, :open => :Open)
rename!(df1, :close => :Close)
rename!(df1, :high => :High)
rename!(df1, :low => :Low)
rename!(df1, :volume => :Volume)
CSV.write("SPY.csv", df1)

df2 = time_series_daily_adjusted("AGG", outputsize = "full", datatype = "csv", parser = x -> CSV.File(IOBuffer(x.body))) |> DataFrame
sort!(df2, :timestamp)
rename!(df2, :timestamp => :Date)
rename!(df2, :adjusted_close => "Adj Close")
rename!(df2, :open => :Open)
rename!(df2, :close => :Close)
rename!(df2, :high => :High)
rename!(df2, :low => :Low)
rename!(df2, :volume => :Volume)
CSV.write("AGG.csv", df2)


text:
function [fig] = Plot2(c,Trans,Param)
OrdPoly=Param.OrdPoly;
Subs=Param.OrdPoly;
fig=Param.fig;
wQX=Param.w;
xQX=Param.xw;
wQZ=Param.w;
xQZ=Param.xw;
figure(fig);
hold on
dS=2/Subs;
for iF=1:size(Param.Grid.Faces,2)
  ksi=-1;
  for is=1:Subs
    eta=-1;
    for js=1:Subs
      cP=0;
      for jx=1:OrdPoly+1
        for jz=1:OrdPoly+1
          cP=cP+Lagrange(ksi+0.5*dS,xQX,jx)*Lagrange(eta+0.5*dS,xQZ,jz)...
            *c(1,jx,jz,iF);
        end
      end
      [X1]=Trans(ksi,eta,Param.Grid.Faces(iF),Param.Grid);
      [X2]=Trans(ksi+dS,eta,Param.Grid.Faces(iF),Param.Grid);
      [X3]=Trans(ksi+dS,eta+dS,Param.Grid.Faces(iF),Param.Grid);
      [X4]=Trans(ksi,eta+dS,Param.Grid.Faces(iF),Param.Grid);

      if strcmp(Param.GridType,'Radial')
        r1=sqrt(X1(1,1,1)^2+X1(1,1,2)^2);
        X1=(Param.RadIP+(r1-Param.RadI)/(Param.RadO-Param.RadI))/r1*X1;
        r2=sqrt(X2(1,1,1)^2+X2(1,1,2)^2);
        X2=(Param.RadIP+(r2-Param.RadI)/(Param.RadO-Param.RadI))/r2*X2;
        r3=sqrt(X3(1,1,1)^2+X3(1,1,2)^2);
        X3=(Param.RadIP+(r3-Param.RadI)/(Param.RadO-Param.RadI))/r3*X3;
        r4=sqrt(X4(1,1,1)^2+X4(1,1,2)^2);
        X4=(Param.RadIP+(r4-Param.RadI)/(Param.RadO-Param.RadI))/r4*X4;
      end
      x=[X1(1) X2(1) X3(1) X4(1) X1(1)];
      y=[X1(2) X2(2) X3(2) X4(2) X1(2)];

      fill(x,y,cP);
      eta=eta+dS;
    end
    ksi=ksi+dS;
  end
end
hold off
fig=fig+1;
end


text:
using Distributed, Pkg
Pkg.activate(joinpath(@__DIR__, ".."))
if haskey(ENV, "FOOD_CACHING_PARALLEL")
    nprocs() < Sys.CPU_THREADS && addprocs(20,
                                           exeflags = "--project",
                                           dir = joinpath(@__DIR__, ".."))
end
using FoodCachingModels
import FoodCachingFitting: resume_fit

resume_fit(ARGS[1])


text:
module ProjectPath

using Pkg

function project(m::Module, xs...)
    path = pathof(m)
    path === nothing && return dirname(Pkg.project().path)
    return normpath(joinpath(dirname(dirname(path)), xs...))
end

src(m::Module, xs...) = project(m, "src", xs...)
deps(m::Module, xs...) = project(m, "deps", xs...)
test(m::Module, xs...) = project(m, "test", xs...)

macro include(xs...)
    return esc(:(Base.include($project($__module__, "src", $(xs...)))))
end

macro project(xs...)
    return esc(:($project($__module__, $(xs...))))
end

end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %%
function collatz(n)
    print(n)
    while n ∉ (1, 0, -1, -5, -17)
        n = iseven(n) ? n ÷ 2 : 3n + 1
        print("→", n)
    end
end

# %%
collatz(73)

# %%
collatz(-15)

# %%
collatz(-7)

# %%
collatz(-50)

# %%
collatz(big(2)^73 + 1)

# %%
@code_warntype collatz(73)

# %%
@code_warntype collatz(big(2)^73 + 1)

# %%
function collatz_seq(n, maxiters = 10^10)
    s = [n]
    l = 0
    while n > 1 && l < maxiters
        n = iseven(n) ? n ÷ 2 : 3n + 1
        push!(s, n)
        l += 1
    end
    s
end

collatz_seq(73)

# %%
function collatz_length(n, maxiters = 10^10)
    l = 0
    while n > 1 && l < maxiters
        n = iseven(n) ? n ÷ 2 : 3n + 1
        l += 1
    end
    l
end

collatz_length(73)

# %%
@code_warntype collatz_length(73, 10^10)

# %%
@code_native debuginfo=:none collatz_length(73, 10^10)

# %%
using Plots
n = 1:10^4
scatter(n, collatz_length.(n); label="", ms=1.5, msc=:auto, ma=0.5)

# %%

max_stars_repo_path: 0019/Collatz problem.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0019/Collatz problem.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0019/Collatz problem.jl

max_forks_repo_name: genkuroki/public


text:
using LACosmic
using Documenter

DocMeta.setdocmeta!(LACosmic, :DocTestSetup, :(using LACosmic); recursive=true)

makedocs(;
    modules=[LACosmic],
    authors="Miles Lucas <mdlucas@hawaii.edu> and contributors",
    repo="https://github.com/JuliaAstro/LACosmic.jl/blob/{commit}{path}#{line}",
    sitename="LACosmic.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://juliaastro.github.io/LACosmic.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Examples" => "examples.md",
        "API/Reference" => "api.md"
    ],
)

deploydocs(;
    repo="github.com/JuliaAstro/LACosmic.jl",
    devbranch="main",
    push_preview=true
)


text:
module SpeedTests

using Base.Test
using DataArrays, DataFrames
using DataFramesMeta
using Devectorize

srand(1)
const n = 5_000_000
a = rand(n)
b = rand(n)
da = data(a)
db = data(b)
df = DataFrame(a = da, b = db)
df2 = DataFrame(Any[a, b])
names!(df2, [:a, :b])

Base.values(da::DataArray) = da.data

function dot1(a::Vector, b::Vector)
    x = 0.0
    for i in 1:length(a)
        x += a[i] * b[i]
    end
    return x
end

function dot2(da::DataVector, db::DataVector)
    T = eltype(da)
    x = 0.0
    for i in 1:length(da)
        x += da[i]::T * db[i]::T
    end
    return x
end

function dot3(df::DataFrame)
    da, db = df[:a], df[:b]
    T = eltype(da)
    x = 0.0
    for i in 1:length(da)
        x += da[i]::T * db[i]::T
    end
    return x
end

function dot4(df::DataFrame)
    da, db = df[:a], df[:b]
    return dot2(da, db)
end

function dot5(da::DataVector, db::DataVector)
    x = 0.0
    for i in 1:length(da)
        x += da.data[i] * db.data[i]
    end
    return x
end

function dot6(da::DataVector, db::DataVector)
    x = 0.0
    for i in 1:length(da)
        x += values(da)[i] * values(db)[i]
    end
    return x
end

function dot7(da::DataVector, db::DataVector)
    x = 0.0
    for i in 1:length(da)
        if !(isna(da, i) || isna(da, i))
            x += values(da)[i] * values(db)[i]
        end
    end
    return x
end

function dot8(a::Vector, b::Vector)
    x = 0.0
    for i in 1:length(a)
        if !(isnan(a[i]) || isnan(a[i]))
            x += a[i] * b[i]
        end
    end
    return x
end

function dot9(df::DataFrame)
    @with df begin
        x = 0.0
        for i in 1:length(:a)
            x += values(:a)[i] * values(:b)[i]
        end
        x
    end
end

## function dot10(df::DataFrame)
##     @with df begin
##         @devec x = sum(:a .* :b)
##         x
##     end
## end
## t10 = @elapsed dot10(df)

t1 = @elapsed dot1(a, b)
t2 = @elapsed dot2(da, db)
t3 = @elapsed dot3(df)
t4 = @elapsed dot4(df)
t5 = @elapsed dot5(da, db)
t6 = @elapsed dot6(da, db)
t7 = @elapsed dot7(da, db)
t8 = @elapsed dot8(a, b)
t9 = @elapsed dot9(df)

end


text:
if isdefined(Base, :bodyfunction)
    __lookup_kwbody__(m) = Base.bodyfunction(m)
else
    __lookup_kwbody__(m) = missing
end

function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    # widgets
    ## slider
    precompile(Tuple{typeof(slider),UnitRange{Int}})
    precompile(Tuple{Core.kwftype(typeof(slider)),NamedTuple{(:signal, :orientation),Tuple{Signal{Int},Char}},typeof(slider),UnitRange{Int}})
    precompile(Tuple{typeof(push!),Slider{Int},UnitRange{Int},Int})
    precompile(Tuple{typeof(push!),Slider{Float32},StepRangeLen{Float32,Float64,Float64},Float32})
    precompile(Tuple{typeof(push!),Slider{Float64},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},Float64})
    precompile(Tuple{Type{Slider},Signal{Int},Gtk.GtkScaleLeaf,UInt,Array{Any,1}})
    precompile(Tuple{Type{Slider},Signal{Float64},GtkScaleLeaf,UInt,Array{Any,1}})
    precompile(Tuple{Type{Slider},Signal{Float32},GtkScaleLeaf,UInt,Array{Any,1}})
    let fbody = __lookup_kwbody__(which(slider, (StepRangeLen{Float32,Float64,Float64},)))
        if !ismissing(fbody)
            precompile(fbody, (Nothing,Nothing,Signal{Float32},String,Bool,Nothing,typeof(slider),StepRangeLen{Float32,Float64,Float64},))
            precompile(fbody, (Nothing,Nothing,Signal{Float64},String,Bool,Nothing,typeof(slider),StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},))
        end
    end
    ## checkbox
    precompile(Tuple{typeof(checkbox)})
    precompile(Tuple{typeof(checkbox),Bool})
    precompile(Tuple{Core.kwftype(typeof(checkbox)),NamedTuple{(:label,),Tuple{String}},typeof(checkbox)})
    precompile(Tuple{typeof(push!),Checkbox,Bool})
    ## togglebutton
    precompile(Tuple{typeof(togglebutton)})
    precompile(Tuple{typeof(togglebutton),Bool})
    precompile(Tuple{Core.kwftype(typeof(togglebutton)),NamedTuple{(:label,),Tuple{String}},typeof(togglebutton)})
    ## button
    precompile(Tuple{typeof(button)})
    precompile(Tuple{typeof(button),String})
    precompile(Tuple{Core.kwftype(typeof(button)),NamedTuple{(:widget,),Tuple{Gtk.GtkToolButtonLeaf}},typeof(button)})
    precompile(Tuple{Core.kwftype(typeof(button)),NamedTuple{(:widget,),Tuple{Gtk.GtkButtonLeaf}},typeof(button)})
    ## spinbutton
    precompile(Tuple{typeof(spinbutton),UnitRange{Int}})
    precompile(Tuple{Type{SpinButton},Signal{Int},Gtk.GtkSpinButtonLeaf,UInt,Array{Any,1}})
    precompile(Tuple{Core.kwftype(typeof(spinbutton)),NamedTuple{(:widget, :signal),Tuple{Gtk.GtkSpinButtonLeaf,Signal{Int}}},typeof(spinbutton),UnitRange{Int}})
    precompile(Tuple{typeof(push!),SpinButton{Int},UnitRange{Int},Int})
    ## cyclicspinbutton
    precompile(Tuple{typeof(cyclicspinbutton),UnitRange{Int},Signal{Bool}})
    precompile(Tuple{Type{CyclicSpinButton},Signal{Int},Gtk.GtkSpinButtonLeaf,UInt,Array{Any,1}})
    ## textbox
    precompile(Tuple{typeof(textbox),String})
    precompile(Tuple{Type{Textbox},Signal{String},Gtk.GtkEntryLeaf,UInt,Array{Any,1},Nothing})
    precompile(Tuple{Type{Textbox},Signal{Int},Gtk.GtkEntryLeaf,UInt,Array{Any,1},UnitRange{Int}})
    precompile(Tuple{Type{Textbox},Signal{Float64},GtkEntryLeaf,UInt,Array{Any,1},Nothing})
    precompile(Tuple{Type{Textbox},Signal{Float32},GtkEntryLeaf,UInt,Array{Any,1},Nothing})
    precompile(Tuple{Core.kwftype(typeof(textbox)),NamedTuple{(:signal,),Tuple{Signal{Int}}},typeof(textbox),Type{Int}})
    precompile(Tuple{Core.kwftype(typeof(textbox)),NamedTuple{(:range,),Tuple{UnitRange{Int}}},typeof(textbox),Int})
    precompile(Tuple{Core.kwftype(typeof(textbox)),NamedTuple{(:widget, :signal, :range),Tuple{GtkEntryLeaf,Signal{Int},Base.OneTo{Int}}},typeof(textbox),Int})
    let fbody = __lookup_kwbody__(which(textbox, (Type{Float32},)))
        if !ismissing(fbody)
            precompile(fbody, (GtkEntryLeaf,Nothing,Nothing,Signal{Float32},Bool,Nothing,Symbol,typeof(textbox),Type{Float32},))
        end
    end
    ## textarea
    precompile(Tuple{typeof(textarea),String})
    ## label
    precompile(Tuple{typeof(label),String})
    ## dropdown
    precompile(Tuple{typeof(dropdown),Array{Pair{String,Function},1}})
    precompile(Tuple{Core.kwftype(typeof(dropdown)),NamedTuple{(:label,),Tuple{String}},typeof(dropdown),Array{Pair{String,Function},1}})
    ## progressbar
    precompile(Tuple{typeof(progressbar),Interval{:closed,:closed,Int}})
    ## player
    precompile(Tuple{typeof(player),UnitRange{Int}})
    precompile(Tuple{typeof(player),Signal{Int},UnitRange{Int}})
    let fbody = __lookup_kwbody__(which(player, (Signal{Int},Base.OneTo{Int},)))
        if !ismissing(fbody)
            precompile(fbody, (String,Int,typeof(player),Signal{Int},Base.OneTo{Int},))
        end
    end
    ## timewidget
    precompile(Tuple{Core.kwftype(typeof(timewidget)),NamedTuple{(:signal,),Tuple{Signal{Time}}},typeof(timewidget),Time})
    precompile(Tuple{typeof(push!),TimeWidget{Time},Time})
    precompile(Tuple{typeof(push!),TimeWidget{DateTime},DateTime})
    ## datetimewidget
    precompile(Tuple{Core.kwftype(typeof(datetimewidget)),NamedTuple{(:signal,),Tuple{Signal{DateTime}}},typeof(datetimewidget),DateTime})

    # canvas
    for unit in (DeviceUnit, UserUnit)
        precompile(Tuple{typeof(canvas),Type{unit}})
        precompile(Tuple{typeof(mousedown_cb),  Ptr{GObject},Ptr{Gtk.GdkEventButton},MouseHandler{unit}})
        precompile(Tuple{typeof(mouseup_cb),    Ptr{GObject},Ptr{Gtk.GdkEventButton},MouseHandler{unit}})
        precompile(Tuple{typeof(mousemove_cb),  Ptr{GObject},Ptr{Gtk.GdkEventScroll},MouseHandler{unit}})
        precompile(Tuple{typeof(mousescroll_cb),Ptr{GObject},Ptr{Gtk.GdkEventScroll},MouseHandler{unit}})
        precompile(Tuple{typeof(init_pan_drag),Canvas{unit},Signal{ZoomRegion{RInt}}})
        precompile(Tuple{typeof(init_zoom_rubberband),Canvas{unit},Signal{ZoomRegion{RInt}}})
        precompile(Tuple{typeof(init_zoom_scroll),Canvas{unit},Signal{ZoomRegion{RInt}}})
        precompile(Tuple{typeof(init_pan_scroll),Canvas{unit},Signal{ZoomRegion{RInt}}})
        precompile(Tuple{typeof(Reactive.run_push),Signal{MouseScroll{unit}},MouseScroll{unit},Function})
        precompile(Tuple{typeof(draw),Function,Canvas{unit},Signal{Int},Vararg{Signal{Int},N} where N})
    end
    precompile(Tuple{Type{ZoomRegion},Tuple{Base.OneTo{Int},Base.OneTo{Int}}})
    precompile(Tuple{typeof(push!),Signal{ZoomRegion{RInt}},Tuple{UnitRange{Int},UnitRange{Int}}})
    # image_surface
    precompile(Tuple{typeof(image_surface),Matrix{RGB{N0f8}}})
    precompile(Tuple{typeof(image_surface),Matrix{RGBA{N0f8}}})
    precompile(Tuple{typeof(image_surface),Matrix{Gray{N0f8}}})
    # map
    # Note: @asserts on some of these cause problems on earlier versions of Julia
    precompile(Tuple{typeof(map),Function,Signal{ZoomRegion{RInt}}})
    precompile(Tuple{typeof(map),Function,Textbox{String},Textbox{Int}})
    precompile(Tuple{typeof(map),Function,Button})
    for unit in (DeviceUnit, UserUnit)
        precompile(Tuple{typeof(map),Function,Signal{MouseButton{unit}}})
        precompile(Tuple{typeof(map),Function,Signal{MouseScroll{unit}}})
    end
    # gc_preserve
    precompile(Tuple{typeof(gc_preserve),GtkWindowLeaf,Tuple{Dict{String,Any},Signal{Nothing},Signal{Nothing}}})
    precompile(Tuple{typeof(gc_preserve),GtkWindowLeaf,Dict{String,Dict{String,Any}}})
    precompile(Tuple{typeof(gc_preserve),GtkFrameLeaf,Player{PlayerWithTextbox}})
    # init_signal2widget
    precompile(Tuple{typeof(init_signal2widget),GtkScaleLeaf,UInt,Signal{Int}})
    precompile(Tuple{typeof(init_signal2widget),GtkScaleLeaf,UInt,Signal{Float64}})
    precompile(Tuple{typeof(init_signal2widget),GtkScaleLeaf,UInt,Signal{Float32}})
end


text:
"""
    entropy(P)

Computes the entropy of a distribution.
"""
function entropy(P; norm=false)
    if norm
        P /= sum(P)
    end
    return - sum(P[P.>0] .* log.(P[P.>0]))
end

"""
    relative_entropy(P, Q)

Computes the relative entropy of P to Q.
"""
relative_entropy(P, Q) = - sum(P .* log.(Q))

"""
    KL(P, Q)

Computes the Kullback-Leibler divergence of P to Q
"""
KL(P, Q) = relative_entropy(P, Q) - entropy(P)

"""
    utility(P, M)

Computes the average utility.
"""
function utility(P, M)
    return sum(P .* M)
end

function report_utility(P, M, bottom_species=nothing, top_species=nothing)
    @assert size(P) == size(M) "P and M should be of the same size"
    n, m = size(P)
    bottom_species == nothing && bottom_species == 1:n
    top_species == nothing && top_species == 1:m
    U = P .* M
    a, b = marginals(P)
    println("Bottom species")
    for i in 1:n
        println("\t$(bottom_species[i]): utility = $(sum(U[i,:]) / a[i]), entropy = $(
                                entropy(P[i,:] |> x -> x / sum(x)))")
    end
    println()
    println("Top species")
    for j in 1:m
        println("\t$(top_species[j]): utility = $(sum(U[:,j] / b[j])), entropy = $(
                                entropy(P[:,j] |> x -> x / sum(x)))")
    end
    println()
    println("Average utility = $(sum(U))")
    println("Average entropy = $(entropy(P))")
end

max_stars_repo_name: PoisotLab/OTSIN

max_issues_repo_name: PoisotLab/OTSIN

max_forks_repo_name: PoisotLab/OTSIN


text:
module JLDTest

export Object

mutable struct Object
    data::Int16
end

end


text:
using ITensors
using ITensorTDVP

n = 10
s = siteinds("S=1/2", n)

function heisenberg(n)
  os = OpSum()
  for j in 1:(n - 1)
    os += 0.5, "S+", j, "S-", j + 1
    os += 0.5, "S-", j, "S+", j + 1
    os += "Sz", j, "Sz", j + 1
  end
  return os
end

H = MPO(heisenberg(n), s)
ψ = randomMPS(s, "↑"; linkdims=10)

@show inner(ψ', H, ψ) / inner(ψ, ψ)

ϕ = tdvp(
  H,
  -1.0,
  ψ;
  nsweeps=20,
  reverse_step=false,
  normalize=true,
  maxdim=30,
  cutoff=1e-10,
  outputlevel=1,
)

@show inner(ϕ', H, ϕ) / inner(ϕ, ϕ)

e2, ϕ2 = dmrg(H, ψ; nsweeps=10, maxdim=20, cutoff=1e-10)

@show inner(ϕ2', H, ϕ2) / inner(ϕ2, ϕ2)

ϕ3 = ITensorTDVP.dmrg(H, ψ; nsweeps=10, maxdim=20, cutoff=1e-10, outputlevel=1)

@show inner(ϕ3', H, ϕ3) / inner(ϕ3, ϕ3)


text:
const block_signs = [0b1000 0b0010;
                     0b0100 0b0001]

const n_block = 16
const block_decode = Vector{Char}(undef, typemax(UInt16))
block_decode[0b0000 + 1] = ' '
block_decode[0b0001 + 1] = '▗'
block_decode[0b0010 + 1] = '▖'
block_decode[0b0011 + 1] = '▄'
block_decode[0b0100 + 1] = '▝'
block_decode[0b0101 + 1] = '▐'
block_decode[0b0110 + 1] = '▞'
block_decode[0b0111 + 1] = '▟'
block_decode[0b1000 + 1] = '▘'
block_decode[0b1001 + 1] = '▚'
block_decode[0b1010 + 1] = '▌'
block_decode[0b1011 + 1] = '▙'
block_decode[0b1100 + 1] = '▀'
block_decode[0b1101 + 1] = '▜'
block_decode[0b1110 + 1] = '▛'
block_decode[0b1111 + 1] = '█'
block_decode[(n_block+1):typemax(UInt16)] = unicode_table[1:(typemax(UInt16)-n_block)]

"""
The `BlockCanvas` is also Unicode-based.
It has half the resolution of the `BrailleCanvas`.
In contrast to BrailleCanvas, the pixels don't
have visible spacing between them.
This canvas effectively turns every character
into 4 pixels that can individually be manipulated
using binary operations.
"""
struct BlockCanvas <: LookupCanvas
    grid::Array{UInt16,2}
    colors::Array{ColorType,2}
    visible::Bool
    pixel_width::Int
    pixel_height::Int
    origin_x::Float64
    origin_y::Float64
    width::Float64
    height::Float64
    xscale::Union{Symbol,Function}
    yscale::Union{Symbol,Function}
end

@inline x_pixel_per_char(::Type{BlockCanvas}) = 2
@inline y_pixel_per_char(::Type{BlockCanvas}) = 2

@inline lookup_encode(c::BlockCanvas) = block_signs
@inline lookup_decode(c::BlockCanvas) = block_decode

BlockCanvas(args...; nargs...) = CreateLookupCanvas(BlockCanvas, args...; nargs...)

function char_point!(c::BlockCanvas, char_x::Int, char_y::Int, char::Char, color::UserColorType)
    if checkbounds(Bool, c.grid, char_x, char_y)
        c.grid[char_x,char_y] = n_block + char
        set_color!(c.colors, char_x, char_y, crayon_256_color(color))
    end
    c
end


text:
struct Stack{T, I, A <: AbstractAllocator{T, I}} # <: AbstractVector{T}
    alloc::A
    Stack{T, I, A}(alloc::A) where {T, I, A <: AbstractAllocator{T, I}} =
        new{T, I, A}(alloc)
end

Base.push!(s::Stack{T, I, A}, t::T) where {T, I, A <: AbstractAllocator{T, I}} =
    allocateend!(s.alloc, t)

Base.isempty(s::Stack{T, I, A}) where {T, I, A <: AbstractAllocator{T, I}} =
    isempty(s.alloc)

function Base.pop!(s::Stack{T, I, A}) where {T, I, A <: AbstractAllocator{T, I}}
    t = s.alloc[s.alloc.last]
    deallocateend!(s.alloc)
    t
end


text:
import Base: getindex, setindex!, get, delete!,
    length, iterate,
    getproperty, setproperty!,
    haskey, isempty,
    union,
    copy, push!

mutable struct OrdDict <: AbstractDict{Symbol,String}
    items::Vector{Pair{Symbol, String}}
end

OrdDict()  = OrdDict(Pair{Symbol, String}[])
OrdDict(items...) = OrdDict(collect(items))

# iteration definitions
length(a::OrdDict)     = length(a.items)
isempty(a::OrdDict)    = isempty(a.items)
iterate(a::OrdDict)    = iterate(a.items)
iterate(a::OrdDict, i) = iterate(a.items, i)

# getindex needed for dictionary-like key lookup
function getindex(a::OrdDict, key::Symbol)
    for item in a.items
        if first(item) == key
            return last(item)
        end
    end
    throw(KeyError(key))
end

# alternative multiple dispatch implementation by number
getindex(a::OrdDict, index::Integer) = a.items[index]

function setindex!(a::OrdDict, value::AbstractString, key::Symbol)
    for (i, item) in enumerate(a.items)
        if first(item) == key
            a.items[i] = key => value
            return
        end
    end
    push!(a.items, key => value)
end

push!(a::OrdDict, x::Pair{Symbol, String}) = push!(a.items, x)

function haskey(a::OrdDict, key::Symbol)
    for item in a.items
        if first(item) == key
            return true
        end
    end
    false
end

function delete!(a::OrdDict, key::Symbol)
    for (i, item) in enumerate(a.items)
        if first(item) == key
            deleteat!(a.items, i)
        end
    end
    a
end

# make it act like a dictionary
# need to keep the OrdDict.items to behave as usual
function getproperty(a::OrdDict, key::Symbol)
    if key == :items
        getfield(a, :items)
    else
        a[key]
    end
end

function setproperty!(a::OrdDict, key::Symbol, value::String)
    if key == :items
        setfield!(a, :item, value)
    else
        a[key] = value
    end
end

# julia> d = OrdDict(:one => "1", :two => "2")
# OrdDict with 2 entries:
#   :one => "1"
#   :two => "2"

# julia> d.one
# "1"

# julia> d.items
# 2-element Vector{Pair{Symbol, String}}:
#  :one => "1"
#  :two => "2"

# julia> d.two
# "2"

# julia> d.three = "3"
# "3"

# julia> d.three
# "3"

# julia> d.items
# 3-element Vector{Pair{Symbol, String}}:
#    :one => "1"
#    :two => "2"
#  :three => "3"

# julia> for (k,v) in d
#            println("key: ", k, " value: ", v)
#        end
# key: one value: 1
# key: two value: 2
# key: three value: 3

max_stars_repo_name: FourMInfo/JuliaforBeginners

max_issues_repo_name: FourMInfo/JuliaforBeginners

max_forks_repo_name: FourMInfo/JuliaforBeginners


text:
## General purpose functions

#User is expected to pass functions
#x_k(mhatk, hyperparam)  : provides a solution with data mhatk using hyperparam (e.g., p0, alpha)
#c_ik(x)
#These are often collected into arrays of functions
# xs			: [x_k for k = 1:K]
# c_k  			: [c_ik for i = 1:d]
# cs			: (d, K) Matrix with elements c_ik

#General purpose functions
###
function shrink(phat_k, p0, alpha, Nhat_k; alphaMax = 1e10)
	if Nhat_k == 0 || alpha > alphaMax
		return p0
	end
	@. alpha * p0 / (Nhat_k + alpha) + Nhat_k * phat_k / (Nhat_k + alpha)
end

#takes un-noramlized counts
#assumes out is correctly sized no bounds checking!
function _shrink!(mhat_k, p0, alpha, Nhat_k, out; alphaMax = 1e10)
	if Nhat_k == 0 || alpha > alphaMax
		return p0
	end
	for i = 1:length(p0)
		@inbounds out[i] = alpha * p0[i] / (Nhat_k + alpha) + mhat_k[i] / (Nhat_k + alpha)
	end
	out
end

#need to be careful because some problems may have zero data
function get_GM_anchor(mhats)
	Nhats = vec(sum(mhats, dims=1))
	non_zero_indx = Nhats .> 0
	vec(mean(mhats[:, non_zero_indx] ./ Nhats[non_zero_indx]', dims=2))
end

function sim_path(p_k, N::Int)
    mhats_k = zeros(length(p_k))

    if N == 0 #allow simulating zero data points
    	return mhats_k
    end

    for j in rand(Categorical(p_k), N)
        mhats_k[j] += 1
    end
    mhats_k
end

#ps is a matrix of pk
function sim_path(ps, Nhats)
    K = size(ps, 2)
    out = zeros(size(ps))
    for k = 1:K
        out[:, k] = sim_path(ps[:, k], Nhats[k])
    end
    out
end

#empirical bayes moment-matched estimates.
#returns p0, alpha0
#VG this method does not seem to consistenty yield alpha in [0, 1]
function eb_mm_estimates(mhats)
	Nhats = sum(mhats, dims=1)
	p0 = mean(mhats ./ Nhats, dims=2)
	K = size(mhats, 2)

	C0 = sum(mhats.^2)/K
	sq_norm_p0 = norm(p0)^2
	Nbar = mean(Nhats)
	var_n = mean( @.(Nhats^2 - Nhats) )

	alpha0 = C0 - sq_norm_p0 * var_n - Nbar
	@assert alpha0 > 0 "Moment matching fails"
	alpha0 /= (1 - sq_norm_p0) / var_n
	alpha0 = 1 / alpha0 - 1

	return p0, alpha0
end

#solves an approximation to the Apriori MSE equation for alpha
#for large K this is like AlphaOR for MSE
#for now, super lazy and just search a grid.
#returns AlphaEstimate, minimizingAlphaIndex
function mse_estimates(mhats, supps, p0, alpha_grid)
	K = size(mhats, 2)
	Nhats = sum(mhats, dims=1)
	phats = mhats ./ Nhats
	mu0s = vec(p0' * supps)

	#summary statistics
	muhats = Vector{Float64}(undef, K)
	sigmas = Vector{Float64}(undef, K)
	mse0 = Vector{Float64}(undef, K)
	for k = 1:K
		muhats[k] = dot(phats[:, k], supps[:, k])
		sigmas[k] = dot(phats[:, k], supps[:, k].^2) - muhats[k]^2
		mse0[k] = dot(phats[:, k], @.((supps[:, k] - mu0s[k])^2) )
	end

	#worker function
	function mse(a)
		out = 0.
		for k = 1:K
			out += (Nhats[k] - a^2) / (Nhats[k] + a)^2 * sigmas[k]
			out += ( a / (a + Nhats[k]) )^2 * mse0[k]
		end
		out/K
	end
	out = map(mse, alpha_grid)
	alpha_grid[argmin(out)], argmin(out)
end

#kth element of true perf
function z_k(x_k, c_k, mhat_k, ps_k, lam_k, lamavg, hyper_param)
	x = x_k(mhat_k, hyper_param)
	cs = map(c_ik -> c_ik(x), c_k)
	lam_k/lamavg * dot(ps_k, cs)
end

#actually returns zLOO_k * N * lamavg
function zLOO_k_unsc(x_k, c_k, mhat_k, hyper_param)
	#only compute for terms where mhat_k[i] > 0
	out = 0.
	for i = 1:length(mhat_k)
		if mhat_k[i] > 0
			mhat_k[i] -= 1
		    x = x_k(mhat_k, hyper_param)
		    mhat_k[i] += 1
		    out += mhat_k[i] * c_k[i](x)
		end
	end
	out
end

#average true performance
#hyper_param = p0, alpha
function zbar(xs, cs, mhats, ps, lams, hyper_param)
	lamavg = mean(lams)
	K = size(cs, 2)
	out = 0.
	#VG maybe change this to a more numerically stable way to do avg?
	for k = 1:K
		out += z_k(xs[k], cs[:, k], mhats[:, k], ps[:, k], lams[k], lamavg, hyper_param)
	end
	out/K
end

#full-info for scaling/comparison
#uses the fact that data scale doesn't matter.
#Only works for the regular xs (which is lame)
zstar(xs, cs, ps, lams) = zbar(xs, cs, ps, ps, lams, (ps[:, 1], 0.))

function zLOObar_unsc(xs, cs, mhats, hyper_param)
	K = size(cs, 2)
	out = 0.
	for k = 1:K
		@inbounds @views out += zLOO_k_unsc(xs[k], cs[:, k], mhats[:, k], hyper_param)
	end
	out /K
end

#actually returns N lam_avg * ZCVbar
function zCVbar_unsc(xs, cs, cv_data, hyper_param)
	out = 0.
	numFolds = length(cv_data)
	lams = ones(length(xs))

	for fold = 1:numFolds
		#form the training set and testing distribution
		#be careful about empty problems in testing
		train_data = sum(cv_data[setdiff(1:numFolds, fold)])

		#dispatch to the zbar method, again using scaling not important
		out += zbar(xs, cs, train_data, cv_data[fold], lams, hyper_param)
	end
	return out
end

#returns alphaOR, minimizingAlphaIndex, curveInAlpha
function oracle_alpha(xs, cs, mhats, ps, lams, p0, alpha_grid)
	alphaOR = 0.
	jstar = -1
	best_val = Inf
	out = zeros(length(alpha_grid))
	for (j, alpha) in enumerate(alpha_grid)
		out[j] = zbar(xs, cs, mhats, ps, lams, (p0, alpha))
		if out[j] < best_val
			jstar = j
			alphaOR = alpha
			best_val = out[j]
		end
	end
	return alphaOR, jstar, out
end

#returns alphaLOO, jstar, and cuveInAlpha_unsc
function loo_alpha(xs, cs, mhats, p0, alpha_grid)
	alphaLOO = 0.
	jstar = -1
	best_val = Inf
	out = zeros(length(alpha_grid))

	for (j, alpha) in enumerate(alpha_grid)
		out[j] = zLOObar_unsc(xs, cs, mhats, (p0, alpha))
		if out[j] < best_val
			jstar = j
			alphaLOO = alpha
			best_val = out[j]
		end
	end
	return alphaLOO, jstar, out
end

function cv_alpha(xs, cs, mhats, p0, alpha_grid, numFolds)
	#divy up the data
	cv_data = split_cv(mhats, numFolds)

	alphaCV = 0.
	jstar = -1
	best_val = Inf
	out = zeros(length(alpha_grid))
	for (j, alpha) in enumerate(alpha_grid)
		out[j] = zCVbar_unsc(xs, cs, cv_data, (p0, alpha))
		if out[j] < best_val
			jstar = j
			alphaCV = alpha
			best_val = out[j]
		end
	end
	return alphaCV, jstar, out
end

#########
#Functions to optimize anchor among beta distributions
########
#perform a naive discretization into d+1 points (d bins)
#p_i = Prob((i-1)/(d) <= Beta i <= i /d ) i = 1:d
#assumes p0_out is of size d
function formBetaAnchor!(theta1, theta2, d, p0_out)
	@assert length(p0_out) == d "p0_out must be of length $d"
	dist = Distributions.Beta(theta1, theta2)
	for i = 1:d
		p0_out[i] = cdf(dist, i/d) - cdf(dist, (i-1)/d)
	end
	p0_out
end

#Fits a beta(theta1, theta2) distribution for anchor (discretized)
#returns alphaLOO, p0, best_loo that optimizes loo error
function loo_betaAnchor(xs, cs, mhats, alpha_grid, theta2_grid, mu_grid; info=false)
	alphaLOO = 0.
	theta1LOO = -Inf
	theta2LOO = -Inf
	best_val = Inf
	d = size(mhats, 1)
	p0 = zeros(d)

	for theta2 in theta2_grid
		#theta1 = mu/(1-mu) * theta2
		#since mu in [0, 1] we search a scaled grid
		for mu in mu_grid
			theta1 = mu / (1-mu) * theta2

			#form p0
			formBetaAnchor!(theta1, theta2, d, p0)

			for alpha in alpha_grid
				out = zLOObar_unsc(xs, cs, mhats, (p0, alpha))
				if out < best_val
					alphaLOO = alpha
					theta1LOO = theta1
					theta2LOO = theta2
					best_val = out
				end
			end #over alpha
		end #over mu  (aka theta1)
	end#over theta 2

    #confirm it does better than mm
    theta1, theta2 = mm_beta_param(get_GM_anchor(mhats))
   	formBetaAnchor!(theta1, theta2, d, p0)
   	alpha, jstar, curve = loo_alpha(xs, cs, mhats, p0, alpha_grid)
   	out = curve[jstar]
   	if out < best_val
   		return alpha, p0, out
   	end

	info && println("theta1Loo:\t", theta1LOO, "\t theta2LOO:\t", theta2LOO, "\tMean:\t", theta1LOO/(theta1LOO + theta2LOO))

	#compute the p0 one more time (hopefully fast)
	formBetaAnchor!(theta1LOO, theta2LOO, d, p0)
	return alphaLOO, p0, best_val
end

#Fits a beta(theta1, theta2) distribution for anchor (discretized)
#returns alphaOR, p0, perf that optimizes oracle error
function oracle_betaAnchor(xs, cs, mhats, ps, lams, alpha_grid, theta2_grid, mu_grid; info=false)
	alphaOR = 0.
	theta1OR = -Inf
	theta2OR = -Inf
	best_val = Inf
	d = size(mhats, 1)
	p0 = zeros(d)

	for theta2 in theta2_grid
		#theta1 = mu/(1-mu) * theta2
		#since mu in [0, 1] we search a scaled grid
		for mu in mu_grid
			theta1 = mu / (1-mu) * theta2

			#form p0
			formBetaAnchor!(theta1, theta2, d, p0)

			for alpha in alpha_grid
				out = zbar(xs, cs, mhats, ps, lams, (p0, alpha))
				if out < best_val
					alphaOR = alpha
					theta1OR = theta1
					theta2OR = theta2
					best_val = out

					info && println("alpha\t", alphaOR, "\t theta1Loo:\t", theta1OR, "\t theta2OR:\t", theta2OR, "\tMean:\t", theta1OR/(theta1OR + theta2OR))
				end
			end #over alpha
		end #over mu  (aka theta1)
	end#over theta 2

    #confirm it does better than mm
    theta1, theta2 = mm_beta_param(get_GM_anchor(mhats))
   	formBetaAnchor!(theta1, theta2, d, p0)
   	alpha, jstar, curve = oracle_alpha(xs, cs, mhats, ps, lams, p0, alpha_grid)
   	out = curve[jstar]
   	if out < best_val
   		return alpha, p0, out
   	end

	#compute the p0
	info && println("theta1Loo:\t", theta1OR, "\t theta2OR:\t", theta2OR, "\tMean:\t", theta1OR/(theta1OR + theta2OR))
	formBetaAnchor!(theta1OR, theta2OR, d, p0)
	return alphaOR, p0, best_val
end

#compute the beta dist param by mm to mean and variance
function mm_beta_param(pGM)
	supp = range(0, stop=1, length=length(pGM))
	mu = dot(pGM, supp)
	var = dot(pGM, supp.^2) - mu^2
	if var > mu * (1 - mu)
		#fail silently
		return (1., 1.)
	end
	v = mu * (1 - mu) / var
	return mu * v, (1 - mu) * v
end

#Fits a beta(theta1, theta2) distribution for anchor (discretized)
#returns alphaLOO, p0, best_loo that optimizes loo error
#uses a black-box optimization routine
function loo_betaAnchor(xs, cs, mhats; info=false, time_limit = 45, iterations=1000,
						init_alpha=1., init_theta1= 1., init_theta2=1.)

	d = size(mhats, 1)
	p0 = zeros(d) #used as shared storage

	#ys = [alpha, theta1, theta2]
	function f(ys)
		alpha = ys[1]^2
		theta1 = ys[2]^2
		theta2 = ys[3]^2
		#form p0
		formBetaAnchor!(theta1, theta2, d, p0)
		zLOObar_unsc(xs, cs, mhats, (p0, alpha))
	end
	ys1 = [init_alpha, init_theta1, init_theta2]
	info && println("Initial Value:\t", f(ys1))
    t= @elapsed results = optimize(f, ys1,
                    ParticleSwarm(n_particles=40),
                    Optim.Options(time_limit=time_limit, iterations=iterations))

    z1 = Optim.minimum(results)
    ystar = Optim.minimizer(results)
    info && println("Theta1:\t", ystar[2]^2, "\t Theta2:\t", ystar[3]^2, "\t mu:\t", (ystar[2]^2)/(ystar[2]^2 + ystar[3]^2))
	info && println("Time:\t", t, "\tTime Per Iteration:\t", t/iterations(results))
	formBetaAnchor!(ystar[2]^2, ystar[3]^2, d, p0)

	return ystar[1]^2, p0, z1
end

#### Helper function for binning data
#dat_vec is vector of actual realizations
#computes a grid of size d
#returns
	#d vector of lefthand binpoints
	#d vector of counts
	#raw histogram object for safety
#assumes NA's already filtered from calcualtion
function bin_data(dat_vec, d; TOL = .001)
	#Explicitly make bins to prevent histogram from rounding
	l = minimum(dat_vec)
	u = maximum(dat_vec)
	bin_width = (u - l)/d
    bin_edges = collect(range(l, stop=u, length=d + 1))
    bin_edges[end] += TOL * bin_width  #inflate the last one a little...
    bin_edges[1] -= TOL * bin_width
    dat_hist = fit(Histogram, dat_vec, bin_edges, closed =:left)
    bin_edges[1:d], dat_hist.weights, dat_hist
end

###Helper function for splitting data for k-fold cross-validation
#This is not a memory efficient implementation
#Just lazy split/iterations
function split_cv(mhats, numFolds)
	d, K = size(mhats)
    Ntot = convert(Int, sum(mhats))
    splits = mod.(1:Ntot, numFolds) .+ 1
    shuffle!(splits)
    cv_data = [zero(mhats) for fold = 1:numFolds]

    #iterate through each k in order
    ix_split = 1

    for k = 1:K
        for i = 1:d
            for j = 1:mhats[i, k]
                cv_data[splits[ix_split]][i, k] += 1
                ix_split += 1
            end
        end
    end
    return cv_data
end