text:
module MultRoot
using Polynomials
import PolynomialRoots
using LinearAlgebra
include("../utils.jl")
## The main function here is `MultRoot.multroot`

## Polynomial root finder for polynomials with multiple roots
##
## Based on "Computing multiple roots of inexact polynomials"
## https://doi.org/10.1090/S0025-5718-04-01692-8
## Author: Zhonggang Zeng
## Journal: Math. Comp. 74 (2005), 869-903
##
## Zeng has a MATLAB package `multroot`, from which this name is derived.
## Basic idea is
## 1) for polynomial p we do gcd decomposition p = u * v; p' = u * w. Then roots(v) are the roots without multiplicities.
## 2) can repeat with u to get multiplicities.
##
## This is from Gauss, as explained in paper. Zeng shows how to get u,v,w when the polynomials
## are inexact due to floating point approximations or even model error. This is done in his
## algorithm II.
## 3) Zeng's algorithm I (pejroot) uses the pejorative manifold of Kahan and Gauss-Newton to
## improve the root estimates from algorithm II (roots(v)). The pejorative manifold is defined by
## the multiplicities l and is operationalized in evalG and evalJ from Zeng's paper.

using ..AGCD
monic(p) = p/p[end]
rcoeffs(p) = reverse(p.a)
function proots(zs)
    rs = PolynomialRoots.roots(zs)
    if all(iszero.(imag(rs)))
        return real.(rs)
    else
        return rs
    end
end

## map monic(p) to a point in C^n
## p = 1x^n + a1x^n-1 + ... + an_1 x + an -> (a1,a2,...,an)
function p2a(p::Poly)
    p = monic(p)
    rcoeffs(p)[2:end]
end
function p2a(p::Vector{T}) where {T}
    a,pn = reverse(p[1:end-1]),p[end]
    a .* (1/p[end])
end

## get value of gl(z). From p16
function evalG(zs::Vector, ls::Vector)
    length(zs) == length(ls) || throw("Length mismatch")

    s = prod([poly([z])^l for (z,l) in zip(zs, ls)])  # \prod (x-z_i)^l_i
    p2a(s)
#    rcoeffs(s)[2:end]
end

## get jacobian J_l(z), p16

## get jacobian J_l(z), p16
function evalJ!(J, zs::Vector{T}, ls::Vector) where {T}
    length(zs) == length(ls) || throw("Length mismatch")

    n, m = sum(ls), length(zs)
    u = evalG(zs, ls .- 1)
    pushfirst!(u, one(T))

#    x = Polynomials.variable(T)
#    u = prod((x-z)^(l-1) for (z,l) in zip(zs, ls))  # \prod (x-z_i)^l_i

    for j in 1:m
        s = -ls[j] * u

        for (l, zl) in zip(1:m, zs)
            l == j && continue
            s = AGCD._polymul(s, (one(T), -zl))
        end
    J[:,j] = s#rceffs(s)
    end
    J
end
function evalJ(zs::Vector{T}, ls) where {T}
    n, m = sum(ls), length(zs)
    J = zeros(T, n, m)
    evalJ!(J, zs, ls)
    J
end

## Gauss-Newton iteration to solve weighted least squares problem
## G_l(z) = a, where a is related to monic version of polynomial p
## l is known multiplicity structure of polynomial p = (x-z1)^l1 * (x-z2)^l2 * ... * (x-zn)^ln
## Algorithm I, p17
pejroot(p::Poly, z0, ls; kwargs...) = pejroot(Polynomials.coeffs(p), z0, ls; kwargs...)
function pejroot(p::Vector{T}, z0::Vector{S}, l::Vector{Int};
                 wts::Union{Vector, Nothing}=nothing, # weight vector
                 τ = sqrt(eps(real(T))),
                 maxsteps = 100
                 ) where {T, S <: Union{T, Complex{T}}}

    a = p2a(p)
    λ = min(sqrt(eps(real(T))), norm(p,2) * eps(real(T))^(2/3))
    if wts == nothing
        wts = map(aj -> min(1, 1/abs(aj)), a)
    end

    ## Solve WJ Δz = W(Gl(z) - a) in algorithm I
    zk = copy(z0);
    J = evalJ(zk, l)

    deltak = AGCD.weighted_least_square(J, evalG(zk,l).-a, wts)
    zk .-=  deltak
    δ0 = norm(deltak, 2)

    cvg = false

    for ctr in 1:maxsteps

        AGCD.weighted_least_square!(deltak, evalJ!(J, zk,l), evalG(zk,l).-a, wts)
        zk .-= deltak
        δ1 = norm(deltak, 2)

        Δ = δ0 - δ1

        if Δ < 0 && ctr > 2
            @debug "Growing delta. Best guess is being returned."
            break
        end

        ## add extra abs(delta) < 100*eps() condition
        if δ1^2 < Δ * τ# || δ1 < λ
            cvg = true
            break
        end

        δ0 = δ1
    end

    if !cvg @info ("""The multiplicity count may be in error--the initial guess for the roots failed to improve when refined along the pejorative manifold.""")
        return(z0)
    end
    return(zk)
end

"""
        identify_z0s_ls(p; ...)

Step one of the algorithm identifies initial guesses for the roots and
identifies the pejorative manifold.

"""
function identify_z0s_ls(p::Vector{T};
                         θ::Real=sqrt(eps(real(T))),
                         ρ::Real=cbrt(eps(real(T)))*θ, # initial residual tolerance
                         ϕ::Real=100.0,          # residual tolerance growth factor
                         precondition = false
                         ) where {T}

    q = AGCD._polyder(p)
    if precondition
        p,q, phi, alpha = AGCD.precondition(p,AGCD._polyder(p))
    else
        phi = one(T)
    end

    u_j, v_j, w_j, residual = AGCD.agcd(p, q, θ=θ,  ρ=ρ)

    ## bookkeeping
    ρ = max(ρ, ϕ * norm(residual))
    rs = PolynomialRoots.roots(v_j)
    zs = proots(v_j)
    N = length(zs)
    ls = ones(Int, N)

    p0 = u_j

    while AGCD._degree(p0) > 0

        if AGCD._degree(p0) == 1
            z = proots(p0)[1]
            tmp, ind = findmin(abs.(zs .- z))
            ls[ind] = ls[ind] + 1
            break
        end

        u_j, v_j, w_j, residual= AGCD.agcd(p0, θ=θ, ρ=ρ, maxk=N+1)

        ## need to worry about residual between
        ## u0 * v0 - monic(p0) and u0 * w0 - monic(Polynomials.polyder(p0))
        ## resiudal tolerance grows with m, here it depends on
        ## initial value and previous residual times a growth tolerance, ϕ
        ρ = max(ρ, ϕ * abs(residual))
        ## update multiplicities
        for z in proots(v_j)
            tmp, ind = findmin(abs.(zs .- z))
            ls[ind] = ls[ind] + 1
        end

        ## rename
        p0 = u_j
    end

    # remove preconditioning from roots
    phi * zs, ls, ρ

end

"""
    multroot(p; [θ, ρ, ϕ, δ])

Find roots of polynomial `p`,

The `multroot` function returns the roots and their multiplicities
for `Poly` objects. It performs better than `roots` when the
polynomial has multiplicities.

Based on "Computing multiple roots of inexact polynomials"
Zhonggang Zeng
Journal: Math. Comp. 74 (2005), 869-903
https://doi.org/10.1090/S0025-5718-04-01692-8

Zeng has a MATLAB package `multroot`, from which this name is derived.

Basic idea is:
* for polynomial p we do gcd decomposition p = u * v; p' = u * w. Then roots(v) are the roots without multiplicities.
* can repeat with u to get multiplicities.

The basic idea is from Gauss, as explained in paper. Zeng shows how to get u,v,w when the polynomials
are inexact due to floating point approximations or even model error. This is done in his
algorithm II.

Zeng's algorithm I (pejroot) uses the pejorative manifold of Kahan and Gauss-Newton to
improve the root estimates from algorithm II (roots(v)). The pejorative manifold is defined by
the multiplicities l and is operationalized in `evalG` and `evalJ` from Zeng's paper.

Examples:
```
x = poly([0.0]);
p = (x-1)^4 * (x-2)^3 * (x-3)^3 * (x-4)
multroot(p) # ([4.0, 3.0, 2.0, 1.0], [1, 3, 3, 4])

## For "prettier" printing, results can be coerced to a dict
Dict(k=>v for (k, v) in zip(multroot(p)...))
## Dict{Float64,Int64} with 4 entries:
##   4.0 => 1
##   2.0 => 3
##   1.0 => 4
##   3.0 => 3

## compare to
roots(p)
# 11-element Array{Complex{Float64},1}:
#   4.000000000049711 + 0.0im
#   3.000340992482669 + 0.0005901104690775108im
#   3.000340992482669 - 0.0005901104690775108im
#  2.9993180137669726 + 0.0im
#  2.0006129094631833 + 0.0im
#  1.9996935459268157 + 0.0005303204210266187im
#  1.9996935459268157 - 0.0005303204210266187im
#  1.0006651061397798 + 0.00066713538542472im
#  1.0006651061397798 - 0.00066713538542472im
#  0.9993348938107887 + 0.0006630950464016094im
#  0.9993348938107887 - 0.0006630950464016094im

## Large order polynomials prove difficult. We can't match the claims in Zeng's paper
## as we don't get the pejorative manifold structure right.
p = poly(1.0:7.0);
multroot(p^2) ## should be 1,2,3,4,...,7 all with multplicity 2, but
## ([7.00028, 6.99972, 6.00088, 5.99912, 5.00102, 4.99898, 4.00055, 3.99945, 3.00014, 2.99986, 2.00002, 1.99998, 1.0, 0.999999], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1])

## nearby roots can be an issue
delta = 0.00001  ## delta = 0.0001 works as desired.
p = (x-1 - delta)*(x-1)*(x-1 + delta)
multroot(p)
## ([1.0], [3])
```
"""
function multroot(ps::Vector{T};
                  θ::Real=1e-2 * sqrt(eps(real(T))),      # 1e-8 in paper
                  ρ::Real=(eps(real(T)))^(5/6), # 1e-10 in paper
                  ϕ::Real=100.0,             # residual tolerance growth factor
                  τ::Real= θ,       # passed to solve y sigma
                  precondition=true
                  ) where {T}

    p = float.(ps[1:findlast(!iszero,ps)])

    # simple cases
    AGCD._degree(p) == 0 && error("Degree of `p` must be atleast 1")
    if AGCD._degree(p) == 1
        return (-p[1:1]/p[end], [1])
    end

    # two steps
    zs, ls, rho = identify_z0s_ls(p, θ=θ, ρ=ρ, ϕ=ϕ, precondition=precondition)

    if maximum(ls) > 1
        zs = pejroot(p, zs, ls, τ=τ)
    end

    return (zs, ls)

end

## Different interfaces

## can pass in Poly too
multroot(p::Poly; kwargs...) = multroot(Polynomials.coeffs(p); kwargs...)
## Can pass in function
multroot(f::Function; kwargs...) = multroot(as_poly(Float64, f); kwargs...)

end


text:
"""
BatchedFeaturedGraph

A batched representation of the FeaturedGraph, to enable parallel computation.

It is deliberately more restrictive to prevent incorrect usage.
"""
struct BatchedFeaturedGraph{T <: Real} <: AbstractFeaturedGraph
    graph::AbstractArray{T, 3}
    nf::AbstractArray{T, 3}
    ef::AbstractArray{T, 4}
    gf::AbstractMatrix{T}

    function BatchedFeaturedGraph{T}(graph, nf, ef, gf) where T <: Real
        check_dimensions(graph, nf, ef, gf)
        return new{T}(graph, nf, ef, gf)
    end

    function BatchedFeaturedGraph{T}(
        graph;
        nf=zeros(0, size(graph, 1), size(graph, 3)),
        ef=zeros(0, size(graph, 1), size(graph, 1), size(graph, 3)),
        gf=zeros(0, size(graph, 3))
    ) where T <: Real
        check_dimensions(graph, nf, ef, gf)
        return new{T}(graph, nf, ef, gf)
    end
end

BatchedFeaturedGraph(graph, nf, ef, gf) = BatchedFeaturedGraph{Float32}(graph, nf, ef, gf)

function BatchedFeaturedGraph{T}(fgs::Vector{FG}) where {T <: Real, FG <: FeaturedGraph}
    ngraphs = length(fgs)
    maxNodes = Base.maximum(nv, fgs)
    nfLength = size(fgs[1].nf, 1)
    efLength = size(fgs[1].ef, 1)
    gfLength = size(fgs[1].gf, 1)

    graph = zeros(T, maxNodes, maxNodes, ngraphs)
    nf = zeros(T, nfLength, maxNodes, ngraphs)
    ef = zeros(T, efLength, maxNodes, maxNodes, ngraphs)
    gf = zeros(T, gfLength, ngraphs)

    for (i, fg) in enumerate(fgs)
        graph[1:nv(fg),1:nv(fg),i] = fg.graph
        nf[:, 1:nv(fg), i] = fg.nf
        ef[:, 1:nv(fg), 1:nv(fg), i] = fg.ef
        gf[:, i] = fg.gf
    end

    return BatchedFeaturedGraph{T}(graph, nf, ef, gf)
end

BatchedFeaturedGraph(fgs::Vector{FG}) where {FG <: FeaturedGraph} = BatchedFeaturedGraph{Float32}(fgs)

# ========== Accessing ==========
# code from GraphSignals.jl

"""
    graph(::BatchedFeaturedGraph)
Get referenced graph.
"""
graph(fgs::BatchedFeaturedGraph) = fgs.graph

"""
    node_feature(::BatchedFeaturedGraph)
Get node feature attached to graph.
"""
node_feature(fgs::BatchedFeaturedGraph) = fgs.nf

"""
    edge_feature(::BatchedFeaturedGraph)
Get edge feature attached to graph.
"""
edge_feature(fgs::BatchedFeaturedGraph) = fgs.ef

"""
    global_feature(::BatchedFeaturedGraph)
Get global feature attached to graph.
"""
global_feature(fgs::BatchedFeaturedGraph) = fgs.gf

"""
    has_graph(::BatchedFeaturedGraph)
Check if graph is available or not.
"""
has_graph(fgs::BatchedFeaturedGraph) = fgs.graph != zeros(0,0,0)

"""
    has_node_feature(::BatchedFeaturedGraph)
Check if node feature is available or not.
"""
has_node_feature(fgs::BatchedFeaturedGraph) = !isempty(fgs.nf)

"""
    has_edge_feature(::BatchedFeaturedGraph)
Check if edge feature is available or not.
"""
has_edge_feature(fgs::BatchedFeaturedGraph) = !isempty(fgs.ef)

"""
    has_global_feature(::BatchedFeaturedGraph)
Check if global feature is available or not.
"""
has_global_feature(fgs::BatchedFeaturedGraph) = !isempty(fgs.gf)

# ========== LightGraphs compatibility ==========

LightGraphs.nv(fgs::BatchedFeaturedGraph) = nv(graph(fgs))
LightGraphs.nv(g::AbstractArray{T, 3}) where T<:Real = size(g,1)

Flux.@functor BatchedFeaturedGraph

ext: jl


text:
props(so::AbstractSceneObject) = so._props
function prop(so::AbstractSceneObject, p::Symbol, default_val = nothing)
    return get(props(so), p, default_val)
end
function prop!(so::AbstractSceneObject, p::Symbol, val::Any)
    props(so)[p] = val
end

name(so::AbstractSceneObject) = prop(so, :name)

tr(so::AbstractSceneObject) = prop(so, :tr)
function tr!(so::AbstractSceneObject, tr::AffineMap)
    prop!(so, :tr, tr)
end

local_tr(so::AbstractSceneObject) = prop(so, :tr)
function local_tr!(so::AbstractSceneObject, local_tr::AffineMap)
    prop!(so, :tr, local_tr)
end

vis(so::AbstractSceneObject) = prop(so, :vis)
function vis!(so::AbstractSceneObject, vis::MeshCat.Visualizer)
    prop!(so, :vis, vis)
end

scene(so::AbstractSceneObject) = prop(so, :scene)
function scene!(so::AbstractSceneObject, s::AbstractScene)
    prop!(so, :scene, s)
end

parent(so::AbstractSceneObject) = prop(so, :parent)
function parent!(so::AbstractSceneObject, parent_so::AbstractSceneObject)
    prop!(so, :parent, parent_so)

    parent_children = prop(parent_so, :children)
    index = findfirst(x -> x===so, parent_children)
    if (index === nothing)
        push!(parent_children, so)
    end

    update!(so)
end

material(so::AbstractSceneObject) = prop(so, :material)
function material!(so::AbstractSceneObject, mat::Material)
    prop!(so, :material, mat)

    update!(so)
end

function update!(so::AbstractSceneObject)
    v = vis(parent(so))[name(so)]
    vis!(so, v)

    render!(so)
end

function Base.delete!(so::AbstractSceneObject)
    parent_so = prop(so, :parent)
    parent_children = prop(parent_so, :children)
    index = findfirst(x -> x===so, parent_children)
    if (index === nothing)
        @error "Can't find the child $(name(so)) to delete"
    end
    deleteat!(parent_children, index)

    v = vis(so)
    delete!(v)

    # render!(so)
end

function display(so::AbstractSceneObject, level=0)
    print(' '^(3*level))
    println(name(so))
    children = prop(so, :children)
    for child in children
        display(child, level+1)
    end
end

function build_props(;kwargs...)
    res = Dict{Symbol, Any}(
        :tr => identityTransform(),
        :children => Vector{AbstractSceneObject}(undef, 0),
    )
    for a in kwargs
        res[a[1]] = a[2]
    end
    return res
end

#-----------------------------------------------------------------
# Empty Scene Object
#-----------------------------------------------------------------
mutable struct EmptySceneObject <: AbstractSceneObject
    _props::Dict{Symbol, Any}

    function EmptySceneObject(;kwargs...)
        props = build_props(;kwargs...)
        return new(props)
    end
end

function render!(e::EmptySceneObject)
    v = vis(e)

    # @info local_tr(e), name(e)
    MeshCat.settransform!(v, local_tr(e))
end

#-----------------------------------------------------------------
# Box
#-----------------------------------------------------------------
mutable struct PointCloud <: AbstractSceneObject
    _props::Dict{Symbol, Any}

    function PointCloud(;name="defaultPoints", points::Vector{SVector{3, Float64}} = [], kwargs...)
        props = build_props(;name=name, points=points, kwargs...)

        return new(props)
    end
end

points(pc::PointCloud) = prop(pc, :points)

function render!(pc::PointCloud)
    v = vis(pc)

    points = prop(pc, :points)

    mat = material(pc)
    if (mat !== nothing)
        mat_color = color(mat)
        colors = repeat([mat_color], length(points))
        pointcloud = MeshCat.PointCloud(points, colors)
    else
        pointcloud = MeshCat.PointCloud(points)
    end

    if (mat !== nothing)
        MeshCat.setobject!(v, pointcloud, points_material(mat))
        # MeshCat.setobject!(v, pointcloud)
    else
        MeshCat.setobject!(v, pointcloud)
    end

    MeshCat.settransform!(v, local_tr(pc))
end

#-----------------------------------------------------------------
# Box
#-----------------------------------------------------------------
mutable struct Box <: AbstractSceneObject
    _props::Dict{Symbol, Any}

    function Box(;name="defaultBox", size::SVector{3, Float64} = SVector(1.0, 1.0, 1.0), kwargs...)
        props = build_props(;name=name, size=size, kwargs...)

        return new(props)
    end
end

Base.size(b::Box) = prop(b, :size)

function render!(b::Box)
    v = vis(b)
    min_point = size(b) * -0.5
    mat = material(b)
    box = GeometryBasics.HyperRectangle(GeometryBasics.Vec(min_point...), GeometryBasics.Vec(size(b)...))
    if (mat !== nothing)
        MeshCat.setobject!(v, box, material(mat))
    else
        MeshCat.setobject!(v, box)
    end

    affine_map = local_tr(b)
    MeshCat.settransform!(v, affine_map)
end

#-----------------------------------------------------------------
# MeshCat.Mesh
#-----------------------------------------------------------------
mutable struct Mesh <: AbstractSceneObject
    _props::Dict{Symbol, Any}

    function Mesh(; name="defaultMesh", mesh::GeometryBasics.AbstractMesh, kwargs...)
        props = build_props(; mesh=mesh, name=name, kwargs...)

        return new(props)
    end
end

mesh(m::Mesh) = prop(m, :mesh)

function render!(m::Mesh)
    v = vis(m)
    mat = material(m)
    the_mesh = mesh(m)

    vert_colors = prop(m, :vert_colors)
    if (vert_colors !== nothing)
        the_mesh = MeshCat.meta(the_mesh, vertexColors=vert_colors)
        mat = MeshCat.MeshLambertMaterial(vertexColors=true)
        MeshCat.setobject!(v, the_mesh, mat)
    else
        if (mat !== nothing)
            MeshCat.setobject!(v, the_mesh, material(mat))
        else
            MeshCat.setobject!(v, the_mesh)
        end
    end

    MeshCat.settransform!(v, local_tr(m))
end

# create a primitive mesh: rect, box, sphere, cone
function primitive(type::Symbol = :sphere; name="primitive", kwargs...)
    mesh = nothing
    if (type == :box || type == :cube)
        box = GeometryBasics.HyperRectangle(GeometryBasics.Vec(-0.5, -0.5, -0.5), GeometryBasics.Vec(1.0, 1.0, 1.0))
        mesh = GeometryBasics.mesh(box, uv=GeometryBasics.Vec2f0)
    elseif (type == :rect || type == :rectangle)
        box = GeometryBasics.HyperRectangle(GeometryBasics.Vec(-0.5, -0.5, 0), GeometryBasics.Vec(1.0, 1.0, 0))
        mesh = GeometryBasics.mesh(box, uv=GeometryBasics.Vec2f0)
    else (type == :sphere)
        sphere = GeometryBasics.HyperSphere(GeometryBasics.Point(0.0, 0.0, 0.0), 0.5)
        mesh = GeometryBasics.mesh(sphere, uv=GeometryBasics.Vec2f0)
    end
    res = Mesh(; name=name, mesh=mesh, kwargs...)
    return res
end

#-----------------------------------------------------------------

#-----------------------------------------------------------------
# Arrow
#-----------------------------------------------------------------
mutable struct Arrow <: AbstractSceneObject
    _props::Dict{Symbol, Any}

    function Arrow(; name="defaultArrow", from::SVector{3, Float64}, to::SVector{3, Float64}, kwargs...)
        props = build_props(; from=from, to=to, name=name, kwargs...)

        return new(props)
    end
end

from(a::Arrow) = prop(a, :from)
to(a::Arrow) = prop(a, :to)

function render!(a::Arrow)
    v = vis(a)
    mat = material(a)
    if (mat === nothing)
        mat = Material()
    end

    f = from(a)
    t = to(a)
    vec = t - f
    point = GeometryBasics.Point(f...)
    vec = GeometryBasics.Vec((t - f)...)

    arrow = MeshCat.ArrowVisualizer(v)
    MeshCat.setobject!(arrow, material(mat))
    MeshCat.settransform!(arrow, point, vec)
end
#-----------------------------------------------------------------

#-----------------------------------------------------------------
# Axes
#-----------------------------------------------------------------
mutable struct Axes <: AbstractSceneObject
    _props::Dict{Symbol, Any}

    function Axes(;
        name="defaultAxes",
        x_axis::SVector{3, Float64} = unitX3(),
        y_axis::SVector{3, Float64} = unitY3(),
        z_axis::SVector{3, Float64} = unitZ3(),
        kwargs...
    )
        props = build_props(; x_axis=x_axis, y_axis=y_axis, z_axis=z_axis, name=name, kwargs...)

        return new(props)
    end
end

x_axis(a::Axes) = prop(a, :x_axis)
y_axis(a::Axes) = prop(a, :y_axis)
z_axis(a::Axes) = prop(a, :z_axis)

function render!(a::Axes)
    v = vis(a)
    mat = material(a)
    if (mat === nothing)
        mat = Material()
    end

    point = GeometryBasics.Point(0.0, 0.0, 0.0)

    shaft_scale = prop(a, :shaft_scale, 0.01)
    axes_scale = prop(a, :axes_scale, 1.0)

    color!(mat, Colors.RGBA(0.9, 0.1, 0.1, 1.0))
    vec = GeometryBasics.Vec(prop(a, :x_axis)) * axes_scale
    arrow = MeshCat.ArrowVisualizer(v["X"])
    MeshCat.setobject!(arrow, material(mat))
    MeshCat.settransform!(arrow, point, vec, shaft_radius = shaft_scale)

    color!(mat, Colors.RGBA(0.1, 0.9, 0.1, 1.0))
    vec = GeometryBasics.Vec(prop(a, :y_axis)) * axes_scale
    arrow = MeshCat.ArrowVisualizer(v["Y"])
    MeshCat.setobject!(arrow, material(mat))
    MeshCat.settransform!(arrow, point, vec, shaft_radius = shaft_scale)

    color!(mat, Colors.RGBA(0.1, 0.1, 0.9, 1.0))
    vec = GeometryBasics.Vec(prop(a, :z_axis)) * axes_scale
    arrow = MeshCat.ArrowVisualizer(v["Z"])
    MeshCat.setobject!(arrow, material(mat))
    MeshCat.settransform!(arrow, point, vec, shaft_radius = shaft_scale)

    MeshCat.settransform!(v, local_tr(a))
end
#-----------------------------------------------------------------

#-----------------------------------------------------------------
# Line Segments
#-----------------------------------------------------------------
mutable struct LineSegments <: AbstractSceneObject
    _props::Dict{Symbol, Any}

    function LineSegments(; name="defaultAxes", points::Vector{SVector{3, Float64}} = [], kwargs...)
        props = build_props(; points=points, name=name, kwargs...)

        return new(props)
    end
end

points(ls::LineSegments) = prop(ls, :points)

function render!(ls::LineSegments)
    v = vis(ls)
    mat = material(ls)

    pts = points(ls)
    pts2 = GeometryBasics.Point.(pts)

    if (mat !== nothing)
        MeshCat.setobject!(v, MeshCat.LineSegments(GeometryBasics.Point.(points(ls)), lines_material(mat)))
    else
        MeshCat.setobject!(v, MeshCat.LineSegments(GeometryBasics.Point.(points(ls))))
    end
    MeshCat.settransform!(v, local_tr(ls))
end

#-----------------------------------------------------------------
# Line Segments
#-----------------------------------------------------------------
mutable struct PointCloud <: AbstractSceneObject
    _props::Dict{Symbol, Any}

    function PointCloud(; name="defaultPointCloud", points::Vector{SVector{3, Float64}} = [], kwargs...)
        props = build_props(; points=points, name=name, kwargs...)

        return new(props)
    end
end

points(pc::PointCloud) = prop(pc, :points)

function render!(pc::PointCloud)
    v = vis(pc)
    mat = material(pc)

    pts = points(pc)
    pts2 = GeometryBasics.Point.(pts)

    if (mat !== nothing)
        # create a list with the same material
        colors = repeat([mat._color], length(points(pc)))
        MeshCat.setobject!(v, MeshCat.PointCloud(GeometryBasics.Point.(points(pc)), colors))
    else
        MeshCat.setobject!(v, MeshCat.PointCloud(GeometryBasics.Point.(points(pc))))
    end
    MeshCat.settransform!(v, local_tr(pc))
end

#-----------------------------------------------------------------


text:
using  Test
import RecursiveFactorization
import LinearAlgebra
using  LinearAlgebra: norm
using  Random

Random.seed!(12)

const baselu = LinearAlgebra.lu
const mylu = RecursiveFactorization.lu

function testlu(A, MF, BF)
    @test MF.info == BF.info
    @test norm(MF.L*MF.U - A[MF.p, :], Inf) < 100sqrt(eps(real(one(float(first(A))))))
    nothing
end

@testset "Test LU factorization" begin
    for _p in (true, false), T in (Float64, Float32, ComplexF64, ComplexF32, Real)
        p = Val(_p)
        for s in [1:10; 50:80:200; 300]
            siz = (s, s+2)
            @info("size: $(siz[1]) × $(siz[2]), T = $T, p = $_p")
            if isconcretetype(T)
                A = rand(T, siz...)
            else
                _A = rand(siz...)
                A  = Matrix{T}(undef, siz...)
                copyto!(A, _A)
            end
            MF = mylu(A, p)
            BF = baselu(A, p)
            testlu(A, MF, BF)
            i = rand(1:s) # test `MF.info`
            A[:, i] .= 0
            MF = mylu(A, p, check=false)
            BF = baselu(A, p, check=false)
            testlu(A, MF, BF)
        end
    end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 629c1314-fdb8-11ea-0810-73f40ec50097
using Pkg, DrWatson

# ╔═╡ 629c6026-fdb8-11ea-222c-01752dee3679
begin
	#@quickactivate "StatisticalRethinkingStan"
	using StanSample
	using StatisticalRethinking
end

# ╔═╡ 62a90e5c-fdb8-11ea-2748-a31802543587
for i in 5:7
	include(projectdir("models", "05", "m5.$(i)s.jl"))
end

# ╔═╡ 4ce0616e-fdb7-11ea-0474-63cb801384f5
md"## Clip-05-40.1s.jl"

# ╔═╡ 629ccafc-fdb8-11ea-1012-d11648fe79d0
md"### snippet 5.39"

# ╔═╡ ed7b94dc-8346-11eb-0644-69bf473d42e2
md"
!!! note
	Restart this notebook to re-execute below `include()`.
"

# ╔═╡ 62a988b4-fdb8-11ea-0cdd-1759b8567402
if success(rc5_5s)
	post5_5s_df = read_samples(m5_5s, :dataframe)
	title5 = "Kcal_per_g vs. neocortex_perc" * "\n89% predicted and mean range"
	fig1 = plotbounds(
		df, :neocortex_perc, :kcal_per_g,
		post5_5s_df, [:a, :bN, :sigma];
		title=title5
	)
end

# ╔═╡ 62b4472e-fdb8-11ea-3bf0-8b453a2872a9
if success(rc5_6s)
	post5_6s_df = read_samples(m5_6s, :dataframe)
	title6 = "Kcal_per_g vs. log mass" * "\n89% predicted and mean range"
	fig2 = plotbounds(
		df, :lmass, :kcal_per_g,
		post5_6s_df, [:a, :bM, :sigma];
		title=title6
	)
end

# ╔═╡ 62b4cd64-fdb8-11ea-38d9-8d6810b1a4d4
if success(rc5_7s)
	post5_7s_df = read_samples(m5_7s, :dataframe)
	title7 = "Counterfactual,\nholding M=0.0"
	fig3 = plotbounds(
		df, :neocortex_perc, :kcal_per_g,
		post5_7s_df, [:a, :bN, :sigma];
		title=title7
	)
end

# ╔═╡ 62be59ce-fdb8-11ea-2278-8732aa168e50
if success(rc5_7s)
	title8 = "Counterfactual,\nholding N=0.0"
	fig4 = plotbounds(
		df, :lmass, :kcal_per_g,
		post5_7s_df, [:a, :bM, :sigma];
		title=title8,
		xlab="log(mass)"
	)
end

# ╔═╡ 62beccba-fdb8-11ea-32e6-0fff022ac3fe
plot(fig1, fig2, fig3, fig4, layout=(2, 2))

# ╔═╡ 62c8c1fc-fdb8-11ea-32c1-b9d03576a104
md"## End of clip-05-40.1s.jl"

# ╔═╡ Cell order:
# ╟─4ce0616e-fdb7-11ea-0474-63cb801384f5
# ╠═629c1314-fdb8-11ea-0810-73f40ec50097
# ╠═629c6026-fdb8-11ea-222c-01752dee3679
# ╟─629ccafc-fdb8-11ea-1012-d11648fe79d0
# ╟─ed7b94dc-8346-11eb-0644-69bf473d42e2
# ╠═62a90e5c-fdb8-11ea-2748-a31802543587
# ╠═62a988b4-fdb8-11ea-0cdd-1759b8567402
# ╠═62b4472e-fdb8-11ea-3bf0-8b453a2872a9
# ╠═62b4cd64-fdb8-11ea-38d9-8d6810b1a4d4
# ╠═62be59ce-fdb8-11ea-2278-8732aa168e50
# ╠═62beccba-fdb8-11ea-32e6-0fff022ac3fe
# ╟─62c8c1fc-fdb8-11ea-32c1-b9d03576a104


text:
#__precompile__()
"Various functions regarding Dynamic Binary Networks and fitness models"
module DynNetsW
#export sample_bern_mat, sample_fitTV, sampleDgpNetsTv

using Distributions, StatsBase,Optim, LineSearches, StatsFuns, Roots
using StaticNets,Utilities
using ForwardDiff

abstract type SdErgm end
abstract type SdErgmW <: SdErgm end
abstract type SdErgmWcount <: SdErgmW end
abstract type SdErgm <: SdErgm end

#constants

#-----------------WEIGHTED DIRECTED NETWORKS------------

# Directed NEtworks
struct  SdErgmDirW1 <: SdErgmW
     """ A Logistic Gas model for Directed Weighted networks and probability depending only
             on time varying parameters.
             ''\p^{(t)}_{ij}  = logistic( \theta^{(t)}_i + \theta^{(t)}_j)''

             # Model Definition and parameters conversion functions
         """
     # node - specific W and group Specific (A B)
     obsT:: Array{<:Real,2} # In strengths degrees and out strenghts for each t . The first dimension is always time. Tx2N
     Par::Array{Array{<:Real,1},1} #[[InW;outW], Agroups , Bgroups]  Group specific parameters
     #Nodes in a group have the same GAS parameters, but different fitnesses

     # each node belongs to a W group, each W group is assigned 2 AB groups, in order
     #to allow different parameters for in and our. The AB groups In and Out can be the same

     groupsInds::Array{Array{Int,1},1} #[[Nodes_2_groupsW, Wgroups_2_groupsABI , Wgroup_2_groupsABO ]  (_2_ stands for "assignement to" )
     # In and out W parameters are assumed to be always different, while In and Out AB
     # can be equal. Moreover, I assume that nodes with the same Wgroup have the same ABgroups
     scoreScalingType::String # String that specifies the rescaling of the score. For a list of possible
     # choices see function scalingMatGas
  end

#inizializza senza specificare assenza di scaling
SdErgmDirW1( obsT , Par  , groupsInds) =
        SdErgmDirW1( obsT,Par, groupsInds,"FISHER-DIAG")
#Initialize by observations only
SdErgmDirW1(obsT:: Array{<:Real,2}) =   (N = round(length(obsT[1,:])/2);
                                                SdErgmDirW1(obsT,
                                                    [ zeros(2N) ,0.9  * ones(1), 0.01 * ones(1)],
                                                    [Int.(1:N),Int.(ones(N)),Int.(ones(N))],"FISHER-DIAG") )
fooSdErgmDirW1 = SdErgmDirW1(ones(100,20))

# Relations between Static and Dynamical models: conventions on storage for
# parameters and observations
StaModType(Model::SdErgmDirW1) = StaticNets.fooErgmDirW1# to be substituted with a conversion mechanism
# linSpacedPar(Model::SdErgmDirBin1,Nnodes::Int;NgroupsW = Nnodes, deltaN::Int=3,graphConstr = true) =
#         StaticNets.linSpacedPar(StaModType(Model),Nnodes;Ngroups = NgroupsW,deltaN=deltaN,graphConstr =graphConstr);
#
# options and conversions of parameters for optimization

function NumberOfGroups(Model::SdErgmDirW1,groupsInds::Array{Array{Int,1},1})
     NGW = length(unique(groupsInds[1]))
     GBA = length( unique( [groupsInds[2][.!(groupsInds[2].==0)] ; groupsInds[3][.!(groupsInds[3].==0)] ]) )
      return NGW, GBA
 end
function NumberOfGroupsAndABindNodes(Model::SdErgmDirW1,groupsInds::Array{Array{Int,1},1})
     NGW = length(unique(groupsInds[1]))
     GBA = length( unique( [groupsInds[2][.!(groupsInds[2].==0)] ; groupsInds[3][.!(groupsInds[3].==0)] ]) )
     ABgroupsIndNodesIO = [groupsInds[2][groupsInds[1]] ; groupsInds[3][groupsInds[1]] ]
     indTvNodesIO =   (.!(ABgroupsIndNodesIO .==0))#if at least one group ind is zero then both in and out tv par are constant
     return NGW, GBA, ABgroupsIndNodesIO, indTvNodesIO
 end
function vec_2_array_all_par(Model::SdErgmDirW1,VecGasPar::Array{<:Real,1})
     groupsInds = Model.groupsInds
     NGW, GBA = NumberOfGroups(Model,groupsInds)
     ArrayGasPar =  [VecGasPar[1:2NGW],ones(GBA).*VecGasPar[2NGW+1:2NGW+GBA],ones(GBA).*VecGasPar[2NGW+GBA+1:2NGW + 2GBA]]
     return ArrayGasPar
     end
function restrictGasPar(Model::SdErgmDirW1,vecUnGasPar::Array{<:Real,1})
     "From the Unrestricted values of the parameters return the restricted ones
     takes as inputs a vector of parameters that can take any value in R and returns
     the appropriately contrained values, e.g. the coefficient of autoregressive
     component in the GAS update rule has to be positive and <1."
     groupsInds = Model.groupsInds
     NGW, GBA= NumberOfGroups(Model,groupsInds)
     W_un = vecUnGasPar[1:2NGW]
     diag_B_Un = vecUnGasPar[2NGW+1:2NGW+GBA]
     diag_A_Un = vecUnGasPar[2NGW+ GBA+1:2NGW+2GBA]
     W_Re = W_un
     diag_A_Re = exp.(diag_A_Un)
     diag_B_Re = 1./(1+exp.(.-diag_B_Un))
     vecRePar =  [W_Re; diag_B_Re; diag_A_Re]
     return vecRePar
     end
function unRestrictGasPar( Model::SdErgmDirW1,vecReGasPar::Array{<:Real,1})
     "From the restricted values of the parameters return the unrestricted ones
     takes as inputs a vector of parameters that can take any value in R and returns
     the appropriately contrained values, i.e. the coefficient of autoregressive
     component in the GAS update rule has to be positive and <1."
     groupsInds = Model.groupsInds
     NGW, GBA= NumberOfGroups(Model,groupsInds)
     W_Re = vecReGasPar[1:2NGW]
     diag_B_Re = vecReGasPar[2NGW+1:2NGW+GBA]
     diag_A_Re = vecReGasPar[2NGW+GBA+1:end]

     W_Un = W_Re
     diag_B_Un = log.(diag_B_Re./(1-diag_B_Re ))
     diag_A_Un = log.(diag_A_Re)

     vecUnPar =  [W_Re;diag_B_Un; diag_A_Un]
     return vecUnPar
     end

function uBndPar2bndPar_T(Model::SdErgmDirW1, uBndParNodesIO_T::Array{<:Real,2}  )
    T,N = size(uBndParNodesIO_T)
    outBndPar_T = zeros(uBndParNodesIO_T)
    for t = 1:T
        uBndParNodesIO_t = uBndParNodesIO_T[t,:]
        outBndPar_T[t,:] = StaticNets.uBndPar2bndPar(StaModType(Model), uBndParNodesIO_t )
    end
    return outBndPar_T
end
#Gas Filter Functions

function scalingMatGas(Model::SdErgmDirW1,expMat::Array{<:Real,2},I_tm1::Array{<:Real,2})
    "Return the matrix required for the scaling of the score, given the expected
     matrix and the Scaling matrix at previous time. "
    if uppercase(Model.scoreScalingType) == "FISHER-EWMA"
         error()
        # λ = 0.8
        #
        # I = expMat.*(1-expMat)
        # diagI = sum(I,2)
        # [I[i,i] = diagI[i] for i=1:length(diagI) ]
        # I_t =  λ*I + (1-λ) *I_tm1
        # scaling = sqrtm(I) ##
    else
        error()
    end
    return scaling
 end

function scaleScoreGas(Model::SdErgmDirW1,  ftotIO_t::Array{<:Real,1},fBndTot_t::Array{<:Real,1},gradIO_t::Array{<:Real,1}, expfIO_t::Array{<:Real,1},exp_Mat_t::Array{<:Real,2},
                        indMin::Int64,  WgasIO::Array{<:Real,1},BgasIO::Array{<:Real,1},AgasIO::Array{<:Real,1},strIO_t::Array{<:Real,1},
                        exp_strIO_t::Array{<:Real,1},t::Int; I_tm1::Array{<:Real,2} = zeros(exp_Mat_t))
    "Return the rescaled score according to the scalingtype of the model"
    if uppercase(Model.scoreScalingType) == ""
        s_t = gradIO_t
    elseif uppercase(Model.scoreScalingType) == "FISHER-DIAG"
        #Rescale the score in order to set its variance to one
        tmpMat = (exp_Mat_t.^2)
        N = round(Int,length(gradIO_t)/2)
        diagI_t = ([sumSq(tmpMat,2) ; sumSq(tmpMat,1)])
        tmpDiagEl =   sum(diagI_t[1:N]) - diagI_t[N+indMin] #
        print_diagI_t = copy(diagI_t)
        if false # unhortodox scaling that uses the observed values
            tmpDiagEl =  ( (sum(diagI_t[1:N]) - diagI_t[N+indMin]) + sum(strIO_t[1:N]) )/2#
            diagI_t = (diagI_t + strIO_t.^2)./2
        end
        diagI_t[N + indMin] = tmpDiagEl
        diagI_t = diagI_t.*( expfIO_t.^2)
        rescVec = sqrt.(diagI_t)
        s_t =  gradIO_t./rescVec
        s_t[gradIO_t.==0] = 0

        if false
            println(t)
        printInd = 45
    display((ftotIO_t[1:N][printInd],fBndTot_t[1:N][printInd],strIO_t[1:N][printInd],exp_strIO_t[1:N][printInd],s_t[1:N][printInd]))
    display((ftotIO_t[N+1:end][printInd],fBndTot_t[N+1:end][printInd],strIO_t[N+1:end][printInd],exp_strIO_t[N+1:end][printInd],s_t[1+N:2N][printInd]))
    display((ftotIO_t[N+1:end][indMin],fBndTot_t[N+1:end][indMin],strIO_t[N+1:end][indMin],exp_strIO_t[N+1:end][indMin],s_t[1+N:end][indMin]))
        end
        if false# maximum(abs.(s_t[isfinite.(s_t)])) > 1e4

            indOprint = indmax(abs.(s_t[1+N:2N]))
            @show indOprint
             println("str")
             println(strIO_t[indOprint])
            # display([strIO_t[1:N] strIO_t[N+1:2N] ] )
             println("exp Str")
             println((indmax(exp_strIO_t[1:N]),exp_strIO_t[1:N][indmax(exp_strIO_t[1:N])] ,indmax(exp_strIO_t[N+1:2N]),exp_strIO_t[N+1:2N][indmax(exp_strIO_t[N+1:2N])]  ) )
             #display([exp_strIO_t[1:N] exp_strIO_t[N+1:2N] ] )
            println("fBndTot_t")
            println(fBndTot_t[N+1:end][indOprint])
            #display([fBndTot_t[1:N] fBndTot_t[N+1:2N] ] )
            println("exp_mat_t")
            println(exp_Mat_t[:,indOprint])
            #display(exp_Mat_t)
            println("grad")
             #display([gradIO_t[1:N] gradIO_t[N+1:2N] ] )
           println([maximum(abs.(gradIO_t[2:N])) indmax(abs.(gradIO_t[2:N])); maximum(abs.(gradIO_t[1+N:2N])) indmax(abs.(gradIO_t[1+N:2N]))])
             println(gradIO_t)
             println("diagI_t")
              #display([print_diagI_t[1:N] print_diagI_t[N+1:2N] ] )
              println(" expfIO_t.^2")
               #display([ (expfIO_t.^2)[1:N]  (expfIO_t.^2)[N+1:2N] ] )
              println("tmpDiagEl")
              #display(tmpDiagEl)
            println("resc vect")
            println(rescVec)
            #display([rescVec[1:N] rescVec[N+1:2N]])
           println("score")
           println(s_t[2:end])
           println([maximum(abs.(s_t[2:N])) indmax(abs.(s_t[2:N])); maximum(abs.(s_t[1+N:2N])) indmax(abs.(s_t[1+N:2N]))])
            # println(s_t[1:N])
            #println(s_t[N+1:2N]  )
            #display([s_t[1:N] s_t[N+1:2N] ] )
            println("A")
            # display([AgasIO[1:N] AgasIO[N+1:2N] ] )
             println("B")
             # display([BgasIO[1:N] BgasIO[N+1:2N] ] )
              println("W")
              # display([WgasIO[1:N] WgasIO[N+1:2N] ] )
            maximum(abs.(s_t[isfinite.(s_t)])) > 1e8?  error():()
        end

    elseif uppercase(Model.scoreScalingType) == "FISHER-EWMA"
         error()
        # λ = 0.8
        #
        # I = scalingMatGas(Model,grad_t,exp_mat_t,I_tm1)
        # diagI = sum(I,2)
        # [I[i,i] = diagI[i] for i=1:length(diagI) ]
        # I_t =  λ*I + (1-λ) *I_tm1
        # scalingMat = sqrtm(I) ##
    end

    return s_t, I_tm1
 end
function predict_score_driven_par( Model::SdErgmDirW1,strIO_t::Array{<:Real,1},
                         ftotIO_t::Array{<:Real,1},I_tm1::Array{<:Real,2},indTvNodesIO::BitArray{1},
                         WgasIO::Array{<:Real,1},BgasIO::Array{<:Real,1},AgasIO::Array{<:Real,1},t::Int ;
                         fBndTot_t::Array{<:Real,1} =  StaticNets.uBndPar2bndPar(StaModType(Model),ftotIO_t) )
     #= likelihood and gradients depend on all the parameters (ftot_t), but
     only the time vaying ones (f_t) are to be updated=#
     fIO_t = ftotIO_t[indTvNodesIO] #Time varying fitnesses
     N,ftotI_t,ftotO_t = splitVec(ftotIO_t)
     indTvNodesI = indTvNodesIO[1:N]

     #find the minimum out parameter
     indMin = indmin(ftotO_t)
     # indMin2 = indmax(strIO_t[N+1:end])
     #
     indTvNodesIO[N+1:end][indMin] ? (scoreAdd = true; indMin = indmin(ftotO_t[indTvNodesIO[N+1:end]]) ) : scoreAdd = false
     #display((BgasIO[1],AgasIO[1]))
     sumParMat, exp_mat_t = StaticNets.expMatrix2(StaModType(Model),fBndTot_t )

     exp_strIO_t = [sumSq(exp_mat_t,2);sumSq(exp_mat_t,1)]
     #The following part is a faster version of logLikelihood_t
     loglike_t =   ((sum(log.(sumParMat)) - sum(fBndTot_t.*strIO_t) ))

     # compute the score of observations at time t wrt parameters at tm1
     expfIO_t = exp.(fIO_t)

     strDiffIO_t =   (exp_strIO_t[indTvNodesIO] -strIO_t[indTvNodesIO])

     gradIO_t = strDiffIO_t.*expfIO_t # only for TV paramaters
    # scoreAdd ?  gradIO_t[N + indMin] +=  -expfIO_t[N + indMin] * sum(strDiffIO_t[1:N]) :()

     sIO_t,I_tm1 =  scaleScoreGas(Model,  ftotIO_t,fBndTot_t,gradIO_t, expfIO_t,exp_mat_t,indMin, WgasIO,BgasIO,AgasIO,strIO_t,exp_strIO_t,t)
     sIO_t[1] = 0 #identification of the update equation
     # sembra che cosi facendo la strength In del primo nodo non contriubisca all'update
     # ma non \'e' cosi. Infatti l'info sul primo in degree è ridondante se considerata
     # insieme a tutti gli altri in degrees e tutti gli out degrees

     if false# indMin_tp1 == indMin
         #update the out parameters
         fO_tp1 = WgasIO[N+1:2N] + BgasIO[N+1:2N].* fIO_t[N+1:2N] + AgasIO[N+1:2N].*sIO_t[N+1:2N]
         #if the minimum out par has changed than transform all the in (static and dynamic)
         indMin_tp1 = indmin(fO_tp1)
         @show (indMin,indMin_tp1)
        # ftotIO_t[1:N] = log.(exp.(ftotI_t) + exp(fO_tp1[indMin_tp1]) - exp(ftotO_t[indMin]) )
        fI_t = ftotIO_t[1:N][indTvNodesI]
        # update the in parameters
        fI_tp1 = WgasIO[1:N] + BgasIO[1:N].* fI_t + AgasIO[1:N].*sIO_t[1:N]
        fIO_tp1 = [fI_tp1 fO_tp1]
    elseif false # unhortodox modification that limits exploding scores
        score_thr_pos = 1
        score_thr_neg = -1
        shockod = AgasIO.*sIO_t
        shockod[(shockod).>score_thr_pos] = score_thr_pos.*( 1+ log.(abs.(shockod[(shockod).>score_thr_pos])))
        shockod[(shockod).<score_thr_neg] = score_thr_neg.*( 1+ log.(abs.(shockod[(shockod).<score_thr_neg])))
        #s_t[abs.(s_t).>score_thr] = sign.(s_t[abs.(s_t).>score_thr]).*score_thr.*( 1+ log.(abs.(s_t[abs.(s_t).>score_thr])))
        fIO_tp1 = WgasIO + BgasIO.* fIO_t + shockod
    elseif true # unhortodox modification that limits exploding scores
        score_thr_pos = 1000
        score_thr_neg = -35
        shockod = sIO_t
        shockod[(shockod).>score_thr_pos] = score_thr_pos.*( 1+ log10.(abs.(shockod[(shockod).>score_thr_pos])))
        shockod[(shockod).<score_thr_neg] = score_thr_neg.*( 1+ log10.(abs.(shockod[(shockod).<score_thr_neg])))
        #s_t[abs.(s_t).>score_thr] = sign.(s_t[abs.(s_t).>score_thr]).*score_thr.*( 1+ log.(abs.(s_t[abs.(s_t).>score_thr])))
        fIO_tp1 = WgasIO + BgasIO.* fIO_t + AgasIO.*shockod
    else
        fIO_tp1 = WgasIO + BgasIO.* fIO_t + AgasIO.*sIO_t

    end

     ftotIO_t[indTvNodesIO] = fIO_tp1 #of all parameters udate dynamic ones with GAS
     ftotIO_tp1 = ftotIO_t
     return ftotIO_tp1,loglike_t,I_tm1
  end
function score_driven_filter_or_dgp( Model::SdErgmDirW1,
                                vResGasPar::Array{<:Real,1} ;
                                obsT::Array{<:Real,2}=Model.obsT,
                                groupsInds:: Array{Array{Int,1},1}=Model.groupsInds,
                                dgpTN = (0,0))
    """GAS Filter the Dynamic Fitnesses from the Observed degrees, given the GAS parameters
     given T observations for the degrees in TxN vector degsT
     """
     #se parametri per dgp sono dati allora campiona il dgp, altrimenti filtra
    sum(dgpTN) == 0 ? dgp = false : dgp = true
    if dgp
        T,N = dgpTN
     Y_T = zeros(Float64,T,N,N)
    else
     T,N2 = size(obsT);N = round(Int,N2/2)
    end

    NGW,GBA,ABgroupsIndNodesIO,indTvNodesIO = NumberOfGroupsAndABindNodes(Model,groupsInds)
    # ABgroups Inds for each par (In and Out). IO means that we have a 2N vector
    # of indices: one index for each In par and below one index for each Out par

    # Organize parameters of the GAS update equation
    WGroupsI = vResGasPar[1:NGW]
    WGroupsO = vResGasPar[NGW+1:2NGW]
    W_allIO = [WGroupsI[groupsInds[1]] ; WGroupsO[groupsInds[1]] ]

    BgasGroups  = vResGasPar[2NGW+1:2NGW+GBA]
    AgasGroups  = vResGasPar[2NGW+GBA+1:2NGW+2GBA]
    #distribute nodes among groups with equal parameters
    AgasIO = AgasGroups[ABgroupsIndNodesIO[indTvNodesIO]]
    BgasIO = BgasGroups[ABgroupsIndNodesIO[indTvNodesIO]]
    WgasIO   = W_allIO[indTvNodesIO]# identify!(Model,W_allIO[indTvNodesIO])

    # start values equal the unconditional mean, but  constant ones remain equal to the unconditional mean, hence initialize as:
    UMallNodesIO = W_allIO
    UMallNodesIO[indTvNodesIO] =  WgasIO./(1-BgasIO)
    if false
        println((vResGasPar[2NGW+1:2NGW+2GBA]))
        println("W")
        display([WgasIO[1:N] WgasIO[N+1:2N] ] )
        println("Unc Means")
         display([UMallNodesIO[1:N] UMallNodesIO[N+1:2N] ] )
     end
    fIOVecT = ones(Real,T,2N)

    ftotIO_t = UMallNodesIO# identify!(Model,UMallNodesIO)
    I_tm1 = eye(Float64,2N)
    loglike = 0

    for t=1:T
        #println(t)
        #display(ftotIO_t)
        fBndTot_t =  StaticNets.uBndPar2bndPar(StaModType(Model),ftotIO_t)

        if dgp
            expMat_t = StaticNets.expMatrix(StaModType(Model),fBndTot_t )# exp.(ftot_t)
            Y_T[t,:,:] = StaticNets.samplSingMatCan(StaModType(Model),expMat_t)
            #display(expMat_t)
            strIO_t =  [ sumSq(Y_T[t,:,:],2); sumSq(Y_T[t,:,:],1)]
        else
            strIO_t = obsT[t,:]  # vector of in and out degrees
        end
        # println("str_t")
        #  display([strIO_t[1:N] strIO_t[N+1:2N] ] )
        #140<=t<=150? display(fBndTot_t'):()
        ftotIO_t,loglike_t,I_tm1 = predict_score_driven_par(Model,strIO_t,ftotIO_t,I_tm1,indTvNodesIO,WgasIO,BgasIO,AgasIO,t ; fBndTot_t = fBndTot_t)
        fIOVecT[t,:] = ftotIO_t #store the filtered parameters from previous iteration
        loglike += loglike_t

    end
    if dgp
        return Y_T , fIOVecT
    else
        return fIOVecT, loglike
    end
    end
score_driven_filter_or_dgp(Model::SdErgmDirW1) = score_driven_filter_or_dgp(Model,array_2_vec_all_par(Model,Model.Par))
sampl(Mod::SdErgmDirW1,T::Int)=( N = length(Mod.groupsInds[1]) ;
                                        tmpTuple = score_driven_filter_or_dgp(Mod,[Mod.Par[1];Mod.Par[2];Mod.Par[3]];  dgpTN = (T,N) );
                                        degsIO_T = [sumSq(tmpTuple[1],3) sumSq(tmpTuple[1],2)];
                                        (SdErgmDirBin1(degsIO_T,Mod.Par,Mod.groupsInds,Mod.scoreScalingType),tmpTuple[1],tmpTuple[2]) )

# # Estimation
#
function estSingSnap(Model::SdErgmDirW1, str_t::Array{<:Real,1}; groupsInds = Model.groupsInds, targetErr::Real=targetErrValDynNets)
    hatBndPar,~,~ = StaticNets.estimate(StaModType(Model); strIO = str_t ,groupsInds = groupsInds[1], targetErr =  targetErr)
    hatUbndPar = StaticNets.bndPar2uBndPar(StaModType(Model),hatBndPar)
    return hatUbndPar
 end
function estimateSnapSeq(Model::SdErgmDirW1; strT::Array{<:Real,2}=Model.obsT,
                     groupsInds = Model.groupsInds,targetErr::Real=targetErrValDynNets,print_t = false)
    T,N2 = size(Model.obsT);N = round(Int,N2/2)
    NGW = length(unique(groupsInds[1]))
    UnParT = zeros(T,2NGW)
    prog = 0
    for t = 1:T
        print_t?println(t):()
        str_t = strT[t,:]
        UnParT[t,:] = estSingSnap(Model,str_t; groupsInds = groupsInds ,targetErr =  targetErr)
        round(t/T,2)>prog ? (prog=round(t/T,2);println(prog) ):()
    end
    return UnParT
 end

function estimate(Model::SdErgmDirW1;start_values = [zeros(Real,10),zeros(Real,3),zeros(Real,3) ],
                    targeting = false, meanEstSS::Array{Float64,1} = zeros(1) )
    T,N2 = size(Model.obsT);N = round(Int,N2/2)
    groupsInds = Model.groupsInds
    NGW,GBA,ABgroupsIndNodesIO,indTvNodesIO = NumberOfGroupsAndABindNodes(Model,groupsInds)
    optims_opt, algo = setOptionsOptim(Model)

    #set the starting points for the optimizations if not given
    if sum([sum(start_values[i]) for i =1:3 ])==0

         B0_ReGroups =  0.8 * ones(GBA)
         A0_ReGroups = 0.01 * ones(GBA)

         BgasNodes_0   = B0_ReGroups[ABgroupsIndNodesIO[indTvNodesIO]]
         AgasNodes_0   = A0_ReGroups[ABgroupsIndNodesIO[indTvNodesIO]]
         # choose the starting points for W in order to fix the Unconditional Means
         # equal to the static estimates for the average degrees(over time)
         MeanStrT =  meanSq(Model.obsT,1)
         if targeting
             if sum(meanEstSS) == 0
                 estSS =  DynNets.estimateSnapSeq(Model)
                 meanEstSS = meanSq(estSS,1)
             end
             UMstaticGroups = meanEstSS
         else
             error("This model doe not work well on real data without a targeting")
             UMstaticGroups = StaticNets.bndPar2uBndPar(StaModType(Model), estSingSnap(Model ,MeanStrT; groupsInds = groupsInds ))
         end
        # display([UMstaticGroups[1:N] UMstaticGroups[N+1:2N] ])
         conv_flag1 = true

         GBA>0?( W0_Groups = UMstaticGroups.*(1-B0_ReGroups[[groupsInds[2];groupsInds[3]]])):error()#W0_Groups = UMstaticGroups
         #W0_Groups = W_prob

    else
         W0_Groups = start_values[1]
         B0_ReGroups = start_values[2]
         A0_ReGroups = start_values[3]
         conv_flag1 = true
    end
    if GBA == 0  return [UMstaticGroups, zeros(1), zeros(1)] ,  conv_flag1 end;
    # Estimate  jointly all  Gas parameters
    vGasParAll0_Un = unRestrictGasPar(Model,[W0_Groups ;B0_ReGroups;A0_ReGroups ] )[2:end] # vGasParGroups0_Un#
    # objective function for the optimization
    function objfunGas(vparUn::Array{<:Real,1})# a function of the groups parameters
        #vParRes = restrictGasPar(Model,[-Inf;vparUn])
        GBA,Bun,Aun = splitVec(vparUn) #ottimizza solo A e B
        Bre = 1./(1+exp.(-Bun))
        Are = exp.(Aun)
        W_new = UMstaticGroups.*(1-Bre)
        vParRes = [W_new;Bre;Are]
        #display((vParRes[2],vParRes[end-1],vParRes[end]))
            foo,loglikelValue = score_driven_filter_or_dgp( Model,vParRes)

        return  -loglikelValue
    end

    #Run the optimization without AD if a matrix sqrt is required in the filter
    if uppercase(Model.scoreScalingType) == "FISHER-EWMA"
        ADobjfunGas = objfunGas
    else
        ADobjfunGas = TwiceDifferentiable(objfunGas, vGasParAll0_Un[end-1:end]; autodiff = :forward);
    end
    display(vGasParAll0_Un)

    display(objfunGas(vGasParAll0_Un[end-1:end]))

    #display(ForwardDiff.gradient(objfunGas, vGasParAll0_Un[end-1:end]))
    #error()
    optim_out2  = optimize(objfunGas,vGasParAll0_Un[end-1:end],algo,optims_opt)
    BAestUn = Optim.minimizer(optim_out2)
    BAestRe = [1./(1+exp.(-BAestUn[1:GBA]));exp.(BAestUn[1+GBA:end])]
    arrayGasParHat_Re  = [UMstaticGroups.*(1 - BAestRe[1:GBA]) , BAestRe[1:GBA],BAestRe[GBA+1:end]]
    conv_flag = conv_flag1*Optim.converged(optim_out2)
    println(optim_out2)
    return  arrayGasParHat_Re, conv_flag,Optim.minimum(optim_out2)
    end

end


text:
export find_clusters;
export KModel;

struct KModel{D,T <: Number}
    k::Int
    dims::Int
    centroids::Vector{NTuple{D,T}}
    data::Vector{NTuple{D,T}}
    assignm::Vector{Int}

    function KModel{D,T}(centroids, data, assignm) where T <: Number where D
        @assert eltype(eltype(centroids)) == eltype(eltype(centroids)) == T
        return new(length(centroids),D,totuples(centroids),totuples(data),assignm);
    end
    function KModel{D,T}(centroids::Vector{NTuple{D,T}},
        data::Vector{NTuple{D,T}},
        assignm::Vector{Int}) where T <: Number where D
        return new(length(centroids),D,centroids,data,assignm);
    end
    function KModel{D,T}(k, data) where T <: Number where D
        data = totuples(data);
        dc = center_coords(data) .* 2;
        centroids = [dc .* c for c in rtuples(Float64,D,k)];
        new(k,D,centroids, data, assign(data,centroids))
    end
end
function KModel(k, data)
    @assert length(data) > 0
    T = eltype(eltype(data));
    D = length(first(data));
    return KModel{D,T}(k,data)
end

data(km) = km.data;
data(km,c) = km.data[findall(isequal(c),km.assignm)];
ave_distance(km) = sum([ave_distance(km,c) for c in 1:km.k])/km.k;
ave_distance(km::KModel, c) = ave_distance(km.centroids[c],data(km,c));
non_empty_centroids(km::KModel) = findall((c) -> length(data(km,c)) > 0, 1:km.k);
empty_centroids(km::KModel) = findall((c) -> length(data(km,c)) == 0, 1:km.k);
KModel(km::KModel) = KModel{km.dims,eltype(eltype(km.data))}(copy(km.centroids),copy(km.data),copy(km.assignm));

function update_kmodel!(km , tol::Float64)
    alltolr = true;
    update_assigment!(km);

    for c in 1:km.k

        cd = data(km,c);
        if length(cd) == 0
            ccoords = center_coords(km.centroids)
        else
            ccoords = center_coords(cd);
        end
        if distance(ccoords, km.centroids[c]) > tol alltolr = false; end
        km.centroids[c] = ccoords;

    end

    return alltolr;
end
update_kmodel!(km; tol::Float64 = 1e-4) = update_kmodel!(km , tol);
function update_kmodel_many_times!(km , tol, itrs)
    for i in 1:itrs
        update_kmodel!(km, tol) && return true;
        i == itrs && return false;
    end
end
update_kmodel_many_times!(km; tol = 1e-4, itrs = 100) = update_kmodel_many_times!(km , tol, itrs)
update_assigment!(km) = km.assignm .= tuple(assign(km.data, km.centroids)...);
update_k!(km) = km.k = length(km.centroids);

function find_clusters(k, data; itrs = 100, uptol = 1e-5, upitrs = 100)
    @assert length(data) > 0
    D = length(first(data));
    T = eltype(first(data))

    bestkm = KModel{D,T}(k,data);
    update_kmodel_many_times!(bestkm, uptol, upitrs);
    for i in 1:itrs
        newkm = KModel{D,T}(k,data);
        update_kmodel_many_times!(newkm, uptol, upitrs);
        if ave_distance(bestkm) > ave_distance(newkm)
            bestkm = newkm;
        end
    end
    return bestkm;
end

function get_clean(km)
    nec = non_empty_centroids(km);
    if length(nec) == 0
        return KModel(km);
    else
        cd = km.centroids[nec];
        newkm = KModel{km.dims, eltype(eltype(km.data))}(cd, copy(km.data), copy(km.assignm));
        update_assigment!(newkm);
        return newkm;
    end
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# tests the output of the embedding example is correct
using Test

if Sys.iswindows()
    # libjulia needs to be in the same directory as the embedding executable or in path
    ENV["PATH"] = string(Sys.BINDIR, ";", ENV["PATH"])
end

function checknum(s, rx, cond)
    m = match(rx, s)
    if m === nothing
        return false
    else
        num = m[1]
        return cond(parse(UInt, num))
    end
end

@test length(ARGS) == 1
@testset "gcext example" begin
    out = Pipe()
    err = Pipe()
    p = run(pipeline(Cmd(ARGS), stdin=devnull, stdout=out, stderr=err), wait=false)
    close(out.in)
    close(err.in)
    out_task = @async readlines(out)
    err_task = @async readlines(err)
    # @test success(p)
    errlines = fetch(err_task)
    lines = fetch(out_task)
    @test length(errlines) == 0
    @test length(lines) == 6
    @test checknum(lines[2], r"([0-9]+) full collections", n -> n >= 10)
    @test checknum(lines[3], r"([0-9]+) partial collections", n -> n > 0)
    @test checknum(lines[4], r"([0-9]+) object sweeps", n -> n > 0)
    @test checknum(lines[5], r"([0-9]+) internal object scan failures",
        n -> n == 0)
    @test checknum(lines[6], r"([0-9]+) corrupted auxiliary roots",
        n -> n == 0)
end

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia


text:
module AdventOfCode2020
# Write your package code here.
export Day
export read_infile, line_parser, solve

abstract type Day end

function solve(::Type{T}, path::String) where T <: Day
    path |> p->read_infile(T,p) |> lines -> solve(T, lines)
end
solve(::Type{T}) where T <: Day = solve(T, "$(@__DIR__)/assets/day$(String(Symbol(T))[4:end-2]).txt")

line_parser(::Type{<:Day}, line) = line |> strip

function read_infile(::Type{T}, path)::Vector where T <: Day
    open(path,"r") do datafile
        [line_parser(T, line) for line in eachline(datafile)]
    end
end

function split_on_empty_line(lines)
    if !("" == lines[end])
        push!(lines, "")
    end
    res = Vector{Vector{String}}()
    last_line = 1
    for i in 1:length(lines)
        if lines[i] == ""

            push!(res, lines[last_line:i-1])
            last_line = i + 1
        end
    end
    res
end

include("day1.jl")
include("day2.jl")
include("day3.jl")
include("day4.jl")
include("day5.jl")
include("day6.jl")
include("day7.jl")
include("day8.jl")
include("day9.jl")

end


text:
type StateResult
    starttime::DateTime
    endtime::DateTime
    startstate::Vector
    result::Vector
    completed::Bool
    # function StateResult(starttime::DateTime,state::Vector)
    #     new(starttime,starttime,state,Float64[])
    # end
end

function processPairs(stateResults::Vector{StateResult}, e::MktEvent)

    completed = StateResult[]
    i = 0
    while length(stateResults) - i > 0
        if e.dt >= stateResults[i+1].endtime # time's up!
            push!(completed,  shift!(stateResults)) # extract the state/result pair
            completed[end].completed =  true
        elseif e.isExec == 1 # is execution not rate update
            # if X milliseconds (tm) of results have not been accumulated,
            # update result
            stateResults[i+1].result = [round((e.bid + e.ask)/2,3)]
        end
        i += 1
    end
    return completed
end

max_stars_repo_name: Andy-P/PriceDist

max_issues_repo_name: Andy-P/PriceDist

max_forks_repo_name: Andy-P/PriceDist


text:
# This file was generated, do not modify it. # hide
lgb = LGBMRegressor() #initialised a model with default params
lgbm = machine(lgb, features[train, :], targets[train, 1])
learning_range = range(lgb, :learning_rate, lower=1e-3, upper=1, scale=:log)
curve = learning_curve!(lgbm, resampling=CV(nfolds=5),
                        range=learning_range, resolution=100,
                        measure=rms)

figure(figsize=(8,6))
plot(curve.parameter_values, curve.measurements)
xscale("log")
xlabel("Learning rate (log scale)", fontsize=14)
ylabel("RMSE", fontsize=14)

plt.savefig(joinpath(@OUTPUT, "lgbm_hp2.svg")) # hide


text:
module  TestDataExploration

using PopGen
using DataFrames
using Test

cats = @nancycats;

@testset "missing()" begin
    @test typeof(missing_data(cats)) == DataFrame
    @test typeof(missing_data(cats, by = "population")) == DataFrame
    @test typeof(missing_data(cats, by = "locus")) == DataFrame
    @test typeof(missing_data(cats, by = "full")) == DataFrame
end

@testset "pairwise_identical" begin
    pw_i = pairwise_identical(cats)
    @test typeof(pw_i) == DataFrame
    @test size(pw_i) == (27966, 4)
    pw_i_2 = pairwise_identical(cats, samples(cats)[1:10])
    @test typeof(pw_i_2) == DataFrame
    @test size(pw_i_2) == (45, 4)
end

@testset "frequency tables" begin
    @test typeof(geno_freqtable(cats)) == DataFrame
    @test typeof(geno_freqtable(cats, by = "population")) == DataFrame
    @test typeof(allele_freqtable(cats)) == DataFrame
    @test typeof(allele_freqtable(cats, by = "population")) == DataFrame
end

end # module


text:
using FluxExtensions, Test, Flux

# we want to verify that threaded and non-threaded versions of
# learning returns the same function value and gradient
@testset "support for threaded learning" begin
	datas = [randn(3, 10) for i in 1:3];
	model = Dense(3, 1);
	loss(model, data) = sum(model(data));
	f1 = 0.0
	for i in 1:length(datas)
		f = loss(model, datas[i])
		@show f
		f1 += Flux.data(f)
		Flux.back!(f)
	end
	f1 /= length(datas)
	FluxExtensions.scalegrads!(model, 1/length(datas));

	models = [deepcopy(model) for _ in 1:length(datas)];
	f2 = FluxExtensions.∇loss!(loss, models, datas, nothing);
	@test f1 ≈ f2
	g1 = map( p -> Flux.Tracker.grad(p), params(model))
	g2 = map( p -> Flux.Tracker.grad(p), params(models[1]))
	@test g1 ≈ g2
end

max_stars_repo_name: pevnak/FluxExtensions

max_issues_repo_name: pevnak/FluxExtensions

max_forks_repo_name: pevnak/FluxExtensions


text:
"""
    mutable struct MCMCSchedule
        num_mcmc_steps::Int64
        num_updates::Int64
        start::NamedTuple{(:mcmciter, :pidx),Tuple{Int64,Int64}}
        exclude_updates::DefaultDict{Int64,OrdinalRange{Int64,Int64},UnitRange{Int64}}
    end

An object used for iterating over the steps of MCMC samplers. `num_mcmc_steps`
is the total number of MCMC steps, whereas `num_updates` is a total number of
separate updates that can be performed in a single MCMC step. Not all of those
updates have to be performed and the pattern of excluding the updates is stored
in `exclude_updates`. The `num_updates` and `exclude_updates` can change in the
midst of iterating through the `MCMCSchedule` and the iterator will register
that and act accordingly.
"""
mutable struct MCMCSchedule{BE,S,EI}
    num_mcmc_steps::Int64
    num_updates::Int64
    start::S
    exclude_updates::DefaultDict{Int64,OrdinalRange{Int64,Int64},UnitRange{Int64}}
    extra_info::EI

    function MCMCSchedule(
            num_mcmc_steps,
            num_updates,
            exclude_updates=[];
            start::S=(prev_mcmciter=nothing, prev_pidx=nothing, mcmciter=1, pidx=1),
            backend::BE=GenericMCMCBackend(),
            extra_info::EI=nothing
        ) where {S,BE,EI}
        tuples_excluded_updates = DefaultDict{Int64,OrdinalRange{Int64,Int64}}(0:0)
        for exclude_update in exclude_updates
            for idx in exclude_update[1]
                tuples_excluded_updates[idx] = exclude_update[2]
            end
        end
        new{BE,S,EI}(
            num_mcmc_steps,
            num_updates,
            start,
            tuples_excluded_updates,
            extra_info
        )
    end
end

"""
    Base.iterate(iter::MCMCSchedule, state=iter.start)

Iterate through the `MCMCSchedule`, outputting the named tuples
`(mcmciter=..., pidx=...)` along the way, indicating the current index of the
mcmc sampler (`mcmciter`), as well as the index of update that is to be
performed.
"""
function Base.iterate(iter::MCMCSchedule, state=iter.start)
    state.mcmciter > iter.num_mcmc_steps && return nothing
    new_state_tmp = (
        prev_mcmciter = state.mcmciter,
        prev_pidx = state.pidx,
        mcmciter = state.mcmciter,
        pidx = state.pidx,
    )
    new_state = transition(iter, new_state_tmp)
    return (state, extra_transitions(iter, new_state))
end

extra_transitions(iter::MCMCSchedule, new_state) = new_state

"""
    transition(schedule::MCMCSchedule, state)

Determine the next state of the `MCMCSchedule` iterator, by skipping through
all updates that are to be excluded, as per `schedule.exclude_updates`.
"""
function transition(schedule::MCMCSchedule, state)
    pidx_to_reset = (state.pidx == schedule.num_updates)
    new_state = (
        prev_mcmciter = state.prev_mcmciter,
        prev_pidx = state.prev_pidx,
        mcmciter = state.mcmciter + pidx_to_reset,
        pidx = ( pidx_to_reset ? 1 : state.pidx + 1 ),
    )
    iterations_to_skip = schedule.exclude_updates[new_state.pidx]
    skip_state = new_state.mcmciter in iterations_to_skip
    skip_state && return transition(schedule, new_state)
    new_state
end

"""
    reschedule!(
        schedule::MCMCSchedule,
        num_new_updates=0,
        idxes_to_remove=[],
        idxes_to_add=[]
    )

Make changes to the `schedule` (possibly in the midst of iterating through it).
Add `num_new_updtes`-many allowable updates that can be performed at each MCMC
step, `idxes_to_remove` lists all parameter udpates that are to be completely
removed from the MCMC sampler, `idxes_to_add` lists entries that are to be
added to `schedule.exclude_updates`.
"""
function reschedule!(
        schedule::MCMCSchedule,
        num_new_updates=0,
        idxes_to_remove=[],
        idxes_to_add=[]
    )
    schedule.num_updates += num_new_updates
    for idx in idxes_to_remove
        schedule.exclude_updates[idx] = 1:schedule.num_mcmc_steps
    end
    for idx_to_add in idxes_to_add
        schedule.exclude_updates[idx_to_add[1]] = idx_to_add[2]
    end
end


text:
# generate the table found in
# https://arxiv.org/abs/2201.05967

using DyadicKDE

# specify parameters
n_data = 500
n_evals = 50
n_repeats = 2000
degeneracies = ["total", "partial", "none"]
kernel_names = ["epanechnikov_order_2", "epanechnikov_order_4"]
evals = collect(range(-2.0, stop=2.0, length=n_evals))
sdp_solver = "cosmo"
n_resample = 10000
significance_level = 0.05
results = []
ps = Dict(
    "total" => [0.5, 0.0, 0.5],
    "partial" => [0.25, 0.0, 0.75],
    "none" => [0.2, 0.2, 0.6],
)

println("Running experiments")
for degen in degeneracies
    for kernel_name in kernel_names

        global counter = 0
        Threads.@threads for rep in 1:n_repeats

            global counter += 1
            println("Degeneracy: $degen, Kernel: $kernel_name, Repeat $counter/$n_repeats")
            data = make_data(n_data, ps[degen])
            h_ROT = estimate_ROT_bandwidth(data, "epanechnikov_order_2")

            est = DyadicKernelDensityEstimator(
                kernel_name, h_ROT, significance_level,
                n_resample, sdp_solver, evals, data,
                Dict("degen" => degen, "p" => ps[degen]))

            fit(est)
            push!(results, est)

        end
    end
end

println("\nResults:")
for degen in degeneracies
    for kernel_name in kernel_names

        results_reduced = [est for est in results if
                               est.meta["degen"] == degen
                               && est.kernel_name == kernel_name]

        p = ps[degen]
        h_ROT = mean([est.bandwidth for est in results_reduced])
        RIMSE = mean([get_RIMSE(est) for est in results_reduced])
        ucb_coverage = mean([get_ucb_coverage(est) for est in results_reduced])
        pci_coverage = mean([get_pci_coverage(est) for est in results_reduced])
        ucb_width = mean([get_ucb_average_width(est) for est in results_reduced])
        pci_width = mean([get_pci_average_width(est) for est in results_reduced])

        println("π = $p")
        println("Degeneracy: $degen")
        println("h_ROT = $h_ROT")
        println("Kernel: $kernel_name")
        println("RIMSE: $RIMSE")
        println("UCB coverage: $ucb_coverage")
        println("UCB width: $ucb_width")
        println("PCI coverage: $pci_coverage")
        println("PCI width: $pci_width")
        println()

    end
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    VariablePartition(var)

A method for partitioning spatial data into subsets of
constant value for variable `var`.
"""
struct VariablePartition <: PartitionMethod
  var::Symbol
end

function partition(sdata, method::VariablePartition)
  var = method.var

  @assert var ∈ name.(variables(sdata)) "invalid variable name"

  # data for variable
  vdata = sdata[var]

  # partition function with missings
  function f(i, j)
    vi, vj = vdata[i], vdata[j]
    mi, mj = ismissing(vi), ismissing(vj)
    (mi && mj) || ((!mi && !mj) && (vi == vj))
  end

  # partition function without missings
  g(i, j) = vdata[i] == vdata[j]

  # select the appropriate predicate function
  pred = Missing <: eltype(vdata) ? f : g

  # perform partition
  p = partition(sdata, PredicatePartition(pred))

  # retrieve value from each subset
  vals = [vdata[s[1]] for s in subsets(p)]

  # save metadata
  metadata = Dict(:values => vals)

  SpatialPartition(sdata, subsets(p), metadata)
end


text:
module CuNESSie

using Adapt
using CUDA
using IterativeSolvers
using LinearAlgebra
using NESSie
using NESSie: yukawa
using NESSie.BEM: BEMResult, LocalBEMResult, NonlocalBEMResult
using Preconditioners

# device code
include("device/model.jl")
include("device/math.jl")
include("device/matrix.jl")
include("device/laplace.jl")
include("device/yukawa.jl")
include("device/nonlocal.jl")
#export CuPosition, CuPositionVector, CuTriangle, CuTriangleVector, NoAccessArray, ReadOnlyArray,
#    laplacepot_single, laplacepot_double, regularyukawapot_single, regularyukawapot_double

# host code
include("host/common.jl")
include("host/model.jl")
include("host/matrix.jl")
include("host/local.jl")
include("host/nonlocal.jl")
include("host/post.jl")
#export LaplacePotentialMatrix, LocalSystem, LocalSystemMatrix, NonlocalSystem,
#    NonlocalSystemMatrix, NonlocalSystemOutputs, PositionVector, PostProcessor,
#    ReYukawaPotentialMatrix, TriangleVector, φ⃰Ω
export rfenergy, solve, φΩ, φΣ

end # module


text:
mutable struct TennisGame1
    m_score1::Int
	m_score2::Int
	player1Name::String
	player2Name::String

    TennisGame1(player1Name, player2Name) = new(0, 0, player1Name, player2Name)
end

function Tenniskata.won_point(game::TennisGame1, playerName::String)
    if playerName == "player1"
        game.m_score1 += 1
    else
        game.m_score2 += 1
    end
end

function Tenniskata.get_score(game::TennisGame1)
    score = ""
    tempScore = 0
    if game.m_score1 == game.m_score2
        if game.m_score1 == 0
			score = "Love-All"
        elseif game.m_score1 == 1
            score = "Fifteen-All"
        elseif game.m_score1 == 2
            score = "Thirty-All"
        else
            score = "Deuce"
        end
    elseif game.m_score1 >= 4 || game.m_score2 >= 4
        minusResult = game.m_score1 - game.m_score2
        if minusResult == 1
			score = "Advantage player1"
		elseif minusResult == -1
			score = "Advantage player2"
		elseif minusResult >= 2
			score = "Win for player1"
		else
			score = "Win for player2"
        end
    else
        for i = 1:2
            if i == 1
				tempScore = game.m_score1
			else
				score *= "-"
				tempScore = game.m_score2
            end
            if tempScore == 0
                score *= "Love"
            elseif tempScore == 1
                score *= "Fifteen"
            elseif tempScore == 2
                score *= "Thirty"
            elseif tempScore == 3
                score *= "Forty"
            end
        end
    end
    return score
end

max_stars_repo_name: jacob-alt-del/Tennis-Refactoring-Kata

max_issues_repo_name: jacob-alt-del/Tennis-Refactoring-Kata

max_forks_repo_name: jacob-alt-del/Tennis-Refactoring-Kata


text:
export BrickletDMXReadFrameLowLevel
struct BrickletDMXReadFrameLowLevel
    frame_length::Integer
    frame_chunk_offset::Integer
    frame_chunk_data::Vector{Integer}
    frame_number::Integer
end

export BrickletDMXFrameErrorCount
struct BrickletDMXFrameErrorCount
    overrun_error_count::Integer
    framing_error_count::Integer
end

export BrickletDMXFrameCallbackConfig
struct BrickletDMXFrameCallbackConfig
    frame_started_callback_enabled::Bool
    frame_available_callback_enabled::Bool
    frame_callback_enabled::Bool
    frame_error_count_callback_enabled::Bool
end

export BrickletDMXSPITFPErrorCount
struct BrickletDMXSPITFPErrorCount
    error_count_ack_checksum::Integer
    error_count_message_checksum::Integer
    error_count_frame::Integer
    error_count_overflow::Integer
end

export BrickletDMXIdentity
struct BrickletDMXIdentity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletDMXReadFrame
struct BrickletDMXReadFrame
    frame::Vector{Integer}
    frame_number::Integer
end

export BrickletDMX
"""
DMX master and slave
"""
mutable struct BrickletDMX <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletDMX(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_dmx")
        deviceInternal = package.BrickletDMX(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export set_dmx_mode
"""
    $(SIGNATURES)

Sets the DMX mode to either master or slave.

Calling this function sets frame number to 0.
"""
function set_dmx_mode(device::BrickletDMX, dmx_mode)
    device.deviceInternal.set_dmx_mode(dmx_mode)
end

export get_dmx_mode
"""
    $(SIGNATURES)

Returns the DMX mode, as set by :func:`Set DMX Mode`.
"""
function get_dmx_mode(device::BrickletDMX)
    return device.deviceInternal.get_dmx_mode()
end

export write_frame_low_level
"""
    $(SIGNATURES)

Writes a DMX frame. The maximum frame size is 512 byte. Each byte represents one channel.

The next frame can be written after the :cb:`Frame Started` callback was called. The frame
is double buffered, so a new frame can be written as soon as the writing of the prior frame
starts.

The data will be transfered when the next frame duration ends, see :func:`Set Frame Duration`.

Generic approach:

* Set the frame duration to a value that represents the number of frames per second you want to achieve.
* Set channels for first frame.
* Wait for the :cb:`Frame Started` callback.
* Set channels for next frame.
* Wait for the :cb:`Frame Started` callback.
* and so on.

This approach ensures that you can set new DMX data with a fixed frame rate.

This function can only be called in master mode.
"""
function write_frame_low_level(device::BrickletDMX, frame_length, frame_chunk_offset, frame_chunk_data)
    device.deviceInternal.write_frame_low_level(frame_length, frame_chunk_offset, frame_chunk_data)
end

export read_frame_low_level
"""
    $(SIGNATURES)

Returns the last frame that was written by the DMX master. The size of the array
is equivalent to the number of channels in the frame. Each byte represents one channel.

The next frame is available after the :cb:`Frame Available` callback was called.

Generic approach:

* Call :func:`Read Frame` to get first frame.
* Wait for the :cb:`Frame Available` callback.
* Call :func:`Read Frame` to get second frame.
* Wait for the :cb:`Frame Available` callback.
* and so on.

Instead of polling this function you can also use the :cb:`Frame` callback.
You can enable it with :func:`Set Frame Callback Config`.

The frame number starts at 0 and it is increased by one with each received frame.

This function can only be called in slave mode.
"""
function read_frame_low_level(device::BrickletDMX)
    return device.deviceInternal.read_frame_low_level()
end

export set_frame_duration
"""
    $(SIGNATURES)

Sets the duration of a frame.

Example: If you want to achieve 20 frames per second, you should
set the frame duration to 50ms (50ms * 20 = 1 second).

If you always want to send a frame as fast as possible you can set
this value to 0.

This setting is only used in master mode.
"""
function set_frame_duration(device::BrickletDMX, frame_duration)
    device.deviceInternal.set_frame_duration(frame_duration)
end

export get_frame_duration
"""
    $(SIGNATURES)

Returns the frame duration as set by :func:`Set Frame Duration`.
"""
function get_frame_duration(device::BrickletDMX)
    return device.deviceInternal.get_frame_duration()
end

export get_frame_error_count
"""
    $(SIGNATURES)

Returns the current number of overrun and framing errors.
"""
function get_frame_error_count(device::BrickletDMX)
    return device.deviceInternal.get_frame_error_count()
end

export set_communication_led_config
"""
    $(SIGNATURES)

Sets the communication LED configuration. By default the LED shows
communication traffic, it flickers once for every 10 received data packets.

You can also turn the LED permanently on/off or show a heartbeat.

If the Bricklet is in bootloader mode, the LED is off.
"""
function set_communication_led_config(device::BrickletDMX, config)
    device.deviceInternal.set_communication_led_config(config)
end

export get_communication_led_config
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Communication LED Config`
"""
function get_communication_led_config(device::BrickletDMX)
    return device.deviceInternal.get_communication_led_config()
end

export set_error_led_config
"""
    $(SIGNATURES)

Sets the error LED configuration.

By default the error LED turns on if there is any error (see :cb:`Frame Error Count`
callback). If you call this function with the Show-Error option again, the LED
will turn off until the next error occurs.

You can also turn the LED permanently on/off or show a heartbeat.

If the Bricklet is in bootloader mode, the LED is off.
"""
function set_error_led_config(device::BrickletDMX, config)
    device.deviceInternal.set_error_led_config(config)
end

export get_error_led_config
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Error LED Config`.
"""
function get_error_led_config(device::BrickletDMX)
    return device.deviceInternal.get_error_led_config()
end

export set_frame_callback_config
"""
    $(SIGNATURES)

Enables/Disables the different callbacks. By default the
:cb:`Frame Started` callback and :cb:`Frame Available` callback are enabled while
the :cb:`Frame` callback and :cb:`Frame Error Count` callback are disabled.

If you want to use the :cb:`Frame` callback you can enable it and disable
the cb:`Frame Available` callback at the same time. It becomes redundant in
this case.
"""
function set_frame_callback_config(device::BrickletDMX, frame_started_callback_enabled, frame_available_callback_enabled, frame_callback_enabled, frame_error_count_callback_enabled)
    device.deviceInternal.set_frame_callback_config(frame_started_callback_enabled, frame_available_callback_enabled, frame_callback_enabled, frame_error_count_callback_enabled)
end

export get_frame_callback_config
"""
    $(SIGNATURES)

Returns the frame callback config as set by :func:`Set Frame Callback Config`.
"""
function get_frame_callback_config(device::BrickletDMX)
    return device.deviceInternal.get_frame_callback_config()
end

export get_spitfp_error_count
"""
    $(SIGNATURES)

Returns the error count for the communication between Brick and Bricklet.

The errors are divided into

* ACK checksum errors,
* message checksum errors,
* framing errors and
* overflow errors.

The errors counts are for errors that occur on the Bricklet side. All
Bricks have a similar function that returns the errors on the Brick side.
"""
function get_spitfp_error_count(device::BrickletDMX)
    return device.deviceInternal.get_spitfp_error_count()
end

export set_bootloader_mode
"""
    $(SIGNATURES)

Sets the bootloader mode and returns the status after the requested
mode change was instigated.

You can change from bootloader mode to firmware mode and vice versa. A change
from bootloader mode to firmware mode will only take place if the entry function,
device identifier and CRC are present and correct.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_bootloader_mode(device::BrickletDMX, mode)
    return device.deviceInternal.set_bootloader_mode(mode)
end

export get_bootloader_mode
"""
    $(SIGNATURES)

Returns the current bootloader mode, see :func:`Set Bootloader Mode`.
"""
function get_bootloader_mode(device::BrickletDMX)
    return device.deviceInternal.get_bootloader_mode()
end

export set_write_firmware_pointer
"""
    $(SIGNATURES)

Sets the firmware pointer for :func:`Write Firmware`. The pointer has
to be increased by chunks of size 64. The data is written to flash
every 4 chunks (which equals to one page of size 256).

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_write_firmware_pointer(device::BrickletDMX, pointer)
    device.deviceInternal.set_write_firmware_pointer(pointer)
end

export write_firmware
"""
    $(SIGNATURES)

Writes 64 Bytes of firmware at the position as written by
:func:`Set Write Firmware Pointer` before. The firmware is written
to flash every 4 chunks.

You can only write firmware in bootloader mode.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function write_firmware(device::BrickletDMX, data)
    return device.deviceInternal.write_firmware(data)
end

export set_status_led_config
"""
    $(SIGNATURES)

Sets the status LED configuration. By default the LED shows
communication traffic between Brick and Bricklet, it flickers once
for every 10 received data packets.

You can also turn the LED permanently on/off or show a heartbeat.

If the Bricklet is in bootloader mode, the LED is will show heartbeat by default.
"""
function set_status_led_config(device::BrickletDMX, config)
    device.deviceInternal.set_status_led_config(config)
end

export get_status_led_config
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Status LED Config`
"""
function get_status_led_config(device::BrickletDMX)
    return device.deviceInternal.get_status_led_config()
end

export get_chip_temperature
"""
    $(SIGNATURES)

Returns the temperature as measured inside the microcontroller. The
value returned is not the ambient temperature!

The temperature is only proportional to the real temperature and it has bad
accuracy. Practically it is only useful as an indicator for
temperature changes.
"""
function get_chip_temperature(device::BrickletDMX)
    return device.deviceInternal.get_chip_temperature()
end

export reset
"""
    $(SIGNATURES)

Calling this function will reset the Bricklet. All configurations
will be lost.

After a reset you have to create new device objects,
calling functions on the existing ones will result in
undefined behavior!
"""
function reset(device::BrickletDMX)
    device.deviceInternal.reset()
end

export write_uid
"""
    $(SIGNATURES)

Writes a new UID into flash. If you want to set a new UID
you have to decode the Base58 encoded UID string into an
integer first.

We recommend that you use Brick Viewer to change the UID.
"""
function write_uid(device::BrickletDMX, uid)
    device.deviceInternal.write_uid(uid)
end

export read_uid
"""
    $(SIGNATURES)

Returns the current UID as an integer. Encode as
Base58 to get the usual string version.
"""
function read_uid(device::BrickletDMX)
    return device.deviceInternal.read_uid()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletDMX)
    return device.deviceInternal.get_identity()
end

export write_frame
"""
Writes a DMX frame. The maximum frame size is 512 byte. Each byte represents one channel.

The next frame can be written after the :cb:`Frame Started` callback was called. The frame
is double buffered, so a new frame can be written as soon as the writing of the prior frame
starts.

The data will be transfered when the next frame duration ends, see :func:`Set Frame Duration`.

Generic approach:

* Set the frame duration to a value that represents the number of frames per second you want to achieve.
* Set channels for first frame.
* Wait for the :cb:`Frame Started` callback.
* Set channels for next frame.
* Wait for the :cb:`Frame Started` callback.
* and so on.

This approach ensures that you can set new DMX data with a fixed frame rate.

This function can only be called in master mode.
"""
function write_frame(device::BrickletDMX, frame)

        frame = convert(Vector{UInt8}, frame)

    if length(frame) > 65535
        throw(TinkerforgeInvalidParameterError("Frame can be at most 65535 items long"))
    end

    frame_length = length(frame)
    frame_chunk_offset = 0

    if frame_length == 0
        frame_chunk_data = [0] * 60
        ret = write_frame_low_level(device, frame_length, frame_chunk_offset, frame_chunk_data)
    else
        lock(device.stream_lock) do
            while frame_chunk_offset < frame_length
                frame_chunk_data = create_chunk_data(frame, frame_chunk_offset, 60, 0)
                ret = write_frame_low_level(device, frame_length, frame_chunk_offset, frame_chunk_data)
                frame_chunk_offset += 60
            end
        end
    end

        return ret
end

export read_frame
"""
Returns the last frame that was written by the DMX master. The size of the array
is equivalent to the number of channels in the frame. Each byte represents one channel.

The next frame is available after the :cb:`Frame Available` callback was called.

Generic approach:

* Call :func:`Read Frame` to get first frame.
* Wait for the :cb:`Frame Available` callback.
* Call :func:`Read Frame` to get second frame.
* Wait for the :cb:`Frame Available` callback.
* and so on.

Instead of polling this function you can also use the :cb:`Frame` callback.
You can enable it with :func:`Set Frame Callback Config`.

The frame number starts at 0 and it is increased by one with each received frame.

This function can only be called in slave mode.
"""
function read_frame(device::BrickletDMX)

    lock(device.stream_lock) do
        ret = read_frame_low_level(device, )
            frame_length = ret.frame_length
        frame_out_of_sync = ret.frame_chunk_offset != 0

        frame_data = ret.frame_chunk_data

        while !frame_out_of_sync && length(frame_data) < frame_length
            ret = read_frame_low_level(device, )
                frame_length = ret.frame_length
            frame_out_of_sync = ret.frame_chunk_offset != length(frame_data)
            frame_data += ret.frame_chunk_data
        end

        if frame_out_of_sync # discard remaining stream to bring it back in-sync
            while ret.frame_chunk_offset + 56 < frame_length
                ret = read_frame_low_level(device, )
                    frame_length = ret.frame_length
            end

            throw(TinkerforgeStreamOutOfSyncError("Frame stream is out-of-sync"))
        end
    end

        return ReadFrame(frame_data[:frame_length], ret.frame_number)
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletDMX, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
using SparseArrays

@testset "right_mul!" begin

    # NOTE
    # currently julia/sparsearrays does not implement a
    # mul!(C::StridedMatrix, X::StridedMatrix, A::SparseMatrixCSC)
    # instead it falls back to some generic matrix multiplciation in LinearAlgebra
    #
    # If this throws:
    # - check if mul!() is now specialized for the above type signature
    #   - if it is, remove right_mul! in helpers.jl & change right_muL! to mul! in
    #     slice_matrices.jl
    #   - otherwise, update the type signature
    #
    # see:
    # Julia PR: https://github.com/JuliaLang/julia/pull/24045
    # Julia Issue: https://github.com/JuliaLang/julia/issues/29956#issuecomment-440867024
    # MC PR: https://github.com/crstnbr/MonteCarlo.jl/pull/9
    # discourse: https://discourse.julialang.org/t/asymmetric-speed-of-in-place-sparse-dense-matrix-product/10256

    mul_dense_sparse_is_missing = !occursin("SparseArrays", string(which(mul!, (Matrix, Matrix, SparseMatrixCSC)).file))
    @test mul_dense_sparse_is_missing

    # check that right_mul! is correct
    for Lx in (4, 8)
        for Ly in (4, 8)
            M = rand(Lx, Ly)
            S = sparse(rand(Ly, Lx))
            X = rand(Lx, Lx)
            @test M * S == MonteCarlo.mul!(X, M, S)
        end
    end
end

@testset "Slice Matrices" begin
    m = HubbardModelAttractive(8, 1)
    dqmc = DQMC(m, beta=5.0)
    eT = dqmc.stack.hopping_matrix_exp
    eV = similar(dqmc.stack.eV)
    A = similar(eT)

    # No checkerboard
    for slice in rand(1:50, 2)
        MonteCarlo.interaction_matrix_exp!(dqmc, m, eV, dqmc.conf, slice, 1.)

        # MonteCarlo.slice_matrix
        @test MonteCarlo.slice_matrix(dqmc, m, slice, 1.) ≈ eT * eT * eV
        @test MonteCarlo.slice_matrix(dqmc, m, slice, -1.) ≈ inv(eT * eT * eV)
        # MonteCarlo.slice_matrix!(dqmc, m, slice, 1., A)
        # @test A ≈ eT * eT * eV
        # MonteCarlo.slice_matrix!(dqmc, m, slice, -1., A)
        # @test A ≈ inv(eT * eT * eV)

        # MonteCarlo.multiply_slice_matrix...
        A = eT * eT * eV
        input = rand(size(A)...)
        result = A * input
        MonteCarlo.multiply_slice_matrix_left!(dqmc, m, slice, input)
        @test input ≈ result
        input = rand(size(A)...)
        result = input * A
        MonteCarlo.multiply_slice_matrix_right!(dqmc, m, slice, input)
        @test input ≈ result

        A = inv(eT * eT * eV)
        input = rand(size(A)...)
        result = A * input
        MonteCarlo.multiply_slice_matrix_inv_left!(dqmc, m, slice, input)
        @test input ≈ result
        input = rand(size(A)...)
        result = input * A
        MonteCarlo.multiply_slice_matrix_inv_right!(dqmc, m, slice, input)
        @test input ≈ result

        A = adjoint(eT * eT * eV)
        input = rand(size(A)...)
        result = A * input
        MonteCarlo.multiply_daggered_slice_matrix_left!(dqmc, m, slice, input)
        @test input ≈ result
    end

    # Checkerboard
    dqmc = DQMC(m, beta=5.0, checkerboard=true)

    for slice in rand(1:50, 2)
        MonteCarlo.interaction_matrix_exp!(dqmc, m, eV, dqmc.conf, slice, 1.)

        # MonteCarlo.slice_matrix
        @test maximum(abs,
            MonteCarlo.slice_matrix(dqmc, m, slice, 1.) .- eT * eT * eV
        ) < 2dqmc.parameters.delta_tau
        @test maximum(abs,
            MonteCarlo.slice_matrix(dqmc, m, slice, -1.) .- inv(eT * eT * eV)
        ) < 2dqmc.parameters.delta_tau
        MonteCarlo.slice_matrix!(dqmc, m, slice, 1., A)
        @test maximum(abs, A .- eT * eT * eV) < 2dqmc.parameters.delta_tau
        MonteCarlo.slice_matrix!(dqmc, m, slice, -1., A)
        @test maximum(abs, A .- inv(eT * eT * eV)) < 2dqmc.parameters.delta_tau

        # MonteCarlo.multiply_slice_matrix...
        A = eT * eT * eV
        input = rand(size(A)...)
        result = A * input
        MonteCarlo.multiply_slice_matrix_left!(dqmc, m, slice, input)
        @test maximum(abs, input .- result) < 2dqmc.parameters.delta_tau
        input = rand(size(A)...)
        result = input * A
        MonteCarlo.multiply_slice_matrix_right!(dqmc, m, slice, input)
        @test maximum(abs, input .- result) < 2dqmc.parameters.delta_tau

        A = inv(eT * eT * eV)
        input = rand(size(A)...)
        result = A * input
        MonteCarlo.multiply_slice_matrix_inv_left!(dqmc, m, slice, input)
        @test maximum(abs, input .- result) < 2dqmc.parameters.delta_tau
        input = rand(size(A)...)
        result = input * A
        MonteCarlo.multiply_slice_matrix_inv_right!(dqmc, m, slice, input)
        @test maximum(abs, input .- result) < 2dqmc.parameters.delta_tau

        A = adjoint(eT * eT * eV)
        input = rand(size(A)...)
        result = A * input
        MonteCarlo.multiply_daggered_slice_matrix_left!(dqmc, m, slice, input)
        @test maximum(abs, input .- result) < 2dqmc.parameters.delta_tau
    end
end

using MonteCarlo: vmul!, lvmul!, rvmul!, rdivp!, udt_AVX_pivot!, rvadd!, vsub!
using MonteCarlo: vmin!, vmininv!, vmax!, vmaxinv!, vinv!
using MonteCarlo: BlockDiagonal#, CMat64, CVec64, StructArray

@testset "Custom Linear Algebra " begin
    for type in (Float64, ComplexF64)
        @testset "avx multiplications ($type)" begin
            A = rand(type, 16, 16)
            B = rand(type, 16, 16)
            C = rand(type, 16, 16)
            atol = 100eps(Float64)
            type == ComplexF64 && (atol *= sqrt(2))

            vmul!(C, A, B)
            @test A * B ≈ C atol = atol

            vmul!(C, A, B')
            @test A * B' ≈ C atol = atol

            vmul!(C, A', B)
            @test A' * B ≈ C atol = atol

            vmul!(C, A', B')
            @test A' * B' ≈ C atol = atol

            D = Diagonal(rand(16))
            vmul!(C, A, D)
            @test A * D ≈ C atol = atol

            copyto!(C, A)
            rvmul!(C, D)
            @test A * D ≈ C atol = atol

            copyto!(C, A)
            lvmul!(D, C)
            @test D * A ≈ C atol = atol

            copyto!(C, A)
            rvadd!(C, D)
            @test A + D ≈ C atol = atol

            copyto!(C, A)
            rvadd!(C, B)
            @test A + B ≈ C atol = atol

            vsub!(C, A, I)
            @test A - I ≈ C atol = atol

            if type == Float64
                v = rand(16) .+ 0.5
                w = copy(v)

                vmin!(v, w)
                @test v ≈ min.(1.0, w) atol = atol

                vmininv!(v, w)
                @test v ≈ 1.0 ./ min.(1.0, w) atol = atol

                vmax!(v, w)
                @test v ≈ max.(1.0, w) atol = atol

                vmaxinv!(v, w)
                @test v ≈ 1.0 ./ max.(1.0, w) atol = atol

                v = copy(w)
                vinv!(w)
                @test w ≈ 1.0 ./ v atol = atol
            end
        end

        @testset "UDT transformations + rdivp! ($type)" begin
            U = Matrix{Float64}(undef, 16, 16)
            D = Vector{Float64}(undef, 16)
            T = rand(16, 16)
            X = copy(T)
            MonteCarlo.udt_AVX!(U, D, T)
            @test U * Diagonal(D) * T ≈ X

            U = Matrix{type}(undef, 16, 16)
            T = rand(type, 16, 16)
            X = copy(T)
            pivot = Vector{Int64}(undef, 16)
            tempv = Vector{type}(undef, 16)
            udt_AVX_pivot!(U, D, T)
            @test U * Diagonal(D) * T ≈ X

            copyto!(T, X)
            pivot = Vector{Int64}(undef, 16)
            tempv = Vector{type}(undef, 16)
            udt_AVX_pivot!(U, D, T, pivot, tempv, Val(false))
            # pivoting matrix
            P = zeros(length(pivot), length(pivot))
            for (i, j) in enumerate(pivot)
                P[i, j] = 1.0
            end
            @test U * Diagonal(D) * UpperTriangular(T) * P ≈ X

            u = copy(U)
            t = copy(T)
            tmp = similar(T)
            rdivp!(u, t, tmp, pivot)
            @test u ≈ U * P' / UpperTriangular(T)
        end
    end

    #=
    @testset "Complex StructArray" begin
        M1 = rand(ComplexF64, 8, 8)
        C1 = StructArray(M1)
        M2 = rand(ComplexF64, 8, 8)
        C2 = StructArray(M2)
        M3 = rand(ComplexF64, 8, 8)
        C3 = StructArray(M3)

        @test C1 isa CMat64
        @test M1 == C1

        # Test avx multiplications
        vmul!(C1, C2, C3)
        vmul!(M1, M2, M3)
        @test M1 == C1

        vmul!(C1, C2, adjoint(C3))
        vmul!(M1, M2, adjoint(M3))
        @test M1 ≈ C1 atol=1e-14 # check

        vmul!(C1, adjoint(C2), C3)
        vmul!(M1, adjoint(M2), M3)
        @test M1 ≈ C1 atol=1e-14 #check

        D = Diagonal(rand(8))
        vmul!(C1, C2, D)
        vmul!(M1, M2, D)
        @test M1 == C1

        rvmul!(C1, D)
        rvmul!(M1, D)
        @test M1 == C1

        lvmul!(D, C1)
        lvmul!(D, M1)
        @test M1 == C1

        # Test UDT and rdivp!
        M2 = Matrix(C2)
        D = rand(8)
        pivot = Vector{Int64}(undef, 8)
        tempv = StructArray(Vector{ComplexF64}(undef, 8))
        @test tempv isa CVec64
        udt_AVX_pivot!(C1, D, C2, pivot, tempv, Val(false))
        P = zeros(length(pivot), length(pivot))
        for (i, j) in enumerate(pivot)
            P[i, j] = 1.0
        end
        @test Matrix(C1) * Diagonal(D) * UpperTriangular(Matrix(C2)) * P ≈ M2 #check

        M1 = Matrix(C1)
        M2 = Matrix(C2)
        rdivp!(C1, C2, C3, pivot)
        @test C1 ≈ M1 * P' / UpperTriangular(M2) # check
    end
    =#

    @testset "BlockDiagonal" begin
        b1 = rand(4, 4)
        b2 = rand(4, 4)
        atol = 100eps(Float64)

        B = BlockDiagonal(b1, b2)
        @test B isa BlockDiagonal{Float64, 2, Matrix{Float64}}

        # Check values/indexing
        for i in 1:4, j in 1:4
            @test B[i, j] == b1[i, j]
            @test B[4+i, 4+j] == b2[i, j]
            @test B[4+i, j] == 0.0
            @test B[i, 4+j] == 0.0
        end

        B1 = copy(B)
        M1 = Matrix(B1)
        @test M1 == B1

        B2 = BlockDiagonal(rand(4, 4), rand(4, 4))
        M2 = Matrix(B2)
        B3 = BlockDiagonal(rand(4, 4), rand(4, 4))
        M3 = Matrix(B3)

        # Test avx multiplications
        vmul!(B1, B2, B3)
        vmul!(M1, M2, M3)
        @test M1 ≈ B1 atol = atol

        vmul!(B1, B2, adjoint(B3))
        vmul!(M1, M2, adjoint(M3))
        @test M1 ≈ B1

        vmul!(B1, adjoint(B2), B3)
        vmul!(M1, adjoint(M2), M3)
        @test M1 ≈ B1 atol = atol

        D = Diagonal(rand(8))
        vmul!(B1, B2, D)
        vmul!(M1, M2, D)
        @test M1 ≈ B1 atol = atol

        rvmul!(B1, D)
        rvmul!(M1, D)
        @test M1 ≈ B1 atol = atol

        lvmul!(D, B1)
        lvmul!(D, M1)
        @test M1 ≈ B1 atol = atol

        rvadd!(B1, D)
        rvadd!(M1, D)
        @test M1 ≈ B1 atol = atol

        rvadd!(B1, B2)
        rvadd!(M1, M2)
        @test M1 ≈ B1 atol = atol

        vsub!(B1, B2, I)
        @test M2 - I ≈ B1 atol = atol

        x = rand()
        @test B2 * x ≈ M2 * x atol = atol

        @test log(B2) ≈ log(M2) atol = atol

        @test det(B2) ≈ det(M2) atol = atol

        transpose!(B1, B2)
        @test B1 ≈ M2' atol = atol

        vmul!(B1, B2', B3')
        @test B1 ≈ M2' * M3' atol = atol

        # Test UDT and rdivp!
        M2 = Matrix(B2)
        D = rand(8)
        pivot = Vector{Int64}(undef, 8)
        tempv = Vector{Float64}(undef, 8)
        udt_AVX_pivot!(B1, D, B2, pivot, tempv, Val(false))

        P = zeros(length(pivot), length(pivot))
        # pivot is per block, temporary make it global (in the sense of matrix indices)
        n = size(B1.blocks[1], 1); N = length(B.blocks)
        for i in 1:N; pivot[(i-1)*n+1 : i*n] .+= (i-1)*n; end
        for (i, j) in enumerate(pivot)
            P[i, j] = 1.0
        end
        for i in 1:N; pivot[(i-1)*n+1 : i*n] .-= (i-1)*n; end

        @test Matrix(B1) * Diagonal(D) * UpperTriangular(Matrix(B2)) * P ≈ M2

        M1 = Matrix(B1)
        M2 = Matrix(B2)
        rdivp!(B1, B2, B3, pivot)
        @test B1 ≈ M1 * P' / UpperTriangular(M2)
    end
end


text:
# This file is part of Kpax3. License is MIT.

function mutation!(o::KOffspring,
                   mrate::Float64)
  n = length(o.R)

  w = zeros(Int, n)
  found = false
  t = 0
  i = 0

  for a in 1:n
    if rand() <= mrate
      # remove a from its cluster
      o.v[o.R[a]] -= 1
      t = o.v[o.R[a]]

      copyto!(w, o.v)

      # do not put i back into its cluster
      if t > 0
        # candidate empty cluster to be filled
        found = false
        i = 0
        while !found && i < n
          i += 1
          if w[i] == 0
            w[i] = 1
            found = true
          end
        end

        w[o.R[a]] = 0
      end

      o.R[a] = StatsBase.sample(StatsBase.ProbabilityWeights(w, n - t))
      o.v[o.R[a]] += 1
    end
  end

  nothing
end