text:
"""
    CircularTurnBuffer{names, types, Tbs} <: AbstractTurnBuffer{names, types}
    CircularTurnBuffer{names, types}(capacities::NTuple{N, Int}, sizes::NTuple{N, NTuple{M, Int} where M}) where {names, types, N}

Using `CircularArrayBuffer` to store each element specified in `names` and `types`.. The memory of the buffer will be pre-allocated.
In RL problems, three of the most common `CircularArrayBuffer` based buffers are: [`CircularSARDBuffer`](@ref), [`CircularSARDSBuffer`](@ref), [`CircularSARDSABuffer`](@ref)
"""
struct CircularTurnBuffer{names, types, Tbs, Tm} <: AbstractTurnBuffer{names, types}
    buffers::Tbs
    meta::Tm
    function CircularTurnBuffer{names, types}(
        capacities::NTuple{N, Int},
        sizes::NTuple{N, NTuple{M, Int} where M};
        args...) where {names, types, N}
        buffers = merge(NamedTuple(),
                        (names[i], CircularArrayBuffer{types.parameters[i]}(capacities[i], sizes[i]...)) for i in 1:N)
        meta = args.data
        new{names, types, typeof(buffers), typeof(meta)}(buffers, meta)
    end
end

length(b::CircularTurnBuffer) = length(b.isdone)
capacity(b::CircularTurnBuffer) = capacity(b.isdone)
isfull(b::CircularTurnBuffer) = isfull(b.isdone)

##############################
# CircularSARDBuffer
##############################

"`CircularSARDBuffer` is just an alias for `CircularTurnBuffer{(:state, :action, :reward, :isdone)}`."
const CircularSARDBuffer = CircularTurnBuffer{SARD}

"""
    CircularSARDBuffer(capacity; state_type::Type=Int, action_type::Type=Int, state_size=(), action_size=())

`capacity` specifies how many latest turns the buffer will store at most.

!!! note
    Note that, the length of state and action is 1 step longer in oder to store the state and action in the next step. This is the supposed behavior of **SARD** buffers.

See also: [`sample`](@ref)
"""
function CircularSARDBuffer(
    capacity;
    state_type::Type=Int,
    action_type::Type=Int,
    state_size=(),
    action_size=())
    CircularSARDBuffer{Tuple{state_type, action_type, Float64, Bool}}(
        (capacity+1, capacity+1, capacity, capacity),
        (state_size, action_size, (), ()))
end

getindex(b::CircularSARDBuffer, i::Int) = (state=b.state[i], action=b.action[i], reward=b.reward[i], isdone=b.isdone[i], nextstate=b.state[i+1], nextaction=b.action[i+1])
function eltype(b::CircularSARDBuffer{names, types}) where {names, types}
    ts, ta, tr, td = types.parameters[2].parameters
    NamedTuple{SARDSA, Tuple{ts, ta, tr, td, ts, ta}}
end

function push!(b::CircularSARDBuffer, s, a)
    if isempty(b)
        push!(b.state, s)
        push!(b.action, a)
    else
        b.state[end] = s
        b.action[end] = a
    end
end

function push!(b::CircularSARDBuffer, r, d, ns, na)
    push!(b.reward, r)
    push!(b.isdone, d)
    push!(b.state, ns)
    push!(b.action, na)
end

"""
    sample(b::CircularSARDBuffer, batch_size::Int)

Sample a random batch of **S**tates, **A**ctions, **R**ewards, is**D**one,
next**S**tates, next**A**ctions without replacement of `batch_size`.
"""
function sample(b::CircularSARDBuffer, batch_size::Int)
    inds = rand(1:length(b), batch_size)
    batch_view(b, inds), inds
end

function batch_view(b::CircularSARDBuffer, inds)
    view(b.state, inds),
    view(b.action, inds),
    view(b.reward, inds),
    view(b.isdone, inds),
    view(b.state, inds .+ 1)
end

##############################
# CircularSARDSBuffer
##############################
"`CircularSARDSBuffer` is just an alias for `CircularTurnBuffer{(:state, :action, :reward, :isdone, :nextstate)}`."
const CircularSARDSBuffer = CircularTurnBuffer{SARDS}

function CircularSARDSBuffer(
    capacity;
    state_type::Type=Int,
    action_type::Type=Int,
    state_size=(),
    action_size=())
    CircularSARDSBuffer{Tuple{state_type, action_type, Float64, Bool, state_type}}(
        Tuple(capacity for i in 1:5),
        (state_size, action_size, (), (), state_size))
end

function push!(b::CircularSARDSBuffer{Tuple{Ts, Ta, Float64, Bool, Ts}}, s::Ts, a::Ta, r::Float64, d::Bool, ns::Ts) where {Ts, Ta}
    push!(b.state, s)
    push!(b.action, a)
    push!(b.reward, r)
    push!(b.isdone, d)
    push!(b.nextstate, ns)
end

##############################
# CircularSARDSABuffer
##############################

"`CircularSARDSABuffer` is just an alias for `CircularTurnBuffer{(:state, :action, :reward, :isdone, :nextstate, :nextaction)}`."
const CircularSARDSABuffer = CircularTurnBuffer{SARDSA}

function CircularSARDSABuffer(
    capacity;
    state_type::Type=Int,
    action_type::Type=Int,
    state_size=(),
    action_size=())
    CircularSARDSABuffer{Tuple{state_type, action_type, Float64, Bool, state_type, action_type}}(
        Tuple(capacity for i in 1:6),
        (state_size, action_size, (), (), state_size, action_size))
end

function push!(b::CircularSARDSABuffer{Tuple{Ts, Ta, Float64, Bool, Ts, Ta}}, s::Ts, a::Ta, r::Float64, d::Bool, ns::Ts, na::Ta) where {Ts, Ta}
    push!(b.state, s)
    push!(b.action, a)
    push!(b.reward, r)
    push!(b.isdone, d)
    push!(b.nextstate, ns)
    push!(b.nextaction, na)
end


text:
using DelaunayTriangulations

@testset "Delaunay" begin
    vertices = [(0.0, 0.0), (0.5, 0.0), (1.0, 0.0),
                (0.0, 0.5), (0.5, 0.5), (1.0, 0.5),
                (0.0, 1.0), (0.5, 1.0), (1.0, 1.0)]
    triangulation = Triangulation(vertices)
    add_triangle!(triangulation, (1, 2, 5))
    add_triangle!(triangulation, (1, 5, 4))
    add_triangle!(triangulation, (2, 3, 6))
    add_triangle!(triangulation, (2, 6, 5))
    add_triangle!(triangulation, (4, 5, 7))
    add_triangle!(triangulation, (5, 8, 7))
    add_triangle!(triangulation, (5, 6, 8))
    add_triangle!(triangulation, (8, 6, 9))
    # Ghost triangles
    add_triangle!(triangulation, (2, 1, -1))
    add_triangle!(triangulation, (1, 4, -1))
    add_triangle!(triangulation, (4, 7, -1))
    add_triangle!(triangulation, (7, 8, -1))
    add_triangle!(triangulation, (8, 9, -1))
    add_triangle!(triangulation, (9, 6, -1))
    add_triangle!(triangulation, (6, 3, -1))
    add_triangle!(triangulation, (3, 2, -1))
    @testset "Insert Vertex" begin
        @testset "Inside Triangulation" begin
            # New vertex has index 10
            insert_vertex!(triangulation, (0.25, 0.01), (1, 2, 5))
            # Triangle (1, 2, 5) and (1, 5, 4) are split into
            # (1, 10, 4), (1, 2, 10), (10, 2, 5) and (10, 5, 4)
            @test adjacent(triangulation, (1, 2)) == 10
            @test adjacent(triangulation, (2, 10)) == 1
            @test adjacent(triangulation, (10, 1)) == 2
            @test adjacent(triangulation, (10, 2)) == 5
            @test adjacent(triangulation, (2, 5)) == 10
            @test adjacent(triangulation, (5, 10)) == 2
            @test adjacent(triangulation, (10, 5)) == 4
            @test adjacent(triangulation, (5, 4)) == 10
            @test adjacent(triangulation, (4, 10)) == 5
            @test adjacent(triangulation, (10, 4)) == 1
            @test adjacent(triangulation, (4, 1)) == 10
            @test adjacent(triangulation, (1, 10)) == 4
            # All other triangles are unchanged
            @test adjacent(triangulation, (2, 3)) == 6
            @test adjacent(triangulation, (3, 6)) == 2
            @test adjacent(triangulation, (6, 2)) == 3
            @test adjacent(triangulation, (6, 5)) == 2
            @test adjacent(triangulation, (5, 2)) == 6
            @test adjacent(triangulation, (2, 6)) == 5
            @test adjacent(triangulation, (4, 5)) == 7
            @test adjacent(triangulation, (5, 7)) == 4
            @test adjacent(triangulation, (7, 4)) == 5
            @test adjacent(triangulation, (5, 6)) == 8
            @test adjacent(triangulation, (6, 8)) == 5
            @test adjacent(triangulation, (8, 5)) == 6
            @test adjacent(triangulation, (6, 9)) == 8
            @test adjacent(triangulation, (9, 8)) == 6
            @test adjacent(triangulation, (8, 6)) == 9
            @test adjacent(triangulation, (2, 1)) == -1
            @test adjacent(triangulation, (1, 4)) == -1
            @test adjacent(triangulation, (4, 7)) == -1
            @test adjacent(triangulation, (7, 8)) == -1
            @test adjacent(triangulation, (8, 9)) == -1
            @test adjacent(triangulation, (9, 6)) == -1
            @test adjacent(triangulation, (6, 3)) == -1
            @test adjacent(triangulation, (3, 2)) == -1
        end
        @testset "Outside Triangulation" begin
            # New vertex has index 11
            insert_vertex!(triangulation, (1.1, 0.6), (9, 6, -1))
            # Ghost triangles (9, 6, -1) and (6, 3, -1) are split into
            # (6, 11, 9), (6, 3, 11), (9, 11, -1) and (11, 3, -1)
            @test adjacent(triangulation, (9, 6)) == 11
            @test adjacent(triangulation, (6, 11)) == 9
            @test adjacent(triangulation, (11, 9)) == 6
            @test adjacent(triangulation, (6, 3)) == 11
            @test adjacent(triangulation, (3, 11)) == 6
            @test adjacent(triangulation, (11, 6)) == 3
            @test adjacent(triangulation, (9, 11)) == -1
            @test adjacent(triangulation, (11, -1)) == 9
            @test adjacent(triangulation, (-1, 9)) == 11
        end
        @testset "Triangulation" begin
            vertices = [(rand(), rand()) for _ in 1:100]
            triangulation = delaunay_triangulation(vertices)
            for (u, v) in keys(triangulation.edges)
                w = adjacent(triangulation, (u, v))
                if u != -1 && v != -1 && w != -1
                    for q in 1:length(triangulation.vertices)
                        if q != u && q != v && q != w
                            @test incircle(triangulation, u, v, w, q) < 0.0
                        end
                    end
                end
            end
        end
    end
end

ext: jl

max_stars_repo_name: PaulMatlashewski/DelaunayTriangulations

max_issues_repo_name: PaulMatlashewski/DelaunayTriangulations

max_forks_repo_name: PaulMatlashewski/DelaunayTriangulations


text:
const NAME = "TestDir"
const test_model = Model(NAME);
const n = 10
const λs = sort(n * rand(n))
const ∇²f = Symmetric(begin
    Q = orth(randn(n, n))
    Q' * diagm(λs) * Q
end)
const ∇²f⁻¹ = inv(∇²f)
const f = x -> 0.5 * x'∇²f * x
const LaTex = "\$\\frac{1}{2} x^\\Top \\nabla^2 f x\$"
const ∇f! = (g, x) -> (g .= ∇²f * x)
const x₀ = randn(n)
const f₀ = f(x₀)
const ∇f₀ = ∇f!(similar(x₀), x₀)
const ∇f₀_norm = norm(∇f₀)
const h₀ = ∇²f * ∇f₀
const test_options = DriverOptions();
const test_driver = Driver(options = test_options);
const Δₘ = 100.0
const δ = 1.0e-12
const ϵ = 10.0e-6
const w = 3
const N = 2000
const WEAVE_LEVEL = ALL
const LOG_LEVEL = INFO
const S_UPDATE = S_update_c
const QN_UPDATE = SR1
const PFLAG = false

function restricted_type_test(type)
    @test_throws AccessError type.afield

    @test propertynames(type) ≡ ()

    return nothing
end

# TODO: Define memory manipulation test here

function rosenbrock_model(n::Int)
    rosen_model = Model("rosenbrock")

    objective!(
        rosen_model,
        x -> begin
            N = lastindex(x)
            return 100sum((x[i+1] - x[i]^2)^2 for i = 1:N-1) +
                   sum((x[i] - 1)^2 for i = 1:N-1)
        end,
    )

    gradient!(
        rosen_model,
        (g, x) -> begin
            N = lastindex(x)
            g[1] = -2 * (1 - x[1]) - 400 * x[1] * (-x[1]^2 + x[2])
            for i = 2:N-1
                g[i] =
                    -2 * (1 - x[i]) + 200 * (-x[i-1]^2 + x[i]) -
                    400 * x[i] * (-x[i]^2 + x[1+i])
            end
            g[N] = 200 * (x[N] - x[N-1]^2)
            return g
        end,
    )

    initial_iterate!(rosen_model, randn(100))

    formula!(
        rosen_model,
        "\$\\sum_{i=1}^{N-1} \\left[100(x_{i+1}^2 - x_i^2)^2 + (1 - x_i)^2\\right]\$",
    )

    return rosen_model
end

# function rosen(x)
#    N = lastindex(x)
#    return 100sum((x[i+1] - x[i]^2)^2 for i = 1:N-1) + sum((x[i] - 1)^2 for i = 1:N-1)
# end

# function ∇rosen!(g, x)
#     N = lastindex(x)
#     g[1] = -2*(1 - x[1]) - 400*x[1]*(-x[1]^2 + x[2])
#     for i in 2:N-1
#         g[i] = -2*(1 - x[i]) + 200*(-x[i-1]^2 + x[i]) - 400*x[i]*(-x[i]^2 + x[1 + i])
#     end
#     g[N] = 200 * (x[N] - x[N-1]^2)
#     return g
# end


text:
#
println("Fc) Symbolic simplification of recoupling coefficients by means of sum rules.")
#
using SymEngine
j1  = Basic(:j1);    j2  = Basic(:j2);    j3  = Basic(:j3);    j4  = Basic(:j4);    j5  = Basic(:j5)
j6  = Basic(:j6);    j7  = Basic(:j7);    j8  = Basic(:j8);    j9  = Basic(:j9);    j10 = Basic(:j10)
j11 = Basic(:j11);   j12 = Basic(:j12);   j13 = Basic(:j13);   j14 = Basic(:j14);   j15 = Basic(:j15)
j16 = Basic(:j16);   j17 = Basic(:j17);   j18 = Basic(:j18);   j19 = Basic(:j19);   j20 = Basic(:j20)
j21 = Basic(:j21);   j22 = Basic(:j22);   j23 = Basic(:j23);   j24 = Basic(:j24);   j25 = Basic(:j25)

L1  = Basic(:L1);    L2  = Basic(:L2);    L3  = Basic(:L3);    L12 = Basic(:L12);    L123 = Basic(:L123);
S1  = Basic(:S1);    S2  = Basic(:S2);    S3  = Basic(:S3);    S12 = Basic(:S12);    S123 = Basic(:S123);
T1  = Basic(:T1);    T2  = Basic(:T2);    T3  = Basic(:T3);    T12 = Basic(:T12)
Jm1 = Basic(:Jm1);   Jp1 = Basic(:Jp1);   J1  = Basic(:J1);    Jm2 = Basic(:Jm2);   Jp2 = Basic(:Jp2);   J2  = Basic(:J2)
Jm3 = Basic(:Jm3);   Jp3 = Basic(:Jp3);   J3  = Basic(:J3);    J12 = Basic(:J12)
Xp1 = Basic(:Xp1);   Xm2 = Basic(:Xm2);   Xp2 = Basic(:Xp2);   Xm3 = Basic(:Xm3);   Xp3 = Basic(:Xp3);
J   = Basic(:J)

if  false
    #
    leftCsq  = RacahAlgebra.Csq( j1, j2, j3)
    wa       = RacahAlgebra.evaluate(leftCsq, leftCsq)
elseif false
    #
    leftCsq  = RacahAlgebra.Csq( RacahAlgebra.Csq( j1, j2, j12), j3, J)
    rightCsq = RacahAlgebra.Csq( j1, RacahAlgebra.Csq( j2, j3, j23), J)
    wa       = RacahAlgebra.evaluate(leftCsq, rightCsq)
elseif false
    # LSJ-recoupling, two open shells
    leftCsq  = RacahAlgebra.Csq( RacahAlgebra.Csq( L1, L2, L12), RacahAlgebra.Csq( S1, S2, S12), J)
    rightCsq = RacahAlgebra.Csq( RacahAlgebra.Csq( L1, S1, T1),  RacahAlgebra.Csq( L2, S2, T2),  J)
    @time  wa       = RacahAlgebra.evaluate(leftCsq, rightCsq)
elseif false
    # jjJ-recoupling, two open shells
    leftCsq  = RacahAlgebra.Csq( RacahAlgebra.Csq( Jm1, Jp1, J1),   RacahAlgebra.Csq( Jm2, Jp2, J2), J)
    rightCsq = RacahAlgebra.Csq( RacahAlgebra.Csq( RacahAlgebra.Csq( Jm1, Jp1, Xp1), Jm2, Xm2), Jp2, J)
    wa       = RacahAlgebra.evaluate(leftCsq, rightCsq)
elseif false
    # LSJ-recoupling, three open shells
    leftCsq  = RacahAlgebra.Csq( RacahAlgebra.Csq( RacahAlgebra.Csq( L1, L2, L12), L3, L123),
                                 RacahAlgebra.Csq( RacahAlgebra.Csq( S1, S2, S12), S3, S123), J)
    rightCsq = RacahAlgebra.Csq( RacahAlgebra.Csq( RacahAlgebra.Csq( L1, S1, T1),  RacahAlgebra.Csq( L2, S2, T2), T12),
                                                   RacahAlgebra.Csq( L3, S3, T3),   J)
    wa       = RacahAlgebra.evaluate(leftCsq, rightCsq)
elseif true
    # jjJ-recoupling, three open shells
    leftCsq  = RacahAlgebra.Csq( RacahAlgebra.Csq( RacahAlgebra.Csq( Jm1, Jp1, J1),   RacahAlgebra.Csq( Jm2, Jp2, J2), J12),
                                                   RacahAlgebra.Csq( Jm3, Jp3, J3), J)
    rightCsq = RacahAlgebra.Csq( RacahAlgebra.Csq( RacahAlgebra.Csq( RacahAlgebra.Csq(
                                 RacahAlgebra.Csq( Jm1, Jp1, Xp1), Jm2, Xm2), Jp2, Xp2), Jm3, Xm3), Jp3, J)
    wa       = RacahAlgebra.evaluate(leftCsq, rightCsq)
elseif true
    #
    leftCsq  = RacahAlgebra.Csq( RacahAlgebra.Csq(j1,j2,j5), RacahAlgebra.Csq(j3,j4,j6), j7 )
    rightCsq = RacahAlgebra.Csq( j1, RacahAlgebra.Csq( RacahAlgebra.Csq(j2,j3,j8), j4, j9), j7 )
    wa       = RacahAlgebra.evaluate(leftCsq, rightCsq)
elseif true
    #
    leftCsq  = RacahAlgebra.Csq( RacahAlgebra.Csq( RacahAlgebra.Csq(j1,j2,j9),j3,j10),
                                 RacahAlgebra.Csq( RacahAlgebra.Csq( RacahAlgebra.Csq(j4,j5,j11),j6,j12), RacahAlgebra.Csq(j7,j8,j13),j14), j15)
    rightCsq = RacahAlgebra.Csq( RacahAlgebra.Csq( RacahAlgebra.Csq(j1,j4,j16), j7, j17),
                                 RacahAlgebra.Csq( RacahAlgebra.Csq( RacahAlgebra.Csq(j2,j5,j18), RacahAlgebra.Csq(j8,j3,j19), j20), j6,j21), j15)
    wa       = RacahAlgebra.evaluate(leftCsq, rightCsq)
end


text:
const Theme = Attributes

Base.broadcastable(x::AbstractScene) = Ref(x)
Base.broadcastable(x::AbstractPlot) = Ref(x)
Base.broadcastable(x::Attributes) = Ref(x)

# The rules that we use to convert values to a Observable in Attributes
value_convert(x::Observables.AbstractObservable) = Observables.observe(x)
value_convert(@nospecialize(x)) = x

# We transform a tuple of observables into a Observable(tuple(values...))
function value_convert(x::NTuple{N, Union{Any, Observables.AbstractObservable}}) where N
    result = Observable(to_value.(x))
    onany((args...)-> args, x...)
    return result
end

value_convert(x::NamedTuple) = Attributes(x)

# Version of `convert(Observable{Any}, obj)` that doesn't require runtime dispatch
node_any(@nospecialize(obj)) = isa(obj, Observable{Any}) ? obj :
                               isa(obj, Observable) ? convert(Observable{Any}, obj) : Observable{Any}(obj)

node_pairs(pair::Union{Pair, Tuple{Any, Any}}) = (pair[1] => node_any(value_convert(pair[2])))
node_pairs(pairs) = (node_pairs(pair) for pair in pairs)

Attributes(; kw_args...) = Attributes(Dict{Symbol, Observable}(node_pairs(kw_args)))
Attributes(pairs::Pair...) = Attributes(Dict{Symbol, Observable}(node_pairs(pairs)))
Attributes(pairs::AbstractVector) = Attributes(Dict{Symbol, Observable}(node_pairs.(pairs)))
Attributes(pairs::Iterators.Pairs) = Attributes(collect(pairs))
Attributes(nt::NamedTuple) = Attributes(; nt...)
attributes(x::Attributes) = getfield(x, :attributes)
Base.keys(x::Attributes) = keys(x.attributes)
Base.values(x::Attributes) = values(x.attributes)
function Base.iterate(x::Attributes, state...)
    s = iterate(keys(x), state...)
    s === nothing && return nothing
    return (s[1] => x[s[1]], s[2])
end

function Base.copy(attributes::Attributes)
    result = Attributes()
    for (k, v) in attributes
        # We need to create a new Signal to have a real copy
        result[k] = copy(v)
    end
    return result
end

function Base.deepcopy(obs::Observable)
    return Observable{Any}(to_value(obs))
end

function Base.deepcopy(attributes::Attributes)
    result = Attributes()
    for (k, v) in attributes
        # We need to create a new Signal to have a real copy
        result[k] = deepcopy(v)
    end
    return result
end

Base.filter(f, x::Attributes) = Attributes(filter(f, attributes(x)))
Base.empty!(x::Attributes) = (empty!(attributes(x)); x)
Base.length(x::Attributes) = length(attributes(x))

function Base.merge!(target::Attributes, args::Attributes...)
    for elem in args
        merge_attributes!(target, elem)
    end
    return target
end

Base.merge(target::Attributes, args::Attributes...) = merge!(copy(target), args...)

@generated hasfield(x::T, ::Val{key}) where {T, key} = :($(key in fieldnames(T)))

@inline function Base.getproperty(x::T, key::Symbol) where T <: Union{Attributes, Transformable}
    if hasfield(x, Val(key))
        getfield(x, key)
    else
        getindex(x, key)
    end
end

@inline function Base.setproperty!(x::T, key::Symbol, value) where T <: Union{Attributes, Transformable}
    if hasfield(x, Val(key))
        setfield!(x, key, value)
    else
        setindex!(x, value, key)
    end
end

function Base.getindex(x::Attributes, key::Symbol)
    x = attributes(x)[key]
    # We unpack Attributes, even though, for consistency, we store them as nodes
    # this makes it easier to create nested attributes
    return x[] isa Attributes ? x[] : x
end

function Base.setindex!(x::Attributes, value, key::Symbol)
    if haskey(x, key)
        x.attributes[key][] = value
    else
        x.attributes[key] = node_any(value)
    end
end

function Base.setindex!(x::Attributes, value::Observable, key::Symbol)
    if haskey(x, key)
        # error("You're trying to update an attribute Observable with a new Observable. This is not supported right now.
        # You can do this manually like this:
        # lift(val-> attributes[$key] = val, Observable::$(typeof(value)))
        # ")
        return x.attributes[key] = node_any(value)
    else
        #TODO make this error. Attributes should be sort of immutable
        return x.attributes[key] = node_any(value)
    end
    return x
end

_indent_attrs(s, n) = join(split(s, '\n'), "\n" * " "^n)

function Base.show(io::IO,::MIME"text/plain", attr::Attributes)

    io = IOContext(io, :compact => true)

    d = Dict()
    print(io, """Attributes with $(length(attr)) $(length(attr) != 1 ? "entries" : "entry")""")

    if length(attr) < 1
        return
    end

    print(io, ":")

    ks = sort(collect(keys(attr)), by = lowercase ∘ String)
    maxlength = maximum(length ∘ String, ks)

    for k in ks
        print(io, "\n  ")
        print(io, k)
        print(io, " => ")
        v = to_value(attr[k])
        if v isa Attributes
            print(io, _indent_attrs(repr(v), 2))
        else
            print(io, to_value(attr[k]))
        end
    end
end

Base.show(io::IO, attr::Attributes) = show(io, MIME"text/plain"(), attr)
theme(x::AbstractPlot) = x.attributes
isvisible(x) = haskey(x, :visible) && to_value(x[:visible])

#dict interface
const AttributeOrPlot = Union{AbstractPlot, Attributes}
Base.pop!(x::AttributeOrPlot, args...) = pop!(x.attributes, args...)
Base.haskey(x::AttributeOrPlot, key) = haskey(x.attributes, key)
Base.delete!(x::AttributeOrPlot, key) = delete!(x.attributes, key)
function Base.get!(f::Function, x::AttributeOrPlot, key::Symbol)
    if haskey(x, key)
        return x[key]
    else
        val = f()
        x[key] = val
        return x[key]
    end
end
Base.get!(x::AttributeOrPlot, key::Symbol, default) = get!(()-> default, x, key)
Base.get(f::Function, x::AttributeOrPlot, key::Symbol) = haskey(x, key) ? x[key] : f()
Base.get(x::AttributeOrPlot, key::Symbol, default) = get(()-> default, x, key)

# This is a bit confusing, since for a plot it returns the attribute from the arguments
# and not a plot for integer indexing. But, we want to treat plots as "atomic"
# so from an interface point of view, one should assume that a plot doesn't contain subplots
# Combined plots break this assumption in some way, but the way to look at it is,
# that the plots contained in a Combined plot are not subplots, but _are_ actually
# the plot itself.
Base.getindex(plot::AbstractPlot, idx::Integer) = plot.converted[idx]
Base.getindex(plot::AbstractPlot, idx::UnitRange{<:Integer}) = plot.converted[idx]
Base.setindex!(plot::AbstractPlot, value, idx::Integer) = (plot.input_args[idx][] = value)
Base.length(plot::AbstractPlot) = length(plot.converted)

function Base.getindex(x::AbstractPlot, key::Symbol)
    argnames = argument_names(typeof(x), length(x.converted))
    idx = findfirst(isequal(key), argnames)
    if idx === nothing
        return x.attributes[key]
    else
        x.converted[idx]
    end
end

function Base.getindex(x::AttributeOrPlot, key::Symbol, key2::Symbol, rest::Symbol...)
    dict = to_value(x[key])
    dict isa Attributes || error("Trying to access $(typeof(dict)) with multiple keys: $key, $key2, $(rest)")
    dict[key2, rest...]
end

function Base.setindex!(x::AttributeOrPlot, value, key::Symbol, key2::Symbol, rest::Symbol...)
    dict = to_value(x[key])
    dict isa Attributes || error("Trying to access $(typeof(dict)) with multiple keys: $key, $key2, $(rest)")
    dict[key2, rest...] = value
end

function Base.setindex!(x::AbstractPlot, value, key::Symbol)
    argnames = argument_names(typeof(x), length(x.converted))
    idx = findfirst(isequal(key), argnames)
    if idx === nothing && haskey(x.attributes, key)
        return x.attributes[key][] = value
    elseif !haskey(x.attributes, key)
        x.attributes[key] = convert(Observable, value)
    else
        return setindex!(x.converted[idx], value)
    end
end

function Base.setindex!(x::AbstractPlot, value::Observable, key::Symbol)
    argnames = argument_names(typeof(x), length(x.converted))
    idx = findfirst(isequal(key), argnames)
    if idx === nothing
        if haskey(x, key)
            # error("You're trying to update an attribute Observable with a new Observable. This is not supported right now.
            # You can do this manually like this:
            # lift(val-> attributes[$key] = val, Observable::$(typeof(value)))
            # ")
            return x.attributes[key] = value
        else
            return x.attributes[key] = value
        end
    else
        return setindex!(x.converted[idx], value)
    end
end

# a few shortcut functions to make attribute conversion easier
function get_attribute(dict, key)
    convert_attribute(to_value(dict[key]), Key{key}())
end

function merge_attributes!(input::Attributes, theme::Attributes)
    for (key, value) in theme
        if !haskey(input, key)
            input[key] = copy(value)
        else
            current_value = input[key]
            if value isa Attributes && current_value isa Attributes
                # if nested attribute, we merge recursively
                merge_attributes!(current_value, value)
            elseif value isa Attributes || current_value isa Attributes
                error("""
                Type missmatch while merging plot attributes with theme for key: $(key).
                Found $(value) in theme, while attributes contains: $(current_value)
                """)
            else
                # we're good! input already has a value, can ignore theme
            end
        end
    end
    return input
end


text:
"""
    utils.jl

Description:
    Common utility methods for the ClusterValidityIndices.jl package.
"""

"""
    sort_cvi_data(data::RealMatrix, labels::IntegerVector)

Sorts the CVI data by the label index, assuring that clusters are provided incrementally.
"""
function sort_cvi_data(data::RealMatrix, labels::IntegerVector)
    index = sortperm(labels)
    data = data[:, index]
    labels = labels[index]

    return data, labels
end # sort_cvi_data(data::RealMatrix, labels::IntegerVector)

"""
    relabel_cvi_data(labels::IntegerVector)

Relabels the vector to present new labels in incremental order.
"""
function relabel_cvi_data(labels::IntegerVector)
    # Get the unique labels and their order of appearance
    unique_labels = unique(labels)
    n_unique_labels = length(unique_labels)
    n_labels = length(labels)

    # Create a new ordered list of unique labels
    new_unique_labels = [x for x in 1:n_unique_labels]

    # Map the old unique labels to the new ones
    label_mapping = Dict(zip(unique_labels, new_unique_labels))

    # Create a new labels vector with ordered labels
    new_labels = zeros(Int, n_labels)
    for ix = 1:n_labels
        new_labels[ix] = label_mapping[labels[ix]]
    end

    return new_labels
end # relabel_cvi_data(labels::IntegerVector)


text:
#centrallimit.jl
#Algorithm 7.2 of Random Eigenvalues by Alan Edelman

using Winston

function centrallimit(m, t)
    # sums m random matrices and histograms their eigenvalues
    # normalized by 1/sqrt(m)
    # Input :
    #  m : number of matrices to sum
    #  t : number of trials to perform
    n = 100     # matrix size
    dx = .1     # bin size
    v = Float64[]
    # we choose a random diagonal {-1,1} matrx as our prototype
    # A = diagm(sign(randn(1,n)))
    # proto = kron([1 0; 0 -1], eye(n/2))
    for i = 1:t                         # for each trial
        B = zeros(n,n)                  # initialize the accumulator
        for z = 1:m                     # sum up the random matrices
            Q = full(QR(randn(n,n) + im*randn(n,n))[:Q])  # (piecewise) Haar orthogonal
            A = diagm(sign(randn(n,1)))
            B = B + Q'*A*Q              # A and Q'AQ are asymptotically free
        end
        append!(v, real(eigvals(B)/sqrt(m)))  # normalize eigenvalues
    end
    grid, count = hist(v, -2-dx/2:dx:2)

    p = FramedPlot()
    h = Histogram(count*pi/2/(dx*n*t), step(grid))
    h.x0 = first(grid)
    add(p, h)
    if isinteractive()
        Winston.display(p)
    else
        file(p, "centrallimit_m=$(m)_t=$(t).png")
    end
end


text:
using Geodesics
using Test

@testset "Forward" begin
    @test_throws ArgumentError Geodesics.forward(0, 2π, 0, 0, 1, 0)
    @test_throws ArgumentError Geodesics.forward(0, 2π, 0, 0, 0, 0)
    @test_throws ArgumentError Geodesics.inverse(0, 0, 1, 1, 1, 1)

    # Compare to Haversine computation on the unit sphere
    let lon0 = 0, lat0 = 0, azimuth = 1, distance = 1, a = 1, f = 0
        lon1, lat1, baz = Geodesics.forward(lon0, lat0, azimuth, distance, a, f)
        @test lon1 ≈ 0.9189892552937267
        @test lat1 ≈ 0.4719777676633856
        @test baz ≈ mod(-1.9047280216911358, 2π) atol=1e-6
    end

    # Make sure if we don't move, we end up in the same place
    let lon0 = 0.1, lat0 = 0.2, azimuth = 1, distance = 0, a = 1, f = 0
        lon1, lat1, baz = Geodesics.forward(lon0, lat0, azimuth, distance, a, f)
        @test lon1 == lon0
        @test lat1 == lat0
        @test baz ≈ mod(azimuth + π, 2π)
    end
end

@testset "Inverse" begin
    @test_throws ArgumentError Geodesics.inverse(0, 0, 1, 1, 0, 0)
    @test_throws ArgumentError Geodesics.inverse(0, 2π, 0, -3π, 1, 1)
    @test_throws ArgumentError Geodesics.inverse(0, 0, 1, 1, 1, 1)

    # Compare to Haversine computation on the unit sphere
    let lon0 = 0, lat0 = 0, lon1 = 1, lat1 = 1, a = 1, f = 0, atol = 1e-6
        distance, azimuth, backazimuth = Geodesics.inverse(lon0, lat0, lon1, lat1, a, f)
        @test distance ≈ 1.2745557823062943 atol=atol
        @test azimuth ≈ 0.49536728921867335 atol=atol
        @test backazimuth ≈ mod(-2.0661636160135703, 2π) atol=atol
    end
end

@testset "Ang dist" begin
    # Degrees
    let lon0 = 360*rand() - 180, lat0 = 180*rand() - 90, lon1 = 360*rand() - 180,
            lat1 = 180*rand() - 90, atol = 1e-6
        @test Geodesics.angular_distance(lon0, lat0, lon1, lat1, f=0) ≈
            rad2deg(Geodesics.inverse(deg2rad.((lon0, lat0, lon1, lat1))..., 1, 0)[1]) atol=atol
        @test Geodesics.angular_distance(lon0, lat0, lon1, lat1) ≈
            rad2deg(Geodesics.inverse(deg2rad.((lon0, lat0, lon1, lat1))..., 1, Geodesics.F_WGS84)[1]) atol=atol
    end

    # Radians
    let lon0 = 2π*(rand() - 0.5), lat0 = π*(rand() - 0.5), lon1 = 2π*(rand() - 0.5),
            lat1 = π*(rand() - 0.5), atol = 1e-6
        @test Geodesics.angular_distance(lon0, lat0, lon1, lat1, false, f=0) ≈
            Geodesics.inverse(lon0, lat0, lon1, lat1, 1, 0)[1] atol=atol
        @test Geodesics.angular_distance(lon0, lat0, lon1, lat1, false) ≈
            Geodesics.inverse(lon0, lat0, lon1, lat1, 1, Geodesics.F_WGS84)[1] atol=atol
    end
end

@testset "Azimuth" begin
    # Degrees
    let lon0 = 360*rand() - 180, lat0 = 180*rand() - 90, lon1 = 360*rand() - 180,
            lat1 = 180*rand() - 90, atol = 1e-6
        @test Geodesics.azimuth(lon0, lat0, lon1, lat1, f=0) ≈
            rad2deg(Geodesics.inverse(deg2rad.((lon0, lat0, lon1, lat1))..., 1, 0)[2]) atol=atol
        @test Geodesics.azimuth(lon0, lat0, lon1, lat1) ≈
            rad2deg(Geodesics.inverse(deg2rad.((lon0, lat0, lon1, lat1))..., 1, Geodesics.F_WGS84)[2]) atol=atol
    end

    # Radians
    let lon0 = 2π*(rand() - 0.5), lat0 = π*(rand() - 0.5), lon1 = 2π*(rand() - 0.5),
            lat1 = π*(rand() - 0.5), atol = 1e-6
        @test Geodesics.azimuth(lon0, lat0, lon1, lat1, false, f=0) ≈
            Geodesics.inverse(lon0, lat0, lon1, lat1, 1, 0)[2] atol=atol
        @test Geodesics.azimuth(lon0, lat0, lon1, lat1, false) ≈
            Geodesics.inverse(lon0, lat0, lon1, lat1, 1, Geodesics.F_WGS84)[2] atol=atol
    end

    # One point at pole
    let lon0 = 0, lat0 = -90, lon1 = 0, lat1 = 0
        for T in (Int16, Int32, Int64, Float16, Float32, Float64)
            @test Geodesics.azimuth(T(lon0), T(lat0), T(lon1), T(lat1), true) == 0.0
        end
    end
    let lon0 = 90, lat0 = 0, lon1 = 90, lat1 = 90
        for T in (Int16, Int32, Int64, Float16, Float32, Float64)
            @test Geodesics.azimuth(T(lon0), T(lat0), T(lon1), T(lat1), true) == 0.0
        end
    end
end

@testset "Ang step" begin
    # Degrees
    let lon0 = 360*rand() - 180, lat0 = 180*rand() - 90, dist = 360*rand(),
            az = 360*rand() - 90, atol = 1e-6
        @test all(isapprox.(Geodesics.angular_step(lon0, lat0, az, dist, f=0),
            rad2deg.(Geodesics.forward(deg2rad.((lon0, lat0, az, dist))..., 1, 0)), atol=atol))
        @test all(isapprox.(Geodesics.angular_step(lon0, lat0, az, dist),
            rad2deg.(Geodesics.forward(deg2rad.((lon0, lat0, az, dist))..., 1, Geodesics.F_WGS84)), atol=atol))
    end

    # Radians
    let lon0 = 2π*(rand() - 0.5), lat0 = π*(rand() - 0.5), dist = 2π*rand(),
            az = 2π*(rand() - 0.5), atol = 1e-6
        @test all(isapprox.(Geodesics.angular_step(lon0, lat0, az, dist, false, f=0),
            Geodesics.forward(lon0, lat0, az, dist, 1, 0), atol=atol))
        @test all(isapprox.(Geodesics.angular_step(lon0, lat0, az, dist, false),
            Geodesics.forward(lon0, lat0, az, dist, 1, Geodesics.F_WGS84), atol=atol))
    end
end

@testset "Surf dist" begin
    # Degrees
    let lon0 = 360*rand() - 180, lat0 = 180*rand() - 90, lon1 = 360*rand() - 180,
            lat1 = 180*rand() - 90, a = 10_000_000*rand(), atol = 1e-6
        @test Geodesics.surface_distance(lon0, lat0, lon1, lat1, a, f=0) ≈
            Geodesics.inverse(deg2rad.((lon0, lat0, lon1, lat1))..., a, 0)[1] atol=atol
        @test Geodesics.surface_distance(lon0, lat0, lon1, lat1, a) ≈
            Geodesics.inverse(deg2rad.((lon0, lat0, lon1, lat1))..., a, Geodesics.F_WGS84)[1] atol=atol
    end

    # Radians
    let lon0 = 2π*(rand() - 0.5), lat0 = π*(rand() - 0.5), lon1 = 2π*(rand() - 0.5),
            lat1 = π*(rand() - 0.5), a = 10_000_000*rand(), atol = 1e-6
        @test Geodesics.surface_distance(lon0, lat0, lon1, lat1, a, false, f=0) ≈
            Geodesics.inverse(lon0, lat0, lon1, lat1, a, 0)[1] atol=atol
        @test Geodesics.surface_distance(lon0, lat0, lon1, lat1, a, false) ≈
            Geodesics.inverse(lon0, lat0, lon1, lat1, a, Geodesics.F_WGS84)[1] atol=atol
    end
end


text:
using Mitosis
using Test
const TEST = true
include("gauss.jl")
include("statespace.jl")

@testset "Example" begin
    include("../example/example.jl")
end


text:
using Cryptoworld
using Test

@testset "Cryptoworld.jl" begin
    # Write your tests here.
end


text:
# Cross Validation procedures
using Dates
using LinearAlgebra
using DelimitedFiles
using Random

struct OneFoldSpec
    in_sample_ids::Set
    out_of_sample_ids::Set
end

struct CV_multiple_models_spec
    all_models::Any
    folds_spec::Array{OneFoldSpec, 1}
    nr_burnin::Int
    nr_samples::Int
    inter_collect_interval::Int
    top_n_fraction::Float64
    nr_repetitions::Int
    log_path::String
end

struct CV_result
    model::Any  # the algorith that is used to produce this result
    in_sample_LL_list::Array{Float64,2}
    out_of_sample_LL_list::Array{Float64,2}
    per_fold_avgMAPs::Array{Float64,2}
    per_fold_MAPs::Array{Float64,2}
    start_time::DateTime
    end_time::DateTime
    duration::Dates.Millisecond
    LL_out_subjectwise::Array{Float64,2}
end

struct CV_multiple_models_result
    all_results::Array{CV_result,1}
    start_time::DateTime
    end_time::DateTime
    duration::Dates.Millisecond
end

function get_k_folds(all_subj_ids::Set, out_of_sample_subj_sets::Array)
    k_folds = length(out_of_sample_subj_sets)
    folds_spec = Array{OneFoldSpec, 1}(undef, k_folds)
    println(typeof(folds_spec))
    for k in 1:k_folds
        in_sample = copy(all_subj_ids)
        out_of_sample = out_of_sample_subj_sets[k]
        for i in out_of_sample delete!(in_sample, i) end
        out_of_sample = Set(out_of_sample)
        oneFspec = OneFoldSpec(in_sample,out_of_sample)
        folds_spec[k] = oneFspec
    end
    folds_spec
end

function get_k_folds_all_out(all_subj_ids::Set, k_folds::Int)
    nr_subj = length(all_subj_ids)
    if k_folds > nr_subj
        error("Error in get_k_folds_all_out: k_folds = $k_folds > nr_subj=$nr_subj")
    end
    if k_folds <= 1
        error("Error in get_k_folds_all_out: k_folds = $k_folds <=1")
    end

    out_of_sample_sets = Array{Set, 1}(undef, k_folds)
    for k = 1:k_folds
        out_of_sample_sets[k] = Set()
    end

    s1=randperm(nr_subj)
    out_of_sample_subj_ids = collect(all_subj_ids)[s1]
    println(out_of_sample_subj_ids)
    for i = 1:nr_subj
        i_fold = 1+((i-1) % k_folds)
        push!(out_of_sample_sets[i_fold], out_of_sample_subj_ids[i])
    end
    get_k_folds(all_subj_ids,  out_of_sample_sets)
end

function get_k_folds(all_subj_ids::Set, k_folds::Int, out_of_sample_N::Int)
    # produces k folds, each having the same number (out_of_sample_N ) of subjects
    # in the out of sample set. Note: It is possible that not ALL subjects are in
    # one of the out-of-sample sets.
    # Use get_k_folds_all_out to get a specification which has ALL subjects out
    # of sample. In that case, not all sets are of the same size.
    nr_subj = length(all_subj_ids)
    max_i = k_folds * out_of_sample_N
    if max_i > nr_subj
        error("Error in get_k_folds: k_folds * out_of_sample_N = $max_i > nr_subj=$nr_subj")
    end
    s1=randperm(nr_subj)[1:max_i]
    out_of_sample_subj_ids = collect(all_subj_ids)[s1]
    out_of_sample_subj_pairs = reshape(out_of_sample_subj_ids,(k_folds,out_of_sample_N))
    # println(out_of_sample_subj_pairs)
    out_of_sample_sets = Array{Set, 1}(undef, k_folds)
    for k = 1:k_folds
        out_of_sample_sets[k] = Set(out_of_sample_subj_pairs[k,:])
    end
    get_k_folds(all_subj_ids, k_folds, out_of_sample_sets)
end

function run_cross_validation(description::String, data::Array{Int,2},
                        cv_spec::CV_multiple_models_spec, rl_module_path::String)
    cv_start_time = now()
    log_buffer = IOBuffer()
    log_print_ln!(log_buffer, "\nCROSS VALIDATION: start: $cv_start_time")
    log_print_ln!(log_buffer, "description: $description")
    data_string = RL_Fit.to_short_string(RL_Fit.get_data_stats(data))
    log_print_ln!(log_buffer, "Data:\n$data_string")
    nr_models = length(cv_spec.all_models)
    log_print_ln!(log_buffer, "nr of models: $nr_models")

    for i_model = 1:nr_models
        current_model = cv_spec.all_models[i_model]
        log_print_ln!(log_buffer, "********************************************************************************")
        log_print_ln!(log_buffer, "*********       Start cross validation of model $(current_model.model_name)       ************")
        log_print_ln!(log_buffer, "********************************************************************************")
        cv_result_model = run_CV_single_model(
            data,
            current_model,
            cv_spec.folds_spec,
            cv_spec.nr_burnin,
            cv_spec.nr_samples,
            cv_spec.inter_collect_interval,
            cv_spec.top_n_fraction,
            cv_spec.nr_repetitions,
            log_buffer,
            rl_module_path
        )
        write_CV_result(cv_result_model, cv_spec.folds_spec)
        log_print_ln!(log_buffer, "summed out-of-sample LL: $(sum(cv_result_model.out_of_sample_LL_list))")
        log_print_ln!(log_buffer, "END cross validation of model $(current_model.model_name), duration: $(cv_result_model.duration)")
        log_print_ln!(log_buffer, "================================================================================")
    end

     log_file_name = "$(cv_spec.log_path)/CV_log_$(Dates.format(now(), "dd_u_yyyy_HH_MM")).txt"

     log_str = String(take!(log_buffer))
     open(log_file_name, "w") do f
         write(f, log_str)
     end
end

function run_CV_single_model(data::Array{Int,2}, model_spec::Any,
                            folds_spec::Array{OneFoldSpec,1},
                            nr_burnin::Int,
                            nr_samples::Int,
                            inter_collect_interval::Int,
                            top_n_fraction::Float64,
                            nr_repetitions::Int,
                            log_buffer::IOBuffer,
                            rl_module_path::String
                            )
    start_time=now()
    nr_of_folds = length(folds_spec)
    nr_of_params = length(model_spec.model_params)
    in_sample_LL_list = zeros(Float64,nr_of_folds, 2)
    out_of_sample_LL_list = zeros(Float64, nr_of_folds, 2)
    per_fold_avgMAPs = Array{Float64, 2}(undef, (nr_of_folds, nr_of_params+2) )  # +2: space for LL and NR of subj
    per_fold_MAPs = Array{Float64, 2}(undef, (nr_of_folds, nr_of_params+2) )  # +2: space for LL and NR of subj

    allSubj = sort(unique(data[:, Int(c_subject)]))
    nr_subjects = length(allSubj)
    LL_out_subjectwise = zeros(Float64, nr_of_folds, nr_subjects)

    for k in 1:nr_of_folds
       log_print_ln!(log_buffer, "   START FOLD #$k")
       in_sample_data = RL_Fit.filter_by_person(data, folds_spec[k].in_sample_ids)
       out_of_sample_data = RL_Fit.filter_by_person(data, folds_spec[k].out_of_sample_ids)
        log_print_ln!(log_buffer, "in-sample_subjects: $(folds_spec[k].in_sample_ids)")
        log_print_ln!(log_buffer, "out-of-sample_subjects: $(folds_spec[k].out_of_sample_ids)")
        mcmc_spec = Metropolis_Spec(
            model_spec.learner_function_V,
            model_spec.model_params,
            nr_burnin,
            nr_samples,
            inter_collect_interval)

        mcmc_samples = RL_Fit.run_parallel_Metropolis(in_sample_data, mcmc_spec, rl_module_path; nr_of_repetitions=nr_repetitions)
        nr_map_samples = Int(ceil(size(mcmc_samples)[1] * top_n_fraction))
        log_print_ln!(log_buffer, "sampling data size: $(size(mcmc_samples)), top_N: $nr_map_samples")
        top_MAP, top_n_avgMAP, top_n_min, top_n_max, top_n_stdMAP, top_n_samples = RL_Fit.get_MAP(mcmc_samples, nr_map_samples)

        in_sample_params = top_MAP[1:end-1]
        per_fold_avgMAPs[k,:] = hcat(top_n_avgMAP, length(folds_spec[k].in_sample_ids))

        top_MAP = reshape(top_MAP, 1, :) # Convert vector (nparam,) to array (nparam x 1)
        per_fold_MAPs[k,:] = hcat(top_MAP, length(folds_spec[k].in_sample_ids))

        LL_in = model_spec.learner_function_V(in_sample_data, in_sample_params)[1]

        in_sample_LL_list[k,:] = [LL_in length(folds_spec[k].in_sample_ids)]

        # get out of sample LL
        LL_out = model_spec.learner_function_V(out_of_sample_data, in_sample_params)[1]
        out_of_sample_LL_list[k,:] = [LL_out  length(folds_spec[k].out_of_sample_ids)]

        # Save the LL of each of the out-of-sample subjects (to be used for Bayesian Model Selection)
        for out_subj in folds_spec[k].out_of_sample_ids
            out_of_sample_data_subject = RL_Fit.filter_by_person(data, out_subj)
            LL_out_subjectwise[k, out_subj] = model_spec.learner_function_V(out_of_sample_data_subject, in_sample_params)[1]
        end

        log_print_ln!(log_buffer, "fold #$k result: top_MAP=$(round.(top_MAP, digits=3)), out of sample LL=$(round(LL_out, digits=3))")
        log_print_ln!(log_buffer, "END of fold #$k. timestamp: $(now())")
    end
    end_time = now()
    results =  CV_result(
        model_spec,
        in_sample_LL_list,
        out_of_sample_LL_list,
        per_fold_avgMAPs,
        per_fold_MAPs,
        start_time,
        end_time,
        end_time-start_time,
        LL_out_subjectwise)
end

function write_CV_result(cv_res::CV_result, f_spec)
    println("Write results to file")

    cv_path = "$(cv_res.model.result_folder)/cross_validation"
    if !ispath(cv_path)
        println("mkpath($(cv_path))")
        mkpath(cv_path)
    end
    path = "$(cv_path)/per_fold_avgMAPs.csv"
    writedlm(path, cv_res.per_fold_avgMAPs)

    path = "$(cv_path)/per_fold_MAPs.csv"
    writedlm(path, cv_res.per_fold_MAPs)

    path = "$(cv_path)/in_sample_LL_list.csv"
    writedlm(path, cv_res.in_sample_LL_list)

    path = "$(cv_path)/out_of_sample_LL_list.csv"
    writedlm(path, cv_res.out_of_sample_LL_list)

    path = "$(cv_path)/out_of_sample_LL_SUM.csv"
    writedlm(path, sum(cv_res.out_of_sample_LL_list[:,1]))

    # keep track of the subject ID's
    for i_fold = 1:length(f_spec)
        path = "$(cv_path)/subj_id_In_$(i_fold).csv"
        writedlm(path, f_spec[i_fold].in_sample_ids)
        path = "$(cv_path)/subj_id_Out_$(i_fold).csv"
        writedlm(path, f_spec[i_fold].out_of_sample_ids)
    end

    path = "$(cv_path)/LL_out_subjectwise.csv"
    writedlm(path, cv_res.LL_out_subjectwise)

end


text:
using Echidna
using Test

print("Testing NSGAII on ZDT1", "\n")
nobjs = 2
nvars = 30
pop_size = 100
n_iters = 10
zdt1_problem = Problem(
    2, [false for i in 1:nobjs], nvars,
    [MOGA_Real(0.0, 1.0) for i in 1:30],
    ZDT1
)

config = NSGAII(zdt1_problem, ZDT1, pop_size, n_iters)
result = garun(config)
@test length(result) == pop_size

max_stars_repo_name: phawthorne/echidna

max_issues_repo_name: phawthorne/echidna

max_forks_repo_name: phawthorne/echidna


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule raptor_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("raptor")
JLLWrappers.@generate_main_file("raptor", UUID("a0baa3b2-8ece-51b8-9952-c0b61b0787fb"))
end  # module raptor_jll


text:
function rvrp_solution_unit_tests()

end


text:
"""
Function calculates decision rules for given prices and Utils
# Arguments
- `mdl::Model`: Ntuple containing coordinates
# Examples
```julia-repl
julia> iterate_decission_rules(mdl)
```
"""
function iterate_decission_rules(mdl::Model)

    # get number of state variables
    l_dim = length( mdl.GridDef.Grids ) ÷ 2

    # initialize DR
    mdl.DR = zeros( Int64, *( length.( mdl.GridDef.Grids[1:l_dim] )... ), length( mdl.GridDef.Grids[end] ) )

    #precompute VF for given prices and etc:
    precomputeVF(mdl)

    #initialize loop parameters
    iter, no_convergeance = (0, true)

    while (no_convergeance)

        #there must be a more efficient way to do the below:
        old_vf, old_dr = copy(mdl.VF), copy(mdl.DR)

        #find optimal decission rules/ vf matrix given old DR / old VF
        dr_CPU(mdl)

        #optional afterprocessing
        mdl.VF_afterprocess(mdl)

        # consitent strategy rule
        check1 = old_dr != mdl.DR

        #might concentrate with above in the future. Right now it is divided in case more condition would arrive.
        no_convergeance = (check1)

        #iteration increment
        iter += 1

        if (!no_convergeance)

            printstyled("\e[s\e[5G✓\e[u",bold = true, color = :green)

        end
        if (iter == vf_iterator_maxiter)

            no_convergeance = false
            printstyled("\e[s\e[5G!\e[u",bold = true, color = :red)
            global DR_punish = true

        end

    end

end


text:
@testset "Riemann matrix" begin
    τ = random_siegel(2);
    R = RiemannMatrix(τ);
    @test R.τ == τ
    @test R.g == 2
    @test R.X == real(τ)
    @test R.Y == imag(τ)
    @test size(R.ellipsoid)[1] == 5
    @test R.T[2,1] == 0
end


text:
# coding: utf-8

# In[1]:

#a simple day trading game
#day trader is only allowed to make at maximum two trades
#the strategy is long only
#lets find out the maximum profit

#more details can be found in the following link
# https://www.geeksforgeeks.org/maximum-profit-by-buying-and-selling-a-share-at-most-twice/

#an alternative version in recursion exists
#its done by using a different approach
#strongly recommend you to take a look
# https://github.com/je-suis-tm/recursion-and-dynamic-programming/blob/master/stock%20trading%20recursion.jl

# In[2]:

#there are two scenarios to maximize the profit
#one trade or two trades
#first we run a reverse iteration
#to obtain the maximum profit from one trade
#then we run a normal iteration
#to obtain the maximum profit
#from one trade plus the result from reverse iteration
function stock_trading(prices)

    #initialize the profit at zero
    profit=[0 for _ in 1:length(prices)]

    #initialize maximum price with the close price
    max_price=prices[end]

    #reverse order iteration
    for i in length(prices)-1:-1:1

        #update the maximum price to compute the maximum profit
        if prices[i]>max_price

            max_price=prices[i]

        end

        #two scenarios to get the maximum profit
        #either the previous iteration is larger
        #or this round of iteration
        profit[i]=max(profit[i+1],max_price-prices[i])

    end

    #initialize minimum price with the open price
    min_price=prices[1]

    #second round of iteration
    for i in 2:length(prices)

        #update the minimum price to compute the maximum profit
        if prices[i]<min_price

            min_price=prices[i]

        end

        #two scenarios to get the maximum profit
        #either the previous iteration is larger
        #or this round of iteration plus the result from single transaction
        profit[i]=max(profit[i-1],profit[i]+prices[i]-min_price)

    end

    return profit[end]

end

# In[3]:

stock_trading([10,22,5,75,65,80])

# In[4]:

stock_trading([2,30,15,10,8,25,80])

# In[5]:

stock_trading([100,30,15,10,8,25,80])

# In[6]:

stock_trading([90,70,35,11,5])

max_stars_repo_path: stock trading dynamic programming.jl

max_issues_repo_path: stock trading dynamic programming.jl

max_forks_repo_path: stock trading dynamic programming.jl


text:
alogab(a,b) = vifelse(andmask(a>0, b>0), a*log(a/b)-a+b, vifelse(andmask(iszero(a), iszero(b)), one(eltype(b)), convert(eltype(b), Inf)))
blogab(a,b) = vifelse(andmask(a>0, b>0), -b*log(a/b)+a-b, convert(eltype(b), Inf))

aloga(a) = vifelse(a>0, a*log(a)-a+one(eltype(a)), convert(eltype(a), Inf))
loga(a) = vifelse(a>0, -log(a)+a-one(eltype(a)), convert(eltype(a), Inf))

andmask(x::Bool, y::Bool) = x && y

γ(::KullbackLeibler, a::T, b::T) where T = alogab(a,b)
γ(::ReverseKullbackLeibler, a::T, b::T) where T = blogab(a,b)
γ(::Hellinger, a::T, b::T) where T = vifelse(andmask(a >= 0, b >= 0), 2*a+(2-4*sqrt(a/b))*b, convert(T, Inf))
γ(::ChiSquared, a::T, b::T) where T = half(T)*abs2(a - b)/b

function γ(d::CressieRead{D}, a::T, b::T) where {T, D}
    α = d.α
    u = a/b
    vifelse(andmask(a>0, b>0), (u^(1+α)+α-u*(1+α))*b/(α*(1+α)),
            α > 0 ? zero(eltype(a)) : convert(eltype(a), Inf))
end

γ(::KullbackLeibler, a::T) where T = aloga(a)
γ(::ReverseKullbackLeibler, a::T) where T = loga(a)
γ(::Hellinger, a::T) where T = vifelse(a>=0, 2*a - 4*sqrt(a) + 2, convert(T, Inf))
γ(::ChiSquared, a::T) where T = abs2(a - one(T))*half(T)

function γ(d::CressieRead{D}, a::T) where {T, D}
    α = d.α
    vifelse(a>=0, (a^(1 + α) + α - a*(1 + α))/(α*(1 + α)),
            α > 0 ? zero(eltype(a)) : convert(eltype(a), Inf))
end

∇ᵧ(::KullbackLeibler, a::T, b::T) where T = vifelse(andmask(a>0, b>0), log(a/b), convert(T, -Inf))
∇ᵧ(::ReverseKullbackLeibler, a::T, b::T) where T = vifelse(andmask(a>0, b>0), -b/a + one(T), convert(T, -Inf))
∇ᵧ(d::CressieRead, a::T, b::T) where T = vifelse(andmask(a>0, b>0), ((a/b)^d.α - one(T))/d.α, convert(T, sign(d.α)*Inf))
∇ᵧ(d::Hellinger, a::T, b::T) where T = vifelse(andmask(a>0, b>0), 2(one(T)-one(T)/sqrt(a/b)), convert(T, -Inf))
∇ᵧ(d::ChiSquared, a::T, b::T) where T = a/b - one(T)

∇ᵧ(::KullbackLeibler, a::T) where T = vifelse(a > 0, log(a), convert(T, -Inf))
∇ᵧ(::ReverseKullbackLeibler, a::T) where T = vifelse(a > 0, -1/a + one(T), convert(T, -Inf))
∇ᵧ(d::CressieRead, a::T) where T = vifelse(a >= 0, (a^d.α - 1)/d.α, convert(T, sign(d.α)*Inf))
∇ᵧ(d::Hellinger, a::T) where T = vifelse(a > 0, 2(one(T)-one(T)/sqrt(a)), convert(T, -Inf))
∇ᵧ(d::ChiSquared, a::T) where T = a - one(T)

Hᵧ(::KullbackLeibler, a::T, b::T) where T = vifelse(andmask(a > 0, b > 0), one(T)/a, convert(T, Inf))
Hᵧ(::ReverseKullbackLeibler, a::T, b::T) where T = vifelse(andmask(a > 0, b > 0), b/a^2, convert(T, Inf))
Hᵧ(d::CressieRead, a::T, b::T) where T = vifelse(andmask(a > 0, b > 0), a^(d.α-1)*b^(-d.α), convert(T, Inf))
Hᵧ(d::Hellinger, a::T, b::T) where T = vifelse(andmask(a > 0, b > 0), sqrt(b)/sqrt(a^3), convert(T, Inf))
Hᵧ(d::ChiSquared, a::T, b::T) where T = vifelse(b >= 0, 1/b, convert(T, Inf))

Hᵧ(::KullbackLeibler, a::T) where T = vifelse(a > 0, one(T)/a, convert(T, Inf))
Hᵧ(::ReverseKullbackLeibler, a::T) where T = vifelse(a > 0, one(T)/a^2, convert(T, Inf))
Hᵧ(d::CressieRead, a::T) where T = vifelse(a > 0, a^(d.α-1), convert(T, Inf))
Hᵧ(d::Hellinger, a::T) where T = vifelse(a > 0, one(T)/sqrt(a^(3)), convert(T, Inf))
Hᵧ(d::ChiSquared, a::T) where T = one(T)

## Modified Divergences Function
function γᵤ(d::D, a::T, b::T) where {T, D<:AbstractModifiedDivergence}
    @unpack γ₀, γ₁, γ₂, ρ  = d.m
    (γ₀ + γ₁*((a/b)-ρ) + half(T)*γ₂*(a/b-ρ)^2)*b
end

function γₗ(d::D, a::T, b::T) where {T, D<:AbstractModifiedDivergence}
    @unpack g₀, g₁, g₂, φ  = d.m
    (g₀ + g₁*((a/b)-φ) + half(T)*g₂*(a/b-φ)^2)*b
end

function γᵤ(d::D, a::T) where {T, D<:AbstractModifiedDivergence}
    @unpack γ₀, γ₁, γ₂, ρ  = d.m
    (γ₀ + γ₁*(a-ρ) + half(T)*γ₂*(a-ρ)^2)
end

function γₗ(d::D, a::T) where {T, D<:AbstractModifiedDivergence}
    @unpack g₀, g₁, g₂, φ  = d.m
    (g₀ + g₁*(a-φ) + half(T)*g₂*(a-φ)^2)
end

function ∇ᵤ(d::D, a::T, b::T) where {T, D<:AbstractModifiedDivergence}
    @unpack γ₀, γ₁, γ₂, ρ  = d.m
    γ₁ + (a/b)*γ₂ - γ₂*ρ
end

function ∇ₗ(d::D, a::T, b::T) where {T, D<:AbstractModifiedDivergence}
    @unpack g₀, g₁, g₂, φ  = d.m
    g₁ + (a/b)*g₂ - g₂*φ
end

function ∇ᵤ(d::D, a::T) where {T, D<:AbstractModifiedDivergence}
    @unpack γ₀, γ₁, γ₂, ρ  = d.m
    (γ₁ + γ₂*(a-ρ))
end

function ∇ₗ(d::D, a::T) where {T, D<:AbstractModifiedDivergence}
    @unpack g₀, g₁, g₂, φ  = d.m
    (g₁ + g₂*(a-φ))
end

function Hᵤ(d::D, a::T, b::T) where {T, D<:AbstractModifiedDivergence}
    @unpack γ₀, γ₁, γ₂, ρ  = d.m
    γ₂/b
end

function Hₗ(d::D, a::T, b::T) where {T, D<:AbstractModifiedDivergence}
    @unpack g₀, g₁, g₂, φ  = d.m
     g₂/b
end

function Hᵤ(d::D, a::T) where {T, D<:AbstractModifiedDivergence}
    @unpack γ₀, γ₁, γ₂, ρ  = d.m
    γ₂
end

function Hₗ(d::D, a::T) where {T, D<:AbstractModifiedDivergence}
    @unpack g₀, g₁, g₂, φ  = d.m
    g₂
end

##
function γ(d::ModifiedDivergence, a::T, b::T) where T
    @unpack ρ  = d.m
    div = d.d
    vifelse(a>ρ*b, γᵤ(d, a, b), γ(div, a, b))
end

function γ(d::ModifiedDivergence, a::T) where T
    @unpack ρ  = d.m
    div = d.d
    vifelse(a>ρ, γᵤ(d, a), γ(div, a))
end

function γ(d::FullyModifiedDivergence, a::T, b::T) where T
    @unpack ρ, φ = d.m
    div = d.d
    vifelse(a>ρ*b, γᵤ(d, a, b), vifelse(a<φ*b, γₗ(d, a, b),  γ(div, a, b)))
end

function γ(d::FullyModifiedDivergence, a::T) where T
    @unpack ρ, φ = d.m
    div = d.d
    vifelse(a>ρ, γᵤ(d, a), vifelse(a<φ, γₗ(d, a), γ(div, a)))
end

function ∇ᵧ(d::ModifiedDivergence, a::T, b::T) where T
    @unpack ρ  = d.m
    div = d.d
    vifelse(a>ρ*b, ∇ᵤ(d, a, b), ∇ᵧ(div, a, b))
end

function ∇ᵧ(d::ModifiedDivergence, a::T) where T
    @unpack ρ  = d.m
    div = d.d
    vifelse(a>ρ, ∇ᵤ(d, a), ∇ᵧ(div, a))
end

function ∇ᵧ(d::FullyModifiedDivergence, a::T, b::T) where T
    @unpack ρ, φ  = d.m
    div = d.d
    vifelse(a>ρ*b, ∇ᵤ(d, a, b),  vifelse(a<φ*b, ∇ₗ(d, a, b),  ∇ᵧ(div, a, b)))
end

function ∇ᵧ(d::FullyModifiedDivergence, a::T) where T
    @unpack ρ, φ = d.m
    div = d.d
    vifelse(a>ρ, ∇ᵤ(d, a),  vifelse(a<φ, ∇ₗ(d, a),  ∇ᵧ(div, a)))
end

function Hᵧ(d::ModifiedDivergence, a::T, b::T) where T
    @unpack ρ = d.m
    div = d.d
    vifelse(a>ρ*b, Hᵤ(d, a, b), Hᵧ(div, a, b))
end

function Hᵧ(d::ModifiedDivergence, a::T) where T
    @unpack ρ = d.m
    div = d.d
    vifelse(a>ρ, Hᵤ(d, a), Hᵧ(div, a))
end

function Hᵧ(d::FullyModifiedDivergence, a::T, b::T) where T
    @unpack ρ, φ = d.m
    div = d.d
    vifelse(a>ρ*b, Hᵤ(d, a, b), vifelse(a<φ*b, Hₗ(d, a, b),  Hᵧ(div, a, b)))
end

function Hᵧ(d::FullyModifiedDivergence, a::T) where T
    @unpack ρ, φ = d.m
    div = d.d
    vifelse(a>ρ, Hᵤ(d, a), vifelse(a<φ, Hₗ(d, a),  Hᵧ(div, a)))
end

eval(d::AbstractDivergence, a::T, b::T) where T<:Real = γ(d, a, b)
gradient(d::AbstractDivergence, a::T, b::T) where T<:Real = ∇ᵧ(d, a, b)
hessian(d::AbstractDivergence, a::T, b::T) where T<:Real = Hᵧ(d, a, b)

eval(d::AbstractDivergence, a::T) where T<:Real = γ(d, a)
gradient(d::AbstractDivergence, a::T) where T<:Real = ∇ᵧ(d, a)
hessian(d::AbstractDivergence, a::T) where T<:Real = Hᵧ(d, a)

## eval
function eval(d::AbstractDivergence, a::AbstractArray{T}, b::AbstractArray{T}) where T <: Real
    r = zero(T)
    @avx for i ∈ eachindex(a,b)
        r += γ(d, a[i], b[i])
    end
    return r
end

function eval(d::AbstractDivergence, a::AbstractArray{T}) where T <: Real
    r = zero(T)
    @avx for i ∈ eachindex(a)
        r += γ(d, a[i])
    end
    return r
end

## gradient
function gradient!(u::AbstractVector{T}, d::AbstractDivergence, a::AbstractArray{T}, b::AbstractArray{T}) where T <: Real
    @avx for i ∈ eachindex(a, b, u)
        u[i] = ∇ᵧ(d, a[i], b[i])
    end
    return u
end

function gradient!(u::AbstractVector{T}, d::AbstractDivergence, a::AbstractArray{T}) where T <: Real
    @avx for i ∈ eachindex(a, u)
        u[i] = ∇ᵧ(d, a[i])
    end
    return u
end

function gradient(d::AbstractDivergence, a::AbstractArray{T}, b::AbstractArray{T}) where T <: Real
    u = similar(a)
    gradient!(u, d, a, b)
end

function gradient(d::AbstractDivergence, a::AbstractArray{T}) where T <: Real
    u = similar(a)
    gradient!(u, d, a)
end

function gradient_sum(d::AbstractDivergence, a::AbstractArray{T}, b::AbstractArray{T}) where T <: Real
    r = zero(T)
    @avx for i ∈ eachindex(a,b)
        r += ∇ᵧ(d, a[i], b[i])
    end
    return r
end

function gradient_sum(d::AbstractDivergence, a::AbstractArray{T}) where T <: Real
    r = zero(T)
    @avx for i ∈ eachindex(a)
        r += ∇ᵧ(d, a[i])
    end
    return r
end

## hessian
function hessian!(u::AbstractVector{T}, d::AbstractDivergence, a::AbstractArray{T}, b::AbstractArray{T}) where T <: Real
    @avx for i ∈ eachindex(a, b, u)
        u[i] = Hᵧ(d, a[i], b[i])
    end
    return u
end

function hessian!(u::AbstractVector{T}, d::AbstractDivergence, a::AbstractArray{T}) where T <: Real
    @avx for i ∈ eachindex(a, u)
        u[i] = Hᵧ(d, a[i])
    end
    return u
end

function hessian(d::AbstractDivergence, a::AbstractArray{T}, b::AbstractArray{T}) where T <: Real
    u = similar(a)
    hessian!(u, d, a, b)
end

function hessian(d::AbstractDivergence, a::AbstractArray{T}) where T <: Real
    u = similar(a)
    hessian!(u, d, a)
end

function hessian_sum(d::AbstractDivergence, a::AbstractArray{T}, b::AbstractArray{T}) where T <: Real
    r = zero(T)
    @avx for i ∈ eachindex(a,b)
        r += Hᵧ(d, a[i], b[i])
    end
    return r
end

function hessian_sum(d::AbstractDivergence, a::AbstractArray{T}) where T <: Real
    r = zero(T)
    @avx for i ∈ eachindex(a)
        r += Hᵧ(d, a[i])
    end
    return r
end

half(::Type{T}) where T = convert(T, 0.5)
#checksize(a, b) = length(a) != length(b) ? throw(DimensionMismatch()) : nothing


text:
# This elixir transforms the setup of elixir_advection_basic to a rotated square.
# The nodal values of the initial condition and the exact solution are the same as
# in elixir_advection_basic.
# However, on this rotated mesh, the metric terms are non-trivial.
# The same errors as with elixir_advection_basic are expected (except for rounding errors).

using OrdinaryDiffEq
using Trixi

# Define new structs inside a module to allow re-evaluating the file.
# This module name needs to be unique among all examples, otherwise Julia will throw warnings
# if multiple test cases using the same module name are run in the same session.
module TrixiExtensionAdvectionRotated

using Trixi

# initial_condition_convergence_test transformed to the rotated rectangle
struct InitialConditionConvergenceTestRotated
  sin_alpha::Float64
  cos_alpha::Float64
end

function InitialConditionConvergenceTestRotated(alpha)
  sin_alpha, cos_alpha = sincos(alpha)

  InitialConditionConvergenceTestRotated(sin_alpha, cos_alpha)
end

function (initial_condition::InitialConditionConvergenceTestRotated)(x, t, equation::LinearScalarAdvectionEquation2D)
  sin_ = initial_condition.sin_alpha
  cos_ = initial_condition.cos_alpha

  # Rotate back to unit square

  # Clockwise rotation by α and translation by 1
  # Multiply with [  cos(α)  sin(α);
  #                 -sin(α)  cos(α)]
  x_rot = SVector(cos_ * x[1] + sin_ * x[2], -sin_ * x[1] + cos_ * x[2])
  a = equation.advection_velocity
  a_rot = SVector(cos_ * a[1] + sin_ * a[2], -sin_ * a[1] + cos_ * a[2])

  # Store translated coordinate for easy use of exact solution
  x_trans = x_rot - a_rot * t

  c = 1.0
  A = 0.5
  L = 2
  f = 1/L
  omega = 2 * pi * f
  scalar = c + A * sin(omega * sum(x_trans))

  return SVector(scalar)
end

end # module TrixiExtensionAdvectionRotated

import .TrixiExtensionAdvectionRotated

###############################################################################
# semidiscretization of the linear advection equation

alpha = pi * 0.1
initial_condition = TrixiExtensionAdvectionRotated.InitialConditionConvergenceTestRotated(alpha)
sin_ = initial_condition.sin_alpha
cos_ = initial_condition.cos_alpha
T = [cos_ -sin_; sin_ cos_]

advection_velocity = Tuple(T * [0.2, -0.7])
equations = LinearScalarAdvectionEquation2D(advection_velocity)

# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)

mapping(xi, eta) = T * SVector(xi, eta)

cells_per_dimension = (16, 16)

# Create curved mesh with 16 x 16 elements
mesh = StructuredMesh(cells_per_dimension, mapping)

# A semidiscretization collects data structures and functions for the spatial discretization
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

###############################################################################
# ODE solvers, callbacks etc.

# Create ODE problem with time span from 0.0 to 1.0
ode = semidiscretize(semi, (0.0, 1.0));

# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
# and resets the timers
summary_callback = SummaryCallback()

# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
analysis_callback = AnalysisCallback(semi, interval=100)

# The SaveSolutionCallback allows to save the solution to a file in regular intervals
save_solution = SaveSolutionCallback(interval=100,
                                     solution_variables=cons2prim)

# The StepsizeCallback handles the re-calculcation of the maximum Δt after each time step
stepsize_callback = StepsizeCallback(cfl=1.6)

# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback)

###############################################################################
# run the simulation

# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);

# Print the timer summary
summary_callback()


text:
using Regression
using Base.Test

srand(54321)  # ensure repeatable experiment
              # avoid CI issues

## auxiliary functions

_relerr(x::Array, x0::Array) = sumabs2(x - x0) / sumabs2(x0)
_corrdist(x::Array, x0::Array) = 1.0 - dot(vec(x), vec(x0)) / (vecnorm(x) * vecnorm(x0))

_classify(U::Matrix) = vec(mapslices(indmax, U, 1))

function _errrate(pred::PredictionModel, X::Matrix, θ::Array, θ0::Array)
    u0 = _classify(predict(pred, θ0, X))
    u  = _classify(predict(pred, θ, X))
    countnz(u .!= u0) / length(u0)
end

const _solvers = [GD(), BFGS(), LBFGS(6)]

function verify_solver(title, loss::Loss, data,
                       θg::Array, vcond::Function, thres::Real)
    println("    $title")

    X = data[1]
    dx = size(X, 1)
    dθ = size(θg, ndims(θg))
    b = dθ == dx ? 0.0 :
        dθ == dx + 1 ? 1.0 :
        error("Unmatched dimensions.")

    pb = ndims(θg) == 1 ? UnivariateRegression(loss, data...; bias=b) :
                          MultivariateRegression(loss, data...; bias=b)

    for solv in _solvers
        println("      - with solver $(solv)")
        n = size(X, ndims(X))
        ret = Regression.solve(pb;
                               reg=SqrL2Reg(1.0e-4),
                               solver=solv,
                               options=Regression.Options(grtol=n * 1.0e-6))

        @test isa(ret, Regression.Solution)
        θe = ret.sol
        @test vcond(θe, θg) < thres
    end
end

## Univariate loss

# data
d = 5
n = 2000
w = randn(d)
wa = [w; randn()]

X = randn(d, n)
u = X'w + 0.01 * randn(n)
ua = u .+ wa[d+1]

# test cases

verify_solver("LinearPred + SqrLoss",
    SqrLoss(), (X, u), w, _relerr, 1.0e-4)

verify_solver("AffinePred + SqrLoss",
    SqrLoss(), (X, ua), wa, _relerr, 1.0e-4)

verify_solver("LinearPred + LogisticLoss",
    SqrLoss(), (X, u), w, _corrdist, 1.0e-3)

verify_solver("AffinePred + LogisticLoss",
    SqrLoss(), (X, ua), wa, _corrdist, 1.0e-3)

## Multivariate loss

# data
d = 5
k = 3
W = randn(k, d)
Wa = [W randn(k)]

X = randn(d, n)
U = W * X + 0.01 * randn(k, n)
Ua = U .+ Wa[:,d+1]

# test cases

verify_solver("MvLinearPred + SumSqrLoss",
    SumSqrLoss(), (X, U, k), W, _relerr, 1.0e-4)

verify_solver("MvAffinePred + SumSqrLoss",
    SumSqrLoss(), (X, Ua, k), Wa, _relerr, 1.0e-4)

y = _classify(U)
verify_solver("MvLinearPred + MultiLogisticLoss",
    MultiLogisticLoss(), (X, y, k), W,
        (θ, θ0) -> _errrate(MvLinearPred(d, k), X, θ, θ0), 0.02)

ya = _classify(Ua)
verify_solver("MvAffinePred + MultiLogisticLoss",
    MultiLogisticLoss(), (X, ya, k), Wa,
        (θ, θ0) -> _errrate(MvAffinePred(d, k), X, θ, θ0), 0.02)


text:
# Wilcoxon.jl
# Wilcoxon signed rank tests
#
# Copyright (C) 2012   Simon Kornblith
#
# Permission is hereby granted, free of charge, to any person obtaining
# a copy of this software and associated documentation files (the
# "Software"), to deal in the Software without restriction, including
# without limitation the rights to use, copy, modify, merge, publish,
# distribute, sublicense, and/or sell copies of the Software, and to
# permit persons to whom the Software is furnished to do so, subject to
# the following conditions:
#
# The above copyright notice and this permission notice shall be
# included in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
# LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
# OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
# WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

export SignedRankTest, ExactSignedRankTest, ApproximateSignedRankTest

## COMMON SIGNED RANK

# Automatic exact/normal selection
"""
    SignedRankTest(x::AbstractVector{<:Real})
    SignedRankTest(x::AbstractVector{<:Real}, y::AbstractVector{<:Real})

Perform a Wilcoxon signed rank test of the null hypothesis that the distribution of `x`
(or the difference `x - y` if `y` is provided) has zero median against the alternative
hypothesis that the median is non-zero.

When there are no tied ranks and ≤50 samples, or tied ranks and ≤15 samples,
`SignedRankTest` performs an exact signed rank test. In all other cases,
`SignedRankTest` performs an approximate signed rank test. Behavior may be further
controlled by using [`ExactSignedRankTest`](@ref) or [`ApproximateSignedRankTest`](@ref)
directly.

Implements: [`pvalue`](@ref), [`confint`](@ref)
"""
function SignedRankTest(x::AbstractVector{T}) where T<:Real
    (W, ranks, signs, tie_adjustment, n, median) = signedrankstats(x)
    n_nonzero = length(ranks)
    if n_nonzero <= 15 || (n_nonzero <= 50 && tie_adjustment == 0)
        ExactSignedRankTest(x, W, ranks, signs, tie_adjustment, n, median)
    else
        ApproximateSignedRankTest(x, W, ranks, signs, tie_adjustment, n, median)
    end
end
SignedRankTest(x::AbstractVector{T}, y::AbstractVector{S}) where {T<:Real,S<:Real} = SignedRankTest(x - y)

# Get W and absolute ranks for signed rank test
function signedrankstats(x::AbstractVector{S}) where S<:Real
   nonzero_x = x[x .!= 0]
   (ranks, tieadj) = tiedrank_adj(abs.(nonzero_x))
   W = 0.0
   for i = 1:length(nonzero_x)
       if nonzero_x[i] > 0
           W += ranks[i]
       end
   end
   (W, ranks, nonzero_x .> 0, tieadj, length(x), median(x))
end

## EXACT WILCOXON SIGNED RANK TEST

struct ExactSignedRankTest{T<:Real} <: HypothesisTest
    vals::Vector{T} # original values
    W::Float64              # test statistic: Wilcoxon rank-sum statistic
    ranks::Vector{Float64}           # ranks without ties (zero values)
    signs::BitArray{1}      # signs of input of ranks
    tie_adjustment::Float64 # adjustment for ties
    n::Int                  # number of observations
    median::Float64         # sample median
end
"""
    ExactSignedRankTest(x::AbstractVector{<:Real}[, y::AbstractVector{<:Real}])

Perform a Wilcoxon exact signed rank U test of the null hypothesis that the distribution of
`x` (or the difference `x - y` if `y` is provided) has zero median against the alternative
hypothesis that the median is non-zero.

When there are no tied ranks, the exact p-value is computed using the `psignrank` function
from the `Rmath` package. In the presence of tied ranks, a p-value is computed by exhaustive
enumeration of permutations, which can be very slow for even moderately sized data sets.

Implements: [`pvalue`](@ref), [`confint`](@ref)
"""
ExactSignedRankTest(x::AbstractVector{T}) where {T<:Real} =
    ExactSignedRankTest(x, signedrankstats(x)...)
ExactSignedRankTest(x::AbstractVector{S}, y::AbstractVector{T}) where {S<:Real,T<:Real} =
    ExactSignedRankTest(x - y)

testname(::ExactSignedRankTest) = "Exact Wilcoxon signed rank test"
population_param_of_interest(x::ExactSignedRankTest) = ("Location parameter (pseudomedian)", 0, x.median) # parameter of interest: name, value under h0, point estimate
default_tail(test::ExactSignedRankTest) = :both

function show_params(io::IO, x::ExactSignedRankTest, ident)
    println(io, ident, "number of observations:      ", x.n)
    println(io, ident, "Wilcoxon rank-sum statistic: ", x.W)
    print(io, ident, "rank sums:                   ")
    show(io, [sum(x.ranks[x.signs]), sum(x.ranks[map(!, x.signs)])])
    println(io)
    println(io, ident, "adjustment for ties:         ", x.tie_adjustment)
end

# Enumerate all possible Wilcoxon rank-sum results for a given vector, determining left-
# and right-tailed p values
function signedrankenumerate(x::ExactSignedRankTest)
    le = 0
    gr = 0
    n = length(x.ranks)
    tot = 2^n
    for i = 0:tot-1
        # Interpret bits of i as signs to generate wp for all possible sign combinations
        Wp = 0
        b = i
        j = 1
        while b != 0
            Wp += (b & 1)*x.ranks[j]
            j += 1
            b >>= 1
        end
        le += Wp <= x.W
        gr += Wp >= x.W
    end
    (le/tot, gr/tot)
end

function pvalue(x::ExactSignedRankTest; tail=:both)
    check_tail(tail)

    n = length(x.ranks)
    if n == 0
        1.0
    elseif x.tie_adjustment == 0
        # Compute exact p-value using method from Rmath, which is fast but cannot account for ties
        if tail == :both
            if x.W <= n * (n + 1)/4
                2 * psignrank(x.W, n, true)
            else
                2 * psignrank(x.W - 1, n, false)
            end
        elseif tail == :left
            psignrank(x.W, n, true)
        else
            psignrank(x.W - 1, n, false)
        end
    else
        # Compute exact p-value by enumerating all possible ranks in the tied data
        if tail == :both
            min(1, 2 * minimum(signedrankenumerate(x)))
        elseif tail == :left
            first(signedrankenumerate(x))
        else
            last(signedrankenumerate(x))
        end
    end
end

StatsBase.confint(x::ExactSignedRankTest; level::Real=0.95, tail=:both) = calculate_ci(x.vals, level, tail=tail)

## APPROXIMATE SIGNED RANK TEST

struct ApproximateSignedRankTest{T<:Real} <: HypothesisTest
    vals::Vector{T} # original values
    W::Float64              # test statistic: Wilcoxon rank-sum statistic
    ranks::Vector{Float64} # ranks without ties (zero values)
    signs::BitArray{1}      # signs of input of ranks
    tie_adjustment::Float64 # adjustment for ties
    n::Int                  # number of observations
    median::Float64         # sample median
    mu::Float64             # normal approximation: mean
    sigma::Float64          # normal approximation: std
end
"""
    ApproximateSignedRankTest(x::AbstractVector{<:Real}[, y::AbstractVector{<:Real}])

Perform a Wilcoxon approximate signed rank U test of the null hypothesis that the
distribution of `x` (or the difference `x - y` if `y` is provided) has zero median against
the alternative hypothesis that the median is non-zero.

The p-value is computed using a normal approximation to the distribution of the signed rank
statistic:
```math
    \\begin{align*}
        μ & = \\frac{n(n + 1)}{4}\\\\
        σ & = \\frac{n(n + 1)(2 * n + 1)}{24} - \\frac{a}{48}\\\\
        a & = \\sum_{t \\in \\mathcal{T}} t^3 - t
    \\end{align*}
```
where ``\\mathcal{T}`` is the set of the counts of tied values at each tied position.

Implements: [`pvalue`](@ref), [`confint`](@ref)
"""
function ApproximateSignedRankTest(x::Vector, W::Float64, ranks::Vector{T}, signs::BitArray{1}, tie_adjustment::Float64, n::Int, median::Float64) where T<:Real
    nz = length(ranks) # num non-zeros
    mu = W - nz * (nz + 1)/4
    std = sqrt(nz * (nz + 1) * (2 * nz + 1) / 24 - tie_adjustment / 48)
    ApproximateSignedRankTest(x, W, ranks, signs, tie_adjustment, n, median, mu, std)
end
ApproximateSignedRankTest(x::AbstractVector{T}) where {T<:Real} =
    ApproximateSignedRankTest(x, signedrankstats(x)...)
ApproximateSignedRankTest(x::AbstractVector{S}, y::AbstractVector{T}) where {S<:Real,T<:Real} =
    ApproximateSignedRankTest(x - y)

testname(::ApproximateSignedRankTest) = "Approximate Wilcoxon signed rank test"
population_param_of_interest(x::ApproximateSignedRankTest) = ("Location parameter (pseudomedian)", 0, x.median) # parameter of interest: name, value under h0, point estimate
default_tail(test::ApproximateSignedRankTest) = :both

function show_params(io::IO, x::ApproximateSignedRankTest, ident)
    println(io, ident, "number of observations:      ", x.n)
    println(io, ident, "Wilcoxon rank-sum statistic: ", x.W)
    print(io, ident, "rank sums:                   ")
    show(io, [sum(x.ranks[x.signs]), sum(x.ranks[map(!, x.signs)])])
    println(io)
    println(io, ident, "adjustment for ties:         ", x.tie_adjustment)
    println(io, ident, "normal approximation (μ, σ): ", (x.mu, x.sigma))
end

function pvalue(x::ApproximateSignedRankTest; tail=:both)
    check_tail(tail)

    if x.mu == x.sigma == 0
        1.0
    elseif tail == :both
        2 * ccdf(Normal(), abs(x.mu - 0.5 * sign(x.mu))/x.sigma)
    elseif tail == :left
        cdf(Normal(), (x.mu + 0.5)/x.sigma)
    else # tail == :right
        ccdf(Normal(), (x.mu - 0.5)/x.sigma)
    end
end

StatsBase.confint(x::ApproximateSignedRankTest; level::Real=0.95, tail=:both) = calculate_ci(x.vals, level, tail=tail)

# implementation method inspired by these notes: http://www.stat.umn.edu/geyer/old03/5102/notes/rank.pdf
function calculate_ci(x::AbstractVector, level::Real=0.95; tail=:both)
    check_level(level)
    check_tail(tail)

    if tail == :both
        c = level
    else
        c = 1 - 2 * (1-level)
    end
    n = length(x)
    m = div(n * (n + 1), 2)
    k_range = 1:div(m, 2)
    l = [1 - 2 * psignrank(i, n, true) for i in k_range]
    k = argmin(abs.(l .- c))
    vals = Float64[]
    enumerated = enumerate(x)
    for (outer_index, outer_value) in enumerated
        for (inner_index, inner_value) in enumerated
            if outer_index > inner_index
                continue
            end
            push!(vals, (inner_value + outer_value) / 2)
        end
    end
    sort!(vals)
    left = vals[k + 1]
    right = vals[m - k]
    if tail == :both
        return (left, right)
    elseif tail == :left
        return (left, Inf)
    else # tail == :right
        return (-Inf, right)
    end
end


text:
using Oceananigans.Fields: FunctionField

"""
    QuasiAdamsBashforth2TimeStepper{T, TG} <: AbstractTimeStepper

Holds tendency fields and the parameter `χ` for a modified second-order
Adams-Bashforth timestepping method.
"""
struct QuasiAdamsBashforth2TimeStepper{T, TG} <: AbstractTimeStepper
     χ :: T
    Gⁿ :: TG
    G⁻ :: TG
end

"""
    QuasiAdamsBashforth2TimeStepper(arch, grid, tracers, χ=0.1;
                                    Gⁿ = TendencyFields(arch, grid, tracers),
                                    G⁻ = TendencyFields(arch, grid, tracers))

Return an QuasiAdamsBashforth2TimeStepper object with tendency fields on `arch` and
`grid` with AB2 parameter `χ`. The tendency fields can be specified via optional
kwargs.
"""
function QuasiAdamsBashforth2TimeStepper(arch, grid, tracers, χ=0.1;
                                         Gⁿ = TendencyFields(arch, grid, tracers),
                                         G⁻ = TendencyFields(arch, grid, tracers))

    return QuasiAdamsBashforth2TimeStepper{eltype(grid), typeof(Gⁿ)}(χ, Gⁿ, G⁻)
end

#####
##### Time steppping
#####

"""
    time_step!(model::AbstractModel{<:QuasiAdamsBashforth2TimeStepper}, Δt; euler=false)

Step forward `model` one time step `Δt` with a 2nd-order Adams-Bashforth method and
pressure-correction substep. Setting `euler=true` will take a forward Euler time step.
"""
function time_step!(model::AbstractModel{<:QuasiAdamsBashforth2TimeStepper}, Δt; euler=false)
    Δt == 0 && @warn "Δt == 0 may cause model blowup!"

    χ = ifelse(euler, convert(eltype(model.grid), -0.5), model.timestepper.χ)

    # Be paranoid and update state at iteration 0, in case run! is not used:
    model.clock.iteration == 0 && update_state!(model)

    calculate_tendencies!(model)

    ab2_step!(model, Δt, χ) # full step for tracers, fractional step for velocities.

    calculate_pressure_correction!(model, Δt)
    pressure_correct_velocities!(model, Δt)

    tick!(model.clock, Δt)
    update_state!(model)
    store_tendencies!(model)
    update_particle_properties!(model, Δt)

    return nothing
end

#####
##### Time stepping in each step
#####

function ab2_step!(model, Δt, χ)

    workgroup, worksize = work_layout(model.grid, :xyz)

    barrier = Event(device(model.architecture))

    step_field_kernel! = ab2_step_field!(device(model.architecture), workgroup, worksize)

    model_fields = fields(model)

    events = []

    for (i, field) in enumerate(model_fields)

        field_event = step_field_kernel!(field, Δt, χ,
                                         model.timestepper.Gⁿ[i],
                                         model.timestepper.G⁻[i],
                                         dependencies=Event(device(model.architecture)))

        push!(events, field_event)
    end

    wait(device(model.architecture), MultiEvent(Tuple(events)))

    return nothing
end

"""
Time step via

    `U^{n+1} = U^n + Δt ( (3/2 + χ) * G^{n} - (1/2 + χ) G^{n-1} )`

"""

@kernel function ab2_step_field!(U, Δt, χ::FT, Gⁿ, G⁻) where FT
    i, j, k = @index(Global, NTuple)

    @inbounds begin
        U[i, j, k] += Δt * (  (FT(1.5) + χ) * Gⁿ[i, j, k] - (FT(0.5) + χ) * G⁻[i, j, k] )

    end
end

@kernel ab2_step_field!(ϕ::FunctionField, args...) = nothing


text:
using LightGraphs
using BenchmarkTools
@show Threads.nthreads()

@benchgroup "parallel" begin
    @benchgroup "egonet" begin
        function vertex_function(g::Graph, i::Int)
            a = 0
            for u in neighbors(g, i)
                a += degree(g, u)
            end
            return a
        end

        function twohop(g::Graph, i::Int)
            a = 0
            for u in neighbors(g, i)
                for v in neighbors(g, u)
                    a += degree(g, v)
                end
            end
            return a
        end

        function mapvertices(f, g::Graph)
            n = nv(g)
            a = zeros(Int, n)
            Threads.@threads for i in 1:n
                a[i] = f(g, i)
            end
            return a
        end

        function mapvertices_single(f, g)
            n = nv(g)
            a = zeros(Int, n)
            for i in 1:n
                a[i] = f(g, i)
            end
            return a
        end

        function comparison(f, g)
            println("Mulithreaded on $(Threads.nthreads())")
            b1 =  @benchmark mapvertices($f, $g)
            println(b1)

            println("singlethreaded")
            b2 = @benchmark mapvertices_single($f, $g)
            println(b2)
            println("done")
        end

        nv_ = 10000
        g = Graph(nv_, 64 * nv_)
        f = vertex_function
        println(g)

        comparison(vertex_function, g)
        comparison(twohop, g)
    end
end


text:
## Forest Fire Simulation
using Random
using PyPlot
pygui(true)

using Distributed

function showforest(forest)
    PyPlot.imshow(forest)
end

# 0 no tree
# 1 tree
# 2 burning
# 3 burnt

# test if something is burning
function stillburning(forest)
    for position in forest
        if position == 2
            return true
        end
    end
    return false
end

function simulation_step(forest, forest_size)
    temp_forest = copy(forest)
    for x = 1:forest_size, y = 1:forest_size
        if forest[x, y] == 2
            temp_forest[x, y] = 3
            for neighbor in [(-1, 0), (0, -1), (1, 0), (0, 1)]
                position = (x, y) .+ neighbor
                if all(position .> (0, 0)) &&
                   all(position .<= (forest_size, forest_size))
                    xp, yp = position
                    if forest[xp, yp] == 1
                        # put on fire
                        temp_forest[xp, yp] = 2
                    end

                end

            end

        end
    end
    return temp_forest
end

function report_burnt(forest)
    count_burnt = 0
    count_tree = 0
    for position in forest
        if position == 3
            count_burnt = count_burnt + 1
        end
        if position == 1
            count_tree = count_tree + 1
        end
    end

    return count_burnt / (count_burnt + count_tree)
end

function simulation(density = 0.57, forest_size = 251)
    forest = [zero(Int8) for x = 1:forest_size, y = 1:forest_size]

    for x = 1:forest_size, y = 1:forest_size
        if rand() < density
            forest[x, y] = 1
        end
    end

    for y = 1:forest_size
        forest[y, 1] = 2
    end

    while stillburning(forest)
        forest = simulation_step(forest, forest_size)
    end
    print(Threads.threadid())
    print(" ")
    println(report_burnt(forest))

    forest
end

Random.seed!(0)

@time for i in 1:100
    results = simulation(0.57, 500)
    #showforest(results)
end

Threads.nthreads()


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

# %%
using SparseArrays

A = zeros( Integer, 2400, 2400 )
for i in 1:2400
    if rand(1:4) == 1
        for j in i+1:2400
            if rand(1:15) == 1
                A[i,j] = 1
                A[j,i] = 1
            end
        end
    end
end

sum( A)
A = sparse( A )
B = A^2

println( "some stats A: size ", size(A), " min value ", minimum( A ), " max value ", maximum( A ) )
println( "some stats B: size ", size(B), " min value ", minimum( B ), " max value ", maximum( B ) )

# %%
A

# %%
B

# %%
a = spzeros(Integer, 5, 5)
a .= 1
a

# %%
b = a^2

# %%

max_stars_repo_path: 0011/eltype of square of sparse array with Integer eltype is Any.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
"""
    computeslip(mom, δ, force)

Calculates the force due to slippage on the substrate.

Since the model only alows for forces in `x` and `y` direction but not in `z` it is necessary to regularize the velocity at `` h \\rightarrow 0 ``.
The so called no-slip boundary condition can be in theory implemented using `δ` = 0.
However as soon as the film dewets and `height` becomes small, the solution diverges.
Choosing `δ` = 1 satisfies the weak slip condition, while choosing it much larger than 1 yields the strong slip solution.

# Math
One and two dimensional approaches are essentially equal.
Behind the force is the idea of a parabolic shape of the velocity along `z`.
With the boundary conditions

`` \\mathbf{u}|_{h=H} = \\mathbf{u}  ``

`` \\mathbf{u}|_{h=0} = 0 ``

The parameter `δ` allows to push the zero velocity into the substrate, such below te computable region.
Solving this equation with some more assumptions yields

`` \\mathbf{F}_{slip} = \\mu\\alpha_{\\delta}(h) \\mathbf{u} = \\frac{6\\mu h \\mathbf{u}}{2 h^2 + 6\\delta h + 3\\delta^2}. ``

# Example

# References
Here I shamelessly cite the paper from within our group
- [Lattice Boltzmann method for thin-liquid-film hydrodynamics](https://journals.aps.org/pre/abstract/10.1103/PhysRevE.100.033313)

See also: [`Inputconstants`](@ref), [`calculatemoments`](@ref)

"""
function computeslip(mom::JuSwalbe.Macroquant{Vector{T}, Vector{T}}, forces::JuSwalbe.Forces{Vector{T}}, input::JuSwalbe.Inputconstants) where {T<:Number}
    # Measure the length and allocate a dummy array
    len = length(mom.height)
    slippage = zeros(T, len)
    num = zeros(T, len)
    denom = zeros(T, len)
    # Slippage can be calculated assuming a parabolic velocity profile in z
    num .= 6 .* mom.height .* mom.velocity
    denom .= 2 .* mom.height .+ 6 * mom.height * input.δ .+ 3 * input.δ^2
    slippage .= input.μ * (num ./ denom)
    # Write the result into the forcing array
    forces.slip = - slippage
    # Write it out as well
    return slippage
end

function computeslip(mom::JuSwalbe.Macroquant{Matrix{T}, JuSwalbe.Twovector{Matrix{T}}}, forces::JuSwalbe.Forces{JuSwalbe.Twovector{Matrix{T}}}, input::JuSwalbe.Inputconstants) where {T<:Number}
    # Measure the length and allocate a dummy array
    width, thick = size(mom.height)
    δ, μ = T(input.δ), T(input.μ)
    slippagex = zeros(T, (width, thick))
    slippagey = zeros(T, (width, thick))
    numx = zeros(T, (width, thick))
    numy = zeros(T, (width, thick))
    denom = zeros(T, (width, thick))
    # Slippage can be calculated assuming a parabolic velocity profile in z
    numx .= 6 .* mom.height .* mom.velocity.x
    numy .= 6 .* mom.height .* mom.velocity.y
    denom .= 2 .* mom.height .+ 6 * mom.height * δ .+ 3 * δ^2
    slippagex .= μ * ( numx ./ denom )
    slippagey .= μ * ( numy ./ denom )
    # Write the result into the forcing array
    slippage = JuSwalbe.Twovector(x= .- slippagex, y= .- slippagey)
    forces.slip = slippage
    # Write it out as well
    return slippage
end

"""
    computethermalcapillarywaves(mom, force, input)

Computes the forcing arising due to the inclusion of thermal capillary waves that undulate the free surface.

Adds a normal distributed term weighted with the thermal energy kbt to the forcing struct.
This addition makes it possible to simulate not only the thin film equation but also the stochastic thin film equation.
The output of the function should look like (of course dependent on input) this:
![Thermo capillary forcing distribution](https://github.com/Zitzeronion/JuSwalbe.jl/tree/thermalforce-feature/images/thermal_forcing.pdf)

# Math
The stochastic thin film equation (STF) much like the thin film equation is derived from the Landau Lifshitz Navier Stokes equation (LLNS).
Assuming the fluctuations are small and not leading order one can perform an integration of the stochastic stresses, which was done by Grün et al.,

`` \\frac{\\partial h}{\\partial t} = \\nabla\\cdot\\Big[\\frac{h^3}{3\\mu}\\nabla(\\Pi(h) - \\gamma \\Delta h) + \\int_0^h (h-y)\\mathcal{S}_{z||}(y)dy]. ``

Taking into account several assumption it can be shown that the above equation is equal to

`` \\frac{\\partial h}{\\partial t} = \\nabla\\cdot\\Big[\\frac{h^3}{3\\mu}\\nabla(\\Pi(h) - \\gamma \\Delta h) + \\sqrt{\\frac{2k_BT}{3\\mu}h^3}\\mathcal{N}(t)]. ``

Such the integral of `` \\mathcal{S} `` can be expressed with a single multipicative noise term `` \\mathcal{N} `` which is gaussian distributed with zero mean and a variance of one.

This force is adding exactly the term with `` \\mathcal{N} `` to the forcing struct, however it is a little modified due to slippage.

# Example
```jldoctest
julia> using JuSwalbe

julia> mom = simplemoment1d(1000)
JuSwalbe.Macroquant{Array{Float64,1},Array{Float64,1}}
  height: Array{Float64}((1000,)) [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
  velocity: Array{Float64}((1000,)) [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1  …  0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1]
  pressure: Array{Float64}((1000,)) [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
  energy: Array{Float64}((1000,)) [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

julia> force = JuSwalbe.Forces(slip=fill(0.1, 1000), thermal=zeros(1000), h∇p=zeros(1000), bathymetry=zeros(1000))
JuSwalbe.Forces{Array{Float64,1}}
  slip: Array{Float64}((1000,)) [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1  …  0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1]
  h∇p: Array{Float64}((1000,)) [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
  bathymetry: Array{Float64}((1000,)) [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
  thermal: Array{Float64}((1000,)) [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

julia> input = JuSwalbe.Inputconstants()
JuSwalbe.Inputconstants
  lx: Int64 512
  ly: Int64 512
  maxruntime: Int64 100000
  dumping: Int64 1000
  τ: Float64 1.0
  gravity: Float64 0.0
  γ: Float64 0.01
  δ: Float64 1.0
  μ: Float64 0.16666666666666666
  kbt: Float64 0.0

julia> thermal = computethermalcapillarywaves(mom, force, input); # No k_BT -> no force!

julia> input2 = JuSwalbe.Inputconstants(kbt = 0.001)
JuSwalbe.Inputconstants
  lx: Int64 512
  ly: Int64 512
  maxruntime: Int64 100000
  dumping: Int64 1000
  τ: Float64 1.0
  gravity: Float64 0.0
  γ: Float64 0.01
  δ: Float64 1.0
  μ: Float64 0.16666666666666666
  kbt: Float64 0.001

julia> thermal = computethermalcapillarywaves(mom, force, input2); # Now this force is not just zeros
```

# Code
This function is used to compute a force such one needs it for the collision operator, and if Luo forcing is used for the calculation of the moments
See also: [`collisionBGK`](@ref), [`calculatemoments`](@ref)

# References
The theoretical foundation for this was laid by G. Grün and H. Stone,
- [Thin-Film Flow Influenced by Thermal Noise](https://link.springer.com/article/10.1007/s10955-006-9028-8)
- [Spreading of Viscous Fluid Drops on a Solid Substrate Assisted by Thermal Fluctuations](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.95.244505)

Some more recent references are
- [Molecular simulation of thin liquid films: Thermal fluctuations and instability](https://journals.aps.org/pre/abstract/10.1103/PhysRevE.100.023108)
- [Fully nonlinear dynamics of stochastic thin-film dewetting](https://journals.aps.org/pre/abstract/10.1103/PhysRevE.92.061002)

Basically it reduces to drastically simplify the LLNS under various assumptions and show that the result is more or less equal to something simple.
"""
function computethermalcapillarywaves(mom::JuSwalbe.Macroquant{Vector{T},Vector{T}}, forces::JuSwalbe.Forces{Vector{T}}, input::Inputconstants) where {T<:Number}
  if input.kbt == 0.0
    return
  else
    len = length(mom.height)
    thermocap = zeros(T, len)
    slip = deepcopy(forces.slip)
    μ = T(input.μ)
    kbt = T(input.kbt)
    # Generate a Gaussian distribution with zero mean and variance of one
    gaussian = Normal()
    # Fill an array with random numbers distributed according to a Normal
    gaussianvec = rand(gaussian, len)
    # Compute the forces due to thermal capillary waves
    thermocap .= sqrt.(2 * kbt * μ * slip) .* gaussianvec
    forces.thermal = thermocap
    return thermocap
  end
end

function computethermalcapillarywaves(mom::JuSwalbe.Macroquant{Matrix{T},JuSwalbe.Twovector{Matrix{T}}}, forces::JuSwalbe.Forces{JuSwalbe.Twovector{Matrix{T}}}, input::Inputconstants) where {T<:Number}
  if input.kbt == 0.0
    return
  else
    width, thick = size(mom.height)
    thermocapx = zeros(T, (width, thick))
    thermocapy = zeros(T, (width, thick))
    slip = deepcopy(forces.slip)
    μ = T(input.μ)
    kbt = T(input.kbt)
    # Generate a Gaussian distribution with zero mean and a variance of one
    gaussian = Normal()
    # Fill an array with random numbers distributed according to a gaussian
    gaussianvecx = rand(gaussian, (width, thick))
    gaussianvecy = rand(gaussian, (width, thick))
    # Compute the forces due to thermal capillary waves
    thermocapx .= sqrt.(2 * kbt * μ * slip.x) .* gaussianvecx
    thermocapy .= sqrt.(2 * kbt * μ * slip.y) .* gaussianvecy

    forces.thermal = JuSwalbe.Twovector(x=thermocapx, y=thermocapy)
    return thermocapx, thermocapy
  end
end

max_stars_repo_name: Zitzeronion/JuSwalbe

max_issues_repo_name: Zitzeronion/JuSwalbe


text:
function lastMapDir()
    celesteDir = config["celeste_dir"]
    mapsDir = joinpath(celesteDir, "Mods")
    lastDir = dirname(get(persistence, "files_lastfile", mapsDir))
    targetDir = ispath(lastDir) ? lastDir : mapsDir
    targetDir = ispath(targetDir) ? targetDir : pwd()

    return targetDir
end

function showFileOpenDialog(leaf::Ahorn.MenuItemsTypes=MenuItem(), folder::String=lastMapDir(), overrideFilename::Union{String, Nothing}=nothing)
    @Ahorn.catchall begin
        filename = openDialog("Select map binary", window, ["*.bin"], folder=folder)

        # Did we actually select a file?
        if filename != ""
            dialogText = "You might have unsaved changes in your map.\nPlease confirm that you want to load '$filename'."
            res = loadedState.lastSavedHash == hash(loadedState.side) || ask_dialog(dialogText, window)

            if res
                # Clear render cache from last map
                if loadedState.map !== nothing && get(config, "clear_render_cache_on_map_load", true)
                    deleteDrawableRoomCache(loadedState.map)
                end

                # Allow backups to "fake" where a file is actually stored
                loadedState.filename = something(overrideFilename, filename)
                loadedState.side = loadSide(filename)
                loadedState.map = loadedState.side.map
                loadedState.lastSavedHash = hash(loadedState.side)
                loadedState.roomName = ""
                loadedState.room = nothing

                persistence["files_lastroom"] = loadedState.roomName
                persistence["files_lastfile"] = loadedState.filename

                packageName = loadedState.map.package

                EntityIds.updateValidIds(loadedState.map)

                loadXMLMeta()
                Ahorn.loadChangedExternalSprites!()

                updateTreeView!(roomList, getTreeData(loadedState.map), 1)
                handleRoomChanged(loadedState.map, loadedState.room)

                GAccessor.title(window, "$baseTitle - $(Maple.getSideName(loadedState.side))")

                draw(canvas)

                return true
            end
        end
    end

    return false
end

# saveDialog warns on overwrite, but only returns the filename
function showFileSaveDialog(leaf::Ahorn.MenuItemsTypes=MenuItem())
    targetDir = lastMapDir()
    filename = saveDialog("Save as", window, ["*.bin"], folder=targetDir)

    if filename != ""
        loadedState.filename = filename
        persistence["files_lastfile"] = loadedState.filename

        saveFile(loadedState.side, filename)

        return true
    end

    return false
end

function menuFileSave(leaf::Ahorn.MenuItemsTypes=MenuItem())
    if loadedState.filename === nothing || loadedState.filename == ""
        showFileSaveDialog(leaf)

    else
        saveFile(loadedState.side, loadedState.filename)
    end
end

function saveFile(side::Side, filename::String)
    try
        PreSaveSanitizers.sortRoomNames(side, config)
        PreSaveSanitizers.warnMissingDecals(side, config)
        PreSaveSanitizers.warnDuplicateIds(side, config)

        loadedState.lastSavedHash = hash(side)

        fn = hasExt(filename, ".bin") ? filename : filename * ".bin"
        encodeSide(side, fn)

    catch e
        println(Base.stderr, e)
        info_dialog("Failed to save map.", Ahorn.window)
    end
end

max_stars_repo_name: CelestialCartographers/Ahorn

max_issues_repo_name: CelestialCartographers/Ahorn

max_forks_repo_name: CelestialCartographers/Ahorn