text:
"""
Computes entropy of a sequence.
"""
function entropy(X :: AbstractVector)
  counts = Dict{eltype(X), Int}()
  n = length(X)
  for i = 1:length(X)
    counts[X[i]] = get(counts, X[i], 0) + 1
  end
  sum(kv -> log2(n / kv[2]) * kv[2], counts) / n
end

"""
Computes conditional entropy `H(Y | X)`.
"""
function conditional_entropy(X :: AbstractVector, Y :: AbstractVector)
  cnts_for_pref = Dict{eltype(X), Dict{eltype(Y), Int}}()
  n = length(X)
  for i = 1:length(X)
    if !haskey(cnts_for_pref, X[i])
      cnts_for_pref[X[i]] = Dict{eltype(Y), Int}()
    end
    counts = cnts_for_pref[X[i]]
    counts[Y[i]] = get(counts, Y[i], 0) + 1
  end

  # this helper function computes "unnormalized entropy" for given counts dict
  # i.e. entropy multiplied by sum of counts
  function _prefix_helper(counts :: Dict{eltype(Y), Int})
    counts_sum = sum(kv -> kv[2], counts)
    sum(kv -> log2(counts_sum / kv[2]) * kv[2], counts)
  end

  sum(kv -> _prefix_helper(kv[2]), cnts_for_pref) / n
end


text:
"""
    couple_flux_matrix(model, Q_blocks[, symmetrical])

Create global flux matrix with coupling between compartments.
"""
function couple_flux_matrix(model, Q_blocks, symmetrical = false)

    # Extract mesh fields
    (; mesh, ρ, κ) = model
    (; point_map, points, facets) = mesh

    # Sizes
    nboundary = size(facets, 2)
    inds_cmpts = cumsum([0; size.(points, 2)])
    npoint = inds_cmpts[end]
    get_inds(icmpt) = inds_cmpts[icmpt]+1:inds_cmpts[icmpt+1]

    # Assemble sparse block diagonal flux matrix
    Q = spzeros(npoint, npoint)

    # Couple flux matrix at boundaries
    for iboundary = 1:nboundary

        # Check if boundary is an interface between two compartments
        cmpts_touch = findall(.!isempty.(facets[:, iboundary]))
        ntouch = length(cmpts_touch)
        if ntouch == 1
            # Only one compartment touches boundary. It is thus an outer boundary,
            # and may possibly have a boundary relaxation coefficient
            cmpt = cmpts_touch[1]

            # Global indices of the boundary
            inds = get_inds(cmpt)

            # Add boundary contribution to global flux matrix for compartment
            Q[inds, inds] += κ[iboundary] * Q_blocks[cmpt, iboundary]
        elseif ntouch == 2
            # Extract flux matrices from corresponding compartments
            cmpt₁, cmpt₂ = cmpts_touch[1], cmpts_touch[2]
            Q₁₁ = Q_blocks[cmpt₁, iboundary]
            Q₂₂ = Q_blocks[cmpt₂, iboundary]
            Q₁₂ = spzeros(size(Q₁₁, 1), size(Q₂₂, 2))
            Q₂₁ = spzeros(size(Q₂₂, 1), size(Q₁₁, 2))

            # Identify pairs of points using global numbering
            inds₁ = unique(facets[cmpt₁, iboundary])
            inds₂ = unique(facets[cmpt₂, iboundary])

            # Check if local indices are already sorted with a one-to-one correspondence
            if all(point_map[cmpt₁][inds₁] .== point_map[cmpt₂][inds₂])
                indinds₁ = 1:length(inds₁)
                indinds₂ = 1:length(inds₂)
            else
                pairmap = point_map[cmpt₁][inds₁] .== point_map[cmpt₂][inds₂]'
                indinds₁ = 1:length(inds₁)
                indinds₂ = [findfirst(pairmap[i, :]) for i ∈ indinds₁]
            end

            # Create coupling blocks
            Q₁₂[:, inds₂[indinds₂]] = Q₁₁[:, inds₁[indinds₁]]
            Q₂₁[:, inds₁[indinds₁]] = Q₂₂[:, inds₂[indinds₂]]

            # Check whether to use same permeability coefficients on both sides
            if symmetrical
                # Use the same permeability coefficient on each side of the boundary
                c₁₂ = 1
                c₂₁ = 1
            else
                # Weigh permeability coefficients on each side of the boundary with
                # initial spin density equilibrium
                ρ₁ = ρ[cmpt₁]
                ρ₂ = ρ[cmpt₂]
                c₂₁ = 2 * ρ₂ / (ρ₁ + ρ₂)
                c₁₂ = 2 * ρ₁ / (ρ₁ + ρ₂)
            end

            # Adjust permeability coefficients
            κ₁ = c₂₁ * κ[iboundary]
            κ₂ = c₁₂ * κ[iboundary]

            # Global indices of the boundary in each of the compartments
            inds₁ = get_inds(cmpt₁)
            inds₂ = get_inds(cmpt₂)

            # Add interface contribution to global flux matrix for compartment 1
            Q[inds₁, inds₁] += κ₁ * Q₁₁
            Q[inds₁, inds₂] -= κ₂ * Q₁₂

            # Add interface contribution to global flux matrix for compartment 2
            Q[inds₂, inds₁] -= κ₁ * Q₂₁
            Q[inds₂, inds₂] += κ₂ * Q₂₂
        else
            error("Each interface only touches 1 or 2 compartments")
        end
    end

    Q
end

ext: jl


text:
function sc(x, m, M)
    y = (x .- minimum(x))./(maximum(x)-minimum(x))
    y .*= M-m
    y .+= m
    return y
end

plot(frame=:none, dpi=400, size=(800,600), legend=false)
n_input = size(params(m)[1],2)
lay1 = sc(collect(1:(n_input+0.5n_input-1)), -1.2, 1.2)
y_in = lay1[vcat(collect(1:3:length(lay1)), collect(2:3:length(lay1)))]

# Scaling
n2 = size(params(m)[1],1)
y_2 = sc(collect(1:n2), -n2/n_input, n2/n_input)
for i in eachindex(y_in)
    for j in eachindex(y_2)
        plot!([1, 2], [y_in[i], y_2[j]], lab=false, line_z = params(m)[1][j,i], alpha=abs(params(m)[1][j,i])^1.2, clim=(-1,1), c=:BrBG)
    end
end

n3 = size(params(m)[3],1)
y_3 = sc(collect(1:n3), -n3/n_input, n3/n_input)
for i in eachindex(y_2)
    for j in eachindex(y_3)
        plot!([2, 3], [y_2[i], y_3[j]], lab=false, line_z = params(m)[3][j,i], alpha=abs(params(m)[3][j,i])^1.2, clim=(-1,1), c=:BrBG)
    end
end

n4 = size(params(m)[5],1)
y_4 = sc(collect(1:n4), -n4/n_input, n4/n_input)
for i in eachindex(y_3)
    for j in eachindex(y_4)
        plot!([3, 4], [y_3[i], y_4[j]], lab=false, line_z = params(m)[5][j,i], alpha=abs(params(m)[5][j,i])^1.2, clim=(-1,1), c=:BrBG)
    end
end

n5 = size(params(m)[7],1)
y_5 = sc(collect(1:n5), -n5/n_input, n5/n_input)
for i in eachindex(y_4)
    for j in eachindex(y_5)
        plot!([4, 5], [y_4[i], y_5[j]], lab=false, line_z = params(m)[7][j,i], alpha=abs(params(m)[7][j,i])^1.2, clim=(-1,1), c=:BrBG)
    end
end

scatter!(fill(1, length(y_in)), y_in, m=:circle, msw=0.0, c=:grey, lab="", ms=1.4)
scatter!(fill(2, length(params(m)[2])), y_2, marker_z = params(m)[2], c=:BrBG, msw=0.2, lab="", clim=(-1,1))
scatter!(fill(3, length(params(m)[4])), y_3, marker_z = params(m)[4], c=:BrBG, msw=0.2, lab="", clim=(-1,1))
scatter!(fill(4, length(params(m)[6])), y_4, marker_z = params(m)[6], c=:BrBG, msw=0.2, lab="", clim=(-1,1))
scatter!(fill(5, length(params(m)[8])), y_5, marker_z = params(m)[8], c=:BrBG, msw=0.2, lab="", clim=(-1,1))

max_stars_repo_name: PoisotLab/NetworkLearningIllustration

max_issues_repo_name: PoisotLab/NetworkLearningIllustration

max_forks_repo_name: PoisotLab/NetworkLearningIllustration


text:
using LinearAlgebra
using GaussianDistributions
using CovarianceFunctions
const Kernel = CovarianceFunctions
using Kernel
using Kernel: gramian
using GaussianDistributions: ConditionalMean, ConditionalKernel
using WoodburyIdentity
using SARA: random_sampling, uncertainty_sampling, inner_sampling, integrated_uncertainty_sampling, ucb, integrated_uncertainty
using ForwardDiff

# C = C | (xi, yi, σ^2) # TODO: make efficient by getting rid of hierarchical structure
function integrated_deviation(f, g, x::AbstractVector, w::AbstractVector = fill(1/length(x), length(x)))
    h(x) = (f(x) - g(x))^2
    return sqrt(integral(h, x, w))
end
function integral(f, x::AbstractVector, w::AbstractVector = fill(1/length(x), length(x)))
    dot(w, f.(x))
end

function sampling_path(policy, G::Gaussian, z::AbstractVector, σ::Real, nsample::Int, ground_truth)
    xs = collect(z) # potential points and integration nodes (z)
    C = G
    x, y = zeros(nsample), zeros(nsample)
    d, dd, u, du = zeros(nsample), zeros(nsample), zeros(nsample), zeros(nsample)
    for i in 1:nsample # random sampling
        xi = policy(C, xs)
        println(xi)
        x[i] = xi
        y[i] = g(xi) # noisy observation
        xs = deleteat!(xs, findfirst(==(xi), xs))
        C = G | @views (x[1:i], y[1:i], σ^2)

        # d[i] = mean((mean(C, z) .- ground_truth).^2)
        u[i] = integrated_uncertainty(C, z)
        D = GaussianDistributions.gradient(C) # multi-ouput process
        d[i] = mean((mean(D[2], z) .- ground_truth).^2)
        du[i] = integrated_uncertainty(D[2], z)
    end
    return C, d, u, du
end

# synthetic example for inner loop acquisition
l = 1/4
k = Kernel.Lengthscale(Kernel.EQ(), l)
G = Gaussian(k)

f(x::Real) = sin(2π*x)
σ = 1e-1
g(x::Real) = f(x) + σ*randn() # noisy observations

# grid of potential points
n = 512
z = range(-1, 1, length = n)
# ground_truth = f.(z)
ground_truth = (x->ForwardDiff.derivative(f, x)).(z)

nsample = 32
sample_helper(policy) = sampling_path(policy, G, z, σ, nsample, ground_truth)
println("random")
Cr, dr, ur, dur = sample_helper(random_sampling)

println("uncertainty")
Cu, du, uu, duu = sample_helper(uncertainty_sampling)

println("integrated uncertainty")
policy = (G, x) -> integrated_uncertainty_sampling(G, x, σ)
Ciu, diu, uiu, duiu = sample_helper(policy)

println("gradient uncertainty")
policy = (G, x) -> inner_sampling(G, x, σ)
Ci, dgu, ugu, dugu = sample_helper(policy)

deviations = (ur, du, diu)
uncertainties = (ur, uu, uiu)

using Plots
plotly()
plot(scale = :log10, legend = :outerbottomleft)
plot!(dr, label = "random")
plot!(du, label = "uncertainty")
plot!(diu, label = "integrated uncertainty")
plot!(dgu, label = "gradient uncertainty")
gui()

# plot(yscale = :log10)
# plot!(ur, label = "random")
# plot!(uu, label = "uncertainty")
# plot!(uiu, label = "integrated uncertainty")
# plot!(ui, label = "gradient uncertainty")
# gui()
#
plot(yscale = :log10)
plot!(dur, label = "random")
plot!(duu, label = "uncertainty")
plot!(duiu, label = "integrated uncertainty")
plot!(dugu, label = "gradient uncertainty")
gui()


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

types = [
    Bool,
    Float16,
    Float32,
    Float64,
    Int128,
    Int16,
    Int32,
    Int64,
    Int8,
    UInt16,
    UInt32,
    UInt64,
    UInt8,
]

# Nullable{T}() = new(true)
for T in types
    x = Nullable{T}()
    @test x.hasvalue === false
    @test isa(x.value, T)
    @test eltype(Nullable{T}) === T
    @test eltype(x) === T
end

# Nullable{T}(value::T) = new(false, value)
for T in types
    x = Nullable{T}(zero(T))
    @test x.hasvalue === true
    @test isa(x.value, T)
    @test x.value === zero(T)
    @test eltype(x) === T

    x = Nullable{T}(one(T))
    @test x.hasvalue === true
    @test isa(x.value, T)
    @test x.value === one(T)
    @test eltype(x) === T
end

# Nullable{T}(value::T, hasvalue::Bool) = new(hasvalue, value)
for T in types
    x = Nullable{T}(zero(T), true)
    @test x.hasvalue === true
    @test isa(x.value, T)
    @test x.value === zero(T)
    @test eltype(x) === T

    x = Nullable{T}(zero(T), false)
    @test x.hasvalue === false
    @test isa(x.value, T)
    @test eltype(Nullable{T}) === T
    @test eltype(x) === T
end

# immutable NullException <: Exception
@test isa(NullException(), NullException)
@test_throws NullException throw(NullException())

# Nullable{T}(value::T) = Nullable{T}(value)
for T in types
    v = zero(T)
    x = Nullable(v)
    @test x.hasvalue === true
    @test isa(x.value, T)
    @test x.value === v

    v = one(T)
    x = Nullable(v)
    @test x.hasvalue === true
    @test isa(x.value, T)
    @test x.value === v
end

# show{T}(io::IO, x::Nullable{T})
io1 = IOBuffer()
io2 = IOBuffer()
for (i, T) in enumerate(types)
    x1 = Nullable{T}()
    x2 = Nullable(zero(T))
    x3 = Nullable(one(T))
    show(io1, x1)
    @test String(take!(io1)) == @sprintf("Nullable{%s}()", T)
    show(io1, x2)
    showcompact(io2, get(x2))
    @test String(take!(io1)) == @sprintf("Nullable{%s}(%s)", T, String(take!(io2)))
    show(io1, x3)
    showcompact(io2, get(x3))
    @test String(take!(io1)) == @sprintf("Nullable{%s}(%s)", T, String(take!(io2)))

    a1 = [x2]
    show(IOContext(io1, compact=false), a1)
    show(IOContext(io2, compact=false), x2)
    @test String(take!(io1)) ==
        @sprintf("Nullable{%s}[%s]", string(T), String(take!(io2)))

    show(io1, a1)
    show(IOContext(io2, compact=true), x2)
    @test String(take!(io1)) ==
        @sprintf("Nullable{%s}[%s]", string(T), String(take!(io2)))
end

module NullableTestEnum
const curmod = current_module()
const curmod_name = fullname(curmod)
const curmod_prefix = "$(["$m." for m in curmod_name]...)"
io = IOBuffer()
@enum TestEnum a b
show(io, Nullable(a))
Base.Test.@test String(take!(io)) == "Nullable{$(curmod_prefix)TestEnum}(a)"
end

# showcompact(io::IO, x::Nullable)
io1 = IOBuffer()
io2 = IOBuffer()
for (i, T) in enumerate(types)
    x1 = Nullable{T}()
    x2 = Nullable(zero(T))
    x3 = Nullable(one(T))
    showcompact(io1, x1)
    @test String(take!(io1)) == "#NULL"
    showcompact(io1, x2)
    showcompact(io2, get(x2))
    @test String(take!(io1)) == String(take!(io2))
    showcompact(io1, x3)
    showcompact(io2, get(x3))
    @test String(take!(io1)) == String(take!(io2))

    a1 = [x2]
    showcompact(io1, a1)
    showcompact(io2, x2)
    @test String(take!(io1)) ==
        @sprintf("Nullable{%s}[%s]", string(T), String(take!(io2)))
end

# get(x::Nullable)
for T in types
    x1 = Nullable{T}()
    x2 = Nullable(zero(T))
    x3 = Nullable(one(T))

    @test_throws NullException get(x1)
    @test get(x2) === zero(T)
    @test get(x3) === one(T)
end

@test_throws NullException get(Nullable())

# get{S, T}(x::Nullable{S}, y::T)
for T in types
    x0 = Nullable()
    x1 = Nullable{T}()
    x2 = Nullable(zero(T))
    x3 = Nullable(one(T))

    @test get(x0, zero(T)) === zero(T)
    @test get(x0, one(T)) === one(T)
    @test get(x1, zero(T)) === zero(T)
    @test get(x1, one(T)) === one(T)
    @test get(x2, one(T)) === zero(T)
    @test get(x3, zero(T)) === one(T)
end

for T in types
    # unsafe_get(x::Nullable)
    x1 = Nullable{T}()
    x2 = Nullable(zero(T))
    x3 = Nullable(one(T))
    a = rand(T)
    x4 = Nullable(a)

    @test isa(unsafe_get(x1), T)
    @test unsafe_get(x2) === zero(T)
    @test unsafe_get(x3) === one(T)
    @test unsafe_get(x4) === a

    # unsafe_get(x)
    x2 = zero(T)
    x3 = one(T)
    x4 = rand(T)

    @test unsafe_get(x2) === zero(T)
    @test unsafe_get(x3) === one(T)
    @test unsafe_get(x4) === x4
end

@test_throws UndefRefError unsafe_get(Nullable())
@test_throws UndefRefError unsafe_get(Nullable{String}())
@test_throws UndefRefError unsafe_get(Nullable{Array}())

for T in types
    # isnull(x::Nullable)
    x1 = Nullable{T}()
    x2 = Nullable(zero(T))
    x3 = Nullable(one(T))

    @test isnull(x1) === true
    @test isnull(x2) === false
    @test isnull(x3) === false

    # isnull(x)
    x1 = zero(T)
    x2 = one(T)
    x3 = rand(T)

    @test isnull(x1) === false
    @test isnull(x2) === false
    @test isnull(x3) === false
end

@test isnull(Nullable())

# function =={S, T}(x::Nullable{S}, y::Nullable{T})
for T in types
    x0 = Nullable()
    x1 = Nullable{T}()
    x2 = Nullable{T}()
    x3 = Nullable(zero(T))
    x4 = Nullable(one(T))

    @test_throws NullException (x0 == x1)
    @test_throws NullException (x0 == x2)
    @test_throws NullException (x0 == x3)
    @test_throws NullException (x0 == x4)

    @test_throws NullException (x1 == x1)
    @test_throws NullException (x1 == x2)
    @test_throws NullException (x1 == x3)
    @test_throws NullException (x1 == x4)

    @test_throws NullException (x2 == x1)
    @test_throws NullException (x2 == x2)
    @test_throws NullException (x2 == x3)
    @test_throws NullException (x2 == x4)

    @test_throws NullException (x3 == x1)
    @test_throws NullException (x3 == x2)
    @test_throws NullException (x3 == x3)
    @test_throws NullException (x3 == x4)

    @test_throws NullException (x4 == x1)
    @test_throws NullException (x4 == x2)
    @test_throws NullException (x4 == x3)
    @test_throws NullException (x4 == x4)
end

# function hash(x::Nullable, h::UInt)
for T in types
    x0 = Nullable()
    x1 = Nullable{T}()
    x2 = Nullable{T}()
    x3 = Nullable(zero(T))
    x4 = Nullable(one(T))

    @test isa(hash(x0), UInt)
    @test isa(hash(x1), UInt)
    @test isa(hash(x2), UInt)
    @test isa(hash(x3), UInt)
    @test isa(hash(x4), UInt)

    @test hash(x0) == hash(x2)
    @test hash(x0) != hash(x3)
    @test hash(x0) != hash(x4)
    @test hash(x1) == hash(x2)
    @test hash(x1) != hash(x3)
    @test hash(x1) != hash(x4)
    @test hash(x2) != hash(x3)
    @test hash(x2) != hash(x4)
    @test hash(x3) != hash(x4)
end

type TestNType{T}
    v::Nullable{T}
end

for T in types
    x1 = TestNType{T}(Nullable{T}())
    @test isnull(x1.v)
    x1.v = one(T)
    @test !isnull(x1.v)
    @test get(x1.v, one(T)) === one(T)
end

# Operators
TestTypes = Union{Base.NullSafeTypes, BigInt, BigFloat,
                  Complex{Int}, Complex{Float64}, Complex{BigFloat},
                  Rational{Int}, Rational{BigInt}}.types
for S in TestTypes, T in TestTypes
    u0 = zero(S)
    u1 = one(S)
    if S <: AbstractFloat
        u2 = S(NaN)
    elseif S <: Complex && S.parameters[1] <: AbstractFloat
        u2 = S(NaN, NaN)
    else
        u2 = u1
    end

    v0 = zero(T)
    v1 = one(T)
    if T <: AbstractFloat
        v2 = T(NaN)
    elseif T <: Complex && T.parameters[1] <: AbstractFloat
        v2 = T(NaN, NaN)
    else
        v2 = v1
    end

    for u in (u0, u1, u2), v in (v0, v1, v2)
        # function isequal(x::Nullable, y::Nullable)
        @test isequal(Nullable(u), Nullable(v)) === isequal(u, v)
        @test isequal(Nullable(u), Nullable(u)) === true
        @test isequal(Nullable(v), Nullable(v)) === true

        @test isequal(Nullable(u), Nullable(v, false)) === false
        @test isequal(Nullable(u, false), Nullable(v)) === false
        @test isequal(Nullable(u, false), Nullable(v, false)) === true

        @test isequal(Nullable(u), Nullable{T}()) === false
        @test isequal(Nullable{S}(), Nullable(v)) === false
        @test isequal(Nullable{S}(), Nullable{T}()) === true

        @test isequal(Nullable(u), Nullable()) === false
        @test isequal(Nullable(), Nullable(v)) === false
        @test isequal(Nullable{S}(), Nullable()) === true
        @test isequal(Nullable(), Nullable{T}()) === true
        @test isequal(Nullable(), Nullable()) === true

        # function isless(x::Nullable, y::Nullable)
        if S <: Real && T <: Real
            @test isless(Nullable(u), Nullable(v)) === isless(u, v)
            @test isless(Nullable(u), Nullable(u)) === false
            @test isless(Nullable(v), Nullable(v)) === false

            @test isless(Nullable(u), Nullable(v, false)) === true
            @test isless(Nullable(u, false), Nullable(v)) === false
            @test isless(Nullable(u, false), Nullable(v, false)) === false

            @test isless(Nullable(u), Nullable{T}()) === true
            @test isless(Nullable{S}(), Nullable(v)) === false
            @test isless(Nullable{S}(), Nullable{T}()) === false

            @test isless(Nullable(u), Nullable()) === true
            @test isless(Nullable(), Nullable(v)) === false
            @test isless(Nullable{S}(), Nullable()) === false
            @test isless(Nullable(), Nullable{T}()) === false
            @test isless(Nullable(), Nullable()) === false
        end
    end
end

# issue #9462
for T in types
    @test isa(convert(Nullable{Number}, Nullable(one(T))), Nullable{Number})
    @test isa(convert(Nullable{Number}, one(T)), Nullable{Number})
    @test isa(convert(Nullable{T}, one(T)), Nullable{T})
    @test isa(convert(Nullable{Any}, Nullable(one(T))), Nullable{Any})
    @test isa(convert(Nullable{Any}, one(T)), Nullable{Any})

    # one(T) is convertible to every type in types
    # let's test that with Nullables
    for S in types
        @test isa(convert(Nullable{T}, one(S)), Nullable{T})
    end
end

@test isnull(convert(Nullable, nothing))
@test isnull(convert(Nullable{Int}, nothing))
@test isa(convert(Nullable{Int}, nothing), Nullable{Int})

@test convert(Nullable, 1) === Nullable(1)
@test convert(Nullable, Nullable(1)) === Nullable(1)
@test isequal(convert(Nullable, "a"), Nullable("a"))
@test isequal(convert(Nullable, Nullable("a")), Nullable("a"))

@test promote_type(Nullable{Int}, Int) === Nullable{Int}
@test promote_type(Nullable{Union{}}, Int) === Nullable{Int}
@test promote_type(Nullable{Float64}, Nullable{Int}) === Nullable{Float64}
@test promote_type(Nullable{Union{}}, Nullable{Int}) === Nullable{Int}
@test promote_type(Nullable{Date}, Nullable{DateTime}) === Nullable{DateTime}

@test Base.promote_op(+, Nullable{Int}, Nullable{Int}) == Nullable{Int}
@test Base.promote_op(-, Nullable{Int}, Nullable{Int}) == Nullable{Int}
@test Base.promote_op(+, Nullable{Float64}, Nullable{Int}) == Nullable{Float64}
@test Base.promote_op(-, Nullable{Float64}, Nullable{Int}) == Nullable{Float64}
@test Base.promote_op(-, Nullable{DateTime}, Nullable{DateTime}) == Nullable{Base.Dates.Millisecond}

# issue #11675
@test repr(Nullable()) == "Nullable{Union{}}()"

# issue #19270
let f19270{S,T}(x::S, y::T) = Base.promote_op(^, S, T)
    @test f19270(Nullable(0.0f0), Nullable(BigInt(0))) == Nullable{Float32}
end

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
module TestMachines

# using Revise
using MLJ
using MLJBase
using Test
using Statistics

@load KNNRegressor

N=50
X = (a=rand(N), b=rand(N), c=rand(N))
y = 2*X.a - X.c + 0.05*rand(N)

train, test = partition(eachindex(y), 0.7);

t = Machine(KNNRegressor(K=4), X, y)
@test_logs (:info, r"Training") fit!(t)
@test_logs (:info, r"Training") fit!(t, rows=train)
@test_logs (:info, r"Not retraining") fit!(t, rows=train)
@test_logs (:info, r"Training") fit!(t)
MLJ.recursive_setproperty!(t, :(model.K),  5)
@test_logs (:info, r"Updating") fit!(t)

predict(t, selectrows(X,test))
@test rms(predict(t, selectrows(X, test)), y[test]) < std(y)

mach = machine(ConstantRegressor(), X, y)
@test_logs (:info, r"Training") fit!(mach)
yhat = predict_mean(mach, X)

n = nrows(X)
rms(yhat, y) ≈ std(y)*sqrt(1 - 1/n)

# test an unsupervised univariate case:
mach = machine(UnivariateStandardizer(), float.(1:5))
@test_logs (:info, r"Training") fit!(mach)
@test isempty(params(mach))

# test a frozen NodalMachine
stand = machine(Standardizer(), source((x1=rand(10),)))
freeze!(stand)
@test_logs (:warn, r"not trained as it is frozen\.$") fit!(stand)

end # module
true


text:
function process_line(line)
    parenthesis_dict = Dict('{'=>'}',
                            '('=>')',
                            '<'=>'>',
                            '['=>']')
    score_dict = Dict(')'=>1,
                      ']'=>2,
                      '}'=>3,
                      '>'=>4)

    open_brackets = []
    for character in line
        if character in keys(parenthesis_dict)
            push!(open_brackets, character)
        elseif character in values(parenthesis_dict)
            if (length(open_brackets)==0) ||
                (character != parenthesis_dict[open_brackets[length(open_brackets)]])
                return -1 # corrupted
            else
                pop!(open_brackets)
            end
        end
    end
    if length(open_brackets) > 0
        # incomplete line - iterate backwards through the list of open brackets
        score = 0
        for bracket in reverse(open_brackets)
            score *= 5
            score += score_dict[parenthesis_dict[bracket]]
        end
        return score
    end
    return 0
end

function test_process_line()
    @assert process_line("[({(<(())[]>[[{[]{<()<>>") == 288957
    @assert process_line("[(()[<>])]({[<{<<[]>>(") == 5566
    @assert process_line("(((({<>}<{<{<>}{[]{[]{}") == 1480781
    @assert process_line("{<[[]]>}<{[{[{[]{()[[[]") == 995444
    @assert process_line("<{([{{}}[<[[[<>{}]]]>[]]") == 294

end

test_process_line()

function process_file(filename)
    input_lines = readlines(filename)
    line_scores = []
    for line in input_lines
        score = process_line(line)
        if score > 0
            push!(line_scores, score)
        end
    end
    line_scores = sort(line_scores)
    return line_scores[Int(ceil(length(line_scores)/2))]
end

score = process_file("input.txt")
println("score is ",score)


text:
using ReachabilityAnalysis

x₀ = Interval(0.9, 1.1)
y₀ = Interval(-0.1, 0.1)
z₀ = Interval(0.9, 1.1)

X₀ = x₀ × y₀ × z₀;

Dict([
    "opts" => Dict(:T=>10.0, :alg=>TMJets(abstol=1e-10, orderT=6, orderQ=2, maxsteps=1_000)),
    "X0"   => X₀,
    "info" => Dict(
        "name" => "steam_governor",
        "dim" => 3,
        "linear" => false,
        "hybrid" => false
    )
])

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
using Revise, Parameters, Logging, Test
using CUDA
includet("C:\\GitHub\\GitHub\\NuclearMedEval\\test\\includeAllUseFullForTest.jl")
using ..CUDAGpuUtils ,..IterationUtils,..ReductionUtils , ..MemoryUtils,..CUDAAtomicUtils
using ..MetadataAnalyzePass,..MetaDataUtils,..WorkQueueUtils,..ProcessMainDataVerB,..HFUtils
using ..MainLoopKernel

##### clear before next dilatation
threads=(32,5);
blocks =1;
mainArrDims= (67,177,90);
dataBdim = (43,21,17);
metaData = MetaDataUtils.allocateMetadata(mainArrDims,dataBdim);
metaDataDims= size(metaData);
resShmem = CUDA.ones(Bool,(dataBdim[1]+2,dataBdim[2]+2,dataBdim[3]+2 ));
sourceShmem = CUDA.ones(Bool,(dataBdim));
resShmemTotalLength = length(resShmem)
sourceShmemTotalLength= length(sourceShmem)
clearIterResShmemLoop= fld(resShmemTotalLength,threads[1]*threads[2])
clearIterSourceShmemLoop= fld(sourceShmemTotalLength,threads[1]*threads[2])
oldBlockCounter= CUDA.zeros(1)
function clearKernel(oldBlockCounter,resShmem,sourceShmem, clearIterResShmemLoop,clearIterSourceShmemLoop,resShmemTotalLength, sourceShmemTotalLength)
    locArr=0
    offsetIter=0
    # localOffset=0

    MainLoopKernel.@clearBeforeNextDilatation( clearIterResShmemLoop,clearIterSourceShmemLoop,resShmemTotalLength, sourceShmemTotalLength)
 return
end

@cuda threads=threads blocks=blocks clearKernel(oldBlockCounter,resShmem,sourceShmem, clearIterResShmemLoop,clearIterSourceShmemLoop,resShmemTotalLength, sourceShmemTotalLength)
@test sum(resShmem)==0
@test sum(sourceShmem)==0

###########iterateOverWorkQueue
using Revise, Parameters, Logging, Test
using CUDA
includet("C:\\GitHub\\GitHub\\NuclearMedEval\\test\\includeAllUseFullForTest.jl")
using ..CUDAGpuUtils ,..IterationUtils,..ReductionUtils , ..MemoryUtils,..CUDAAtomicUtils
using ..MetadataAnalyzePass,..MetaDataUtils,..WorkQueueUtils,..ProcessMainDataVerB,..HFUtils
using ..MainLoopKernel

numberInWorkQueue = 2777
fpTotal=2000
fnTotal=1000
workQueaue = CUDA.zeros(Int64,4,4000)#allocateWorkQueue(fpTotal,fnTotal)
threads=(32,5);
blocks =17;
workQueauecounter= CuArray([numberInWorkQueue])
goldToBeDilatated= CUDA.ones(Bool,1)
segmToBeDilatated= CUDA.ones(Bool,1)
shmemSumLengthMaxDiv4= fld(36*14,4)*4
indexesArr = CUDA.zeros(Int64,1000000)
singleVal = CUDA.zeros(UInt32,1)

for i in 1:length(workQueaue)
    workQueaue[i]= i#using linear index
end
for i in 1 :(fpTotal+fnTotal)
    for j in 1:3
        workQueaue[j,i]= i*j
    end
    workQueaue[4,i]=1 # always gold pass
end
# 70*4
# for i in 1:length(workQueaue)
#     workQueaue[i]= mod(i,255)#using linear index
# end
# for i in 1 :(fpTotal+fnTotal)
#     for j in 1:3
#         workQueaue[j,i]= mod(i*j,255)
#     end
#     workQueaue[4,i]=1 # always gold pass
# end

# for i in 1:length(workQueaue)
#     workQueaue[i]=i#using linear index
# end
# for i in 1 :((fpTotal+fnTotal)*4)
#     workQueaue[i]=i
#     # for j in 1:3
#     #     workQueaue[j,i]=i*j
#     # end
#     # workQueaue[4,i]=1 # always gold pass
# end
Int64.(Array(workQueaue[:,1:8]))

# workQueaue[*4+4]
function iterateOverWorkQueueKernel(indexesArr,workQueauecounter,workQueaue,goldToBeDilatated, segmToBeDilatated,shmemSumLengthMaxDiv4,singleVal)
    shmemSum =  @cuStaticSharedMem(UInt32,(36,14)) # we need this additional spots
    MainLoopKernel.@iterateOverWorkQueue(workQueauecounter,workQueaue,goldToBeDilatated, segmToBeDilatated,shmemSumLengthMaxDiv4,begin

    @ifXY 1 1 CUDA.@atomic singleVal[1]+=1
    end)

 return
end

@cuda threads=threads blocks=blocks iterateOverWorkQueueKernel(indexesArr,workQueauecounter,workQueaue,goldToBeDilatated, segmToBeDilatated,shmemSumLengthMaxDiv4,singleVal)
@test Int64(singleVal[1])==numberInWorkQueue

aa = (cld(33,2))
bb = shmemSumLengthMaxDiv4

fld((cld(33,2)),shmemSumLengthMaxDiv4)

(numberInWorkQueue*4)-sum(indexesArr)

coords = CartesianIndices(indexesArr[1:numberInWorkQueue*4])
aaa = sort(filter(coo-> indexesArr[coo]==0, coords)) #3128


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

using Markdown
using InteractiveUtils

# ╔═╡ c581d17a-c965-11eb-1607-bbeb44933d25
# This file imports an outdated version of PlutoPkgTestA: 0.2.2 (which is stored in the embedded Manifest file).

# It is generated on Julia 1.5 (our oldest supported Julia version, Manifest.toml is not backwards-compatible):

# 1. add our test registry:
#   pkg> registry add https://github.com/JuliaPluto/PlutoPkgTestRegistry

# 2. using Pluto, create this notebook

# 3. use `Pluto.activate_notebook_environment` to change the version of PlutoPkgTestA to 0.2.2

# 4. open the notebook in Pluto again. Add a second package (PlutoPkgTestD) and then remove it again. This adds the auto-generated compat entry for PlutoPkgTestA.

import PlutoPkgTestA

# ╔═╡ aef57966-ea36-478f-8724-e71430f10be9
PlutoPkgTestA.MY_VERSION |> Text

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
PlutoPkgTestA = "419c6f8d-b8cd-4309-abdc-cee491252f94"

[compat]
PlutoPkgTestA = "~0.2.2"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[LibGit2]]
deps = ["Printf"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[Pkg]]
deps = ["Dates", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "UUIDs"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[PlutoPkgTestA]]
deps = ["Pkg"]
git-tree-sha1 = "6c9aa67135641123c559d59ba88e8cb93841773a"
uuid = "419c6f8d-b8cd-4309-abdc-cee491252f94"
version = "0.2.2"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
"""

# ╔═╡ Cell order:
# ╠═c581d17a-c965-11eb-1607-bbeb44933d25
# ╠═aef57966-ea36-478f-8724-e71430f10be9
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
## R code 16.2
m16.1 <- ulam(
    alist(
        w ~ dlnorm( mu , sigma ),
        exp(mu) <- 3.141593 * k * p^2 * h^3,
        p ~ beta( 2 , 18 ),
        k ~ exponential( 0.5 ),
        sigma ~ exponential( 1 )
    ), data=d , chains=4 , cores=4 )

# run the sampler
m16.2 <- stan( model_code=Boxes_model , data=dat_list , chains=3 , cores=3 )

## R code 16.11
dat_list <- list(
    n = as.integer( Panda_nuts$nuts_opened ),
    age = Panda_nuts$age / max(Panda_nuts$age),
    seconds = Panda_nuts$seconds )

m16.4 <- ulam(
    alist(
        n ~ poisson( lambda ),
        lambda <- seconds*phi*(1-exp(-k*age))^theta,
        phi ~ lognormal( log(1) , 0.1 ),
        k ~ lognormal( log(2) , 0.25 ),
        theta ~ lognormal( log(5) , 0.25 )
    ), data=dat_list , chains=4 )

## R code 16.14
sim_lynx_hare <- function( n_steps , init , theta , dt=0.002 ) {
    L <- rep(NA,n_steps)
    H <- rep(NA,n_steps)
    L[1] <- init[1]
    H[1] <- init[2]
    for ( i in 2:n_steps ) {
        H[i] <- H[i-1] + dt*H[i-1]*( theta[1] - theta[2]*L[i-1] )
        L[i] <- L[i-1] + dt*L[i-1]*( theta[3]*H[i-1] - theta[4] )
    }
    return( cbind(L,H) )
}

## R code 16.15
theta <- c( 0.5 , 0.05 , 0.025 , 0.5 )
z <- sim_lynx_hare( 1e4 , as.numeric(Lynx_Hare[1,2:3]) , theta )

## R code 16.18
dat_list <- list(
    N = nrow(Lynx_Hare),
    pelts = Lynx_Hare[,2:3] )

m16.5 <- stan( model_code=Lynx_Hare_model , data=dat_list , chains=3 , cores=3 ,
    control=list( adapt_delta=0.95 ) )


text:
#   This file is part of UnitSystems.jl
#   It is licensed under the MIT license
#   UnitSystems Copyright (C) 2020 Michael Reed
#       _           _                         _
#      | |         | |                       | |
#   ___| |__   __ _| | ___ __ __ ___   ____ _| | __ _
#  / __| '_ \ / _` | |/ / '__/ _` \ \ / / _` | |/ _` |
# | (__| | | | (_| |   <| | | (_| |\ V / (_| | | (_| |
#  \___|_| |_|\__,_|_|\_\_|  \__,_| \_/ \__,_|_|\__,_|
#
#   https://github.com/chakravala
#   https://crucialflow.com

"""
    kilograms(m::Real) = $(lb*g₀/ft)m

Converts mass `m` from slugs to kilogram (kg).
"""
@pure kilograms(m::Number,U::UnitSystem=English) = mass(m,Metric,U)

"""
    slugs(m::Real) = $(ft/lb/g₀)m

Converts mass `m` from kilograms to slugs (slug).
"""
@pure slugs(m::Number,U::UnitSystem=Metric) = mass(m,English,U)

"""
    feet(d) = $(1/ft)d

Converts distance `d` from meters to feet (ft).
"""
@pure feet(d,U::UnitSystem=Metric) = length(d,English,U)

"""
    meters(d) = $(ft)d

Converts distance `d` from feet to meters (m).
"""
@pure meters(d,U::UnitSystem=English) = length(d,Metric,U)

# special exact cases

# IAU to SI
@pure length(U::UnitSystem{kB,ħ,𝘤} where {kB,ħ},S::UnitSystem{kB,ħ,DAY*𝘤/au} where {kB,ħ}) = length(U,S,1/au)
@pure time(U::UnitSystem{kB,ħ,𝘤} where {kB,ħ},S::UnitSystem{kB,ħ,DAY*𝘤/au} where {kB,ħ}) = time(U,S,1/DAY)

# SI to IAU
@pure length(U::UnitSystem{kB,ħ,DAY*𝘤/au} where {kB,ħ},S::UnitSystem{kB,ħ,𝘤} where {kB,ħ}) = length(U,S,au)
@pure time(U::UnitSystem{kB,ħ,DAY*𝘤/au} where {kB,ħ},S::UnitSystem{kB,ħ,𝘤} where {kB,ħ}) = time(U,S,DAY)

# IAU to CGS
@pure length(::UnitSystem{kB,ħ,100𝘤} where {kB,ħ},S::UnitSystem{kB,ħ,DAY*𝘤/au} where {kB,ħ}) = length(U,S,1/au)
@pure time(::UnitSystem{kB,ħ,100𝘤} where {kB,ħ},S::UnitSystem{kB,ħ,DAY*𝘤/au} where {kB,ħ}) = time(U,S,1/DAY)

# CGS to IAU
@pure length(U::UnitSystem{kB,ħ,DAY*𝘤/au} where {kB,ħ},S::UnitSystem{kB,ħ,100𝘤} where {kB,ħ}) = length(U,S,au)
@pure time(U::UnitSystem{kB,ħ,DAY*𝘤/au} where {kB,ħ},S::UnitSystem{kB,ħ,100𝘤} where {kB,ħ}) = time(U,S,DAY)

# IAU to English
@pure length(U::UnitSystem{kB,ħ,𝘤/ft} where {kB,ħ},S::UnitSystem{kB,ħ,DAY*𝘤/au} where {kB,ħ}) = length(U,S,ft/au)
@pure time(U::UnitSystem{kB,ħ,𝘤/ft} where {kB,ħ},S::UnitSystem{kB,ħ,DAY*𝘤/au} where {kB,ħ}) = time(U,S,1/DAY)

# English to IAU
@pure length(U::UnitSystem{kB,ħ,DAY*𝘤/au} where {kB,ħ},S::UnitSystem{kB,ħ,𝘤/ft} where {kB,ħ}) = length(U,S,au/ft)
@pure time(U::UnitSystem{kB,ħ,DAY*𝘤/au} where {kB,ħ},S::UnitSystem{kB,ħ,𝘤/ft} where {kB,ħ}) = time(U,S,DAY)

# base unit conversions

@pure length(U::UnitSystem,S::UnitSystem,l=1) = isquantity(U,S) ? evaldim(length)(U,S) : unit((turn(S)/turn(U))*(planckreduced(S)*electronmass(U)*lightspeed(U)*gravity(S))/(planckreduced(U)*electronmass(S)*lightspeed(S)*gravity(U)),l)
@doc """
$(convertext(:length,"planck(U,S)/mass(U,S)/speed(U,S)"))

Extent of one-dimensional shape or `length` (m), unit conversion factor.

```Julia
julia> length(CGS,Metric) # m⋅cm⁻¹
$(length(CGS,Metric))

julia> length(IAU,Metric) # m⋅au⁻¹
$(length(IAU,Metric))

julia> length(English,Metric) # m⋅ft⁻¹
$(length(English,Metric))

julia> length(Survey,English) # ft⋅ftUS⁻¹
$(length(Survey,English))

julia> length(PlanckGauss,Metric) # m⋅ℓP⁻¹
$(length(PlanckGauss,Metric))
```
""" length, ft, ftUS, L

@pure time(U::UnitSystem,S::UnitSystem,t=1) = isquantity(U,S) ? evaldim(time)(U,S) : unit(length(U,S)/lightspeed(U,S),1)
@doc """
$(convertext(:time,"length(U,S)/speed(U,S)"))

Dimension along which events are ordered or `time` (s), unit conversion factor.

```Julia
julia> time(IAU,Metric) # s⋅day⁻¹
$(time(IAU,Metric))

julia> time(PlanckGauss,Metric) # s⋅tP⁻¹
$(time(PlanckGauss,Metric))
```
""" time(::UnitSystem,::UnitSystem), T

# spacetime

@pure area(U::UnitSystem,S::UnitSystem) = unit(length(U,S)^2)
@pure volume(U::UnitSystem,S::UnitSystem) = unit(length(U,S)^3)
@pure wavenumber(U::UnitSystem,S::UnitSystem) = unit(length(S,U))
@pure angularwavenumber(U::UnitSystem,S::UnitSystem) = unit(angle(U,S)*length(S,U))
@pure fuelefficiency(U::UnitSystem,S::UnitSystem) = area(S,U)
@pure numberdensity(U::UnitSystem,S::UnitSystem) = volume(S,U)
@pure frequency(U::UnitSystem,S::UnitSystem) = time(S,U)
@pure angularfrequency(U::UnitSystem,S::UnitSystem) = unit(angle(U,S)*time(S,U))
@pure frequencydrift(U::UnitSystem,S::UnitSystem) = unit(time(S,U)^2)
@pure speed(U::UnitSystem,S::UnitSystem) = lightspeed(U,S)
@pure acceleration(U::UnitSystem,S::UnitSystem) = unit(speed(U,S)/time(U,S))
@pure jerk(U::UnitSystem,S::UnitSystem) = unit(speed(U,S)/time(U,S)^2)
@pure snap(U::UnitSystem,S::UnitSystem) = unit(speed(U,S)/time(U,S)^3)
@pure crackle(U::UnitSystem,S::UnitSystem) = unit(speed(U,S)/time(U,S)^4)
@pure pop(U::UnitSystem,S::UnitSystem) = unit(speed(U,S)/time(U,S)^5)
@pure volumeflow(U::UnitSystem,S::UnitSystem) = unit(area(U,S)*speed(U,S))
@pure specificenergy(U::UnitSystem,S::UnitSystem) = unit(speed(U,S)^2/gravity(U,S))

# kinematic

@pure inertia(U::UnitSystem,S::UnitSystem) = unit(mass(U,S)/gravity(U,S))
@pure energy(U::UnitSystem,S::UnitSystem) = unit(mass(U,S)*specificenergy(U,S))
@pure power(U::UnitSystem,S::UnitSystem) = unit(energy(U,S)/time(U,S))
@pure force(U::UnitSystem,S::UnitSystem) = unit(inertia(U,S)*acceleration(U,S))
@pure specificforce(U::UnitSystem,S::UnitSystem) = unit(acceleration(U,S)/gravity(U,S))
@pure gravityforce(U::UnitSystem,S::UnitSystem) = unit(gravity(U,S))
@pure pressure(U::UnitSystem,S::UnitSystem) = unit(force(U,S)/area(U,S))

# mechanical

@pure impulse(U::UnitSystem,S::UnitSystem) = unit(force(U,S)*time(U,S))
@pure momentum(U::UnitSystem,S::UnitSystem) = unit(mass(U,S)*speed(U,S))
@pure angularmomentum(U::UnitSystem,S::UnitSystem) = unit(momentum(U,S)*length(U,S)*angle(U,S))
@pure yank(U::UnitSystem,S::UnitSystem) = unit(mass(U,S)*jerk(U,S))
@pure areadensity(U::UnitSystem,S::UnitSystem) = unit(mass(U,S)/area(U,S))
@pure density(U::UnitSystem,S::UnitSystem) = unit(mass(U,S)/volume(U,S))
@pure specificweight(U::UnitSystem,S::UnitSystem) = unit(force(U,S)/volume(U,S))
@pure specificvolume(U::UnitSystem,S::UnitSystem) = unit(volume(U,S)/mass(U,S))
@pure action(U::UnitSystem,S::UnitSystem) = unit(energy(U,S)*time(U,S))
#@pure stiffness(U::UnitSystem,S::UnitSystem) = unit(force(U,S)/length(U,S))
@pure intensity(U::UnitSystem,S::UnitSystem) = unit(power(U,S)/area(U,S))
@pure diffusivity(U::UnitSystem,S::UnitSystem) = unit(speed(U,S)*length(U,S))
@pure viscosity(U::UnitSystem,S::UnitSystem) = unit(force(U,S)/speed(U,S)/length(U,S))
@pure lineardensity(U::UnitSystem,S::UnitSystem) = unit(mass(U,S)/length(U,S))
@pure massflow(U::UnitSystem,S::UnitSystem) = unit(mass(U,S)/time(U,S))
@pure spectralflux(U::UnitSystem,S::UnitSystem) = unit(power(U,S)/length(U,S))
@pure powerdensity(U::UnitSystem,S::UnitSystem) = unit(power(U,S)/volume(U,S))
@pure compressibility(U::UnitSystem,S::UnitSystem) = pressure(S,U)
@pure fluence(U::UnitSystem,S::UnitSystem) = unit(energy(U,S)/area(U,S))
@pure rotationalinertia(U::UnitSystem,S::UnitSystem) = unit(mass(U,S)*area(U,S))

# acoustic

@pure soundexposure(U::UnitSystem,S::UnitSystem) = unit(time(U,S)*pressure(U,S)^2)
@pure specificimpedance(U::UnitSystem,S::UnitSystem) = unit(pressure(U,S)/speed(U,S))
@pure impedance(U::UnitSystem,S::UnitSystem) = unit(specificimpedance(U,S)/area(U,S))
@pure admittance(U::UnitSystem,S::UnitSystem) = unit(area(U,S)/specificimpedance(U,S))
@pure compliance(U::UnitSystem,S::UnitSystem) = unit(time(U,S)^2/mass(U,S))
@pure inertance(U::UnitSystem,S::UnitSystem) = unit(mass(U,S)/length(U,S)^4)

include("kinematicdocs.jl")


text:
using DynamicPolynomials
MP = DynamicPolynomials.MultivariatePolynomials
using GroupsCore
using PermutationGroups

struct VariablePermutation <: SymbolicWedderburn.ByPermutations end

function SymbolicWedderburn.action(a::VariablePermutation, g::PermutationGroups.AbstractPerm, m::Monomial)
    v = variables(m)
    return m(v => SymbolicWedderburn.action(a, g, v))
end

function SymbolicWedderburn.action(::VariablePermutation, g::PermutationGroups.AbstractPerm, v::AbstractVector)
    return map(i -> v[i^g], eachindex(v))
end

abstract type OnMonomials <: SymbolicWedderburn.ByLinearTransformation end

function SymbolicWedderburn.action(
    a::Union{VariablePermutation,OnMonomials},
    el::GroupElement,
    term::MP.AbstractTerm,
)
    return MP.coefficient(term) * SymbolicWedderburn.action(a, el, MP.monomial(term))
end
function SymbolicWedderburn.action(
    a::Union{VariablePermutation,OnMonomials},
    el::GroupElement,
    poly::MP.AbstractPolynomial,
)
    return sum([SymbolicWedderburn.action(a, el, term) for term in MP.terms(poly)])
end

function SymbolicWedderburn.decompose(
    k::MP.AbstractPolynomialLike,
    hom::SymbolicWedderburn.InducedActionHomomorphism,
)
    # correct only if features(hom) == monomials

    indcs = [hom[mono] for mono in MP.monomials(k)]
    coeffs = MP.coefficients(k)

    return indcs, coeffs
end


text:
# Raytracer.jl
# Raytracing for the generation of photorealistic images in Julia
# Copyright (c) 2021 Samuele Colombo, Paolo Galli

# Token implementation

"""
    Keyword

Type wrapping a `Symbol` representing a command or type in a SceneLang script.

# Fields

- `value::Symbol`: the value of the token
"""
struct Keyword
    value::Symbol
end

@enum Command begin
    USING
    SET
    UNSET
    SPAWN
    DUMP
    #construction commands
    LOAD
    ROTATE
    TRANSLATE
    SCALE
    UNITE
    INTERSECT
    DIFF
    FUSE
    TIME
end

@doc """
    Command

Enum type listing all commands of SceneLang.

# Instances

$(join( "- `" .* repr.(instances(Command)) .* "`", "\n"))
""" Command

@enum LiteralType begin
    ColorType
    PointType
    ListType
    TransformationType
    MaterialType
    BrdfType
    PigmentType
    ShapeType
    LightType
    ImageType
    RendererType
    CameraType
    PcgType
    TracerType
end

@doc """
    LiteralType

Enum type listing all main types of SceneLang.

# Instances

$(join( "- `" .* repr.(instances(LiteralType)) .* "`", "\n"))
""" LiteralType

"""
    Identifier

Type wrapping a `Symbol` representing an identifier in a SceneLang script.

# Fields

- `value::Symbol`: the value of the token
"""
struct Identifier
    value::Symbol
end

"""
    LiteralString

Type wrapping a `String` representing a literal string in a SceneLang script.

# Fields

- `value::String`: the value of the token
"""
struct LiteralString
    value::String
end

"""
    LiteralNumber

Type wrapping a `Float32` representing a floating-point number in a SceneLang script.

# Fields

- `value::Float32`: the value of the token
"""
struct LiteralNumber
    value::Float32
end

"""
    LiteralSymbol

Type wrapping a `Symbol` representing a symbol in a SceneLang script.

# Fields

- `value::Symbol`: the value of the token
"""
struct LiteralSymbol
    value::Symbol
end

"""
    MathExpression

Type wrapping a `Expr` representing a mathematical expression in a SceneLang script.

# Fields

- `value::Expr`: the value of the token
"""
struct MathExpression
    value::Expr
end

"""
    StopToken

Convenience empty type to help identify the end of a SceneLang script.

# Fields

- `value::Nothing`: the value of the token
"""
struct StopToken
    value::Nothing

    function StopToken()
        new(nothing)
    end
end

"""
    TokenValue

Union of all types that can be used as token values while interpreting a SceneLang script.

# Types

- [`Keyword`](@ref)
- [`Identifier`](@ref)
- [`LiteralString`](@ref)
- [`LiteralNumber`](@ref)
- [`LiteralSymbol`](@ref)
- [`StopToken`](@ref)

"""
const TokenValue = Union{Keyword, Command, LiteralType, Identifier, MathExpression, LiteralString, LiteralNumber, LiteralSymbol, StopToken}

"""
    Token{T <: TokenValue}

Type representing a language token of a SceneLang script.

# Fields

- `loc::SourceLocation`: a [`SourceLocation`](@ref) representing the position in the script at which the token starts,
- `value::T`: a [`TokenValue`](@ref) representing the value of the token (see [`TokenValue`](@ref))
- `length::Int`: length of the input token.
"""
struct Token{T <: TokenValue}
    loc::SourceLocation
    value::T
    length::Int
end


text:
using Libdl

# Types and global definitions

export Cursor, CXCursor
export CXType

export # TypeKind
    Invalid,
    Unexposed,
    VoidType,
    BoolType,
    Char_U,
    UChar,
    Char16,
    Char32,
    UShort,
    UInt,
    ULong,
    ULongLong,
    #UInt128,
    Char_S,
    SChar,
    WChar,
    Short,
    IntType,
    Long,
    LongLong,
    Int128,
    Float,
    Double,
    LongDouble,
    NullPtr,
    Overload,
    Dependent,
    ObjCId,
    ObjCClass,
    ObjCSel,
    FirstBuiltin,
    LastBuiltin,
    Complex,
    Pointer,
    BlockPointer,
    LValueReference,
    RValueReference,
    Record,
    Enum,
    Typedef,
    ObjCInterface,
    ObjCObjectPointer,
    FunctionNoProto,
    FunctionProto,
    ConstantArray,
    CLVector,
    IncompleteArray

export # CursorKind
    UnexposedDecl,
    StructDecl,
    UnionDecl,
    ClassDecl,
    EnumDecl,
    FieldDecl,
    EnumConstantDecl,
    FunctionDecl,
    VarDecl,
    ParmDecl,
    ObjCInterfaceDecl,
    ObjCCategoryDecl,
    ObjCProtocolDecl,
    ObjCPropertyDecl,
    ObjCIvarDecl,
    ObjCInstanceMethodDecl,
    ObjCClassMethodDecl,
    ObjCImplementationDecl,
    ObjCCategoryImplDecl,
    TypedefDecl,
    CXXMethod,
    Namespace,
    LinkageSpec,
    Constructor,
    Destructor,
    ConversionFunction,
    TemplateTypeParameter,
    NonTypeTemplateParameter,
    TemplateTemplateParameter,
    FunctionTemplate,
    ClassTemplate,
    ClassTemplatePartialSpecialization,
    NamespaceAlias,
    UsingDirective,
    UsingDeclaration,
    TypeAliasDecl,
    ObjCSynthesizeDecl,
    ObjCDynamicDecl,
    CXXAccessSpecifier,
    FirstDecl,
    LastDecl,
    FirstRef,
    ObjCSuperClassRef,
    ObjCProtocolRef,
    ObjCClassRef,
    TypeRef,
    CXXBaseSpecifier,
    TemplateRef,
    NamespaceRef,
    MemberRef,
    LabelRef,
    OverloadedDeclRef,
    VariableRef,
    LastRef,
    FirstInvalid,
    InvalidFile,
    NoDeclFound,
    NotImplemented,
    InvalidCode,
    LastInvalid,
    FirstExpr,
    UnexposedExpr,
    DeclRefExpr,
    MemberRefExpr,
    CallExpr,
    ObjCMessageExpr,
    BlockExpr,
    IntegerLiteral,
    FloatingLiteral,
    ImaginaryLiteral,
    StringLiteral,
    CharacterLiteral,
    ParenExpr,
    UnaryOperator,
    ArraySubscriptExpr,
    BinaryOperator,
    CompoundAssignOperator,
    ConditionalOperator,
    CStyleCastExpr,
    CompoundLiteralExpr,
    InitListExpr,
    AddrLabelExpr,
    StmtExpr,
    GenericSelectionExpr,
    GNUNullExpr,
    CXXStaticCastExpr,
    CXXDynamicCastExpr,
    CXXReinterpretCastExpr,
    CXXConstCastExpr,
    CXXFunctionalCastExpr,
    CXXTypeidExpr,
    CXXBoolLiteralExpr,
    CXXNullPtrLiteralExpr,
    CXXThisExpr,
    CXXThrowExpr,
    CXXNewExpr,
    CXXDeleteExpr,
    UnaryExpr,
    ObjCStringLiteral,
    ObjCEncodeExpr,
    ObjCSelectorExpr,
    ObjCProtocolExpr,
    ObjCBridgedCastExpr,
    PackExpansionExpr,
    SizeOfPackExpr,
    LambdaExpr,
    ObjCBoolLiteralExpr,
    ObjCSelfExpr,
    LastExpr,
    FirstStmt,
    UnexposedStmt,
    LabelStmt,
    CompoundStmt,
    CaseStmt,
    DefaultStmt,
    IfStmt,
    SwitchStmt,
    WhileStmt,
    DoStmt,
    ForStmt,
    GotoStmt,
    IndirectGotoStmt,
    ContinueStmt,
    BreakStmt,
    ReturnStmt,
    GCCAsmStmt,
    AsmStmt,
    ObjCAtTryStmt,
    ObjCAtCatchStmt,
    ObjCAtFinallyStmt,
    ObjCAtThrowStmt,
    ObjCAtSynchronizedStmt,
    ObjCAutoreleasePoolStmt,
    ObjCForCollectionStmt,
    CXXCatchStmt,
    CXXTryStmt,
    CXXForRangeStmt,
    SEHTryStmt,
    SEHExceptStmt,
    SEHFinallyStmt,
    MSAsmStmt,
    NullStmt,
    DeclStmt,
    LastStmt,
    TranslationUnit,
    FirstAttr,
    UnexposedAttr,
    IBActionAttr,
    IBOutletAttr,
    IBOutletCollectionAttr,
    CXXFinalAttr,
    CXXOverrideAttr,
    AnnotateAttr,
    AsmLabelAttr,
    LastAttr,
    PreprocessingDirective,
    MacroDefinition,
    MacroExpansion,
    MacroInstantiation,
    InclusionDirective,
    FirstPreprocessing,
    LastPreprocessing,
    ModuleImportDecl,
    FirstExtraDecl,
    LastExtraDecl

# Type definitions for wrapped types
const CXIndex = Ptr{Cvoid}
const CXUnsavedFile = Ptr{Cvoid}
const CXFile = Ptr{Cvoid}
const CXTypeKind = Int32
const CXCursorKind = Int32
const CXTranslationUnit = Ptr{Cvoid}

###############################################################################
# Container types
#   for now we use these as the element type of the target array
###############################################################################

struct CXSourceLocation
    ptr_data::NTuple{2, Csize_t}
    int_data::Cuint
    CXSourceLocation() = new(0,0)
end

struct CXSourceRange
    ptr_data::NTuple{2, Csize_t}
    begin_int_data::Cuint
    end_int_data::Cuint
    CXSourceRange() = new(0,0,0)
end

struct CXTUResourceUsageEntry
    kind::Cint
    amount::Culong
    CXTUResourceUsageEntry() = new(0,0)
end

struct CXTUResourceUsage
    data::Nothing
    numEntries::Cuint
    entries::Ptr{Cptrdiff_t}
    CXTUResourceUsage() = new(0,0,0)
end

# Generate container types
#for st in Any[
#        :CXTUResourceUsageEntry, :CXTUResourceUsage ]
#    sz_name = Symbol(st, "_size")
#    st_base = Symbol("_", st)
#    @eval begin
#        const $sz_name = get_sz($("$st"))
#        struct $(st)
#            data::Array{$st_base,1}
#            $st(d) = new(d)
#        end
#        $st() = $st(Array($st_base, 1))
#        $st(d::$st_base) = $st([d])
#    end
#end

struct CXCursor
    kind::Cint
    xdata::Cint
    data::NTuple{3, Csize_t}
    CXCursor() = new(0,0,0)
end

struct CXType
    kind::Int32
    data::NTuple{2, Csize_t}
    CXType() = new(0,(0,0))
end

struct CXString
    data::Ptr{UInt8}
    private_flags::Cuint
    CXString() = new(C_NULL,0)
end

function Base.convert(::Type{String}, x::CXString)
    ret = unsafe_string(x.data)
    ccall( (:clang_disposeString, libclang), Cvoid, (CXString,), x)
    return ret
end

###############################################################################
# Set up CXToken wrapping
#   Note that this low-level interface should be used with care.
#   Accessing a CXToken after the originating TokenList has been
#   finalized is undefined and will most likely segfault.
#   The high-level TokenList[] and CLToken interface is preferred.
###############################################################################

struct CXToken
    int_data::NTuple{4, Cuint}
    ptr_data::Csize_t
    CXToken() = new((0,0,0,0),0)
end

struct TokenList
    ptr::Ptr{CXToken}
    size::Cuint
    tunit::CXTranslationUnit
end

abstract type CLToken end
for sym in names(TokenKind, all=true)
    if(sym == :TokenKind) continue end
    @eval begin
        struct $sym <: CLToken
            text::String
        end
    end
end

finalizer(l::TokenList) = cindex.disposeTokens(l)

###############################################################################
# Set up CXCursor wrapping
#
#   Each CXCursorKind enum gets a specific Julia type wrapping
#   so that we can dispatch directly on node kinds.
###############################################################################

abstract type CLCursor end
const CXCursorMap = Dict{Int32,Any}()

# create the typed wrapper cursors
for sym in names(CursorKind, all=true)
    sym == :CursorKind && continue
    rval = getfield(CursorKind, sym)
    @eval begin
        struct $(sym) <: CLCursor
            cursor::CXCursor
        end
        CXCursorMap[Int32($rval)] = $sym
    end
end

CLCursor(c::CXCursor) = CXCursorMap[c.kind](c)
Base.convert(::Type{CXCursor}, x::T) where {T<:CLCursor} = x.cursor
Base.convert(::Type{CLCursor}, x::CXCursor) = CLCursor(x)

###############################################################################
# Set up CXType wrapping
#
#   Each CXTypeKind enum gets a specific Julia type wrapping
#   so that we can dispatch directly on node kinds.
###############################################################################

abstract type CLType end
CLTypeMap = Dict{Int32,Any}()

for sym in names(TypeKind, all=true)
    sym == :TypeKind && continue
    rval = getfield(TypeKind, sym)
    @eval begin
        struct $(sym) <: CLType
            typ::CXType
        end
        CLTypeMap[Int32($rval)] = $sym;
    end
end

CLType(t::CXType) = CLTypeMap[t.kind](t)
Base.convert(::Type{CXType}, x::T) where {T<:CLType} = x.typ
Base.convert(::Type{CLType}, x::CXType) = CLType(x)

# Duplicates
const LastDecl = CXXAccessSpecifier
const ObjCSuperClassRef = FirstRef
const LastRef = VariableRef
const InvalidFile = FirstInvalid
const LastInvalid=InvalidCode
const LastExpr = ObjCSelfExpr
const AsmStmt = GCCAsmStmt
const UnexposedAttr = FirstAttr
const LastPreprocessing = InclusionDirective


text:
"""
    KLocalOperatorZero

A simple KLocalOperator representing zero, but holding informations
about the hilbert basis.
"""
struct KLocalOperatorZero{H<:AbstractHilbert,T} <: AbsLinearOp{T}
    hilb::H
end

KLocalOperatorZero(T::Type{<:Number}, hilb) =
    KLocalOperatorZero{typeof(hilb), T}(hilb)

KLocalOperatorZero(hilb) = KLocalOperatorZero(ComplexF64, hilb)

LocalOperator(hilb) = KLocalOperatorZero(hilb)
LocalOperator(T, hilb) = KLocalOperatorZero(T, hilb)

export LocalOperator

sites(op::KLocalOperatorZero) = []

QuantumOpticsBase.basis(op::KLocalOperatorZero) = op.hilb

# Copy
duplicate(op::KLocalOperatorZero) = copy(op)
Base.copy(op::KLocalOperatorZero) = KLocalOperatorZero(basis(op))

##
function _row_valdiff!(conn::OpConnection, op::KLocalOperatorZero, v) end
function row_valdiff_index!(conn::OpConnectionIndex, op::KLocalOperatorZero, v) end
function accumulate_connections!(acc::AbstractAccumulator, op::KLocalOperatorZero, v) end
function map_connections(fun::Function , op::KLocalOperatorZero, v) end

Base.conj(op::KLocalOperatorZero) = copy(op)
Base.adjoint(op::KLocalOperatorZero) = copy(op)
Base.transpose(op::KLocalOperatorZero) = copy(op)

Base.:+(op_l::KLocalOperatorZero, op_r::KLocalOperator) = duplicate(op_r)
Base.:+(op_l::KLocalOperator, op_r::KLocalOperatorZero) = duplicate(op_l)

Base.:-(op_l::KLocalOperatorZero, op_r::KLocalOperator) = -op_r
Base.:-(op_l::KLocalOperator, op_r::KLocalOperatorZero) = -op_l

_sum_samesite(op_l::KLocalOperatorZero,
                op_r::AbsLinearOperator) = duplicate(opop_r_r)

_sum_samesite(op_l::AbsLinearOperator,
              op_r::KLocalOperatorZero) = duplicate(op_l)

Base.:*(a::Number, b::KLocalOperatorZero) = b

Base.:*(opl::KLocalOperatorZero, opr::KLocalOperator) =
    return opl
Base.:*(opl::KLocalOperator, opr::KLocalOperatorZero) =
    return opr

Base.:(==)(l::KLocalOperatorZero, r::KLocalOperatorZero) =
    l.hilb == r.hilb

Base.isapprox(l::KLocalOperatorZero, r::KLocalOperatorZero; kwargs...) =
    l.hilb == r.hilb

Base.show(io::IO, op::KLocalOperatorZero) = begin
    hilb   = basis(op)

    print(io, "KLocalOperator($hilb)")
end

Base.show(io::IO, m::MIME"text/plain", op::KLocalOperatorZero) = begin
    hilb   = basis(op)

    print(io, "empty KLocalOperator on space:  $hilb")
end


text:
module LibClang

using Clang_jll
export Clang_jll

using CEnum

include("ctypes.jl")
export Ctm, Ctime_t, Cclock_t

include(joinpath(@__DIR__, "..", "gen", "libclang_common.jl"))
include(joinpath(@__DIR__, "..", "gen", "libclang_api.jl"))

foreach(names(@__MODULE__, all=true)) do s
    if startswith(string(s), "CX") || startswith(string(s), "clang_")
        @eval export $s
    end
end

end # module


text:
using Plots, CSV
import DataFrames: DataFrame
using Glob
import MyterialColors: blue_grey_dark, blue_grey
import Colors: HSL
using NaturalSort

import jcontrol: Bicycle, FULLTRACK, naturalsort, Solution
import jcontrol.io: PATHS, load_cached_trials
using jcontrol.visuals
import jcontrol: trimtrial
import jcontrol.forwardmodel: trimsolution

S0 = 0
S1 = 260

trials = map(t -> trimtrial(t, S0, S1), load_cached_trials(; keep_n = nothing))

bike = Bicycle()
globalsolution = DataFrame(CSV.File(joinpath(PATHS["horizons_sims_cache"], "global_solution.csv")))

globalsolution = trimsolution(Solution(globalsolution), S0, S1)

# initialize plots with global solution
xyplot = draw(:arena, legend=:bottomleft)
plot_bike_trajectory!(globalsolution, bike; showbike=false, color=blue_grey_dark, lw=6, alpha=.8, label="global solution")

uplot = plot(;  xlim=[0, 200], ylim=[0, 100], ylabel="long.speed. (cm/s)")
progressplot = plot(; ylabel="time (s)", xlabel="s (cm)", legend=:bottomright)

# plot trials kinematics
for trial in trials
    plot!(uplot, trial.s, trial.u; color="black", alpha=0.015, label=nothing)
end

# plot global solution kinematics
global_kwargs = Dict(:label=>nothing, :color=>blue_grey_dark, :lw=>6, :alpha=>.4)
plot!(uplot, globalsolution.s, globalsolution.u; global_kwargs...)
plot!(progressplot, globalsolution.s, globalsolution.t; global_kwargs...)

# durations histogram
h = histogram(map(t->t.duration, trials), label="trial durations", color="black", alpha=.2, xlim=[0, 15])
plot!(h, [globalsolution.t[end], globalsolution.t[end]], [0, 150], color=blue_grey_dark, lw=6, alpha=1, label="global solution")

# plot simulations
files = sort(glob("multiple_horizons_mtm_horizon_length*.csv", PATHS["horizons_sims_cache"]), lt=natural)
println.(files)
colors = range(HSL(326, .9, .68), stop=HSL(212, .9, .6), length=max(2, length(files)))
alphas = range(.9, .6, length=max(5, length(files)))
for (file, color, alpha) in zip(files, colors, alphas)
    data = trimsolution(Solution(DataFrame(CSV.File(file))), S0, S1)
    name = "horizon length: " * (split(file, "_")[end][1:end-4])

    plot_kwargs = Dict(:lw=>4, :label=>name, :color=>color, :alpha=>1)

    # plot trajectory
    plot!(xyplot, data.x, data.y, lw=4, color=color, label=name)

    # plot kinematics
    plot!(uplot, data.s, data.u; plot_kwargs...)
    scatter!(uplot, [data.s[end]], [data.u[end]]; markersize=10, color=color, label=nothing)
    plot!(progressplot, data.s, data.t; plot_kwargs...)
    plot!(h, [data.t[end], data.t[end]], [0, 100]; plot_kwargs...)
end

# layout and plot
l = @layout [
        b{0.4h}
        c{0.4h}
        d{0.2h}
]

# display(plot(xyplot, size=(800, 1200) ))
display(
    plot(
        uplot,
        progressplot,
        h,
        layout=l,
        size=(1400, 1200),
    ),
)

max_stars_repo_name: FedeClaudi/LocomotionControl

max_issues_repo_name: FedeClaudi/LocomotionControl

max_forks_repo_name: FedeClaudi/LocomotionControl


text:
#  ------------------------------------------------------------------------------------------
#  Copyright (c) Microsoft Corporation. All rights reserved.
#  Licensed under the MIT License (MIT). See LICENSE in the repo root for license information.
#  ------------------------------------------------------------------------------------------

###################################################################################################
# Fetch methods

# Fetch output blobs for given task index (blocking)
"""
    output = fetch(object::Union{BatchController, BlobFuture, BatchFuture}; destroy_blob=false, timeout=60, task_timeout=60)

 Fetch the output from the batch job or future (i.e. the return arguments of the executed function).

    output = fetch(batch_controller::BatchController, idx; destroy_blob=false, timeout=60, wait_for_completion=true)

 Fetch the output of task `idx` only (only supported for inputs of type `BatchController`).

 *Input*:

 - `object` (BatchController, BlobFuture or BatchFuture): Batch control structure or future whose reference will be fetched.

 - `idx` (Integer): Task number

 - `destroy_blob` (Bool): Remove the original blob upon successful fetching (default is `false`)

 - `timeout`: Maximum time to wait for all tasks to finish (in minutes). Default is `60`.

 - `task_timeout`: Maximum runtime per task (in minutes). Default is `60`.

 - `num_restart`: Allowed retries for failed tasks (default is `0`).

 *Output*:

 - `output`: List of return arguments of executed function (one entry per task) or data from blob/batch future.

"""
function fetch(batch_controller::BatchController, idx; destroy_blob=false, timeout=60, wait_for_completion=true, num_restart=0)

    # Wait for specified task to finish
    task_id = batch_controller.task_id[idx]["taskname"]
    pool_no = batch_controller.task_id[idx]["pool"]
    if wait_for_completion
        wait_for_task_to_complete(batch_controller.batch_client[pool_no], batch_controller.job_id[pool_no], task_id, timeout;
            verbose=__verbose__, num_restart=num_restart)
    end

    # Loop over entries in Future for i-th task
    num_files = length(batch_controller.output[idx].blob.name)
    out_files = []
    for blob in batch_controller.output[idx].blob.name

        # Fetch blob and add to collection
        if ~isnothing(batch_controller.blob_client[pool_no])
            try
                val = batch_controller.blob_client[pool_no].get_blob_to_bytes(batch_controller.blobcontainer, blob)
                push!(out_files, deserialize(IOBuffer(val.content)))
                # Delete blob (default is true)
                destroy_blob && batch_controller.blob_client[pool_no].delete_blob(batch_controller.blobcontainer, blob)
                println("correct blob, ", blob)
            catch e
                println("wrong idx is, ", idx)
                println("wrong blob, ", blob)
                println("Batch controller task id taskname, ", batch_controller.task_id[idx]["taskname"])
                println("Batch controller task id pool, ", batch_controller.task_id[idx]["pool"])
                println("Batch controller blob container, ", batch_controller.blobcontainer)
                println("Batch controller output, ", batch_controller.output)
                @warn "Blob does not exist or task has not finished yet. Return nothing."
                println(e)
                push!(out_files, nothing)
            end
        end
    end

    # Replace blob name in Future with value
    if num_files > 1
        return tuple(out_files...)
    elseif length(out_files) == 1
        return out_files[1]
    else
        return nothing
    end
end

"""
    data = fetch!(object::Union{BatchController, BlobFuture, BatchFuture}; destroy_blob=true, timeout=60, task_timeout=60)

 Fetch the output from the batch job or future in place (return output arguments of executed function and overwrite blob future).

    data = fetch!(batch_controller::BatchController, idx; destroy_blob=false, timeout=60, wait_for_completion=true)

 Inplace fetch of the output of task `idx` only (only supported for inputs of type `BatchController`).

 *Input*:

 - `object` (BatchController, BlobFuture or BatchFuture): Batch control structure or future whose reference will be fetched.

 - `idx` (Integer): Task number

 - `destroy_blob` (Bool): Remove the original blob upon successful fetching (default is `true` for inplace fetch)

 - `timeout`: Maximum time to wait for all tasks to finish (in minutes). Default is `60`.

 - `task_timeout`: Maximum runtime per task (in minutes). Default is `60`.

 - `num_restart`: Allowed retries for failed tasks (default is `0`).

 *Output*:

 - `output`: List of return arguments of executed function (one entry per task) or data from blob/batch future.

"""
function fetch!(batch_controller::BatchController, idx; destroy_blob=true, timeout=60, wait_for_completion=true, num_restart=0)

    # Wait for specified task to finish
    task_id = batch_controller.task_id[idx]["taskname"]
    pool_no = batch_controller.task_id[idx]["pool"]
    if wait_for_completion
        wait_for_task_to_complete(batch_controller.batch_client[pool_no], batch_controller.job_id[pool_no], task_id, timeout;
            verbose=__verbose__, num_restart=num_restart)
    end

    # Loop over entries in Future for i-th task
    num_files = length(batch_controller.output[idx].blob.name)
    out_files = []
    for blob in batch_controller.output[idx].blob.name

        # Fetch blob and add to collection
        if ~isnothing(batch_controller.blob_client[pool_no])
            try
                val = batch_controller.blob_client[pool_no].get_blob_to_bytes(batch_controller.blobcontainer, blob)
                push!(out_files, deserialize(IOBuffer(val.content)))

                # Delete blob (default is true)
                destroy_blob && batch_controller.blob_client[pool_no].delete_blob(batch_controller.blobcontainer, blob)
                println("correct blob, ", blob)
            catch e
                println("wrong idx is, ", idx)
                println("wrong blob, ", blob)
                println("Batch controller task id taskname, ", batch_controller.task_id[idx]["taskname"])
                println("Batch controller task id pool, ", batch_controller.task_id[idx]["pool"])
                println("Batch controller blob container, ", batch_controller.blobcontainer)
                println("Batch controller output, ", batch_controller.output)
                @warn "Blob does not exist or task has not finished yet. Return nothing."
                println(e)
                push!(out_files, nothing)
            end
        end
    end

    # Replace blob name in Future with value
    if num_files > 1 && length(out_files) > 1
        batch_controller.output[idx].blob = tuple(out_files...)
    elseif length(out_files) == 1
        batch_controller.output[idx].blob = out_files[1]
    else
        batch_controller.output[idx].blob = nothing
    end

    # Return value
    return batch_controller.output[idx].blob
end

# Fetch output blobs of all tasks in order of completion (blocking operation)
function fetch(batch_controller::BatchController; destroy_blob=false, timeout=60, task_timeout=60, num_restart=0)

    out_files = Array{Any}(undef, length(batch_controller.output))
    remaining_tasks = deepcopy(batch_controller.task_id)
    task_id = nothing
    __verbose__ && print("Monitoring tasks for 'Completed' state, timeout in $timeout minutes ...")
    while true

        # Wait for one task from task list to finish
        try
            task_id = wait_for_one_task_from_multi_pool(batch_controller.batch_client, batch_controller.job_id, remaining_tasks;
                task_timeout=task_timeout, fetch_timeout=timeout, verbose=__verbose__, num_restart=num_restart)[1]
        catch
           throw("Reached timeout for task completion.")
        end
        __verbose__ && print("\nFetch output from task $task_id")

        # Fetch its output
        task_no = findall(i -> i["taskname"] == task_id, batch_controller.task_id)[1]
        if ~isempty(remaining_tasks)
            out_files[task_no] = fetch(batch_controller, task_no; destroy_blob=destroy_blob, timeout=task_timeout, wait_for_completion=false,
                num_restart=0)
            local_id = findall(i -> i["taskname"] == task_id, remaining_tasks)[1]
            popat!(remaining_tasks, local_id)
        end

        # Return it task list is empty
        if isempty(remaining_tasks)
            __verbose__ && print("\n")
            if length(out_files) > 1
                return out_files
            elseif length(out_files) == 1
                return out_files[1]
            else
                return nothing
            end
        end
    end
end

function fetch!(batch_controller::BatchController; destroy_blob=true, timeout=60, task_timeout=60, num_restart=0)

    out_files = Array{Any}(undef, length(batch_controller.output))
    remaining_tasks = deepcopy(batch_controller.task_id)
    task_id = nothing
    __verbose__ && print("Monitoring tasks for 'Completed' state, timeout in $timeout minutes ...")
    while true

        # Wait for one task from task list to finish
        try
            task_id = wait_for_one_task_from_multi_pool(batch_controller.batch_client, batch_controller.job_id, remaining_tasks;
                task_timeout=task_timeout, fetch_timeout=timeout, verbose=__verbose__, num_restart=num_restart)[1]
        catch
           throw("Reached timeout for task completion.")
        end
        __verbose__ && print("\nFetch output from task $task_id")

        # Fetch its output
        task_no = findall(i -> i["taskname"] == task_id, batch_controller.task_id)[1]
        if ~isempty(remaining_tasks)
            out_files[task_no] = fetch!(batch_controller, task_no; destroy_blob=destroy_blob, timeout=timeout, wait_for_completion=false)
            local_id = findall(i -> i["taskname"] == task_id, remaining_tasks)[1]
            popat!(remaining_tasks, local_id)
        end

        # Return it task list is empty
        if isempty(remaining_tasks)
            __verbose__ && print("\n")
            if length(out_files) > 1
                return out_files
            else
                return out_files[1]
            end
        end
    end
end

# Reduction code for remote execution
reduction_code = quote
    @batchdef function remote_reduction(_x, _y; op=+)
        x = fetch(_x)
        y = fetch(_y)
        output = broadcast(op, x, y)
        return output
    end
end

function check_for_reduction_code_in_ast(expr, count)

    # Reached AST leaf
    typeof(expr) != Expr && return

    if expr.head == :call && expr.args[1] == :remote_reduction
        count[1] += 1
    else
        # Step through AST
        for i=1:length(expr.args)
            if typeof(expr.args[i]) == Expr
                check_for_reduction_code_in_ast(expr.args[i], count)
            end
        end
    end
end

"""
    output = fetchreduce(batch_controller::BatchController; op=+, destroy_blob=false, timeout=60, remote=false)

 Fetch the output from the batch job and apply the specified reduction operation to it (across tasks).

 *Input*:

 - `batch_controller` (BatchController): Batch control structure

 - `op`: Algebraic operation to apply to the output from different tasks.

 - `destroy_blob` (Bool): Remove the original blob upon successful fetching (default is `false`)

 - `timeout`: Timeout in minutes.

 - `num_restart`: Allowed retries for failed tasks (default is `0`).

 - `remote` (Bool): If `true`, execute the reduction operation as additional batch tasks. Otherwise, the reduction happens locally.

 *Output*:

 - `output`: Return argument(s) of executed function after application of reduction operation (along tasks)

"""
function fetchreduce(batch_controller::BatchController; op=+, destroy_blob=false, timeout=60, task_timeout=60,
    remote=false, reduction_code=reduction_code, num_restart=0)

    if remote == false
        return fetchreduce_local(batch_controller; op=op, destroy_blob=destroy_blob, timeout=timeout,
            task_timeout=task_timeout, num_restart=num_restart)
    else
        if batch_controller.num_tasks == 1
            return fetch(batch_controller; destroy_blob=destroy_blob, timeout=timeout)
        else
            num_pools = length(batch_controller.batch_client)
            job_ids = Array{Any}(undef, num_pools)
            task_ids = Array{Any}(undef, num_pools)
            output = Array{Any}(undef, num_pools)
            orig_blobs = Array{Any}(undef, num_pools)
            temp_blobs = Array{Any}(undef, num_pools)

            count_code = [0]
            check_for_reduction_code_in_ast(__expressions__, count_code)
            count_code[1] == 0 && eval(reduction_code)

            # First, reduce all outputs within each storage account
            @sync begin
                for i=1:num_pools
                    @async job_ids[i], task_ids[i], output[i], orig_blobs[i], temp_blobs[i] =
                        fetchreduce_remote(
                            batch_controller.batch_client[i],
                            batch_controller.blob_client[i],
                            batch_controller.job_id[i],
                            batch_controller.task_id[findall(k -> k["pool"] == i, batch_controller.task_id)],
                            batch_controller.output[findall(k -> k.client_index == i, batch_controller.output)],
                            batch_controller.blobcontainer,
                            i;
                            op=op, destroy_blob=destroy_blob, timeout=timeout, task_timeout=task_timeout,
                            num_restart=num_restart)
                end
            end

            # Batch controller for remaining tasks (1 per storage account)
            reduce_ctrl = BatchController(job_ids, task_ids, num_pools, output)
            output = fetchreduce_local(reduce_ctrl; op=op, destroy_blob=true, timeout=timeout, task_timeout=task_timeout,
                num_restart=num_restart)

            # Clean up
            for i=1:num_pools
                job_ids[i] != batch_controller.job_id[i] && (batch_controller.batch_client[i].job.delete(job_ids[i]))
                # Output blobs
                if destroy_blob
                    for orig_file in orig_blobs[i]
                        batch_controller.blob_client[i].delete_blob(batch_controller.blobcontainer, orig_file)
                    end
                end
                # Temp blobs
                for temp_file in temp_blobs[i]
                    batch_controller.blob_client[i].delete_blob(batch_controller.blobcontainer, temp_file)
                end
            end
            return output
        end
    end
end

# Fetch and reduce operation locally
function fetchreduce_local(batch_controller::BatchController; op=+, destroy_blob=false, timeout=60, task_timeout=60, num_restart=0)

    # Single task: simply fetch output
    if batch_controller.num_tasks == 1
        output = fetch(batch_controller)

    # Multiple tasks: fetch and sum output
    else
        remaining_tasks = deepcopy(batch_controller.task_id)
        task_id = nothing
        output = nothing

        __verbose__ && print("Monitoring tasks for 'Completed' state, timeout in $timeout minutes ...")
        while true

            # Wait for one task from task list to finish
            try
                task_id = wait_for_one_task_from_multi_pool(batch_controller.batch_client, batch_controller.job_id, remaining_tasks;
                    task_timeout=task_timeout, fetch_timeout=timeout, verbose=__verbose__, num_restart=num_restart)[1]
            catch
                throw("Reached timeout for task completion.")
            end
            __verbose__ && print("\nFetch output from task $task_id")

            # Fetch its output
            task_no = findall(i -> i["taskname"] == task_id, batch_controller.task_id)[1]
            if ~isempty(remaining_tasks)

                # Fetch output
                if isnothing(output)
                    output = fetch(batch_controller, task_no; destroy_blob=destroy_blob, timeout=timeout, wait_for_completion=false)
                else
                    output = broadcast(op, output, fetch(batch_controller, task_no; destroy_blob=destroy_blob,
                        timeout=timeout, wait_for_completion=false))
                end

                # Remove task from task list
                local_id = findall(i -> i["taskname"] == task_id, remaining_tasks)[1]
                popat!(remaining_tasks, local_id)
            end

            # Return it task list is empty
            if isempty(remaining_tasks)
                __verbose__ && print("\n")
                if length(output) > 1
                    return output
                elseif length(output) == 1
                    return output[1]
                else
                    return nothing
                end
            end
        end
    end
end

function wait_for_task_from_job_list!(batch_client, output, job_ids, remaining_tasks, timeout, task_timeout,
    temp_blobs, orig_blobs, original_job_id; num_restart=0)

    # Wait for task to finish
    task_id = wait_for_one_task_from_multi_jobs(batch_client, job_ids, remaining_tasks;
    task_timeout=task_timeout, fetch_timeout=timeout, verbose=__verbose__, num_restart=num_restart)[1]
    local_id = findall(i -> i["taskname"] == task_id, remaining_tasks)[1]

    # Remove completed task from task/job/output list
    popat!(remaining_tasks, local_id)
    completed_job_id = popat!(job_ids, local_id)
    output_ref = popat!(output, local_id)

    # Keep track of temporary blobs
    if completed_job_id != original_job_id
        push!(temp_blobs, output_ref.blob.name[1])
        batch_client.job.delete(completed_job_id)
    else
        push!(orig_blobs, output_ref.blob.name[1])
    end
    return output_ref
end

# Execute reduce function remotely
function fetchreduce_remote(batch_client, blob_client, job_id, task_id, func_output, blobcontainer, pool_no;
    op=+, destroy_blob=false, timeout=60, task_timeout=60, num_restart=0)

    # Single task: simply fetch output
    num_tasks = length(func_output)
    if num_tasks == 1
        return job_id, task_id[1], func_output[1], [], []
    else
        output = deepcopy(func_output)
        remaining_tasks = deepcopy(task_id)
        original_job_id = job_id
        job_ids = Array{Any}(undef, length(remaining_tasks)) .= job_id
        temp_blobs = []; orig_blobs = []

        __verbose__ && print("Monitoring tasks for 'Completed' state, timeout in $timeout minutes ...")
        while true

            # Wait for two tasks to finish
            output_ref1 = wait_for_task_from_job_list!(batch_client, output, job_ids, remaining_tasks,
                timeout, task_timeout, temp_blobs, orig_blobs, original_job_id, num_restart=num_restart)
            output_ref2 = wait_for_task_from_job_list!(batch_client, output, job_ids, remaining_tasks,
                timeout, task_timeout, temp_blobs, orig_blobs, original_job_id, num_restart=num_restart)

            # Submit summation as new batch job (non-blocking) and add to pending tasks/jobs
            jt_id = objectid([output_ref1,output_ref2])
            job_name = join(["batch_reduce_", jt_id])
            task_name = join(["task_", jt_id])
            bctrl = @batchexec(remote_reduction(output_ref1, output_ref2), Options(job_name=job_name,
                task_name=task_name, priority=100, pool=pool_no, reset_mpi=true))
            push!(remaining_tasks, bctrl.task_id[1]); push!(job_ids, bctrl.job_id[1]); push!(output, bctrl.output[1])

            # If only task is left -> fetch it, clean up and return it
            if length(remaining_tasks) == 1
                return job_ids[1], remaining_tasks[1], output[1], orig_blobs, temp_blobs
            end
        end
    end
end

# Inplace Fetch-reduce for output tuple (one temp. copy of output argument)
"""
    output = fetchreduce!(batch_controller::BatchController, output; op=+, destroy_blob=false, timeout=60)

 Fetch the output from the batch job in-place and apply the specified reduction operation to
 the provided `output` (across tasks).

 *Input*:

 - `batch_controller` (BatchController): Batch control structure

 - `output`: Output arguments of executed batch function which will be overwritten.

 - `op`: Algebraic operation to apply to the output from different tasks.

 - `destroy_blob` (Bool): Remove the original blob upon successful fetching (default is `false`)

 - `timeout`: Timeout in minutes.

 - `num_restart`: Allowed retries for failed tasks (default is `0`).

 *Output*:

 - `output`: Return argument(s) of executed function after application of reduction operation (along tasks)

"""
function fetchreduce!(batch_controller::BatchController, output::Tuple; op=+, destroy_blob=true, timeout=60,
    task_timeout=60, num_restart=0)

    # Single task: simply fetch output
    if batch_controller.num_tasks == 1
        temp = fetch!(batch_controller)
        for (i, entry) in enumerate(output)
            if ~isnothing(entry)
                entry[:] = temp[i]
            end
        end

    # Multiple tasks: fetch and sum output
    else

        remaining_tasks = deepcopy(batch_controller.task_id)
        task_id = nothing

        __verbose__ && print("Monitoring tasks for 'Completed' state, timeout in $timeout minutes ...")
        while true

            # Wait for one task from task list to finish
            try
                task_id = wait_for_one_task_from_multi_pool(batch_controller.batch_client, batch_controller.job_id, remaining_tasks;
                    task_timeout=task_timeout, fetch_timeout=timeout, verbose=__verbose__, num_restart=num_restart)[1]
            catch
                throw("Reached timeout for task completion.")
            end
            __verbose__ && print("\nFetch output from task $task_id")

            # Fetch its output
            task_no = findall(i -> i["taskname"] == task_id, batch_controller.task_id)[1]
            if ~isempty(remaining_tasks)

                # Fetch output
                temp = fetch!(batch_controller, task_no; destroy_blob=destroy_blob, timeout=timeout,
                    wait_for_completion=false)
                for (i, entry) in enumerate(output)
                    if ~isnothing(entry) && ~isnothing(temp)
                        entry[:] = broadcast(op, entry, temp[i])
                    end
                end

                # Remove task from task list
                local_id = findall(i -> i["taskname"] == task_id, remaining_tasks)[1]
                popat!(remaining_tasks, local_id)
            end

            # Return it task list is empty
            if isempty(remaining_tasks)
                __verbose__ && print("\n")
                if length(output) > 1
                    return output
                else
                    return output[1]
                end
            end
        end
    end
end

# Inplace fetch-reduce for single output argument
function fetchreduce!(batch_controller::BatchController, output; op=+, destroy_blob=true, timeout=60, task_timeout=60,
    num_restart=0)

    # Single task: simply fetch output
    if batch_controller.num_tasks == 1
        if ~isnothing(output)
            output[:] = fetch!(batch_controller)
        end

    # Multiple tasks: fetch and sum output
    else

        remaining_tasks = deepcopy(batch_controller.task_id)
        task_id = nothing

        __verbose__ && print("Monitoring tasks for 'Completed' state, timeout in $timeout minutes ...")
        while true

            # Wait for one task from task list to finish
            try
                task_id = wait_for_one_task_from_multi_pool(batch_controller.batch_client, batch_controller.job_id, remaining_tasks;
                    task_timeout=task_timeout, fetch_timeout=timeout, verbose=__verbose__, num_restart=num_restart)[1]
            catch
                throw("Reached timeout for task completion.")
            end
            __verbose__ && print("\nFetch output from task $task_id")

            # Fetch its output
            task_no = findall(i -> i["taskname"] == task_id, batch_controller.task_id)[1]
            if ~isempty(remaining_tasks)

                # Fetch output
                temp = fetch!(batch_controller, task_no; destroy_blob=destroy_blob, timeout=timeout,
                    wait_for_completion=false)
                if ~isnothing(output) && ~isnothing(temp)
                    output[:] = broadcast(op, output, temp)
                end

                # Remove task from task list
                local_id = findall(i -> i["taskname"] == task_id, remaining_tasks)[1]
                popat!(remaining_tasks, local_id)
            end

            # Return it task list is empty
            if isempty(remaining_tasks)
                __verbose__ && print("\n")
                if length(output) > 1
                    return output
                else
                    return output[1]
                end
            end
        end
    end
end


text:
# This file contains "tests for humans:" functions you might use to reassure
# yourself that the heuristics are working properly or chose a heuristic for your
# particular problem. The functions in this file are not run by Travis and this file
# requires several packages (such as Gadfly) not listed in the REQUIRE.
using Gadfly
using TravelingSalesmanHeuristics
using Distances

# use Gadfly to plot a TSP visualization
# pts must be a 2 x n matrix and path the vector of column indices
function visTSP(pts, path)
	pts = pts[:,path]
	plot(x = pts[1,:], y = pts[2,:], Geom.path, Geom.point)
end
# visualize a planar TSP in the unit square
# specify number of points and solution method
function visPlanarTSP(n = 10, heuristic = nearest_neighbor; extraArgs...)
	pts = rand(2, n)
	distMat = pairwise(Euclidean(), pts, pts)
	path, _ = heuristic(distMat; extraArgs...)
	visTSP(pts, path)
end

# compare some heuristics
function compareHeuristics(n = 20)
	pts = rand(2, n)
	dm = pairwise(Euclidean(), pts, pts)
	println("Simple lower bound: $(lowerbound(dm))")
	singlenn_no2opt = nearest_neighbor(dm, do2opt = false)
	println("Single start nearest neighbor without 2-opt has cost $(singlenn_no2opt[2])")
	singlenn_2opt = nearest_neighbor(dm, do2opt = true)
	println("Single start nearest neighbor with 2-opt has cost $(singlenn_2opt[2])")
	repnn_no2opt = repetitive_heuristic(dm, nearest_neighbor, do2opt = false)
	println("Multi start nearest neighbor without 2-opt has cost $(repnn_no2opt[2])")
	repnn_2opt = repetitive_heuristic(dm, nearest_neighbor, do2opt = true)
	println("Multi start nearest neighbor with 2-opt has cost $(repnn_2opt[2])")
	singleci_no2opt = cheapest_insertion(dm, do2opt = false)
	println("Single start cheapest insert without 2-opt has cost $(singleci_no2opt[2])")
	singleci_2opt = cheapest_insertion(dm, do2opt = true)
	println("Single start cheapest insert with 2-opt has cost $(singleci_2opt[2])")
	multici_no2opt = repetitive_heuristic(dm, cheapest_insertion, do2opt = false)
	println("Multi start cheapest insert without 2-opt has cost $(multici_no2opt[2])")
	multici_2opt = repetitive_heuristic(dm, cheapest_insertion, do2opt = true)
	println("Multi start cheapest insert with 2-opt has cost $(multici_2opt[2])")
	simanneal = simulated_annealing(dm)
	println("Simulated annealing with default arguments has cost $(simanneal[2])")
end


text:
function mass_matrix_barycentric(v, f)
    f1 = [face[1] for face in f]
    f2 = [face[2] for face in f]
    f3 = [face[3] for face in f]
    # F x 3 matrices of triangle edge vectors, named after opposite vertices
    v1 = [v[f3[i]] - v[f2[i]] for i in eachindex(f)];
    v2 = [v[f1[i]] - v[f3[i]] for i in eachindex(f)];
    v3 = [v[f2[i]] - v[f1[i]] for i in eachindex(f)];
    n = cross.(v1, v2)
    area = norm.(n) / 2
    i = [f1;f2;f3]
    j = [f1;f2;f3]

    diag_v = area / 3
    val = [diag_v;diag_v;diag_v];
    M = sparse(i, j, val, size(v, 1), size(v, 1));
    return M
end

function mass_matrix_barycentric(v::Matrix{<: AbstractFloat}, f::Matrix{<: Signed})
    f1 = @view f[:, 1]
    f2 = @view f[:, 2]
    f3 = @view f[:, 3]

    l1 = sqrt.(sum((view(v, f2, :) - view(v, f3, :)).^2, dims=2))
    l2 = sqrt.(sum((view(v, f3, :) - view(v, f1, :)).^2, dims=2))
    l3 = sqrt.(sum((view(v, f1, :) - view(v, f2, :)).^2, dims=2))
    s = (l1 + l2 + l3) / 2.0
    # @show s
    squared_area = s .* (s - l1) .* (s - l2) .* (s - l3)
    area = sqrt.(squared_area)
    i = [f1;f2;f3]
    j = [f1;f2;f3]

    diag_v = reshape(area / 3, :)
    val = [diag_v;diag_v;diag_v];
    M = sparse(i, j, val, size(v, 1), size(v, 1));
    return M
end

function mass_matrix_barycentric2(v::Matrix{<: AbstractFloat}, f::Matrix{<: Signed})
    f1 = @view f[:,1]
    f2 = @view f[:,2]
    f3 = @view f[:,3]

    v1 = view(v, f3, :) - view(v, f2, :);
    v2 = view(v, f1, :) - view(v, f3, :);
    v3 = view(v, f2, :) - view(v, f1, :);
    n = [cross(v1[i,:], v2[i,:]) for i in 1:size(v1, 1)]
    area = norm.(n) / 2
    i = [f1;f2;f3]
    j = [f1;f2;f3]

    diag_v = area / 3
    val = [diag_v;diag_v;diag_v];
    M = sparse(i, j, val, size(v, 1), size(v, 1));
    return M
end
"""
This function return the cotangent weighted laplacian matrix.
"""
function cotangent_laplacian(v, f)
    nv = length(v)
    nf = length(f)
    f1 = [face[1] for face in f]
    f2 = [face[2] for face in f]
    f3 = [face[3] for face in f]

    cot12 = Vector{Float32}(undef, nf)
    cot23 = Vector{Float32}(undef, nf)
    cot31 = Vector{Float32}(undef, nf)
    diag1 = Vector{Float32}(undef, nf)
    diag2 = Vector{Float32}(undef, nf)
    diag3 = Vector{Float32}(undef, nf)
    for i in 1:nf
        l1 = norm(v[f[i][2]] - v[f[i][3]])
        l2 = norm(v[f[i][3]] - v[f[i][1]])
        l3 = norm(v[f[i][1]] - v[f[i][2]])
        s = (l1 + l2 + l3) / 2.0
        square_area = s * (s - l1) * (s - l2) * (s - l3)
        if square_area <= 0.0; square_area = 1e-12 end
        area = sqrt(square_area)
        cot12_ = (l1^2 + l2^2 - l3^2) / area / 2
        cot23_ = (l2^2 + l3^2 - l1^2) / area / 2
        cot31_ = (l1^2 + l3^2 - l2^2) / area / 2
        cot12[i] = cot12_
        cot23[i] = cot23_
        cot31[i] = cot31_
        diag1[i] = -cot12_ - cot31_
        diag2[i] = -cot12_ - cot23_
        diag3[i] = -cot31_ - cot23_
    end
    II = [f1; f2; f2; f3; f3; f1; f1; f2; f3]
    JJ = [f2; f1; f3; f2; f1; f3; f1; f2; f3]
    V = [cot12; cot12; cot23; cot23; cot31; cot31; diag1; diag2; diag3]
    L = sparse(II, JJ, V, nv, nv)
    return L
end

function cotangent_laplacian(v::Matrix{<:AbstractFloat}, f::Matrix{<:Signed})

    nv = size(v, 1)
    f1 = @view f[:,1]
    f2 = @view f[:,2]
    f3 = @view f[:,3]

    l1 = sqrt.(sum((v[f2,:] - v[f3,:]).^2, dims=2))
    l2 = sqrt.(sum((v[f3,:] - v[f1,:]).^2, dims=2))
    l3 = sqrt.(sum((v[f1,:] - v[f2,:]).^2, dims=2))
    s = (l1 + l2 + l3) / 2.0

    area = sqrt.(s .* (s - l1) .* (s - l2) .* (s - l3))
    cot12 = (l1.^2 + l2.^2 - l3.^2) ./ area / 2
    cot23 = (l2.^2 + l3.^2 - l1.^2) ./ area / 2
    cot31 = (l1.^2 + l3.^2 - l2.^2) ./ area / 2
    diag1 = -cot12 - cot31
    diag2 = -cot12 - cot23
    diag3 = -cot31 - cot23

    II = [f1; f2; f2; f3; f3; f1; f1; f2; f3]
    JJ = [f2; f1; f3; f2; f1; f3; f1; f2; f3]
    V = [cot12; cot12; cot23; cot23; cot31; cot31; diag1; diag2; diag3]
    L = sparse(II, JJ, V, nv, nv)
    return L
end

function vertex_adjacent_area(verts, faces)
    areas = zeros(Float64, length(verts))
    for j in 1:length(faces)
        f = faces[j]
        area = calculate_area(verts[f[1]], verts[f[2]], verts[f[3]])
        for i in 1:3
            areas[f[i]] += area / 3.0
        end
    end
    return areas
end
"""
This function return the gaussian curvature of the vertex
"""
function gaussian_curvature(verts, faces)

    int_angles = internal_angles(verts, faces)
    gauss_K = repeat([2 * π], length(verts))
    for j in 1:length(faces)
        for i in 1:3
            vert = faces[j][i]
            gauss_K[vert] -= int_angles[j,i]
        end
    end
    # return gauss_K ./ vertex_adjacent_area(verts,faces)
    return gauss_K
end
function gaussian_curvature(mesh::Mesh)
    verts = decompose(Point3f0, mesh)
    faces = decompose(TriangleFace{Int}, mesh)
    return gaussian_curvature(verts, faces)
end
"""
This function return the mean curvature of the vertex
"""
function mean_curvature(verts, faces)
    l = cotangent_laplacian(verts, faces)
    m = mass_matrix_barycentric(verts, faces)
    v = MovingFrame.array_of_array_to_matrix(verts)
    minv = spdiagm(0 => (1.0 ./ diag(m)))

    hn = -minv * (l * v )
    h = sqrt.(sum(hn.^2, dims=2))
    # return gauss_K ./ vertex_adjacent_area(verts,faces)
    return reshape(h, :)
end
function mean_curvature(mesh::Mesh)
    verts = decompose(Point3f0, mesh)
    faces = decompose(TriangleFace{Int}, mesh)
    return mean_curvature(verts, faces)
end
function mean_curvature_flow(verts, faces; n_step=10, delta=0.001,min_diff=1e-13)
    new_verts = verts
    mean_coord = reduce(+, verts) / length(verts)
    if ~isapprox(mean_coord, [0f0,0f0,0f0], atol=1e-7)
        new_verts, _ = mesh_centering(new_verts, faces)
    end

    if ~isapprox(mesh_volume(new_verts, faces), 0.0, atol=1e-7)
        new_verts, _ = volume_normalizing(new_verts, faces)
    end

    new_v = MovingFrame.array_of_array_to_matrix(new_verts)
    f = MovingFrame.array_of_array_to_matrix(faces)
    new_v = new_v / sqrt(surface_area2(new_v, f))
    l = cotangent_laplacian(new_verts, faces)

    for i in 1:n_step
        old_v = copy(new_v)
        m = mass_matrix_barycentric(new_v, f)
        s = (m - delta * l)
        b = m * new_v
        new_v = s \ b
        new_v = new_v .- centroid(new_v)
        area = surface_area_list(new_v, f)
        new_v = new_v / sqrt(sum(area))
        diff = tr((new_v - old_v)' * mass_matrix_barycentric(new_v, f) * (new_v - old_v))
        if diff <= min_diff
            # println("Converged")
            break
        end
    end
    new_verts::typeof(verts) = MovingFrame.matrix_to_vertex(new_v)
    return new_verts
end
function spherization(verts, faces;n_step=10, delta=0.001,min_diff=1e-13)
    new_verts =  mean_curvature_flow(verts, faces; n_step=n_step, delta=delta,min_diff=min_diff)
    return new_verts, faces
end
function spherization(mesh::Mesh;n_step=10, delta=0.001,min_diff=1e-13)
    verts = decompose(Point3f0, mesh)
    triangles = decompose(TriangleFace{Int}, mesh)
    new_verts =  mean_curvature_flow(verts, triangles; n_step=n_step, delta=delta,min_diff=min_diff)
    return Mesh(new_verts, faces(mesh))
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export chips

"""
    chips(;kwargs...)

A Chips component.
Alternative to Select and RadioGroup. For more information, see: https://mantine.dev/core/chips/
Keyword arguments:
- `id` (String; optional): The ID of this component, used to identify dash components in callbacks
- `align` (a value equal to: "stretch", "center", "flex-end", "flex-start"; optional): Defines align-items css property
- `className` (String; optional): Often used with CSS to style elements with common properties
- `color` (a value equal to: "dark", "gray", "red", "pink", "grape", "violet", "indigo", "blue", "cyan", "teal", "green", "lime", "yellow", "orange"; optional): Active chip color, defaults to theme.primaryColor
- `data` (optional): Chips. data has the following type: Array of lists containing elements 'label', 'value', 'disabled'.
Those elements have the following types:
  - `label` (String; required): The option's label
  - `value` (String; required): Option's value
  - `disabled` (Bool; optional): If true, this option is disabled and cannot be selecteds
- `direction` (a value equal to: "row", "column"; optional): Defines flex-direction property, row for horizontal, column for vertical
- `grow` (Bool; optional): Defines flex-grow property for each element, true -> 1, false -> 0
- `multiple` (Bool; optional): Allow multiple values to be picked
- `noWrap` (Bool; optional): Defined flex-wrap property
- `position` (a value equal to: "right", "center", "left", "apart"; optional): Defines justify-content property
- `radius` (a value equal to: "xs", "sm", "md", "lg", "xl" | Real; optional): Chip border-radius from theme or number to set value in px
- `size` (a value equal to: "xs", "sm", "md", "lg", "xl"; optional): Predefined chip size
- `spacing` (a value equal to: "xs", "sm", "md", "lg", "xl" | Real; optional): Spacing between chips from theme or number to set value in px
- `style` (Dict; optional): Inline style override
- `value` (String | Array of Strings; optional): Controlled component value
- `variant` (a value equal to: "outline", "filled"; optional): Controls chip appearance, defaults to filled with dark theme and to outline in light theme
"""
function chips(; kwargs...)
        available_props = Symbol[:id, :align, :className, :color, :data, :direction, :grow, :multiple, :noWrap, :position, :radius, :size, :spacing, :style, :value, :variant]
        wild_props = Symbol[]
        return Component("chips", "Chips", "dash_mantine_components", available_props, wild_props; kwargs...)
end


text:
filename = length(ARGS) > 0 ? ARGS[1] : "page1.txt"

cf = Dict()

f = open(filename)
for line in eachline(f)
    for c in line
        cf[c] = get(cf, c, 0) + 1
    end
end
close(f)

cf = sort(collect(cf), by=x->x[2], rev=true)

for (k, v) in cf
    println("$k: $v")
end


text:
using IncrementalInference
# using KernelDensityEstimate, Gadfly # for vstack

fg = emptyFactorGraph()

N=100

doors = [-100.0;0.0;100.0;300.0]'
cov = [3.0]

v1 = addNode!(fg,"x1",doors,N=N)
f1  = addFactor!(fg,[v1], Obsv2(doors, cov', [1.0]))

tem = 2.0*randn(1,N)+getVal(v1)+50.0
v2 = addNode!(fg,"x2", tem, N=N)
addFactor!(fg,[v1;v2],Odo([50.0]',[2.0]',[1.0]))

v3=addNode!(fg,"x3",4.0*randn(1,N)+getVal(v2)+50.0, N=N)
addFactor!(fg,[v2;v3],Odo([50.0]',[4.0]',[1.0]))

v4=addNode!(fg,"x4",2.0*randn(1,N)+getVal(v3)+50.0, N=N)
addFactor!(fg,[v3;v4],Odo([50.0]',[2.0]',[1.0]))

v5=addNode!(fg,"x5",2.0*randn(1,N)+getVal(v4)+50.0, N=N)
addFactor!(fg,[v4;v5],Odo([50.0]',[2.0]',[1.0]))

l1=addNode!(fg, "l1", 0.5*randn(1,N)+getVal(v3)+64.0, N=N)
addFactor!(fg, [v3,l1], Ranged([64.0],[0.5],[1.0]))

l2=addNode!(fg, "l2", 0.5*randn(1,N)+getVal(v3)+64.0, N=N)
addFactor!(fg, [v4,l2], Ranged([16.0],[0.5],[1.0]))

l3=addNode!(fg, "l3", 0.5*randn(1,N)+getVal(v3)+64.0, N=N)
addFactor!(fg, [v5,l3], Ranged([16.0],[0.5],[1.0]))

# writeGraphPdf(fg);

tree = prepBatchTree!(fg,drawpdf=true);

# do belief propagation inference over tree once
# inferOverTreeR!(fg, tree, N=100)


text:
using Test, Shuffle, Documenter, OffsetArrays
using Random: MersenneTwister, default_rng
import Random

@testset "Shuffle.jl" begin

include("shuffles.jl")

include("defaults.jl")

if false
    DocMeta.setdocmeta!(
        Shuffle,
        :DocTestSetup,
        quote
            using Shuffle
            import Random
            using Random: MersenneTwister
        end;
        recursive=true
    )

    doctest(Shuffle)
end

end