text:
"""
    hilbert_x(eb::Eigenbrot)
Return the Hilbert x-transform of `eb`.
"""
function hilbert_x(eb::Eigenbrot)
    nd = copy(eb)
    h, w = size(nd)
    centre_c = 1 + w >> 1
    centre_r = 1 + h >> 1
    for c in 1:w
        for r in 1:h
            sign = if c == centre_c
                r > centre_r ? im : -im
            elseif r == centre_r
                c > centre_c ? im : -im
            else
                c > centre_c ? im : -im
            end
            nd[r, c] *= sign
        end
    end
    return nd
end

"""
    hilbert_y(eb::Eigenbrot)
Return the Hilbert x-transform of `eb`.
"""
function hilbert_y(eb::Eigenbrot)
    nd = copy(eb)
    h, w = size(nd)
    centre_c = 1 + w >> 1
    centre_r = 1 + h >> 1
    for c in 1:w
        for r in 1:h
            sign = if c == centre_c
                r > centre_r ? im : -im
            elseif r == centre_r
                c > centre_c ? im : -im
            else
                r > centre_r ? im : -im
            end
            nd[r, c] *= sign
        end
    end
    return nd
end

"""
    fft_xshear(eb::Eigenbrot, xShear::Real, expand::Bool = true)
Shear `eb` by a factor 'xShear' in the x-direction using the Fourier
shift theorem, and return the result. If `expand` is true the input
is padded so that the shear doesn't wrap the data as a result of
aliasing in the Fourier domain.
"""
function fft_xshear(eb::Eigenbrot, xShear::Real, expand::Bool = true)
    xform = if expand
        extra = ceil(Int, abs(xShear) * height(eb) / 2)
        fftx(pad(eb, left = extra, right = extra), false)
    else
        fftx(eb, false)
    end
    h, w = size(xform)
    halfHeight = h / 2
    lim = iseven(w) ? w >> 1 : (w + 1) >> 1
    for row in 1:h
        δ = xShear * (halfHeight - row + 1)
        for i in 1:lim
            θ = (2.0 * π * δ * i) / w
            cs = cos(θ) + im * sin(θ)
            xform[row, i + 1] *= cs
            xform[row, w - i + 1] *= cs'
        end
        if iseven(w)
            i2 = trunc(Int, abs(δ))
            δ < 0 && (i2 = -i2)
            θ = (2.0 * π * δ * i2) / w
            xform[row, lim] *= cos(θ) + im * sin(θ)
        end
    end

    return fftx(xform, false)
end

"""
    fft_yshear(eb::Eigenbrot, xShear::Real, expand::Bool = true)
Shear `eb` by a factor 'yShear' in the y-direction using the Fourier
shift theorem, and return the result. If `expand` is true the input
is padded so that the shear doesn't wrap the data as a result of
aliasing in the Fourier domain.
"""
function fft_yshear(eb::Eigenbrot, yShear::Real, expand::Bool = true)
    xform = if expand
        extra = ceil(Int, abs(yShear * width(eb) / 2))
        ffty(pad(eb, top = extra, bottom = extra), false)
    else
        ffty(eb, false)
    end
    h, w = size(xform)
    halfWidth = w / 2
    lim = iseven(h) ? h >> 1 : (h + 1) >> 1
    for col in 1:w
        δ = yShear * (halfWidth - col + 1)
        for i in 1:lim
            θ = (2.0 * π * δ * i) / h
            cs = cos(θ) + im * sin(θ)
            xform[i + 1, col] *= cs
            xform[h - i + 1, col] *= cs'
        end
        if iseven(h)
            i2 = trunc(Int, abs(δ))
            δ < 0.0 && (i2 = -i2)
            θ = (2.0 * π * i2 * lim) / h
            xform[lim + 1, col] *= cos(θ) + im * sin(θ)
        end
    end

    return ffty(xform, false)
end

"""
    rotate90(eb::Eigenbrot, n::Integer)
Return the result of rotating `eb` through an angle of
`n` * π /2 (i.e. `n` right angles).
"""
function rotate90(eb::Eigenbrot, n::Integer)
    if n == 0
        return copy(eb)
    elseif n == 1
        h, w = size(eb)
        e2 = Eigenbrot(w, h, eb.fft)
        for r in 1:w
            e2.vals[r, 1:h] = eb.vals[h:-1:1, r]
        end
        return e2
    elseif n == 2
        e2 = Eigenbrot(height(eb), width(eb), eb.fft)
        l = length(eb.vals)
        e2.vals[1:l] = eb.vals[l:-1:1]
        return e2
    elseif n == 3
        h, w = size(eb)
        e2 = Eigenbrot(w, h, eb.fft)
        for r in 1:w
            e2.vals[w - r + 1, h:-1:1] = eb.vals[h:-1:1, r]
        end
        return e2
    else
        return rotate90(eb, mod(n, 4))
    end
end

"""
    fourierRotation(eb::Eigenbrot, θ::Real, expand::Bool = true)
Return the result of rotating `eb` through an angle `θ`,
using the algorithm described by Larkin, Oldfield and Klemm
(Optics Communications 139 (1997) 99-106). If `expand` is true,
the size of the output image is increased to contain the rotated image,
otherwise the aliasing described in Larkin et al. occurs.
"""
function fourierRotation(eb::Eigenbrot, θ::Real, expand::Bool = true)
    h, w = size(eb)
    θ = -θ  # reverse rotation direction

    # Deconstruct θ as the sum of an integral number, n90,
    # of right angles, and a remainder, ϕ, with absolute
    # value less than π/4...
    n90, ϕ = divrem(θ, π / 2)
    if abs(ϕ) > π / 4
        s = sign(ϕ)
        n90 += s
        ϕ = ϕ - s * π / 2
    end
    n90 = trunc(Int, mod(n90, 4))
    result = rotate90(eb, n90)
    ϕ ≈ 0 && return result

    xshear = tan(ϕ / 2)
    yshear = -sin(ϕ)
    if expand
        hh, ww = size(result)
        hm, wm = 0, 0
        msx = [1.0 xshear; 0.0 1.0]
        msy = [1.0 0.0; yshear 1.0]
        cnrs = [0 ww 0 ww; 0 0 hh hh]
        for mm in [msx, msy, msx]
            cnrs = mm * cnrs
            wm = max(wm, maximum(cnrs[1, :]) - minimum(cnrs[1, :]))
            hm = max(hm, maximum(cnrs[2, :]) - minimum(cnrs[2, :]))
        end
        lr = ceil(Int, (wm - ww) / 2)
        tb = ceil(Int, (hm - hh) / 2)
        result = pad(result, left = lr, right = lr, top = tb, bottom = tb)
    end

    result = fft_xshear(result, xshear, false)
    result = fft_yshear(result, yshear, false)
    result = fft_xshear(result, xshear, false)
    return result
end


text:
using AlgebraOfGraphics
using Makie
using DataFrames

function plot_results(results::DataFrame;y=:estimate,color=:coefname,col=:basisname,row=:group,stderror=false,pvalue = DataFrame(:from=>[],:to=>[],:pval=>[]))
    results = deepcopy(results)

    # replace missing/nothing valuess
    results.group = results.group .|> a -> isnothing(a) ? :fixef : a
    results.stderror = results.stderror .|> a -> isnothing(a) ? 0. : a

    results[!,:se_low]  = results[:,y] .- results.stderror
    results[!,:se_high] = results[:,y] .+ results.stderror

    m = mapping(color=color,col=col,row=row)

    m_li = mapping(:time,y)

    basic =  visual(Lines) * m_li

    if stderror
        m_se = mapping(:time,:se_low,:se_high)
        basic = basic + visual(Band,alpha=0.5)*m_se
    end
    basic = basic*data(results)

    # add the pvalues
    if !isempty(pvalue)
        p = deepcopy(pvalue)

        # for now, add them to the fixed effect
        p[!,layout] .= :fixef

        # rename to match the res-dataframe
        p[!,:coefname] = p.coefname

        # probably not needed
        #p[!,:colname_basis] .= results.colname_basis[1]

        #
        un = unique(p.coefname)
        scaleY = [minimum(results.estimate),maximum(results.estimate)]
        stepY = scaleY[2]-scaleY[1]
        posY = stepY*-0.05+scaleY[1]
        Δt = diff(results.colname_basis[1:2])[1]
        Δy = 0.01
        p[!,:color] .=  [findfirst(un .== x) for x in p.coefname]
        p[!,:segments] = [Rect(Vec(x,posY+stepY*(Δy*(n-1))),Vec(y-x+Δt,0.5*Δy*stepY)) for (x,y,n) in zip(p.from,p.to,p.color)]

        basic =  basic + (data(p)*mapping(:segments)*visual(Poly))

    end
    d = basic*m|> draw
    return d

end

ext: jl


text:
using HMMs
using LinearAlgebra
using Distributions

mean = zeros(2)
cov = zeros(2,2) + I
A = [zeros(1,1) ones(1,1)]
b = [0.0]
(x,tol) = HMMs.lenormalcubature( mean, cov, A, b )
@assert( abs(x - 0.5) < tol )

A = zeros(2,2) + I
b = zeros(2)
(x,tol) = HMMs.lenormalcubature( mean, cov, A, b )
@assert( abs(x - 0.25) < tol )

b = [0.0,1.0]
(x,tol) = HMMs.lenormalcubature( mean, cov, A, b )
n1 = cdf( Normal(0,1), 1 )
@assert( abs(x - n1/2) < tol )

b = [1.0,1.0]
(x,tol) = HMMs.lenormalcubature( mean, cov, A, b )
@assert( abs(x - n1^2) < tol )

A = [1.0 0.0;-1.0 -1.0]
b = zeros(2)
@time (x,tol) = HMMs.lenormalcubature( mean, cov, A, b, maxevals=10^7 )
# tol doesn't seem to work for this one
@assert( abs(x - 1/8) < 1e-5 )


text:
# ints_trap = map(hs) do h
#     x = y1:h:y2
#     -(g(x[1]) + g(x[end]))*h/2.0 + sum(g.(x)*h)
# end

function integration_scheme(x::AbstractVector{T}; scheme::Symbol = :trapezoidal, kws...) where T<:AbstractFloat
    if length(collect(x)) == 1
        return [0.]
    end

    if scheme == :trapezoidal
        trap_scheme(x; kws...)
    elseif scheme == :simpson
        simpson_scheme(x; kws...)
    else
        @warn("Integration scheme $scheme unknown, will use trapezoidal")
        trap_scheme(x; kws...)
    end
end

function trap_scheme(x::AbstractVector{T}; x0::T = first(x), xn::T = last(x)) where T<:AbstractFloat

    inds = axes(x,1)

    h = (x[inds[2]]-x[inds[1]])
    σs = similar(x)
    σs .= h

    σs[inds[1]] -= h/T(2)

    σs[end] -= h/T(2)
    # accounts for the ends
    σs[inds[1]] +=   (x[inds[1]] - x0)*(x[inds[2]] - x0 + h)/(2.0h)
    σs[inds[2]] +=  -(x0 - x[inds[1]])^T(2)/(T(2)*h)
    σs[end] +=  (xn - x[end])*(xn - x[end-1] + h)/(2h)
    σs[end-1] += - (xn - x[end])^T(2)/(T(2)*h)

    return σs
end

function simpson_scheme(x::AbstractVector{T}; x0::T = first(x), xn::T = last(x)) where T<:AbstractFloat
    inds = axes(x,1)

    if iseven(length(inds))
        @warn("Simpson's integration scheme is designed for an odd number of mesh points")
    end

    h = (x[inds[2]]-x[inds[1]])
    σs = similar(x).*zero(T)
    for j in 1:length(inds)
        σs[inds[j]] = iseven(j) ? 4.0 : 2.0
    end
    σs = σs.*(h/3.0)

    # accounts for the ends
    σs[inds[1]] = h/3.0
    σs[inds[1]] +=  (x[inds[1]] - x0)*(x[inds[2]] - x0 + h)/(2.0h)
    σs[inds[2]] +=  - (x0-x[inds[1]])^2.0/(2.0h)
    σs[inds[end]] = h/3.0
    σs[inds[end]] +=  (xn-x[inds[end]])*(xn - x[inds[end-1]] + h)/(2h)
    σs[inds[end-1]] += - (xn-x[inds[end]])^2.0/(2h)

    return σs
end


text:
@testset "Scalar operators" begin
    B = FEDVR(1.0:7, 4)
    B̃ = B[:,2:end-1]
    r = axes(B,1)

    V = B'*QuasiDiagonal(identity.(r))*B
    @test V == Diagonal(B.x)

    Ṽ = apply(*, B̃', QuasiDiagonal(identity.(r)), B̃)
    @test Ṽ == Diagonal(B.x[2:end-1])
end

@testset "Mass matrices and inverses" begin
    t = range(0,stop=20,length=5)
    @testset "Order $k" for k ∈ [2,5]
        R = FEDVR(t, k)
        R̃ = R[:,2:end-1]

        m = size(R,2)
        n = size(R̃,2)

        @testset "Mass matrices" begin
            d = R'R
            @test d == Diagonal(Ones(m))

            d̃ = R̃'R̃
            @test d̃ == Diagonal(Ones(n))

            @test R'R̃ == (R̃'R)' == BandedMatrix((-1 => ones(Int,n),), (m,n), (1,-1))
            @test R̃'R[:,1:end-1] == BandedMatrix((1 => ones(Int,n),), (n,m-1), (-1,1))
        end

        @testset "Inverses" begin
            R⁻¹ = pinv(R)
            # R̃⁻¹ = pinv(R̃)

            @test R⁻¹*R === I
            @test R*R⁻¹ === I

            cu = rand(size(R,2))
            cv = rand(size(R,2))
            cuv = [cu cv]

            u = R*cu
            v = R*cv
            uv = R*cuv

            @test R⁻¹*u === cu
            @test R⁻¹*v === cv
            @test R⁻¹*uv === cuv

            ut = u'
            # Fails with: ERROR: MethodError: no method matching axes(::UniformScaling{Float64}, ::Int64)
            # @test ut*R⁻¹' === ut.args[1]

            @warn "Need to implement/test basis inverses for restricted bases"
        end
    end
end


text:
(g::GCNConv)(L̃::AbstractMatrix, X::CuArray) = g(cu(L̃), X)

(g::GCNConv)(L̃::CuArray, X::CuArray) = g.σ.(g.weight * X * L̃ .+ g.bias)

(c::ChebConv)(L̃::AbstractMatrix, X::CuArray) = c(cu(L̃), X)

function (c::ChebConv)(L̃::CuArray, X::CuArray)
    @assert size(X, 1) == c.in_channel "Input feature size must match input channel size."
    @assert size(X, 2) == size(L̃, 1) "Input vertex number must match Laplacian matrix size."

    Z_prev = X
    Z = X * L̃
    Y = view(c.weight,:,:,1) * Z_prev
    Y += view(c.weight,:,:,2) * Z
    for k = 3:c.k
        Z, Z_prev = 2*Z*L̃ - Z_prev, Z
        Y += view(c.weight,:,:,k) * Z
    end
    return Y .+ c.bias
end

# Avoid ambiguity
update_batch_edge(g::GATConv, adj, E::Fill{S,2,Axes}, X::CuMatrix, u) where {S,Axes} = update_batch_edge(g, adj, X)

update_batch_vertex(g::GATConv, M::CuMatrix, X::CuMatrix, u) = update_batch_vertex(g, M)


text:
module CVGSurfaceFile

    include("CVGStruct.jl")
    export cvgseek, save_cvgfile, cvgseek2
    export CVGSurface, Triangle, TrBoundry, Surface
    #将原始的.dat文件转化为自定义类型
    function cvgseek(file)
        filestream = open(file)

        fileline = readline(filestream)

        n_max_vertex, n_vertices, n_triangles = parse.(Int,split(fileline))

        triangles = Triangle[]
        trBoundry = TrBoundry[]

        for i in 1:n_vertices
            pos = position(filestream)
            seek(filestream, pos)
            string_triangleID, string_xyz... = split(readline(filestream),' ')
            triangleID = parse(Int,string_triangleID)
            x, y, z =  parse.(Float64,string_xyz)
            push!(triangles,Triangle(triangleID, x, y, z))
        end

        for i in 1:n_triangles
            pos = position(filestream)
            seek(filestream, pos)
            triangleID1, triangleID2, triangleID3, boundryID = parse.(Int,split(readline(filestream), ' '))
            push!(trBoundry, TrBoundry(triangleID1, triangleID2, triangleID3, boundryID))
        end

        return Surface(n_max_vertex, n_vertices, n_triangles, triangles, trBoundry)

    end

    function cvg_split(ios::IOStream)
        split(readline(ios),' ')
    end

    function cvgseek2(file)
        filestream = open(file)

        fileline = readline(filestream)

        n_max_vertex, n_vertices, n_triangles = parse.(Int,split(fileline))

        triangles = Vector{Triangle}(undef,n_vertices)
        trBoundry = Vector{TrBoundry}(undef,n_triangles)

        for i in 1:n_vertices
            pos = position(filestream)
            seek(filestream, pos)
            string_triangleID, string_xyz... = cvg_split(filestream)
            triangleID = parse(Int,string_triangleID)
            x, y, z =  parse.(Float64,string_xyz)
            triangles[i] = Triangle(triangleID, x, y, z)
        end

        for i in 1:n_triangles
            pos = position(filestream)
            seek(filestream, pos)
            triangleID1, triangleID2, triangleID3, boundryID = parse.(Int, cvg_split(filestream))
            trBoundry[i]  = TrBoundry(triangleID1, triangleID2, triangleID3, boundryID)
        end

        return Surface(n_max_vertex, n_vertices, n_triangles, triangles, trBoundry)

    end

    #给string函数添加方法
    Base.string(cvgsurface::Triangle) = string(cvgsurface.triangleID, " ", cvgsurface.x, " ", cvgsurface.y, " ", cvgsurface.z)
    Base.string(cvgTrBoundry::TrBoundry) = string(cvgTrBoundry.triangleID1, " ", cvgTrBoundry.triangleID2, " ", cvgTrBoundry.triangleID3, " ", cvgTrBoundry.boundryID)

    #将自定义数据转化为字符串数组数据
    function cvg_surface_string(surface::CVGSurface)
        arr_string = String[]
        push!(arr_string,string(surface.n_max_vertex, " ", surface.n_vertices, " ",surface.n_triangles))
        append!(arr_string,string.(surface.triangles))
        append!(arr_string,string.(surface.surfaces))
        return arr_string
    end

    #将自定义数据直接保存为.dat文件
    function save_cvgfile(file::String, surface::CVGSurface)
        open(file, "w") do io
            write(io, string(surface.n_max_vertex, " ", surface.n_vertices, " ",surface.n_triangles, '\n'))
            for i in string.(surface.triangles)
                write(io, i * '\n')
            end

            for i in string.(surface.surfaces)
                write(io, i * '\n')
            end
        end
    end

end

max_stars_repo_name: powerzf/CVGSurfaceFile

max_issues_repo_name: powerzf/CVGSurfaceFile

max_forks_repo_name: powerzf/CVGSurfaceFile


text:
using CachedArrays
using Test
using Random
using BenchmarkTools
using TimerOutputs
import MacroTools
using Distributions
import ProgressMeter

# Helper Functions
inlocal(manager, A; kw...) = inpool(manager, A, CachedArrays.Local; kw...)
inremote(manager, A; kw...) = inpool(manager, A, CachedArrays.Remote; kw...)

function inpool(manager, A, pool; primary_only = true)
    block = CachedArrays.metadata(A)
    id = CachedArrays.getid(block)
    return in(id, CachedArrays.visible_ids(manager, pool; primary_only))
end

include("utils/utils.jl")
include("utils/findnexttree.jl")
include("utils/lru.jl")
include("utils/freebuffer.jl")

include("allocators.jl")

include("memory/block.jl")
include("memory/freelist.jl")
include("memory/compactheap.jl")
include("memcpy.jl")

include("manager.jl")
include("array/array.jl")
include("lib.jl")
include("annotation.jl")

include("integration/corner.jl")
include("integration/stress.jl")


text:
## Grouping machinery

fast_hashed(v::AbstractVector) = isbitstype(eltype(v)) ? PooledArray(v) : v

function permutation_ranges(cols)
    isempty(cols) && return nothing, [:]
    grouping_sa = StructArray(map(fast_hashed, cols))
    gp = GroupPerm(grouping_sa)
    return sortperm(gp), collect(gp)
end

allvariables(e::Entry) = (e.primary..., e.positional..., e.named...)
allvariables(l::Layer) = (l.positional..., l.named...)

function shape(x::Union{Entry, Layer})
    arrays = map(var -> var isa ArrayLike ? var : fill(nothing), allvariables(x))
    return Broadcast.combine_axes(arrays...)
end

function assert_equal(a, b)
    if !isequal(a, b)
        msg = "$a and $b must be equal"
        throw(ArgumentError(msg))
    end
    return a
end

getuniquevalue(v) = reduce(assert_equal, v)

haszerodims(::AbstractArray) = false
haszerodims(::AbstractArray{<:Any, 0}) = true
haszerodims(::Ref) = true
haszerodims(::Tuple) = false

function getnewindex(u, c)
    v = Broadcast.broadcastable(u)
    return v[Broadcast.newindex(v, c)]
end

function subgroups(vs, perm, rgs, axs)
    return map(Iterators.product(rgs, CartesianIndices(axs))) do (rg, c)
        v = getnewindex(vs, c)
        return haszerodims(v) || rg === (:) ? v : view(v, perm[rg])
    end
end

shiftdims(v::AbstractArray) = reshape(v, 1, axes(v)...)
shiftdims(v::Tuple) = shiftdims(collect(v))
shiftdims(v) = v

function group(entry::Entry)
    grouping = foldl(entry.primary, init=()) do acc, v
        return haszerodims(v) ? (acc..., only(v)) : acc
    end
    perm, rgs = permutation_ranges(grouping)
    axs = shape(entry)
    primary, positional, named = map((entry.primary, entry.positional, entry.named)) do tup
        return map(vs -> subgroups(vs, perm, rgs, axs), tup)
    end
    labels = copy(entry.labels)
    map!(shiftdims, values(labels))
    return Entry(entry; primary, positional, named, labels)
end

function separate(nt::NamedTuple)
    primary_keys = filter(key -> hascategoricalentry(nt[key]), keys(nt))
    primary = NamedTuple{primary_keys}(nt)
    return primary, Base.structdiff(nt, primary)
end

function getlabeledarray(layer::Layer, s)
    data, axs = layer.data, shape(layer)
    isdims = s isa DimsSelector || s isa Pair && first(s) isa DimsSelector
    if isdims
        vs, (f, label) = select(data, s)
        d = only(vs) # multiple dims selectors in the same mapping are disallowed
        sz = ntuple(length(axs)) do n
            return n in d.dims ? length(axs[n]) : 1
        end
        arr = map(fill∘f, CartesianIndices(sz))
    elseif isnothing(data)
        vs, (f, label) = select(data, s)
        isprim = any(hascategoricalentry, vs)
        arr = isprim ? map(fill∘f, vs...) : map(x -> map(f, x...), zip(vs...))
    else
        selector = s isa AbstractArray ? s : fill(s)
        labeled_arr = map(selector) do s
            local vs, (f, label) = select(data, s)
            return label, map(f, vs...)
        end
        label, arr = map(first, labeled_arr), map(last, labeled_arr)
    end
    return label, arr
end

function process_mappings(layer::Layer)
    labels = Dict{KeyType, Any}()
    positional, named′ = map((layer.positional, layer.named)) do tup
        return mapkeys(tup) do key
            label, arr = getlabeledarray(layer, tup[key])
            labels[key] = label
            return arr
        end
    end
    primary, named = separate(named′)
    return Entry(; primary, positional, named, labels)
end

"""
    to_entry(layer::Layer)

Convert `layer` to equivalent entry, excluding transformations.
"""
function to_entry(layer::Layer)
    entry = process_mappings(layer)
    grouped_entry = isnothing(layer.data) ? entry : group(entry)
    primary = map(vs -> map(getuniquevalue, vs), grouped_entry.primary)
    return Entry(grouped_entry; primary)
end

function process(layer::Layer)
    init = to_entry(layer)
    return foldl(|>, layer.transformations; init)
end


text:
export encode_icm_cuda, train_lsq_cuda

"Encodes a database with ILS in cuda"
function encode_icm_cuda_single(
  RX::Matrix{Float32},         # in. The data to encode
  B::Matrix{Int16},            # in. Initial list of codes
  C::Vector{Matrix{Float32}},  # in. Codebooks
  ilsiters::Vector{Int64},     # in. ILS iterations to record Bs and obj function. Its max is the total number of iteration we will run.
  icmiter::Integer,            # in. Number of ICM iterations
  npert::Integer,              # in. Number of entries to perturb
  randord::Bool,               # in. Whether to randomize the order in which nodes are visited in ILS
  V::Bool=false)               # in. Whether to print progress

  d, n = size( RX )
  m    = length( C )
  _, h = size( C[1] )

  # Number of results to keep track of
  nr   = length( ilsiters )

  # Make space for the outputs
  Bs   = Vector{Matrix{Int16}}(undef, nr)
  objs = Vector{Float32}(undef, nr)

  # === Compute binary terms (products between all codebook pairs) ===
  binaries, cbi = get_binaries(C)
  _, ncbi       = size(cbi)

  # Create a transposed copy of the binaries for cache-friendliness
  binaries_t = similar(binaries)
  for j = 1:ncbi
    binaries_t[j] = collect(binaries[j]')
  end

  # Create an index from codebook pairs to indices
  cbpair2binaryidx = zeros(Int32, m, m)
  for j = 1:ncbi
    cbpair2binaryidx[ cbi[1,j], cbi[2,j] ] = j
  end

  dev = CuDevice(0)
  ctx = CuContext(dev)

  # Initialize the cuda module, and choose the GPU
  gpuid = 0
  CudaUtilsModule.init( gpuid, cudautilsptx )
  # CUDArt.device( devlist[gpuid] )

  # === Create a state for random number generation ===
  if V; @printf("Creating %d random states... ", n); start_time = time_ns(); end
  d_state = CUDAdrv.Mem.alloc(n*64)
  CudaUtilsModule.setup_kernel(cld(n, 1024), 1024, Cint(n), d_state)

  CUDAdrv.synchronize(ctx)
  if V; @printf("done in %.2f seconds\n", (time_ns() - start_time)/1e9); end

  # Measure time for encoding
  if V; start_time = time_ns(); end

  # Copy X and C to the GPU
  d_RX = CuArrays.CuArray(RX)
  d_C  = CuArrays.CuArray(cat(C..., dims=2))

  # === Get unaries in the gpu ===
  d_prevcost = CuArrays.CuArray{Cfloat}(n)
  d_newcost  = CuArrays.CuArray{Cfloat}(n)

  # CUBLAS.cublasCreate_v2( CUBLAS.cublashandle )
  d_unaries   = Vector{CuArrays.CuArray{Float32}}(undef, m)
  d_codebooks = Vector{CuArrays.CuArray{Float32}}(undef, m)
  for j = 1:m
    d_codebooks[j] = CuArrays.CuArray(C[j])
    # -2 * C' * X
    d_unaries[j] = CuArrays.BLAS.gemm('T', 'N', -2.0f0, d_codebooks[j], d_RX)
    # d_unaries[j] = -2.0f0 * d_codebooks[j]' * d_RX <-- thus runs out of memory real fast

    # Add self-codebook interactions || C_{i,i} ||^2
    CudaUtilsModule.vec_add(n, (1,h), d_unaries[j].buf, CuArrays.CuArray(diag( C[j]' * C[j] )).buf, Cint(n), Cint(h))
  end

  # === Get binaries to the GPU ===
  d_binaries  = Vector{CuArrays.CuArray{Float32}}(undef, ncbi)
  d_binariest = Vector{CuArrays.CuArray{Float32}}(undef, ncbi)
  for j = 1:ncbi
    d_binaries[j]  = CuArrays.CuArray(binaries[j])
    d_binariest[j] = CuArrays.CuArray(binaries_t[j])
  end

  # Allocate space for temporary results
  bbs = Vector{Matrix{Cfloat}}(undef, m-1)

  # Initialize the previous cost
  prevcost = Rayuela.veccost(RX, B, C)
  IDX = 1:n;

  # For codebook i, we have to condition on these codebooks
  to_look      = 1:m
  to_condition = zeros(Int32, m-1, m)
  for j = 1:m
    tmp = collect(1:m)
    splice!( tmp, j )
    to_condition[:,j] = tmp
  end

  # Loop for the number of requested ILS iterations
  for i = 1:maximum(ilsiters)

    # @show i, ilsiters
    # @time begin

    to_look_r      = to_look
    to_condition_r = to_condition

    # Compute the cost of the previous assignments
    # CudaUtilsModule.veccost(n, (1, d), d_RX, d_C, CuArray( convert(Matrix{Cuchar}, (B')-1) ), d_prevcost, Cint(m), Cint(n))
    CudaUtilsModule.veccost2(n, (1, d), d_RX.buf, d_C.buf, CuArrays.CuArray( convert(Matrix{Cuchar}, collect(B').-1) ).buf, d_prevcost.buf, Cint(d), Cint(m), Cint(n))
    CUDAdrv.synchronize(ctx)
    prevcost = Array(d_prevcost)

    newB = copy(B)

    # Randomize the visit order in ICM
    if randord
      to_look_r      = randperm( m )
      to_condition_r = to_condition[:, to_look_r]
    end

    d_newB = CuArrays.CuArray(convert(Matrix{Cuchar}, newB.-1 ))

    # Perturb `npert` entries in each code
    CudaUtilsModule.perturb(n, (1,m), d_state, d_newB.buf, Cint(n), Cint(m), Cint(npert))

    newB = Array( d_newB )
    Bt   = collect(newB')
    d_Bs = CuArrays.CuArray( Bt )

    CUDAdrv.synchronize(ctx)

    # Run the number of ICM iterations requested
    for j = 1:icmiter

      kidx = 1;
      # Loop through each MRF node
      for k = to_look_r

        lidx = 1;
        for l = to_condition_r[:, kidx]
          # Determine the pairwise tables that we'll use in this conditioning
          if k < l
            binariidx = cbpair2binaryidx[ k, l ]
            bbs[lidx] = binaries[ binariidx ]
          else
            binariidx = cbpair2binaryidx[ l, k ]
            bbs[lidx] = binaries_t[ binariidx ]
          end
          lidx = lidx+1;
        end

        # Transfer pairwise tables to the GPU
        d_bbs = CuArrays.CuArray( convert(Matrix{Cfloat}, cat(bbs..., dims=2)) )
        # Sum binaries (condition) and minimize
        CudaUtilsModule.condition_icm3(
          n, (1, h),
          d_unaries[k].buf, d_bbs.buf, d_Bs.buf, Cint(k-1), Cint(m), Cint(n) )
        CUDAdrv.synchronize(ctx)

        kidx = kidx + 1;
      end # for k = to_look_r
    end # for j = 1:icmiter

    newB = Array( d_Bs )
    newB = convert(Matrix{Int16}, newB')
    newB .+= 1

    # Keep only the codes that improved
    CudaUtilsModule.veccost2(n, (1, d), d_RX.buf, d_C.buf, d_Bs.buf, d_newcost.buf, Cint(d), Cint(m), Cint(n))
    CUDAdrv.synchronize(ctx)

    newcost = Array( d_newcost )

    areequal = newcost .== prevcost
    if V; @printf(" ILS iteration %d/%d done. ", i, maximum(ilsiters)); end
    if V; @printf("%5.2f%% new codes are equal. ", 100*sum(areequal)/n ); end

    arebetter = newcost .< prevcost
    if V; @printf("%5.2f%% new codes are better.\n", 100*sum(arebetter)/n ); end

    newB[:, .~arebetter] = B[:, .~arebetter]
    B = copy( newB )
    # end # @time

    # Check if this # of iterations was requested
    if i in ilsiters

      ithidx = findall( i .== ilsiters )[1]

      # Compute and save the objective
      obj = qerror( RX, B, C )
      # @show obj
      objs[ ithidx ] = obj
      # @show size(B)
      Bs[ ithidx ] = B

    end # end if i in ilsiters
  end # end for i=1:max(ilsiters)

  # CUBLAS.cublasDestroy_v2( CUBLAS.cublashandle )
  CudaUtilsModule.finit()

  destroy!(ctx)
  # end # do devlist

  if V; @printf(" Encoding done in %.2f seconds\n", (time_ns() -start_time)/1e9); end

  return Bs, objs
end

function encode_icm_cuda(
  RX::Matrix{Float32},         # in. The data to encode
  B::Matrix{Int16},            # in. Initial list of codes
  C::Vector{Matrix{Float32}},  # in. Codebooks
  ilsiters::Vector{Int64},     # in. ILS iterations to record Bs and obj function. Its max is the total number of iteration we will run.
  icmiter::Integer,            # in. Number of ICM iterations
  npert::Integer,              # in. Number of entries to perturb
  randord::Bool,               # in. Whether to randomize the order in which nodes are visited in ILS
  nsplits::Integer=2,          # in. Number of splits of the data (for limited memory GPUs)
  V::Bool=false)

  # TODO check that splits >= 1
  if nsplits == 1
    Bs, objs =  encode_icm_cuda_single(RX, B, C, ilsiters, icmiter, npert, randord, V)
    GC.gc()
    return Bs, objs
  end

  # Split the data
  d, n = size(RX)
  splits = splitarray(1:n, nsplits)
  nr = length(ilsiters)

  Bs   = Vector{Matrix{Int16}}(undef, nr)
  objs = zeros(Float32, nr)

  # Create storage space for the codes
  for i = 1:nr; Bs[i] = Matrix{Int16}(undef, size(B)); end

  # Run encoding in the GPU for each split
  for i = 1:nsplits
    aaBs, _ = encode_icm_cuda_single(RX[:,splits[i]], B[:,splits[i]], C, ilsiters, icmiter, npert, randord, V)
    GC.gc()
    for j = 1:nr
      # Save the codes
      Bs[j][:,splits[i]] = aaBs[j]
    end
  end

  # Compute the cost again
  for j = 1:nr
    objs[j] = qerror(RX, Bs[j], C)
  end

  GC.gc()
  return Bs, objs
end

"LSQ training but some things happen in the GPU"
function train_lsq_cuda(
  X::Matrix{T},         # d-by-n matrix of data points to train on.
  m::Integer,           # number of codebooks
  h::Integer,           # number of entries per codebook
  R::Matrix{T},         # init rotation
  B::Matrix{Int16},     # init codes
  C::Vector{Matrix{T}}, # init codebooks
  niter::Integer,       # number of optimization iterations
  ilsiter::Integer,     # number of ILS iterations to use during encoding
  icmiter::Integer,     # number of iterations of ICM
  randord::Bool,        # whether to use random order
  npert::Integer,       # The number of codes to perturb
  nsplits::Integer=1,   # The number of splits for icm encoding (for limited memory GPUs)
  V::Bool=false) where T <: AbstractFloat # whether to print progress

  if V
    println()
    println("**********************************************************************************************");
    println("Training LSQ GPU with $m codebooks, $npert perturbations, $icmiter icm iterations and random order = $randord");
    println("**********************************************************************************************");
  end

  d, n = size( X )

  # Initialize C
  RX = R' * X
  C = update_codebooks(RX, B, h, V, "fastbin")

  # Apply the rotation to the codebooks
  for i = 1:m; C[i] = R * C[i]; end
  if V; @printf("%3d %e \n", -2, qerror(X, B, C)); end

  # Initialize B
  B, _ = encode_icm_cuda(X, B, C, [ilsiter], icmiter, npert, randord, nsplits, V)
  B    = B[end]
  if V; @printf("%3d %e \n", -1, qerror(X, B, C)); end

  obj = zeros(T, niter)

  for iter = 1:niter
    obj[iter] = qerror(X, B, C)
    if V; @printf("%3d %e \n", iter, obj[iter]); end

    # Update the codebooks
    C = update_codebooks(X, B, h, V, "fastbin")

    # Update the codes B
    B, _ = encode_icm_cuda(X, B, C, [ilsiter], icmiter, npert, randord, nsplits, V)
    B    = B[end]
  end

  return C, B, obj
end

function experiment_lsq_cuda(
  Xt::Matrix{T}, # d-by-n. Data to learn codebooks from
  B::Matrix{T2}, # codes
  C::Vector{Matrix{T}}, # codebooks
  R::Matrix{T}, # rotation
  Xb::Matrix{T}, # d-by-n. Base set
  Xq::Matrix{T}, # d-by-n. Queries
  gt::Vector{UInt32}, # ground truth
  m::Integer,    # number of codebooks
  h::Integer,    # number of entries per codebook
  niter::Integer=25, # Number of C-B update updates
  ilsiter::Integer=8, # number of ILS iterations
  icmiter::Integer=4, # number of ICM iterations
  randord::Bool=true, # whether to explore the codes in random order
  npert::Integer=4, # numebr of codes to perturb in each ILS iterations
  knn::Integer=1000, # Compute recall @N for this value of N
  nsplits_train::Integer=1, # Number of splits for training data so GPU does not run out of memory
  nsplits_base::Integer=1, # Number of splits for training data so GPU does not run out of memory
  V::Bool=false) where {T <: AbstractFloat, T2 <: Integer} # whether to print progress

  # Train LSQ
  d, _ = size(Xt)
  @printf("Running CUDA LSQ training... ")
  C, B, train_error = Rayuela.train_lsq_cuda(Xt, m, h, R, B, C, niter, ilsiter, icmiter, randord, npert, nsplits_train, V)
  @printf("done\n")

  norms_B, norms_C = get_norms_codebook(B, C)

  # === Encode the base set ===
  B_base = convert(Matrix{Int16}, rand(1:h, m, size(Xb,2)))
  Bs_base, _ = Rayuela.encode_icm_cuda(Xb, B_base, C, [ilsiter * 4], icmiter, npert, randord, nsplits_base, V)
  B_base = Bs_base[end]
  base_error = qerror(Xb, B_base, C)
  if V; @printf("Error in base is %e\n", base_error); end

  # Compute and quantize the database norms
  B_base_norms, db_norms_X = quantize_norms( B_base, C, norms_C )
  db_norms     = vec( norms_C[ B_base_norms ] )

  if V; print("Querying m=$m ... "); end
  @time dists, idx = linscan_lsq(B_base, Xq, C, db_norms, Matrix{Float32}(1.0I, d, d), knn)
  # @time dists, idx = linscan_lsq(B_base, Xq, C, db_norms_X, eye(Float32, d), knn)
  if V; println("done"); end

  recall = eval_recall(gt, idx, knn)
  return C, B, R, train_error, B_base, recall
end

function experiment_lsq_cuda_query_base(
  Xt::Matrix{T}, # d-by-n. Data to learn codebooks from
  B::Matrix{T2}, # codes
  C::Vector{Matrix{T}}, # codebooks
  R::Matrix{T}, # rotation
  Xq::Matrix{T}, # d-by-n. Queries
  gt::Vector{UInt32}, # ground truth
  m::Integer,    # number of codebooks
  h::Integer,    # number of entries per codebook
  niter::Integer=25, # Number of C-B update updates
  ilsiter::Integer=8, # Number of ILS iterations
  icmiter::Integer=4, # Number of ICM iterations
  randord::Bool=true, # Whether to explore the codes in random order
  npert::Integer=4, # Number of codes to perturb in each ILS iterations
  knn::Integer=1000, # Compute recall @N for this value of N
  nsplits_train::Integer=1, # Number of splits for training data so GPU does not run out of memory
  V::Bool=false) where {T <: AbstractFloat, T2 <: Integer} # whether to print progress

  # Train LSQ
  d, _ = size(Xt)
  if V; @printf("Running CUDA LSQ training... "); end
  C, B, train_error = Rayuela.train_lsq_cuda(Xt, m, h, R, B, C, niter, ilsiter, icmiter, randord, npert, nsplits_train, V)
  if V; @printf("done\n"); end

  norms_B, norms_C = get_norms_codebook(B, C)
  db_norms         = vec(norms_C[ norms_B ])

  if V; print("Querying m=$m ... "); end
  @time dists, idx = linscan_lsq(B, Xq, C, db_norms, Matrix{Float32}(1.0*I, d, d), knn)
  if V; println("done"); end

  recall = eval_recall(gt, idx, knn)
  return C, B, R, train_error, recall
end

"Runs an lsq experiment/demo"
function experiment_lsq_cuda(
  Xt::Matrix{T}, # d-by-n. Data to learn codebooks from
  Xb::Matrix{T}, # d-by-n. Base set
  Xq::Matrix{T}, # d-by-n. Queries
  gt::Vector{UInt32}, # ground truth
  m::Integer,    # number of codebooks
  h::Integer,    # number of entries per codebook
  niter::Integer=25, # Number of C-B update updates
  ilsiter::Integer= 8, # number of ILS iterations
  icmiter::Integer=4, # number of ICM iterations
  randord::Bool=true, # whether to explore the codes in random order
  npert::Integer=4, # numebr of codes to perturb in each ILS iterations
  knn::Integer=1000,
  nsplits_train::Integer=1, # Number of splits for training data so GPU does not run out of memory
  nsplits_base::Integer=1, # Number of splits for training data so GPU does not run out of memory
  V::Bool=false) where T <: AbstractFloat # whether to print progress

  # OPQ initialization
  @printf("Running OPQ initialization... ")
  C, B, R, _ = train_opq(Xt, m, h, niter, "natural", V)
  @printf("done\n")

  # ChainQ (second initialization)
  # @printf("Running ChainQ initialization... ")
  # C, B, R, train_error = train_chainq(Xt, m, h, R, B, C, niter, V)
  # @printf("done\n")

  # Actual experiment
  experiment_lsq_cuda(Xt, B, C, R, Xb, Xq, gt, m, h, niter, ilsiter, icmiter, randord, npert, knn, nsplits_train, nsplits_base, V)
end

function experiment_lsq_cuda_query_base(
  Xt::Matrix{T}, # d-by-n. Data to learn codebooks from
  Xq::Matrix{T}, # d-by-n. Queries
  gt::Vector{UInt32}, # ground truth
  m::Integer,    # number of codebooks
  h::Integer,    # number of entries per codebook
  niter::Integer=25, # Number of C-B update updates
  ilsiter::Integer= 8, # number of ILS iterations
  icmiter::Integer= 4, # number of ICM iterations
  randord::Bool=true, # whether to explore the codes in random order
  npert::Integer= 4, # numebr of codes to perturb in each ILS iterations
  init::String="natural",
  niter_opq::Integer=25,
  noter_chainq::Integer=25,
  knn::Integer=1000,
  nsplits_train::Integer=1, # Number of splits for training data so GPU does not run out of memory
  V::Bool=false) where T <: AbstractFloat # whether to print progress

  # OPQ initialization
  if V; @printf("Running OPQ initialization... "); end
  C, B, R, opq_error = train_opq(Xt, m, h, niter_opq, init, V)
  if V; @printf("done\n"); end

  # ChainQ (second initialization)
  if V; @printf("Running ChainQ initialization... "); end
  C, B, R, chainq_error = train_chainq(Xt, m, h, R, B, C, niter_chainq, V)
  if V; @printf("done\n"); end

  # Actual experiment
  experiment_lsq_cuda_query_base(Xt, B, C, R, Xq, gt, m, h, niter, knn, nsplits_train, V), opq_error
end


text:
##############################################################################################################################################
# Implements the Whole History Rating algorithm from Remi Coulom as outlined in the paper: https://www.remi-coulom.fr/WHR/WHR.pdf
#
# The "predict" function is mine, based on my understanding of the algorithm. It gives good results, but I have not checked it with the paper
# and it does not match algorithms given elsewhere (e.g., in the Ruby implementation here: https://github.com/goshrine/whole_history_rating)
#
# This algorithm is also several orders of magnitude faster than the open-source versions (Ruby, Python).
##############################################################################################################################################

import Distributions
import QuadGK
using LinearAlgebra
using OrderedCollections

struct WHR
	playerdayratings::Dict
	playerdaygames::Dict
	default_rating::Float32
	default_score::Float32
	w2::Float32
end

function WHR(;playerdayratings = Dict{Int64, OrderedDict{Int64, Float32}}(), playerdaygames = Dict{Int64, OrderedDict{Int64, Array{Tuple{Int64, Int64, Int64}}}}(), default_rating = 0.0, default_score = 1.0, w2 = 0.000424)
	return WHR(playerdayratings, playerdaygames, default_rating, default_score, w2)
end

#Add games to an WHR object
function add_games!(whr::WHR, original_games::DataFrame; dummy_games::Bool = true, verbose = false)

    #Duplicate whole of new games rating dataframe with reversed results
	games = dupe_for_rating(original_games)
	P1, P2, P1_wins, P2_wins, Day = games.P1, games.P2, games.P1_wins, games.P2_wins, games.Period

	games_added = 0
	players_added = 0

	for row in 1:length(P1)
		#If this row introduces a new player1, create their (empty) entry in the playerdaygames dictionary, create a new playerday and add games
		if !haskey(whr.playerdaygames, P1[row])
			add_player!(whr::WHR, P1[row])
			players_added += 1
		end

		#If this row introduces a new gameday for a player, add that.
		if !haskey(whr.playerdaygames[P1[row]], Day[row])
			add_gameday!(whr::WHR, P1[row], Day[row])
		end

		#Now, add the game in question
		add_each_game!(whr::WHR, P1[row], Day[row], P2[row], P1_wins[row], P2_wins[row])
		games_added += 1

	end
	verbose && println("Added ", games_added ÷ 2, " new games out of ", length(P1) ÷ 2)
	verbose && println("Added ", players_added, " new players.")
end

function add_player!(whr::WHR, player::Int64)
	#Create an empty Dictionary for a new player: day => gamesarray and a new ratings Dictionary day => rating
	whr.playerdaygames[player] = OrderedDict{Int64, Vector{Tuple{Int64, Int64, Int64}}}()
	whr.playerdayratings[player] = OrderedDict{Int64, Float32}()
end

function add_gameday!(whr::WHR, player::Int64, day::Int64)
	#Create an empty Dict for a new gameday and a new rating for that gameday
	#Note that one could come up with a better starting point, e.g., the previous rating

	#Check whether the new gameday is messing up the order of the existing days, if so
	sort_toggle = !isempty(keys(whr.playerdayratings[player])) && day < maximum(keys(whr.playerdayratings[player]))

	whr.playerdaygames[player][day] = Vector{Tuple{Int64, Int64, Int64}}() #Initialise new day, with no games yet
	whr.playerdayratings[player][day] = whr.default_rating #Set a default rating on that gameday

	if sort_toggle
		sort!(whr.playerdaygames[player])
		sort!(whr.playerdayratings[player])
	end
end

function add_each_game!(whr::WHR, P1::Int64, day::Int64, P2::Int64, P1_wins::Int64, P2_wins::Int64)
	#Create a new opponent (P2) and add the scores of the game as a Tuple
	whr.playerdaygames[P1][day] = vcat(whr.playerdaygames[P1][day], (P2, P1_wins, P2_wins))
end

function iterate!(whr::WHR, iterations; exclude_non_ford::Bool = false, delta = 0.001f0, verbose = false, show_ll = false, ll_every = 1, new_invert_method = true)

	#Check and get rid of any players not Ford-connected

	if exclude_non_ford
		non_ford = check_ford(whr)
		connected_players = filter(x -> ∉(first(x), non_ford) , bt.playergames)
		println(length(non_ford), " players excluded due to not being connected. ", length(connected_players), " players being rated.")
	else
		connected_players = whr.playerdaygames
	end

	#Iterate over players selected
	for i = 1:iterations
		for (player, _) in connected_players
			if length(whr.playerdaygames[player]) == 1
				#Need to use a separate update method if we only have one day of games
				update_rating_1dim!(whr, player)
			else
				#If we have several days of games use general update method
				update_rating_ndim!(whr, player, delta = delta, new_invert_method = new_invert_method)
			end
		end
		verbose && println("Iteration on whole WHR: ", i)

		if show_ll
			if i % ll_every  == 0
				println("Log likelihood = ", log_likelihood(whr))
			end
		end
	end
end

function iterate!(whr::WHR, players::Array{Int64}, iterations::Int64; exclude_non_ford::Bool = false, delta::Float32 = 0.001f0)
	#Iterate over players selected
	for i = 1:iterations
		for player in players
			if length(whr.playerdaygames[player]) == 1
				#Need to use a separate update method if we only have one day of games
				update_rating_1dim!(whr, player)
			else
				#If we have several days of games use general update method
				update_rating_ndim!(whr, player, delta = delta)
			end
		end
	end
	#println(iterations, " iterations on added players (", players,") completed.")
end

function iterate_to_convergence!(whr::WHR; step::Int64 = 1, delta::Float32 = 0.001f0, tolerance::Float32 = 0.00001f0)
	ll_percentage = 0.1
	ll = log_likelihood(whr)
	ll_old = ll
	n = 0
	while ll_percentage > tolerance
		n += step
		iterate!(whr, step, delta = delta)
		ll = log_likelihood(whr)
		diff = abs(ll_old) - abs(ll)
		ll_percentage = - (ll - ll_old) / ll_old
		println("iteration: ", n, ". Log-likelihood: ", ll, ". Change: ", ll_percentage)
		ll_old = ll
	end
end

#Use Newton-Raphson method
function update_rating_ndim!(whr::WHR, player::Int64; delta::Float32 = 0.001f0, new_invert_method = true)

	#Get all the players' days played and initial ratings
	#Note: using "get" might allow to set initial rating at average of others - big speed-up

	playerdays = whr.playerdayratings[player].keys #All the days the player played
	playerratings = whr.playerdayratings[player].vals

	#Get the sigma2 vector (one-over it in this case)
	s2inv = sigma2inv(playerdays, whr.w2)

	#Obtain the log likelihood derivative and second derivative at one sweep
	lld, ll2d = llderivatives(whr, player, playerdays)

	if new_invert_method
		change_in_rating = invHG(lld, ll2d, s2inv, playerratings, delta = delta)
		replace!(change_in_rating, NaN => 0)
		whr.playerdayratings[player].vals -= invHG(lld, ll2d, s2inv, playerratings, delta = delta)
	else
		h = hessian(ll2d, s2inv, delta = delta)
		g = gradient(lld, s2inv, playerratings)
		newratings = playerratings - inv(h) * g

		for (i, day) in enumerate(playerdays)
			whr.playerdayratings[player][day] = newratings[i]
		end
	end
end

#Update using Newton-Raphson method for players with only one game
function update_rating_1dim!(whr::WHR, player::Int64)
	playerdays = collect(keys(whr.playerdaygames[player]))
	lld, ll2d = llderivatives(whr, player, playerdays)
	@inbounds dr = lld[1] / ll2d[1]
	whr.playerdayratings[player][playerdays[1]] = whr.playerdayratings[player][playerdays[1]] - dr
end

#Broadcast the element-wise calculation of derivative elements over the playerdays
function llderivatives(whr, player, playerdays)
	days_played = length(playerdays)
	lld = zeros(Float32, days_played)
	ll2d = zeros(Float32, days_played)

	#Add the term for 1 dummy win and 1 dummy loss against opponent of rating 0 on first day
	@fastmath r1 = exp(whr.playerdayratings[player][first(playerdays)])
	@inbounds lld[1] = (1 - r1) / (1 + r1)
	@inbounds ll2d[1] = -2 * r1 / (1 + r1)^2

	#Loop round playerdays and add in loglikelihood derivates to vectors
	@inbounds for (daynum, day) in enumerate(playerdays)
		llde, ll2de = llderivativeselements(whr, player, day)
		@inbounds lld[daynum] += llde
		@inbounds ll2d[daynum] += ll2de
	end
	return lld, ll2d
end

#For each playerday, calculate the llderivates elements (first and second derivatives)
function llderivativeselements(whr, player::Int64, day::Int64)
	lld_tally = 0.0
	ll2d_tally = 0.0
	win_tally = 0.0
	@inbounds @fastmath a = exp(get(whr.playerdayratings[player], day, 0.0))

	@inbounds for (opponent, resultsw, resultsl) in whr.playerdaygames[player][day]
		@inbounds @fastmath b = exp(get(whr.playerdayratings[opponent], day, 0.0))
		@fastmath sumab = 1.0 / (a + b)
		@fastmath lld_tally_add = (resultsw + resultsl) * sumab

		@fastmath win_tally += resultsw
		@fastmath lld_tally += lld_tally_add
		@fastmath ll2d_tally += lld_tally_add * b * sumab
	end

	return (win_tally - (a * lld_tally), -a * ll2d_tally)
end

"""
function llelements(whr, player::Int64, day::Int64)
	ll_tally = 0.0
	own_rating = get(whr.playerdayratings[player], day, 0.0)

	for (opponent, result) in whr.playerdaygames[player][day]
		if opponent == 0
			opponent_rating = 1.0
		else
			opponent_rating = get(whr.playerdayratings[opponent], day, 0.0)
		end
		@fastmath ll_tally += result * own_rating + ((1 - result) * opponent_rating) - log(exp(own_rating) + exp(opponent_rating))
	end

	return ll_tally
end

function log_likelihood(whr::WHR, player)
	#Calculate current log-likelihood for a player history
	player_tally = 0.0
	player_tally = sum(llelements(whr, player, day) for day in keys(whr.playerdaygames[player]))
	return player_tally
end
"""

function log_likelihood(whr::WHR)
	#Calculate current total log-likelihood for all players
	full_tally = sum(log_likelihood(whr::WHR, player) for player in keys(whr.playerdayratings))
	return full_tally
end

function sigma2inv(playerdays, w2::Float32)
	l = length(playerdays) - 1
	s2inv = Array{Float32}(undef, l)
	for i in 1:l
		@inbounds s2inv[i] = 1 / (w2 * (playerdays[i + 1] - playerdays[i]))
	end
	return s2inv
end

function hessian(ll2d, s2inv; delta::Float32 = 0.001f0)
	#Construct a symmetric tridiagonal hessian matrix
	l = length(ll2d)

	#Most of the work goes into the prior (principal diagonal)
	pdiag = Array{Float32}(undef, l)

	@inbounds pdiag[1] = ll2d[1] - s2inv[1] - delta
	for i in 2:l-1
		@inbounds pdiag[i] = ll2d[i] - s2inv[i - 1] - s2inv[i] - delta
	end
	@inbounds pdiag[l] = ll2d[l] - s2inv[l - 1] - delta
	return SymTridiagonal(pdiag, s2inv)
end

function gradient(lld, s2inv, playerratings)
	#Construct gradient vector
	l = length(lld)
	grad = Array{Float32}(undef, l)

	@inbounds grad[1] = lld[1] + (playerratings[2] - playerratings[1]) * s2inv[1]
	for i in 2:(l-1)
		@inbounds grad[i] = lld[i] + (playerratings[i + 1] - playerratings[i]) * s2inv[i] - (playerratings[i] - playerratings[i - 1]) * s2inv[i - 1]
	end
	@inbounds grad[l] = lld[l] - (playerratings[l] - playerratings[l - 1]) * s2inv[l - 1]

	return grad
end

#Create a fast inverse of the hessian multiplied by the gradient using the algorithm from the WHR paper
function invHG(lld, ll2d, s2inv, playerratings; delta = delta)

	#Prepare a, b, d vectors (using vectors, may be faster using SVectors)
	d, a = d_a_vector(ll2d, s2inv, delta)
	b = s2inv

	#Prep y vector (again, could use SVectors)
	g = gradient(lld, s2inv, playerratings)
	y = y_vector(g, a)

	#Prepare x vector
	return x_vector(y, b, d)
end

function d_a_vector(ll2d, s2inv, delta)
	#Construct the d and a vectors from WHR paper (pp.10-11)
	l = length(ll2d)
	d = Array{Float32}(undef, l)
	a = Array{Float32}(undef, l - 1)

	@inbounds d[1] = ll2d[1] - s2inv[1] - delta
	for i in 2:l-1
		@inbounds a[i - 1] = s2inv[i - 1] / d[i - 1]
		@inbounds d[i] = ll2d[i] - s2inv[i - 1] - s2inv[i] - delta - a[i - 1] * s2inv[i - 1]
	end
	@inbounds a[l - 1] = s2inv[l - 1] / d[l - 1]
	@inbounds d[l] = ll2d[l] - s2inv[l - 1] - delta - a[l - 1] * s2inv[l - 1]
	return d, a
end

function y_vector(g, a)
	#Construct the y vector from WHR paper (p.11)
	l = length(g)
	y = Array{Float32}(undef, l)

	@inbounds y[1] = g[1]
	for i in 2:l
		@inbounds y[i] = g[i] - a[i - 1] * y[i - 1]
	end
	return y
end

function x_vector(y, b, d)
	#Construct the x vector from WHR paper (p.11)
	l = length(y)
	x = Array{Float32}(undef, l)

	#Note this one goes in reverse
	@inbounds x[l] = y[l] / d[l]
	for i in reverse(1:l - 1)
		@inbounds x[i] = (y[i] - b[i] * x[i + 1]) / d[i]
	end
	return x
end

function check_ford(whr::WHR)
	#Check for players who have either never won or never lost
	unconnected = []
	for (player, games) in bt.playergames
		if length(filter(y -> y[3] == 0.0, games)) == 0
			unconnected = vcat(unconnected, player)
		end
		if length(filter(y -> y[3] == 1.0, games)) == 0
			unconnected = vcat(unconnected, player)
		end
	end
	return unconnected
end

function covariance(whr::WHR, player::Int64)

	#Following Appendix B in the WHR paper -- there's almost certainly a better way to do this in Julia
	playerdays = collect(keys(whr.playerdaygames[player]))
	playerratings = [whr.playerdayratings[player][day] for day in playerdays]

	s2inv = sigma2inv(playerdays, whr.w2)
	lld, ll2d = llderivatives(whr, player, playerdays)
	h = hessian(ll2d, s2inv)
	g = gradient(lld, s2inv, playerratings)

	n = length(playerdays)

	a = zeros(Float32, n)
	d = zeros(Float32, n)
	b = zeros(Float32, n)

	d[1] = h[1, 1]
	if size(h)[1] > 2
		b[1] = h[1, 2]
	else
		b[1] = 0.0
	end

	for i in range(2, stop = n)
		a[i] = h[i, i - 1] / d[i - 1]
		d[i] = h[i, i] - a[i] * b[i - 1]
		if i < n
			b[i] = h[i, i + 1]
		end
	end

	dp = zeros(Float32, n)
	dp[n] = h[n, n]
	bp = zeros(Float32, n)
	bp[n] = h[n, n - 1]

	ap = zeros(Float32, n)
	for i in range(n - 1, step = -1, stop = 1)
		ap[i] = h[i, i + 1] / dp[i + 1]
		dp[i] = h[i, i] - ap[i] * bp[i + 1]
		if i > 1
			bp[i] = h[i, i-1]
		end
	end

	v = zeros(Float32, n)
	for i in range(1, stop = n - 1)
		v[i] = dp[i + 1]/(b[i] * bp[i + 1] - d[i] * dp[i + 1])
	end
	v[n] = -1 / d[n]

	ev = - a[2:n] .* v[1:n - 1]

	return Bidiagonal(v, ev, :U)
end

function uncertainty(whr::WHR, player::Int64)
	#Create a dictionary of uncertainties around each rating for a player
	if haskey(whr.playerdaygames, player)
		playerdays = collect(keys(get(whr.playerdaygames, player, [0])))
		ndays = length(whr.playerdaygames[player])
	else
		playerdays = [0]
		ndays = 0
	end

	if ndays > 1
		c = covariance(whr, player)
		u = diag(c)
	else
		u = [1.0]
	end

	d = Dict(playerdays .=> u)
	return d
end

function ratings(whr::WHR, player::Int64)
	#Return the two ratings dictionaries - of ratings, and of the uncertainties around them
	return get(whr.playerdayratings, player, whr.default_rating), uncertainty(whr, player)
end

function rating(whr::WHR, player::Int64; rating_day::Int64)
	#Return the player's rating on a day, and the uncertainties it has, projected forward from last rating period

	if haskey(whr.playerdayratings, player)
		rating_dict, var_dict = ratings(whr, player)
	else
		return whr.default_rating, 1.0
	end

	if haskey(rating_dict, rating_day)
		rating = rating_dict[rating_day]
		var = var_dict[rating_day]
	else
		b = collect(keys(rating_dict))
		if isempty(b[b .< rating_day])
			rating = whr.default_rating
			var = 1.0
		else
			last_day_rated = maximum(b[map(x -> x < rating_day, b)])
			rating = rating_dict[last_day_rated]
			var = get(var_dict, last_day_rated, 0)  + abs(last_day_rated - rating_day) * abs(whr.w2)
		end
	end
	return rating, var
end

function rating(whr::WHR, P1::Int64, P2::Int64; rating_day::Int64 = 0)
	#Return a pair of players' ratings as a pair of tuples
	return rating(whr, P1; rating_day = rating_day), rating(whr, P2; rating_day = rating_day)
end

function predict(whr::WHR, P1::Int, P2::Int; rating_day = missing, raw = false, kwargs...)
	#Predict with logitnormal distribution, pulling forward the variance to the rating day (default, present day)

	#If no rating_day provided, use the last day on which either player was rated
	if !haskey(whr.playerdaygames, P1)
		r1 = whr.default_rating
		var1 = 1.0
	elseif length(whr.playerdaygames[P1]) == 0
		r1 = whr.default_rating
		var1 = 1.0
	else
		lastdayP1 = maximum(keys(whr.playerdaygames[P1]))
		r1, var1_lastday = rating(whr, P1, rating_day = lastdayP1)
		var1 = var1_lastday + (abs(rating_day - lastdayP1) * whr.w2)
	end

	if !haskey(whr.playerdaygames, P2)
		r2 = whr.default_rating
		var2 = 1.0
	elseif length(whr.playerdaygames[P2]) == 0
		r2 = 0.0
		var2 = 1.0
	else
		lastdayP2 = maximum(keys(whr.playerdaygames[P2]))
		r2, var2_lastday = rating(whr, P2, rating_day = lastdayP2)
		var2 = var2_lastday + (abs(rating_day - lastdayP2) * whr.w2)
	end
	#Difference between two normally distributed RVs is normally distributed with mean of the difference of the two means, and variance as sum of the variances.
	#For the probability, we are looking for the logistic function of this normal distribution: given by the mean of the logit-normal (I hope).
	#Using approximation here with series of 10 terms

	return mean_logitnormal_approx(r1 - r2, var1 + var2, 10)
end

function mean_logitnormal(mu, var)
	#Returns the mean of the logitnormal. No analytic way, so using Gaussian quadrature
	#Can be made much faster with approximations, but not a performance-critical piece (I hope).
	dist = Distributions.LogitNormal(mu, sqrt(var))
	prob_win = QuadGK.quadgk(x -> x * Distributions.pdf(dist, x), 0, 1)[1]
    return prob_win
end

function mean_logitnormal_approx(mu, var, K)
	#Approximate function for logitnormals - about 1000 times faster and also more stable than numerical quadrature

	cumulative = 0.0
	d = Distributions.Normal(mu, var)
    for i in 1:K-1
        cumulative += std_logistic(Distributions.quantile(d, i/K))
    end
    return cumulative / (K-1)
end

function std_logistic(x)
    return 1 / (1 + exp(-x))
end


text:
using BinaryBuilder

name = "Tokei"
version = v"12.1.2"

sources = [
    ArchiveSource("https://github.com/XAMPPRocky/tokei/archive/refs/tags/v$(version).tar.gz",
                  "81ef14ab8eaa70a68249a299f26f26eba22f342fb8e22fca463b08080f436e50"),
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/tokei*/
cargo build --release
install -D -m 755 "target/${rust_target}/release/tokei${exeext}" "${bindir}/tokei${exeext}"
install_license LICENCE-MIT LICENCE-APACHE
"""

platforms = supported_platforms()
# Our Rust toolchain for i686 Windows is unusable
filter!(p -> !Sys.iswindows(p) || arch(p) != "i686", platforms)

# The products that we will ensure are always built
products = [
    ExecutableProduct("tokei", :tokei),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; compilers=[:c, :rust])

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
function sbar(arg1)
    println("sbar: $arg1")
end


text:
module SolveEig

using LinearAlgebra
using ..Interp:getD
using ..Interp:getCpts
using ..Interp:getCpts1
using ..Interp:baryW1
using ..Interp:barymat
using ..Interp:forward
using ..Interp:rev
using ..EdgeDetect:find
using ..Op:setup

#=
function many_ham(V, N; d2 = -0.5, d1 = 0.0, a = -1.0, b = 1.0, Nsearch=103)

    edges, contin = find(V, Nsearch, a=a, b=b)
#    edges = [0.0 0.0]

    a1=a
    HL = []

    for i in 1:size(edges)[1]
        b1 = edges[i,1]
        h = setup_ham(V, N, a=a1, b = b1, d2=d2, d1=d1)
        push!(HL, h)
        a1=edges[i,2]
    end
    h = setup_ham(V, N, a=a1, b = b, d2=d2, d1=d1)
    push!(HL, h)

    return HL
#    for i in range(edges
#    H1 = setup_ham(V, a=a

end

function total_ham(V, N; d2 = -0.5, d1 = 0.0, a = -1.0, b = 1.0, Nsearch=103)

    DL = []
    VL = []
    Vpts, Dab = setup_ham2(V, N; a = a, b = -1.0 - 1e-30)
    push!(DL, Dab)
    push!(VL, Vpts)
    Vpts, Dab = setup_ham2(V, N; a = -1.0+1e-30, b = 1.0-1e-30)
    push!(DL, Dab)
    push!(VL, Vpts)
    Vpts, Dab = setup_ham2(V, N; a =  1.0+1e-30, b = b)
    push!(DL, Dab)
    push!(VL, Vpts)

#    println("SIZE HL ", size(HL))
    n = 0
    for h in DL
        println( size(h)[1])
        n += size(h)[1]
    end

    nh = n - length(DL) + 1
    #nh = n
    H = zeros(nh, nh)

    Dbig = zeros(nh, nh)
    Vbig = zeros(nh, nh)
    n1 = 0
    n2 = 0

    num_secs = length(DL)

    for (c,(d,v)) in enumerate(zip(DL, VL))
        n2 += size(d)[1]
        Dbig[(n1+1):n2,(n1+1):n2] += d
        Vbig[(n1+1):n2,(n1+1):n2] += v
#        println(n1+1, " ", n2)
        n2 = n2-1
        n1 = n2

    end

    n1 = 0
    n2 = 0
    for (c,(d,v)) in enumerate(zip(DL, VL))
#        println("c $c")
#        println(v)

        n2 += size(d)[1]
        if c != 1
            Vbig[n1+1,n1+1] += v[1,1] / 2.0
        end
        if c != num_secs
            Vbig[n2,n2] = v[end,end] / 2.0
        end
        n2 = n2-1
        n1 = n2
    end

    D2 = Dbig*Dbig

    H = d2 * D2 + Vbig
#    H= Vbig

#=    n1 = 0
    n2 = 0
    for h in HL
        n2 += size(h)[1]
        H[(n1+1):n2,(n1+1):n2] += h
        println(n1+1, " ", n2)
        n2 = n2
        n1 = n2

    end
=#

#=    n1 = 0
    n2 = 0
    for h in HL
        n2 += size(h)[1]
        H[(n1+1):n2,(n1+1):n2] += h
        println(n1+1, " ", n2)
        n2 = n2-1
        n1 = n2

    end
=#

    return H
end

function total_ham_piecewise(V, N, ranges; d2 = -0.5, d1 = 0.0)

    DL = []
    VL = []
    for c in 1:length(ranges)-1
        a = ranges[c]
        b = ranges[c+1]
        Vpts, Dab = setup_ham2(V, N; a = a+1e-15, b = b-1e-15)
        push!(DL, Dab)
        push!(VL, Vpts)
    end

    num_parts = length(DL)
    hsize = 0
    for (c,(d,v)) in enumerate(zip(DL,VL))
        hsize += size(d)[1]
    end

    H = zeros(hsize, hsize)
    S = zeros(hsize, hsize)

    B = hcat(1.0, zeros(1,hsize-1))  #boundary conditions
    B = [B; hcat(zeros(1,hsize-1), 1.0)]

    n1r = 0
    n2r = 0
    n1c = 0
    n2c = 0

    for (c,(d,v)) in enumerate(zip(DL,VL))
        D2 = d2*(d*d)

        N = size(d)[1] - 1

        x = getCpts(N)
        w = baryW1(N)
        y = getCpts1(N-2)
        P = barymat(y, x, w)

        H0 = v + D2

        n2r += size(v)[1] - 2
        n2c += size(v)[1]

 #       println("a ", size(H[(n1r+1):n2r,(n1c+1):n2c]))
 #       println("p ", size(P))
 #       println("H0 ", size(H0))
        H[(n1r+1):n2r,(n1c+1):n2c] += P*H0
        S[(n1r+1):n2r,(n1c+1):n2c] += P*I(N+1)

        if c < num_parts #add continuity of psi, psi'
#            println("B ", size(B))
#            println("hcat ", size(hcat(zeros(1,n2c-1), 1, -1, zeros(1, hsize-(n2c+1)))))
            B = [B; hcat(zeros(1,n2c-1), 1, -1, zeros( 1,hsize-(n2c+1)))] #make interface psi continuous
            st = size(DL[2])[2]
            B = [B; vcat(zeros(n1c,1), DL[c][end,:], -DL[c+1][1,:], zeros(hsize-st-n2c ,1))'] #make deriviates contin
            #            B = [B; hcat(zeros(1,n1c), d[end,:], -DL[c+1][1,:], zeros(1, hsize - (n2c+1))) ]
        end

        n2r = n2r
        n1r = n2r
        n2c = n2c
        n1c = n2c

    end

    sb = size(B)[1]
    H[end-sb+1:end, :] = B
    S[end-sb+1:end, :] .= 0.0

#    println("size Vpts ", size(Vpts))
#    println("size Dab ", size(Dab))
#    println("size x ", size(x))
#    println("size y ", size(y))
#    println("size P ", size(P))

#    println("size PH ", size(PH))

#    B = [ hcat(1.0, zeros(1,N));  hcat(zeros(1,N), 1.0)]

#    println("size B ", size(B))

#    Hbig = [PH; B]
#    s = size(H0)[1]
#    S = collect(I(size(H0)[1]))
#    PS = P * S
#    Sbig = [PS; zeros(2, s)]

    return H, S

    #    Vpts, Dab = setup_ham2(V, N; a = a, b = -1.0 - 1e-30)
#    push!(DL, Dab)
#    push!(VL, Vpts)
#    Vpts, Dab = setup_ham2(V, N; a = -1.0+1e-30, b = 1.0-1e-30)
#    push!(DL, Dab)
#    push!(VL, Vpts)

end
=#

function setup_ham(V, N; d2 = -0.5, d1 = 0.0, a = -1.0, b = 1.0)

    D = getD(N);
    pts = getCpts(N);

    Dab = D/ ((b-a)/2.0)
    Dab2 = Dab*Dab

    pts_f = rev(pts,a,b)

    Vpts = diagm(V.(pts_f))
    H = d2*Dab2 + d1*Dab + Vpts

    return H

end

#=
function setup_ham2(V, N; d2 = -0.5, d1 = 0.0, a = -1.0, b = 1.0)

    D = getD(N);
    pts = getCpts(N);

    Dab = D/ ((b-a)/2.0)
#    Dab2 = Dab*Dab

    pts_f = rev(pts,a,b)

    Vpts = diagm(V.(pts_f))
#    H = d2*Dab2 + d1*Dab + Vpts

    return Vpts, Dab

end

=#

function finite_diff_ham(V, N; d2 = -0.5, d1 = 0.0, a = -1.0, b = 1.0)

    h = (b-a)/N
    grid = a:h:(b+1e-15)

    Vpts = diagm(V.(grid))

    D2 = zeros(length(grid), length(grid))
    for i = 2:(length(grid)-1)
        D2[i,i] = 2.0/h^2*d2
        D2[i,i+1] = -1.0/h^2*d2
        D2[i,i-1] = -1.0/h^2*d2
    end
    D2[1,1] = 1.0/h^2*d2
    D2[1,2] = -1.0/h^2*d2
    D2[end,end] = 1.0/h^2*d2
    D2[end, end-1] = -1.0/h^2*d2

    return Vpts - D2

end

function solve(;N = 30, bc1 =  [:a,0,0.0], bc2 = [:b,0,0.0], ranges=zeros(0,2), contin = Bool[], A0=0.0, A1=0.0, A2=-0.5, B=1.0, a = -1.0, b = 1.0, V = missing, dosplit=true)

    if !ismissing(V)
        A0 = V
    end

    H, S, rhs =  setup(N=N, a=a,b=b,bc1=bc1, bc2=bc2, ranges=ranges, contin=contin, A0=A0, A1=A1,A2=A2,B=B,dosplit=dosplit)

    vals, vects = eigen( H, S)

    return vals, vects

    #vals, vects = eigen( H[2:end-1,2:end-1]);
    #vals, vects = eigen( H);

#    vects_n = normalize_vects(vects)

    #return vals, vects

end

end #end module


text:
"""
**Template Likelihood**

Template file for likelihood creation

```julia
TemplateLikelihood()
```
See all functions you need to implement
---

"""
struct TemplateLikelihood{T<:Real,A<:AbstractVector{T}} <: AbstractLikelihood{T}
    ## Additional parameters can be added
    θ::A
    function TemplateLikelihood{T}() where {T<:Real}
        return new{T,Vector{T}}()
    end
    function TemplateLikelihood{T}(θ::A) where {T<:Real,A<:AbstractVector{T}}
        return new{T,A}(θ)
    end
end

function TemplateLikelihood()
    return TemplateLikelihood{Float64}()
end

function implemented(
    ::TemplateLikelihood, ::Union{<:AnalyticVI,<:QuadratureVI,<:GibbsSampling}
)
    return true
end

function pdf(l::TemplateLikelihood, y::Real, f::Real)
end

function Base.show(io::IO, model::TemplateLikelihood)
    return print(io, "Template Likelihood")
end

function compute_proba(
    l::TemplateLikelihood{T}, μ::AbstractVector, σ²::AbstractVector
) where {T<:Real}
    N = length(μ)
    pred = zeros(T, N)
    sig_pred = zeros(T, N)
    return pred, sig_pred
end

### Local Updates Section ###

function local_updates!(
    l::TemplateLikelihood{T}, y::AbstractVector, μ::AbstractVector, diagΣ::AbstractVector
)
end

function sample_local!(
    l::TemplateLikelihood{T}, y::AbstractVector, f::AbstractVector
) where {T}
    return nothing
end

### Natural Gradient Section ###

∇E_μ(l::TemplateLikelihood, ::AOptimizer, y::AbstractVector) = (nothing,)
∇E_Σ(l::TemplateLikelihood, ::AOptimizer, y::AbstractVector) = (nothing,)


text:
using MPCCsolve

MPCCsolve_success=true
println("On teste le module MPCCsolve")

#Bilan : donner une sortie bilan
if MPCCsolve_success==true
 println("MPCCsolve.jl passes the test !")
else
 println("MPCCsolve.jl contains some error")
end


text:
# Tests for Hölder DEA Models
@testset "HolderDEAModel" begin

    # ------------------
    # Input oriented
    # ------------------
    X = [2 2; 1 4; 4 1; 4 3; 5 5; 6 1; 2 5; 1.6 8]
    Y = [1; 1; 1; 1; 1; 1; 1; 1]

    # Hölder l = 1
    holderio1 = deaholder(X, Y, l = 1, orient = :Input, rts = :VRS)

    @test typeof(holderio1) == HolderL1DEAModel

    @test nobs(holderio1) == 8
    @test ninputs(holderio1) == 2
    @test noutputs(holderio1) == 1
    @test efficiency(holderio1) ≈ [0; 0; 0; 2.0; 4.0; 0.0; 1.0; 0.6]
    @test efficiency(holderio1, :min) == [1; 1; 1; 2; 1; 2; 1; 1]
    @test convert(Matrix, peers(holderio1)) ≈
            [ 1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
            0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0
            0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0
            0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0
            0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0
            0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0
            0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0
            0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0]
    @test slacks(holderio1, :X) ≈ [0.0 0.0; 0.0 0.0; 0.0 0.0; 0.0 0.0; 0.0 1.0; 2.0 0.0; 0.0 1.0; 0.0 4.0]
    @test slacks(holderio1, :Y) ≈ zeros(8,1)

    @test efficiency(deaholder(targets(holderio1, :X), targets(holderio1, :Y), l = 1, orient = :Input, rts = :VRS)) ≈ zeros(8,1)

    # Test no slacks
    holderio1noslacks = deaholder(X, Y, l = 1, orient = :Input, rts = :VRS, slack = false)
    @test efficiency(holderio1noslacks) == efficiency(holderio1)
    @test isempty(slacks(holderio1noslacks, :X)) == 1
    @test isempty(slacks(holderio1noslacks, :Y)) == 1

    # Weighted (weakly)
    holderio1weight = deaholder(X, Y, l = 1, orient = :Input, rts = :VRS, weight = true)
    @test efficiency(holderio1weight) ≈ [0; 0; 0; 0.625; 0.8; 0.0; 0.5; 0.375]

    # Hölder l = 2
    logs, value = Test.collect_test_logs() do
        holderio2 = deaholder(X, Y, l = 2, orient = :Input, rts = :VRS, optimizer = DEAOptimizer(:NLP))
    end

    @test typeof(holderio2) == HolderL2DEAModel

    @test nobs(holderio2) == 8
    @test ninputs(holderio2) == 2
    @test noutputs(holderio2) == 1

    if (Base.find_package("Gurobi") !== nothing)
        using Gurobi
        @info ("Testing Hölder L2 with Gurobi")

        holderio2 = deaholder(X, Y, l = 2, orient = :Input, rts = :VRS, optimizer = DEAOptimizer(Gurobi.Optimizer))
        @test efficiency(holderio2) ≈ [0.0; 0.0; 0.0; 1.788854; 4.0; 0.0; 1.0; 0.6] atol = 1e-5
        @test slacks(holderio2, :X) ≈ [0.0 0.0; 0.0 0.0; 0.0 0.0; 0.0 0.0; 1.0 0.0; 0.0 0.0; 0.0 1.0; 0.0 4.0]
        @test slacks(holderio2, :Y) ≈ zeros(8,1)
        @test convert(Matrix, peers(holderio2)) ≈
            [ 1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
            0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0
            0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0
            0.4  0.0  0.6  0.0  0.0  0.0  0.0  0.0
            0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0
            0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0
            0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0
            0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0]

        @test efficiency(deaholder(targets(holderio2, :X), targets(holderio2, :Y), l = 2, orient = :Input, rts = :VRS, optimizer = DEAOptimizer(Gurobi.Optimizer))) ≈ zeros(8,1)

        # Weighted (weakly)
        holderio2weight = deaholder(X, Y, l = 2, orient = :Input, rts = :VRS, weight = true, optimizer = DEAOptimizer(Gurobi.Optimizer))
        @test efficiency(holderio2weight) ≈ [0.0; 0.0; 0.0; 0.554700; 0.8; 0.0; 0.468521; 0.375] atol = 1e-5
    else
        # Function runs without error
        @test occursin("Model solved with Ipopt.Optimizer is innacuate. Use a solver that supports SOS1 constraints.", string(logs))

        # Weighted (weakly)
        logs, value = Test.collect_test_logs() do
            holderio2weight = deaholder(X, Y, l = 2, orient = :Input, rts = :VRS, weight = true, optimizer = DEAOptimizer(:NLP))
        end
    end

    # Hölder l = Inf
    holderioinf = deaholder(X, Y, l = Inf, orient = :Input, rts = :VRS)

    @test typeof(holderioinf) == HolderLInfDEAModel

    @test nobs(holderioinf) == 8
    @test ninputs(holderioinf) == 2
    @test noutputs(holderioinf) == 1
    @test efficiency(holderioinf) ≈ [0.0; 0.0; 0.0; 4/3; 3.0; 0.0; 1.0; 0.6]
    @test convert(Matrix, peers(holderioinf)) ≈
            [ 1.0       0.0  0.0       0.0  0.0  0.0  0.0  0.0
            0.0       1.0  0.0       0.0  0.0  0.0  0.0  0.0
            0.0       0.0  1.0       0.0  0.0  0.0  0.0  0.0
            2/3       0.0  1/3       0.0  0.0  0.0  0.0  0.0
            1.0       0.0  0.0       0.0  0.0  0.0  0.0  0.0
            0.0       0.0  1.0       0.0  0.0  0.0  0.0  0.0
            0.0       1.0  0.0       0.0  0.0  0.0  0.0  0.0
            0.0       1.0  0.0       0.0  0.0  0.0  0.0  0.0]
    @test slacks(holderioinf, :X) ≈ [0.0 0.0; 0.0 0.0; 0.0 0.0; 0.0 0.0; 0.0 0.0; 2.0 0.0; 0.0 0.0; 0.0 3.4]
    @test slacks(holderioinf, :Y) ≈ zeros(8,1)

    @test efficiency(deaholder(targets(holderioinf, :X), targets(holderioinf, :Y), l = 1, orient = :Input, rts = :VRS)) ≈ zeros(8,1) atol = 1e-5

    # Test no slacks
    holderioinfnoslacks = deaholder(X, Y, l = Inf, orient = :Input, rts = :VRS, slack = false)
    @test efficiency(holderioinfnoslacks) == efficiency(holderioinf)
    @test isempty(slacks(holderioinfnoslacks, :X)) == 1
    @test isempty(slacks(holderioinfnoslacks, :Y)) == 1

    # Weighted (weakly)
    holderioinfweight = deaholder(X, Y, l = Inf, orient = :Input, rts = :VRS, weight = true)
    @test efficiency(holderioinfweight) ≈ [0; 0; 0; 0.4; 0.6; 0.0; 1/3; 0.375]

    # Print
    show(IOBuffer(), holderio1)
    show(IOBuffer(), holderio1noslacks)
    show(IOBuffer(), holderio1weight)
    show(IOBuffer(), holderio2)
    show(IOBuffer(), holderio2weight)
    show(IOBuffer(), holderioinf)
    show(IOBuffer(), holderioinfnoslacks)
    show(IOBuffer(), holderioinfweight)

    # ------------------
    # Output oriented
    # ------------------
    X = [1; 1; 1; 1; 1; 1; 1; 1]
    Y = [7 7; 4 8; 8 4; 3 5; 3 3; 8 2; 6 4; 1.5 5]

    # Hölder l = 1
    holderoo1 = deaholder(X, Y, l = 1, orient = :Output, rts = :VRS)

    @test typeof(holderoo1) == HolderL1DEAModel

    @test nobs(holderoo1) == 8
    @test ninputs(holderoo1) == 1
    @test noutputs(holderoo1) == 2
    @test efficiency(holderoo1) ≈ [0; 0; 0; 3.0; 5.0; 0.0; 2.0; 3.0]
    @test efficiency(holderoo1, :min) == [2; 2; 2; 3; 2; 2; 2; 3]
    @test convert(Matrix, peers(holderoo1)) ≈
            [ 1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
            0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0
            0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0
            0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0
            0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0
            0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0
            0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0
            0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0]
    @test slacks(holderoo1, :X) ≈ zeros(8,1)
    @test slacks(holderoo1, :Y) ≈ [0.0 0.0; 0.0 0.0; 0.0 0.0; 1.0 0.0; 0.0 1.0; 0.0 2.0; 0.0 0.0; 2.5 0.0]

    @test efficiency(deaholder(targets(holderoo1, :X), targets(holderoo1, :Y), l = 1, orient = :Output, rts = :VRS)) ≈ zeros(8,1) atol = 1e-5

    # Test no slacks
    holderoo1noslacks = deaholder(X, Y, l = 1, orient = :Output, rts = :VRS, slack = false)
    @test efficiency(holderoo1noslacks) == efficiency(holderoo1)
    @test isempty(slacks(holderoo1noslacks, :X)) == 1
    @test isempty(slacks(holderoo1noslacks, :Y)) == 1

    # Weighted (weakly)
    holderoo1weight = deaholder(X, Y, l = 1, orient = :Output, rts = :VRS, weight = true)
    @test efficiency(holderoo1weight) ≈ [0; 0; 0; 0.6; 5/3; 0.0; 1/3; 0.6]

    # Hölder l = 2
    logs, value = Test.collect_test_logs() do
        holderoo2 = deaholder(X, Y, l = 2, orient = :Output, rts = :VRS, optimizer = DEAOptimizer(:NLP))
    end

    @test typeof(holderoo2) == HolderL2DEAModel

    @test nobs(holderoo2) == 8
    @test ninputs(holderoo2) == 1
    @test noutputs(holderoo2) == 2

    if (Base.find_package("Gurobi") !== nothing)
        using Gurobi
        @info ("Testing Hölder L2 with Gurobi")

        holderoo2 = deaholder(X, Y, l = 2, orient = :Output, rts = :VRS, optimizer = DEAOptimizer(Gurobi.Optimizer))
        @test efficiency(holderoo2) ≈ [0.0; 0.0; 0.0; 3.0; 5.0; 0.0; 1.897367; 3.0] atol = 1e-5
        @test slacks(holderoo2, :X) ≈ zeros(8,1)
        @test slacks(holderoo2, :Y) ≈ [0.0 0.0; 0.0 0.0; 0.0 0.0; 1.0 0.0; 1.0 0.0; 0.0 0.0; 0.0 0.0; 2.5 0.0]
        @test convert(Matrix, peers(holderoo2)) ≈
            [ 1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
            0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0
            0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0
            0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0
            0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0
            0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0
            0.2  0.0  0.8  0.0  0.0  0.0  0.0  0.0
            0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0]

        @test efficiency(deaholder(targets(holderoo2, :X), targets(holderoo2, :Y), l = 2, orient = :Output, rts = :VRS, optimizer = DEAOptimizer(Gurobi.Optimizer))) ≈ zeros(8,1) atol = 1e-5

        # Weighted (weakly)
        holderoo2weight = deaholder(X, Y, l = 2, orient = :Output, rts = :VRS, weight = true, optimizer = DEAOptimizer(Gurobi.Optimizer))
        @test efficiency(holderoo2weight) ≈ [0.0; 0.0; 0.0; 0.6; 5/3; 0.0; 0.325396; 0.6] atol = 1e-5
    else
        # Function runs without error
        @test occursin("Model solved with Ipopt.Optimizer is innacuate. Use a solver that supports SOS1 constraints.", string(logs))

        # Weighted (weakly)
        logs, value = Test.collect_test_logs() do
            holderoo2weight = deaholder(X, Y, l = 2, orient = :Output, rts = :VRS, weight = true, optimizer = DEAOptimizer(:NLP))
        end
    end

    # Hölder l = Inf
    holderooinf = deaholder(X, Y, l = Inf, orient = :Output, rts = :VRS)

    @test typeof(holderooinf) == HolderLInfDEAModel

    @test nobs(holderooinf) == 8
    @test ninputs(holderooinf) == 1
    @test noutputs(holderooinf) == 2
    @test efficiency(holderooinf) ≈ [0.0; 0.0; 0.0; 2.5; 4.0; 0.0; 1.5; 2.875]
    @test convert(Matrix, peers(holderooinf)) ≈
            [ 1.0    0.0    0.0  0.0  0.0  0.0  0.0  0.0
            0.0    1.0    0.0  0.0  0.0  0.0  0.0  0.0
            0.0    0.0    1.0  0.0  0.0  0.0  0.0  0.0
            0.5    0.5    0.0  0.0  0.0  0.0  0.0  0.0
            1.0    0.0    0.0  0.0  0.0  0.0  0.0  0.0
            0.0    0.0    1.0  0.0  0.0  0.0  0.0  0.0
            0.5    0.0    0.5  0.0  0.0  0.0  0.0  0.0
            0.125  0.875  0.0  0.0  0.0  0.0  0.0  0.0]
    @test slacks(holderooinf, :X) ≈ zeros(8,1)
    @test slacks(holderooinf, :Y) ≈ [0.0 0.0; 0.0 0.0; 0.0 0.0; 0.0 0.0; 0.0 0.0; 0.0 2.0; 0.0 0.0; 0.0 0.0]

    @test efficiency(deaholder(targets(holderooinf, :X), targets(holderooinf, :Y), l = 1, orient = :Output, rts = :VRS)) ≈ zeros(8,1) atol = 1e-5

    # Test no slacks
    holderooinfnoslacks = deaholder(X, Y, l = Inf, orient = :Output, rts = :VRS, slack = false)
    @test efficiency(holderooinfnoslacks) == efficiency(holderooinf)
    @test isempty(slacks(holderooinfnoslacks, :X)) == 1
    @test isempty(slacks(holderooinfnoslacks, :Y)) == 1

    # Weighted (weakly)
    holderooinfweight = deaholder(X, Y, l = Inf, orient = :Output, rts = :VRS, weight = true)
    @test efficiency(holderooinfweight) ≈ [0.0; 0.0; 0.0; 5/9; 4/3; 0.0; 3/11; 0.6]

    # Print
    show(IOBuffer(), holderoo1)
    show(IOBuffer(), holderoo1noslacks)
    show(IOBuffer(), holderoo1weight)
    show(IOBuffer(), holderoo2)
    show(IOBuffer(), holderoo2weight)
    show(IOBuffer(), holderooinf)
    show(IOBuffer(), holderooinfnoslacks)
    show(IOBuffer(), holderooinfweight)

    # ------------------
    # Graph
    # ------------------
    X = [2; 4; 8; 12; 6; 14; 14; 9.412]
    Y = [1; 5; 8; 9; 3; 7; 9; 2.353]

    # Hölder l = 1
    holdergr1 = deaholder(X, Y, l = 1, orient = :Graph, rts = :VRS)

    @test typeof(holdergr1) == HolderL1DEAModel

    @test nobs(holdergr1) == 8
    @test ninputs(holdergr1) == 1
    @test noutputs(holdergr1) == 1
    @test efficiency(holdergr1) ≈ [0; 0; 0; 0.0; 3.0; 2.0; 0.0; 6.0]
    @test efficiency(holdergr1, :min) == [1; 1; 1; 1; 1; 2; 2; 2]
    @test convert(Matrix, peers(holdergr1)) ≈
            [1.0  0.0  0.0    0.0    0.0  0.0  0.0  0.0
            0.0  1.0  0.0    0.0    0.0  0.0  0.0  0.0
            0.0  0.0  1.0    0.0    0.0  0.0  0.0  0.0
            0.0  0.0  0.0    1.0    0.0  0.0  0.0  0.0
            0.5  0.5  0.0    0.0    0.0  0.0  0.0  0.0
            0.0  0.0  0.0    1.0    0.0  0.0  0.0  0.0
            0.0  0.0  0.0    1.0    0.0  0.0  0.0  0.0
            0.0  0.0  0.647  0.353  0.0  0.0  0.0  0.0]
    @test slacks(holdergr1, :X) ≈ [0; 0; 0; 0.0; 0.0; 2.0; 2.0; 0.0]
    @test slacks(holdergr1, :Y) ≈ zeros(8,1) atol = 1e-5

    @test efficiency(deaholder(targets(holdergr1, :X), targets(holdergr1, :Y), l = 1, orient = :Graph, rts = :VRS)) ≈ zeros(8,1) atol = 1e-5

    # Test no slacks
    holdergr1noslacks = deaholder(X, Y, l = 1, orient = :Graph, rts = :VRS, slack = false)
    @test efficiency(holdergr1noslacks) == efficiency(holdergr1)
    @test isempty(slacks(holdergr1noslacks, :X)) == 1
    @test isempty(slacks(holdergr1noslacks, :Y)) == 1

    # Weighted (weakly)
    holdergr1weight = deaholder(X, Y, l = 1, orient = :Graph, rts = :VRS, weight = true)
    @test efficiency(holdergr1weight) ≈ [0; 0; 0; 0.0; 0.5; 2/7; 0.0; 0.715629] atol = 1e-5

    # Hölder l = 2
    logs, value = Test.collect_test_logs() do
        holdergr2 = deaholder(X, Y, l = 2, orient = :Graph, rts = :VRS, optimizer = DEAOptimizer(:NLP))
    end

    @test typeof(holdergr2) == HolderL2DEAModel

    @test nobs(holdergr2) == 8
    @test ninputs(holdergr2) == 1
    @test noutputs(holdergr2) == 1

    if (Base.find_package("Gurobi") !== nothing)
        using Gurobi
        @info ("Testing Hölder L2 with Gurobi")

        holdergr2 = deaholder(X, Y, l = 2, orient = :Graph, rts = :VRS, optimizer = DEAOptimizer(Gurobi.Optimizer))
        @test efficiency(holdergr2) ≈ [0.0; 0.0; 0.0; 0.0; 2.683282; 2.0; 0.0; 5.3648] atol = 1e-5
        @test slacks(holdergr2, :X) ≈ [0.0; 0.0; 0.0; 0.0; 0.0; 2.0; 0.0; 0.0]
        @test slacks(holdergr2, :Y) ≈ [0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0]
        @test convert(Matrix, peers(holdergr2)) ≈
            [1.0  0.0      0.0      0.0  0.0  0.0  0.0  0.0
            0.0  1.0      0.0      0.0  0.0  0.0  0.0  0.0
            0.0  0.0      1.0      0.0  0.0  0.0  0.0  0.0
            0.0  0.0      0.0      1.0  0.0  0.0  0.0  0.0
            0.2  0.8      0.0      0.0  0.0  0.0  0.0  0.0
            0.0  0.0      0.0      1.0  0.0  0.0  0.0  0.0
            0.0  0.0      0.0      0.0  0.0  0.0  1.0  0.0
            0.0  0.45172  0.54828  0.0  0.0  0.0  0.0  0.0]

        @test efficiency(deaholder(targets(holdergr2, :X), targets(holdergr2, :Y), l = 2, orient = :Graph, rts = :VRS, optimizer = DEAOptimizer(Gurobi.Optimizer))) ≈ zeros(8,1) atol = 1e-5

        # Weighted (weakly)
        holdergr2weight = deaholder(X, Y, l = 2, orient = :Graph, rts = :VRS, weight = true, optimizer = DEAOptimizer(Gurobi.Optimizer))
        @test efficiency(holdergr2weight) ≈  [0.0; 0.0; 0.0; 0.0; 0.485071; 2/7; 0.0; 0.710103] atol = 1e-5
    else
        # Function runs without error
        @test occursin("Model solved with Ipopt.Optimizer is innacuate. Use a solver that supports SOS1 constraints.", string(logs))

        # Weighted (weakly)
        logs, value = Test.collect_test_logs() do
            holdergr2weight = deaholder(X, Y, l = 2, orient = :Graph, rts = :VRS, weight = true, optimizer = DEAOptimizer(:NLP))
        end
    end

    # Hölder l = Inf
    holdergrinf = deaholder(X, Y, l = Inf, orient = :Graph, rts = :VRS)

    @test typeof(holdergrinf) == HolderLInfDEAModel

    @test nobs(holdergrinf) == 8
    @test ninputs(holdergrinf) == 1
    @test noutputs(holdergrinf) == 1
    @test efficiency(holdergrinf) ≈ [0.0; 0.0; 0.0; 0.0; 2.0; 2.0; 0.0; 3.832]
    @test convert(Matrix, peers(holdergrinf)) ≈
            [ 1.0  0.0    0.0    0.0  0.0  0.0  0.0  0.0
            0.0  1.0    0.0    0.0  0.0  0.0  0.0  0.0
            0.0  0.0    1.0    0.0  0.0  0.0  0.0  0.0
            0.0  0.0    0.0    1.0  0.0  0.0  0.0  0.0
            0.0  1.0    0.0    0.0  0.0  0.0  0.0  0.0
            0.0  0.0    0.0    1.0  0.0  0.0  0.0  0.0
            0.0  0.0    0.0    1.0  0.0  0.0  0.0  0.0
            0.0  0.605  0.395  0.0  0.0  0.0  0.0  0.0]
    @test slacks(holdergrinf, :X) ≈ [0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 2.0; 0.0] atol = 1e-5
    @test slacks(holdergrinf, :Y) ≈ [0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0] atol = 1e-5

    @test efficiency(deaholder(targets(holdergrinf, :X), targets(holdergrinf, :Y), l = 1, orient = :Graph, rts = :VRS)) ≈ zeros(8,1) atol = 1e-5

    # Test no slacks
    holdergrinfnoslacks = deaholder(X, Y, l = Inf, orient = :Graph, rts = :VRS, slack = false)
    @test efficiency(holdergrinfnoslacks) == efficiency(holdergrinf)
    @test isempty(slacks(holdergrinfnoslacks, :X)) == 1
    @test isempty(slacks(holdergrinfnoslacks, :Y)) == 1

    # Weighted (weakly)
    holdergrinfweight = deaholder(X, Y, l = Inf, orient = :Graph, rts = :VRS, weight = true)
    @test efficiency(holdergrinfweight) ≈ [0.0; 0.0; 0.0; 0.0; 0.4; 5/21; 0.0; 0.636115] atol = 1e-5

    # Print
    show(IOBuffer(), holdergr1)
    show(IOBuffer(), holdergr1noslacks)
    show(IOBuffer(), holdergr1weight)
    show(IOBuffer(), holdergr2)
    show(IOBuffer(), holdergr2weight)
    show(IOBuffer(), holdergrinf)
    show(IOBuffer(), holdergrinfnoslacks)
    show(IOBuffer(), holdergrinfweight)

    # ------------------
    # Test errors
    # ------------------
    @test_throws DimensionMismatch deaholder([1; 2 ; 3], [4 ; 5], l = 1) #  Different number of observations
    @test_throws DimensionMismatch deaholder([1; 2], [4 ; 5], l = 1, Xref = [1; 2; 3; 4]) # Different number of observations in reference sets
    @test_throws DimensionMismatch deaholder([1 1; 2 2], [4 4; 5 5], l = 1, Xref = [1 1 1; 2 2 2]) # Different number of inputs
    @test_throws DimensionMismatch deaholder([1 1; 2 2], [4 4; 5 5], l = 1, Yref = [4 4 4; 5 5 5]) # Different number of inputs
    @test_throws ArgumentError deaholder([1; 2; 3], [4; 5; 6], l = 1, orient = :Error) # Invalid orientation
    @test_throws ArgumentError deaholder([1; 2; 3], [4; 5; 6], l = 1, rts = :Error) # Invalid returns to scale
    @test_throws ArgumentError deaholder([1; 2; 3], [4; 5; 6], l = 2, rts = :Error) # CRS in L2
    @test_throws ArgumentError deaholder([1; 2; 3], [4; 5; 6], l = 5) # Invalid r
    @test_throws ArgumentError deaholder([1; 2; 3], [4; 5; 6], l = 2, rts = :VRS) # No solver
    @test_throws ArgumentError efficiency(holderio1, :Error) # Invalid efficiency type

end


text:
#=
Linear Regression using Flux.jl

pkg> add Flux#master
pkg> add Plots

Reference:
https://qiita.com/cometscome_phys/items/f58174c0dad7ecb811ed

=#

using Random
using Flux
using Plots
#using CuArrays

function vandermonde_matrix(xs,maxdim)
    φ(xs,k) = xs.^k
    ret = zeros(Float32,maxdim,length(xs))
    for k in 0:maxdim-1
        ret[k+1,:]=φ(xs,k)
    end
    ret
end

function target_func(x)
    a₀,a₁,b₀ = 3.0,2.0,1.0
    return @. a₀*x + a₁*x^2 + b₀+ 3*cos(20*x)
end

function train(epochs,batchsize)
    k=4
    xs = range(-2,length=100,stop=2)
    X = vandermonde_matrix(xs,k)
    ys = target_func(xs)

    model = Dense(k,1)
    lossfn(x,y)=Flux.mse(model(x),y)
    optimizer = ADAM(params(model))
    for e in 1:epochs
        shuffle_indices=shuffle(1:length(ys))
        loss=0
        for i in 1:floor(Int,length(ys)/batchsize)
            batch_indice=shuffle_indices[(i-1)*batchsize+1:i*batchsize]
            batch_X=X[:,batch_indice]
            batch_ys=ys[batch_indice]
            data=[(batch_X[:,i], batch_ys[i]) for i in 1:batchsize]
            Flux.train!(lossfn, data, optimizer)
            loss=sum(lossfn(batch_X[:,i],batch_ys[i]) for i in 1:batchsize)/batchsize
        end
        if e%100==0
            @show e
            @show loss
        end
    end

    y_predict = [model(X[:,i]).data[1] for i=1:length(ys)]
    plt = plot(xs, [ys,y_predict], marker=:circle,label=["Data","Estimation"])

end

function main()
    train(1000,16)
end

main()

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
using JetPack, Jets, Test

n1,n2 = 101,201

@testset "JopReshape, linearity test, T=$T" for T in (Float64,Float32)
    A = JopReshape(JetSpace(T,n1,n2), JetSpace(T,n1,n2,1))
    lhs,rhs = linearity_test(A)
    @test lhs ≈ rhs
end

@testset "JopReshape, dot product test, T=$T" for T in (Float64,Float32)
    A = JopReshape(JetSpace(T,n1,n2,1), JetSpace(T,n1,n2))
    lhs,rhs = dot_product_test(A, rand(domain(A)), rand(range(A)))
    @test lhs ≈ rhs rtol=1e-6

    A = JopReshape(JetSpace(T,n1,n2), JetSpace(T,n1,n2,1))
    @test lhs ≈ rhs rtol=1e-6
end

@testset "JopReshape, correctness" begin
    A = JopReshape(JetSpace(Float64,n1*n2), JetSpace(Float64,n1,n2))
    m = rand(domain(A))
    d = A * m
    @test isa(d, Array{Float64,2})
    @test d[:] ≈ m[:]
    a = A' * d
    @test isa(a, Vector{Float64})
    @test a[:] ≈ m[:]
end


text:
"""
    centroidal quadruped
    q = (p, r, f1, f2, f3, f4)
        p - body position
        r - body orientation (modified Rodriques parameters)
        f1 - foot 1 position
        f2 - foot 2 position
        f3 - foot 3 position
        f4 - foot 4 position
"""
mutable struct CentroidalQuadrupedIC{T} <: Model{T}
    # dimensions
	nq::Int # generalized coordinates
    nu::Int # controls
    nw::Int # parameters
    nc::Int # contact points

    # parameters
    mass_body::T
    inertia_body::Matrix{T}
    mass_foot::T

    # environment
    friction_joint::Vector{T}
    friction_body_world::Vector{T}
    friction_foot_world::Vector{T}
    gravity::Vector{T}
end

# function skew(x)
#     return [0.0  -x[3]  x[2];
#             x[3]   0.0 -x[1];
#            -x[2]  x[1]   0.0]
# end

# function L_mult(x)
#     [x[1] -transpose(x[2:4]);
#      x[2:4] x[1] * I(3) + skew(x[2:4])]
# end

# # right quaternion multiply as matrix
# function R_mult(x)
#     [x[1] -transpose(x[2:4]); x[2:4] x[1] * I(3) - skew(x[2:4])]
# end

# # rotation matrix
# function quaternion_rotation_matrix(q)
#     H = [zeros(1, 3); I(3)]
#     transpose(H) * L_mult(q) * transpose(R_mult(q)) * H
# end

# function quaternion_from_mrp(p)
#     """Quaternion (scalar first) from MRP"""
#     return (1.0 / (1.0 + dot(p, p))) * [(1 - dot(p, p)); 2.0 * p]
# end

# function mrp_rotation_matrix(x)
#     quaternion_rotation_matrix(quaternion_from_mrp(x))
# end

# #https://en.wikipedia.org/wiki/Rotation_matrix#In_three_dimensions

# function euler_rotation_matrix(θ)
# 	a = θ[1]
# 	b = θ[2]
# 	c = θ[3]

# 	[cos(a) * cos(b) (cos(a) * sin(b) * sin(c) - sin(a) * cos(c)) (cos(a) * sin(b) * cos(c) + sin(a) * sin(c));
# 	 sin(a) * cos(b) (sin(a) * sin(b) * sin(c) + cos(a) * cos(c)) (sin(a) * sin(b) * cos(c) - cos(a) * sin(c));
# 	 -sin(b) cos(b) * sin(c) cos(b) * cos(c)]
# end

function mass_matrix(model::CentroidalQuadrupedIC, q)
    cat(
        model.mass_body * Diagonal(ones(3)),     # body position
        model.inertia_body,                      # body orienation
        model.mass_foot * Diagonal(ones(3 * 4)), # feet position
        dims=(1, 2)
        )
end

function dynamics_bias(model::CentroidalQuadrupedIC, q, q̇)
    [
        model.mass_body * model.gravity;            # body position
        skew(q̇[4:6]) * model.inertia_body * q̇[4:6]; # body orienation
        model.mass_foot * model.gravity;
        model.mass_foot * model.gravity;
        model.mass_foot * model.gravity;
        model.mass_foot * model.gravity;
    ]
end

function signed_distance(model::CentroidalQuadrupedIC, q)

    position_body = q[1:3]
    orientation_body = q[3 .+ (1:3)]

    position_foot1 = q[6 .+ (1:3)]
    position_foot2 = q[9 .+ (1:3)]
    position_foot3 = q[12 .+ (1:3)]
	position_foot4 = q[15 .+ (1:3)]

    return [position_foot1[3]; position_foot2[3]; position_foot3[3]; position_foot4[3]; position_body[3]]
end

function input_jacobian(model::CentroidalQuadrupedIC, q)
    position_body = q[1:3]
    orientation_body = q[3 .+ (1:3)]
    # R = mrp_rotation_matrix(orientation_body)
    R = euler_rotation_matrix(orientation_body)

    # kinematics in world frame
	r1 = q[6 .+ (1:3)] - position_body
	r2 = q[9 .+ (1:3)] - position_body
	r3 = q[12 .+ (1:3)] - position_body
	r4 = q[15 .+ (1:3)] - position_body

	z3 = zeros(3, 3)

	transpose([
        I(3) z3 I(3) I(3) I(3) I(3);
        z3 I(3) transpose(R) * skew(r1) transpose(R) * skew(r2) transpose(R) * skew(r3) transpose(R) * skew(r4);
        z3 z3 -I(3)    z3    z3   z3;
        z3 z3 z3    -I(3)    z3   z3;
        z3 z3 z3       z3 -I(3)   z3;
        z3 z3 z3       z3    z3 -I(3)
    ])
end

function contact_jacobian(model::CentroidalQuadrupedIC, q)
    z3 = zeros(3, 3)

    [
        z3   z3 I(3)   z3   z3   z3;
        z3   z3   z3 I(3)   z3   z3;
        z3   z3   z3   z3 I(3)   z3;
        z3   z3   z3   z3   z3 I(3);
        I(3) z3   z3   z3   z3   z3;
    ]
end

# nominal configuration
function nominal_configuration(model::CentroidalQuadrupedIC)
    [
        0.0; 0.0; 0.15;
        0.0; 0.0; 0.0;
        0.1; 0.1; 0.0;
        0.1;-0.1; 0.0;
       -0.1; 0.1; 0.0;
       -0.1;-0.1; 0.0;
    ]
end

# friction coefficients
function friction_coefficients(model::CentroidalQuadrupedIC)
    return [model.friction_foot_world; model.friction_body_world]
end

# dimensions
nq = 3 + 3 + 3 * 4       # generalized coordinates
nu = 3 * 4 + 6           # controls
nw = 0                   # parameters
nc = 5                   # contact points

# parameters
gravity = [0.0; 0.0; 9.81]                 # gravity
friction_body_world = [0.5]                # coefficient of friction
friction_foot_world = [0.5; 0.5; 0.5; 0.5] # coefficient of friction

# inertial properties
mass_body = 1.0
inertia_body = Array(Diagonal(ones(3)))
mass_foot = 0.1

centroidal_quadruped_ic = CentroidalQuadrupedIC(nq, nu, nw, nc,
				mass_body,
                inertia_body,
                mass_foot,
                zeros(nq - 6),
                friction_body_world,
                friction_foot_world,
                gravity)

centroidal_quadruped_ic_contact_kinematics = [
    q -> q[6  .+ (1:3)],
    q -> q[9  .+ (1:3)],
    q -> q[12 .+ (1:3)],
    q -> q[15 .+ (1:3)],
    q -> q[ 0 .+ (1:3)],
]

centroidal_quadruped_ic_contact_kinematics_jacobians = [
    q -> [zeros(3, 6) I(3) zeros(3, 9)],
    q -> [zeros(3, 9) I(3) zeros(3, 6)],
    q -> [zeros(3, 12) I(3) zeros(3, 3)],
    q -> [zeros(3, 15) I(3)],
    q -> [I(3) zeros(3, 15)],
]

name(::CentroidalQuadrupedIC) = :centroidal_quadruped_ic
floating_base_dim(::CentroidalQuadrupedIC) = 6


text:
module UnitfulPlots
using Unitful
using RecipesBase

function unit_formatter(T, num)
    string(num)*string(unit(zero(T)))
end

# 2D plots (1 recipe for shorthand form)
@recipe function f(val::AbstractArray{T,N}) where {T<:Quantity,N}
    yformatter := x->unit_formatter(T,x)
    ustrip.(val)
end

# 2D plots (2^4-1 user recipes)
@recipe function f(val::AbstractArray{S,N}, val2::AbstractArray{T,M}
                  ) where {S<:Quantity, T<:Quantity, N, M}
    xformatter := x->unit_formatter(S,x)
    yformatter := x->unit_formatter(T,x)
    ustrip.(val), ustrip.(val2)
end

@recipe function f(val, val2::AbstractArray{T,N}) where {T<:Quantity, N}
    yformatter := x->unit_formatter(T,x)
    val, ustrip.(val2)
end

@recipe function f(val::AbstractArray{T,N}, val2) where {T<:Quantity, N}
    xformatter := x->unit_formatter(T,x)
    ustrip.(val), val2
end

# 3D plots (2^8-1 user recipes)
@recipe function f(val::AbstractArray{S,N}, val2::AbstractArray{T,M}, val3::AbstractArray{U,P}
                  ) where {S<:Quantity, T<:Quantity, U<:Quantity, N, M, P}
    xformatter := x->unit_formatter(S,x)
    yformatter := x->unit_formatter(T,x)
    zformatter := x->unit_formatter(U,x)
    ustrip.(val), ustrip.(val2), ustrip.(val3)
end

@recipe function f(val::AbstractArray{S,N}, val2::AbstractArray{T,M}, val3
                  ) where {S<:Quantity, T<:Quantity, N, M}
    xformatter := x->unit_formatter(S,x)
    yformatter := x->unit_formatter(T,x)
    ustrip.(val), ustrip.(val2), val3
end

@recipe function f(val::AbstractArray{S,N}, val2, val3::AbstractArray{U,P}
                  ) where {S<:Quantity, U<:Quantity, N, P}
    xformatter := x->unit_formatter(S,x)
    zformatter := x->unit_formatter(U,x)
    ustrip.(val), val2, ustrip.(val3)
end

@recipe function f(val, val2::AbstractArray{T,M}, val3::AbstractArray{U,P}
                  ) where {T<:Quantity, U<:Quantity, M, P}
    yformatter := x->unit_formatter(T,x)
    zformatter := x->unit_formatter(U,x)
    val, ustrip.(val2), ustrip.(val3)
end

@recipe function f(val::AbstractArray{T,N}, val2, val3) where {T<:Quantity, N}
    xformatter := x->unit_formatter(T,x)
    ustrip.(val), val2, val3
end

@recipe function f(val, val2::AbstractArray{T,N}, val3) where {T<:Quantity, N}
    yformatter := x->unit_formatter(T,x)
    val, ustrip.(val2), val3
end

@recipe function f(val, val2, val3::AbstractArray{T,N}) where {T<:Quantity, N}
    zformatter := x->unit_formatter(T,x)
    val, val2, ustrip.(val3)
end

end


text:
new_note_header = "======== Input New Note ========"

function notes_handle_query(tags...)
    get_notes(tags...)
    nothing
end

function notes_handle_new(body, increment_id)
    tags = map(x->x[2:end], filter(t->startswith(t, "#") && length(t)>1, split(body)))
    if !("#" in tags) push!(tags, "#") end
    latest_id = increment_id()
    body = body[length(new_note_header)+1 : end]
    body = body[1 : search(body, "#")[1]-1]
    new_note = Note(latest_id, tags, body)
    add_note(new_note)
    nothing
end

function notes_on_done(query, increment_id)
    tags = split(query)

    if isempty(tags)
        tags = [""]
    end

    is_new = tags[1] == "new" || startswith(query, new_note_header)
    if is_new && length(tags) > 1
        body = strip(query)
        notes_handle_new(body, increment_id)
    else
        notes_handle_query(tags)
    end
    nothing
end

function create_on_done_fn(increment_id)
    done_func(query) = notes_on_done(query, increment_id)
    return done_func
end

function update_prompt(cur_mode, is_new, new_prompt, base_prompt)
    if is_new
        cur_mode.prompt = new_prompt
    else
        cur_mode.prompt = base_prompt
    end
end

new_note_header = "======== Input New Note ========"

function notes_is_complete(s, new_prompt, base_prompt)
    cur_mode = LineEdit.mode(s)
    # update_prompt(cur_mode, false, new_prompt, base_prompt)
    str = String(take!(copy(LineEdit.buffer(s))))
    if isempty(str) return true end
    if length(str) < 3 || (length(str) >= 3 && !startswith(str, "new") && !startswith(str, new_note_header)) return true end

    # Deal with the `new` keyword
    if length(str) >= 3 && str[1:3] == "new"
        # reset line and change input
        LineEdit.edit_delete_prev_word(s)
        LineEdit.edit_insert(s, new_note_header)
        # LineEdit.edit_kill_line(s)
        # update_prompt(cur_mode, true, new_prompt, base_prompt)
        return false
     end
    return str[end-1:end] == "##"
end

function create_on_complete_fn(new_prompt, base_prompt)
    return s -> notes_is_complete(s, new_prompt, base_prompt)
end


text:
export sendAnalogSignal, disableAnalogOutput, sendAnalogBurstSignal, trapeziod

### Analog Output Interface ###

function sendAnalogSignal(rp::RedPitaya, chan::Integer, func::String,
                          freq::Number, amplitude::Number; offset::Number=0,AWG_data::Array=[])
  if !in(func,["SINE", "SQUARE", "TRIANGLE", "SAWU", "SAWD", "PWD", "AWG"])
    error("waveform $func not supported")
  end
  source = "SOUR$(chan)"
  send(rp,"GEN:RST")
  send(rp,"$(source):FUNC $(func)") # Set function of output signal
  send(rp,"$(source):FREQ:FIX $(freq)") # Set frequency of output signal
  send(rp,"$(source):VOLT $(amplitude)") # Set amplitude of output signal
  send(rp,"$(source):VOLT:OFFS $(offset)") # Set offset of output signal
  if func=="AWG"
    send(rp, "$(source):TRAC:DATA:DATA $(AWG_data)")
  end
  send(rp,"OUTPUT$(chan):STATE ON") # Set output to ON
end

function sendAnalogBurstSignal(rp::RedPitaya, chan::Integer, func::String,
                          freq::Number, amplitude::Number; offset::Number=0,
                          ncycBurst::Integer=0, repetitions::Integer=0,
                          repetitionCycleTime::Float64=0 , AWG_data::Array=[])
  if !in(func,["SINE", "SQUARE", "TRIANGLE", "SAWU", "SAWD", "PWD", "ARBITRARY"])
    error("waveform $func not supported")
  end
  source = "SOUR$(chan)"
  send(rp,"GEN:RST")
  send(rp,"$(source):FUNC $(func)") # Set function of output signal
  send(rp,"$(source):FREQ:FIX $(freq)") # Set frequency of output signal
  send(rp,"$(source):VOLT $(amplitude)") # Set amplitude of output signal
  send(rp,"$(source):VOLT:OFFS $(offset)") # Set offset of output signal
  if ncycBurst > 0
   println("BurstMode On. Sending $(ncycBurst) Bursts")
   println("$(source):BURS:NCYC $ncycBurst")

   send(rp,"$(source):BURS:STAT ON") #Set burst mode to ON
  #bufferSize(rp)
   send(rp,"$(source):BURS:NCYC $ncycBurst") #Set 1 pulses of sine wave
   if repetitions!=0
     send(rp, "$(source):BURS:NOR $repetitions")
     println("$(repetitionCycleTime)")
     send(rp, "$(source):BURS:INT:PER $(repetitionCycleTime)")
   end
 end
   if func=="ARBITRARY"
    AWG_data=join(AWG_data, ",")
    send(rp, "$(source):TRAC:DATA:DATA $(AWG_data)")
  end
  send(rp,"OUTPUT$(chan):STATE ON") # Set output to ON
end

function disableAnalogOutput(rp::RedPitaya, chan::Integer)
  send(rp,"OUTPUT$(chan):STATE OFF") # Set output to ON
end
function trapeziod(flankratio; len=2^14, offset=0)
    x = range(0, stop=2*pi, length=len)     #generates clean periodic signal

    rise = x * 2/(1-flankratio)      #calc width of rising edge and mirror it ...
    fall = - rise
    edge_point = trunc(Int, (1-flankratio)/2 * len)
    top = trunc(Int, (2^14 - 4*edge_point ) / 2 )

    p1 = collect( rise[1 : edge_point] / (2*pi))
    p2 = collect( ones(top))
    p3 = collect( (fall[1 : 2*edge_point] / (2*pi) .+ 1))
    p4 = collect( -ones(top))
    p5 = collect( rise[1 : edge_point] / (2*pi) .- 1)

    trapez = vcat(p1,p2,p3,p4,p5).-offset

     if length(trapez) != len   #test
        error("ERROR: failed to concatenate array parts correctly, $(len) $(length(trapez))")
    end
return trapez
end