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

struct Set{T} <: AbstractSet{T}
    dict::Dict{T,Nothing}

    Set{T}() where {T} = new(Dict{T,Nothing}())
    Set{T}(s::Set{T}) where {T} = new(Dict{T,Nothing}(s.dict))
end

Set{T}(itr) where {T} = union!(Set{T}(), itr)
Set() = Set{Any}()

"""
    Set([itr])

Construct a [`Set`](@ref) of the values generated by the given iterable object, or an
empty set. Should be used instead of [`BitSet`](@ref) for sparse integer sets, or
for sets of arbitrary objects.
"""
Set(itr) = _Set(itr, IteratorEltype(itr))

_Set(itr, ::HasEltype) = Set{eltype(itr)}(itr)

function _Set(itr, ::EltypeUnknown)
    T = @default_eltype(itr)
    (isconcretetype(T) || T === Union{}) || return grow_to!(Set{T}(), itr)
    return Set{T}(itr)
end

empty(s::AbstractSet{T}, ::Type{U}=T) where {T,U} = Set{U}()

# return an empty set with eltype T, which is mutable (can be grown)
# by default, a Set is returned
emptymutable(s::AbstractSet{T}, ::Type{U}=T) where {T,U} = Set{U}()

_similar_for(c::AbstractSet, ::Type{T}, itr, isz) where {T} = empty(c, T)

function show(io::IO, s::Set)
    print(io, "Set(")
    show_vector(io, s)
    print(io, ')')
end

isempty(s::Set) = isempty(s.dict)
length(s::Set)  = length(s.dict)
in(x, s::Set) = haskey(s.dict, x)
push!(s::Set, x) = (s.dict[x] = nothing; s)
pop!(s::Set, x) = (pop!(s.dict, x); x)
pop!(s::Set, x, default) = (x in s ? pop!(s, x) : default)

function pop!(s::Set)
    isempty(s) && throw(ArgumentError("set must be non-empty"))
    return pop!(s.dict)[1]
end

delete!(s::Set, x) = (delete!(s.dict, x); s)

copy(s::Set) = copymutable(s)

# Set is the default mutable fall-back
copymutable(s::AbstractSet{T}) where {T} = Set{T}(s)

sizehint!(s::Set, newsz) = (sizehint!(s.dict, newsz); s)
empty!(s::Set) = (empty!(s.dict); s)
rehash!(s::Set) = (rehash!(s.dict); s)

iterate(s::Set, i...)       = iterate(KeySet(s.dict), i...)

# In case the size(s) is smaller than size(t) its more efficient to iterate through
# elements of s instead and only delete the ones also contained in t.
# The threshold for this decision boils down to a tradeoff between
# size(s) * cost(in() + delete!()) ≶ size(t) * cost(delete!())
# Empirical observations on Ints point towards a threshold of 0.8.
# To be on the safe side (e.g. cost(in) >>> cost(delete!) ) a
# conservative threshold of 0.5 was chosen.
function setdiff!(s::Set, t::Set)
    if 2 * length(s) < length(t)
        for x in s
            x in t && delete!(s, x)
        end
    else
        for x in t
            delete!(s, x)
        end
    end
    return s
end

"""
    unique(itr)

Return an array containing only the unique elements of collection `itr`,
as determined by [`isequal`](@ref), in the order that the first of each
set of equivalent elements originally appears. The element type of the
input is preserved.

# Examples
```jldoctest
julia> unique([1, 2, 6, 2])
3-element Array{Int64,1}:
 1
 2
 6

julia> unique(Real[1, 1.0, 2])
2-element Array{Real,1}:
 1
 2
```
"""
function unique(itr)
    T = @default_eltype(itr)
    out = Vector{T}()
    seen = Set{T}()
    y = iterate(itr)
    y === nothing && return out
    x, i = y
    if !isconcretetype(T) && IteratorEltype(itr) == EltypeUnknown()
        S = typeof(x)
        return _unique_from(itr, S[x], Set{S}((x,)), i)
    end
    push!(seen, x)
    push!(out, x)
    return unique_from(itr, out, seen, i)
end

_unique_from(itr, out, seen, i) = unique_from(itr, out, seen, i)
@inline function unique_from(itr, out::Vector{T}, seen, i) where T
    while true
        y = iterate(itr, i)
        y === nothing && break
        x, i = y
        S = typeof(x)
        if !(S === T || S <: T)
            R = promote_typejoin(S, T)
            seenR = convert(Set{R}, seen)
            outR = convert(Vector{R}, out)
            if !in(x, seenR)
                push!(seenR, x)
                push!(outR, x)
            end
            return _unique_from(itr, outR, seenR, i)
        end
        if !in(x, seen)
            push!(seen, x)
            push!(out, x)
        end
    end
    return out
end

"""
    unique(f, itr)

Returns an array containing one value from `itr` for each unique value produced by `f`
applied to elements of `itr`.

# Examples
```jldoctest
julia> unique(x -> x^2, [1, -1, 3, -3, 4])
3-element Array{Int64,1}:
 1
 3
 4
```
"""
function unique(f, C)
    out = Vector{eltype(C)}()

    s = iterate(C)
    if s === nothing
        return out
    end
    (x, i) = s
    y = f(x)
    seen = Set{typeof(y)}()
    push!(seen, y)
    push!(out, x)

    return _unique!(f, out, C, seen, i)
end

function _unique!(f, out::AbstractVector, C, seen::Set, i)
    s = iterate(C, i)
    while s !== nothing
        (x, i) = s
        y = f(x)
        if y ∉ seen
            push!(out, x)
            if y isa eltype(seen)
                push!(seen, y)
            else
                seen2 = convert(Set{promote_typejoin(eltype(seen), typeof(y))}, seen)
                push!(seen2, y)
                return _unique!(f, out, C, seen2, i)
            end
        end
        s = iterate(C, i)
    end

    return out
end

"""
    unique!(f, A::AbstractVector)

Selects one value from `A` for each unique value produced by `f` applied to
elements of `A` , then return the modified A.

!!! compat "Julia 1.1"
    This method is available as of Julia 1.1.

# Examples
```jldoctest
julia> unique!(x -> x^2, [1, -1, 3, -3, 4])
3-element Array{Int64,1}:
 1
 3
 4

julia> unique!(n -> n%3, [5, 1, 8, 9, 3, 4, 10, 7, 2, 6])
3-element Array{Int64,1}:
 5
 1
 9

julia> unique!(iseven, [2, 3, 5, 7, 9])
2-element Array{Int64,1}:
 2
 3
```
"""
function unique!(f, A::AbstractVector)
    if length(A) <= 1
        return A
    end

    i = firstindex(A)
    x = @inbounds A[i]
    y = f(x)
    seen = Set{typeof(y)}()
    push!(seen, y)
    return _unique!(f, A, seen, i, i+1)
end

function _unique!(f, A::AbstractVector, seen::Set, current::Integer, i::Integer)
    while i <= lastindex(A)
        x = @inbounds A[i]
        y = f(x)
        if y ∉ seen
            current += 1
            @inbounds A[current] = x
            if y isa eltype(seen)
                push!(seen, y)
            else
                seen2 = convert(Set{promote_typejoin(eltype(seen), typeof(y))}, seen)
                push!(seen2, y)
                return _unique!(f, A, seen2, current, i+1)
            end
        end
        i += 1
    end
    return resize!(A, current - firstindex(A) + 1)
end

# If A is not grouped, then we will need to keep track of all of the elements that we have
# seen so far.
_unique!(A::AbstractVector) = unique!(identity, A::AbstractVector)

# If A is grouped, so that each unique element is in a contiguous group, then we only
# need to keep track of one element at a time. We replace the elements of A with the
# unique elements that we see in the order that we see them. Once we have iterated
# through A, we resize A based on the number of unique elements that we see.
function _groupedunique!(A::AbstractVector)
    isempty(A) && return A
    idxs = eachindex(A)
    y = first(A)
    # We always keep the first element
    it = iterate(idxs, iterate(idxs)[2])
    count = 1
    for x in Iterators.drop(A, 1)
        if !isequal(x, y)
            y = A[it[1]] = x
            count += 1
            it = iterate(idxs, it[2])
        end
    end
    resize!(A, count)
end

"""
    unique!(A::AbstractVector)

Remove duplicate items as determined by [`isequal`](@ref), then return the modified `A`.
`unique!` will return the elements of `A` in the order that they occur. If you do not care
about the order of the returned data, then calling `(sort!(A); unique!(A))` will be much
more efficient as long as the elements of `A` can be sorted.

# Examples
```jldoctest
julia> unique!([1, 1, 1])
1-element Array{Int64,1}:
 1

julia> A = [7, 3, 2, 3, 7, 5];

julia> unique!(A)
4-element Array{Int64,1}:
 7
 3
 2
 5

julia> B = [7, 6, 42, 6, 7, 42];

julia> sort!(B);  # unique! is able to process sorted data much more efficiently.

julia> unique!(B)
3-element Array{Int64,1}:
  6
  7
 42
```
"""
function unique!(A::Union{AbstractVector{<:Real}, AbstractVector{<:AbstractString},
                          AbstractVector{<:Symbol}})
    if isempty(A)
        return A
    elseif issorted(A) || issorted(A, rev=true)
        return _groupedunique!(A)
    else
        return _unique!(A)
    end
end
# issorted fails for some element types, so the method above has to be restricted to
# elements with isless/< defined.
function unique!(A)
    if isempty(A)
        return A
    else
        return _unique!(A)
    end
end

"""
    allunique(itr) -> Bool

Return `true` if all values from `itr` are distinct when compared with [`isequal`](@ref).

# Examples
```jldoctest
julia> a = [1; 2; 3]
3-element Array{Int64,1}:
 1
 2
 3

julia> allunique([a, a])
false
```
"""
function allunique(C)
    seen = Set{eltype(C)}()
    for x in C
        if in(x, seen)
            return false
        else
            push!(seen, x)
        end
    end
    true
end

allunique(::Union{AbstractSet,AbstractDict}) = true

allunique(r::AbstractRange{T}) where {T} = (step(r) != zero(T)) || (length(r) <= 1)
allunique(r::StepRange{T,S}) where {T,S} = (step(r) != zero(S)) || (length(r) <= 1)
allunique(r::StepRangeLen{T,R,S}) where {T,R,S} = (step(r) != zero(S)) || (length(r) <= 1)

filter!(f, s::Set) = unsafe_filter!(f, s)

const hashs_seed = UInt === UInt64 ? 0x852ada37cfe8e0ce : 0xcfe8e0ce
function hash(s::AbstractSet, h::UInt)
    hv = hashs_seed
    for x in s
        hv ⊻= hash(x)
    end
    hash(hv, h)
end

convert(::Type{T}, s::T) where {T<:AbstractSet} = s
convert(::Type{T}, s::AbstractSet) where {T<:AbstractSet} = T(s)

## replace/replace! ##

function check_count(count::Integer)
    count < 0 && throw(DomainError(count, "`count` must not be negative (got $count)"))
    return min(count, typemax(Int)) % Int
end

# TODO: use copy!, which is currently unavailable from here since it is defined in Future
_copy_oftype(x, ::Type{T}) where {T} = copyto!(similar(x, T), x)
# TODO: use similar() once deprecation is removed and it preserves keys
_copy_oftype(x::AbstractDict, ::Type{T}) where {T} = merge!(empty(x, T), x)
_copy_oftype(x::AbstractSet, ::Type{T}) where {T} = union!(empty(x, T), x)

_copy_oftype(x::AbstractArray{T}, ::Type{T}) where {T} = copy(x)
_copy_oftype(x::AbstractDict{K,V}, ::Type{Pair{K,V}}) where {K,V} = copy(x)
_copy_oftype(x::AbstractSet{T}, ::Type{T}) where {T} = copy(x)

_similar_or_copy(x::Any) = similar(x)
_similar_or_copy(x::Any, ::Type{T}) where {T} = similar(x, T)
# Make a copy on construction since it is faster than inserting elements separately
_similar_or_copy(x::Union{AbstractDict,AbstractSet}) = copy(x)
_similar_or_copy(x::Union{AbstractDict,AbstractSet}, ::Type{T}) where {T} = _copy_oftype(x, T)

# to make replace/replace! work for a new container type Cont, only
# _replace!(new::Callable, res::Cont, A::Cont, count::Int)
# has to be implemented

"""
    replace!(A, old_new::Pair...; [count::Integer])

For each pair `old=>new` in `old_new`, replace all occurrences
of `old` in collection `A` by `new`.
Equality is determined using [`isequal`](@ref).
If `count` is specified, then replace at most `count` occurrences in total.
See also [`replace`](@ref replace(A, old_new::Pair...)).

# Examples
```jldoctest
julia> replace!([1, 2, 1, 3], 1=>0, 2=>4, count=2)
4-element Array{Int64,1}:
 0
 4
 1
 3

julia> replace!(Set([1, 2, 3]), 1=>0)
Set{Int64} with 3 elements:
  0
  2
  3
```
"""
replace!(A, old_new::Pair...; count::Integer=typemax(Int)) =
    replace_pairs!(A, A, check_count(count), old_new)

function replace_pairs!(res, A, count::Int, old_new::Tuple{Vararg{Pair}})
    @inline function new(x)
        for o_n in old_new
            isequal(first(o_n), x) && return last(o_n)
        end
        return x # no replace
    end
    _replace!(new, res, A, count)
end

"""
    replace!(new::Function, A; [count::Integer])

Replace each element `x` in collection `A` by `new(x)`.
If `count` is specified, then replace at most `count` values in total
(replacements being defined as `new(x) !== x`).

# Examples
```jldoctest
julia> replace!(x -> isodd(x) ? 2x : x, [1, 2, 3, 4])
4-element Array{Int64,1}:
 2
 2
 6
 4

julia> replace!(Dict(1=>2, 3=>4)) do kv
           first(kv) < 3 ? first(kv)=>3 : kv
       end
Dict{Int64,Int64} with 2 entries:
  3 => 4
  1 => 3

julia> replace!(x->2x, Set([3, 6]))
Set{Int64} with 2 elements:
  6
  12
```
"""
replace!(new::Callable, A; count::Integer=typemax(Int)) =
    _replace!(new, A, A, check_count(count))

"""
    replace(A, old_new::Pair...; [count::Integer])

Return a copy of collection `A` where, for each pair `old=>new` in `old_new`,
all occurrences of `old` are replaced by `new`.
Equality is determined using [`isequal`](@ref).
If `count` is specified, then replace at most `count` occurrences in total.

The element type of the result is chosen using promotion (see [`promote_type`](@ref))
based on the element type of `A` and on the types of the `new` values in pairs.
If `count` is omitted and the element type of `A` is a `Union`, the element type
of the result will not include singleton types which are replaced with values of
a different type: for example, `Union{T,Missing}` will become `T` if `missing` is
replaced.

See also [`replace!`](@ref).

# Examples
```jldoctest
julia> replace([1, 2, 1, 3], 1=>0, 2=>4, count=2)
4-element Array{Int64,1}:
 0
 4
 1
 3

julia> replace([1, missing], missing=>0)
2-element Array{Int64,1}:
 1
 0
```
"""
function replace(A, old_new::Pair...; count::Union{Integer,Nothing}=nothing)
    V = promote_valuetype(old_new...)
    if count isa Nothing
        T = promote_type(subtract_singletontype(eltype(A), old_new...), V)
        replace_pairs!(_similar_or_copy(A, T), A, typemax(Int), old_new)
    else
        U = promote_type(eltype(A), V)
        replace_pairs!(_similar_or_copy(A, U), A, check_count(count), old_new)
    end
end

promote_valuetype(x::Pair{K, V}) where {K, V} = V
promote_valuetype(x::Pair{K, V}, y::Pair...) where {K, V} =
    promote_type(V, promote_valuetype(y...))

# Subtract singleton types which are going to be replaced
function subtract_singletontype(::Type{T}, x::Pair{K}) where {T, K}
    if issingletontype(K)
        Core.Compiler.typesubtract(T, K)
    else
        T
    end
end
subtract_singletontype(::Type{T}, x::Pair{K}, y::Pair...) where {T, K} =
    subtract_singletontype(subtract_singletontype(T, y...), x)

"""
    replace(new::Function, A; [count::Integer])

Return a copy of `A` where each value `x` in `A` is replaced by `new(x)`.
If `count` is specified, then replace at most `count` values in total
(replacements being defined as `new(x) !== x`).

# Examples
```jldoctest
julia> replace(x -> isodd(x) ? 2x : x, [1, 2, 3, 4])
4-element Array{Int64,1}:
 2
 2
 6
 4

julia> replace(Dict(1=>2, 3=>4)) do kv
           first(kv) < 3 ? first(kv)=>3 : kv
       end
Dict{Int64,Int64} with 2 entries:
  3 => 4
  1 => 3
```
"""
replace(new::Callable, A; count::Integer=typemax(Int)) =
    _replace!(new, _similar_or_copy(A), A, check_count(count))

# Handle ambiguities
replace!(a::Callable, b::Pair; count::Integer=-1) = throw(MethodError(replace!, (a, b)))
replace!(a::Callable, b::Pair, c::Pair; count::Integer=-1) = throw(MethodError(replace!, (a, b, c)))
replace(a::Callable, b::Pair; count::Integer=-1) = throw(MethodError(replace, (a, b)))
replace(a::Callable, b::Pair, c::Pair; count::Integer=-1) = throw(MethodError(replace, (a, b, c)))
replace(a::AbstractString, b::Pair, c::Pair) = throw(MethodError(replace, (a, b, c)))

### replace! for AbstractDict/AbstractSet

askey(k, ::AbstractDict) = k.first
askey(k, ::AbstractSet) = k

function _replace!(new::Callable, res::T, A::T,
                   count::Int) where T<:Union{AbstractDict,AbstractSet}
    count == 0 && return res
    c = 0
    if res === A # cannot replace elements while iterating over A
        repl = Pair{eltype(A),eltype(A)}[]
        for x in A
            y = new(x)
            if x !== y
                push!(repl, x => y)
                c += 1
                c == count && break
            end
        end
        for oldnew in repl
            pop!(res, askey(first(oldnew), res))
        end
        for oldnew in repl
            push!(res, last(oldnew))
        end
    else
        for x in A
            y = new(x)
            if x !== y
                pop!(res, askey(x, res))
                push!(res, y)
                c += 1
                c == count && break
            end
        end
    end
    res
end

### replace! for AbstractArray

function _replace!(new::Callable, res::AbstractArray, A::AbstractArray, count::Int)
    c = 0
    if count >= length(A) # simpler loop allows for SIMD
        if res === A # for optimization only
            for i in eachindex(A)
                @inbounds Ai = A[i]
                y = new(Ai)
                @inbounds A[i] = y
            end
        else
            for i in eachindex(A)
                @inbounds Ai = A[i]
                y = new(Ai)
                @inbounds res[i] = y
            end
        end
    else
        for i in eachindex(A)
            @inbounds Ai = A[i]
            if c < count
                y = new(Ai)
                @inbounds res[i] = y
                c += (Ai !== y)
            else
                @inbounds res[i] = Ai
            end
        end
    end
    res
end


text:
using StokesParticles, Test

# no cl test
pList = [Particle2D(Point2D(-0.5,0.3),0.02,1),
         Particle2D(Point2D(-0.5,-0.3),0.01,1),
         Particle2D(Point2D(-0.5,-0.28),0.01,1)]
wList = [LineWall(Point2D(2.0,1.0),Point2D(2.0,0.0),0.1)]
param = sp_params()
data = scratch_data(length(pList))
update_particles_nofluid_nocl!(pList,wList,param,data)

L = 0.5
cl = generate_cell_list(pList,[-3.0,3.0,-3.0,3.0],L)
update_particles_nofluid!(pList,wList,param,data,cl)

ext: jl


text:
using WriteVTK

function vtk_connectivity_map(Nqi, Nqj = 1, Nqk = 1)
    connectivity = Array{Int, 1}(undef, Nqi * Nqj * Nqk)
    L = LinearIndices((1:Nqi, 1:Nqj, 1:Nqk))

    corners = (
        (1, 1, 1),
        (Nqi, 1, 1),
        (Nqi, Nqj, 1),
        (1, Nqj, 1),
        (1, 1, Nqk),
        (Nqi, 1, Nqk),
        (Nqi, Nqj, Nqk),
        (1, Nqj, Nqk),
    )
    edges = (
        (2:(Nqi - 1), 1:1, 1:1),
        (Nqi:Nqi, 2:(Nqj - 1), 1:1),
        (2:(Nqi - 1), Nqj:Nqj, 1:1),
        (1:1, 2:(Nqj - 1), 1:1, 1:1),
        (2:(Nqi - 1), 1:1, Nqk:Nqk),
        (Nqi:Nqi, 2:(Nqj - 1), Nqk:Nqk),
        (2:(Nqi - 1), Nqj:Nqj, Nqk:Nqk),
        (1:1, 2:(Nqj - 1), Nqk:Nqk),
        (1:1, 1:1, 2:(Nqk - 1)),
        (Nqi:Nqi, 1:1, 2:(Nqk - 1)),
        (1:1, Nqj:Nqj, 2:(Nqk - 1)),
        (Nqi:Nqi, Nqj:Nqj, 2:(Nqk - 1)),
    )
    faces = (
        (1:1, 2:(Nqj - 1), 2:(Nqk - 1)),
        (Nqi:Nqi, 2:(Nqj - 1), 2:(Nqk - 1)),
        (2:(Nqi - 1), 1:1, 2:(Nqk - 1)),
        (2:(Nqi - 1), Nqj:Nqj, 2:(Nqk - 1)),
        (2:(Nqi - 1), 2:(Nqj - 1), 1:1),
        (2:(Nqi - 1), 2:(Nqj - 1), Nqk:Nqk),
    )
    if Nqj == Nqk == 1
        @assert Nqi > 1
        corners = (corners[1:2]...,)
        edges = (edges[1],)
        faces = ()
    elseif Nqk == 1
        @assert Nqi > 1
        @assert Nqj > 1
        corners = (corners[1:4]...,)
        edges = (edges[1:4]...,)
        faces = (faces[5],)
    end

    # corners
    iter = 1
    for (i, j, k) in corners
        connectivity[iter] = L[i, j, k]
        iter += 1
    end
    # edges
    for (is, js, ks) in edges
        for k in ks, j in js, i in is
            connectivity[iter] = L[i, j, k]
            iter += 1
        end
    end
    # faces
    for (is, js, ks) in faces
        for k in ks, j in js, i in is
            connectivity[iter] = L[i, j, k]
            iter += 1
        end
    end
    # interior
    for k in 2:(Nqk - 1), j in 2:(Nqj - 1), i in 2:(Nqi - 1)
        connectivity[iter] = L[i, j, k]
        iter += 1
    end
    return connectivity
end

#=
This is the 1D WriteMesh routine
=#
function writemesh(
    base_name,
    x1;
    x2 = nothing,
    x3 = nothing,
    fields = (),
    realelems = 1:size(x1)[end],
)
    (Nqr, _) = size(x1)

    con = vtk_connectivity_map(Nqr)

    M = MeshCell{VTKCellTypes.VTKCellType, Array{Int, 1}}
    cells = Array{M, 1}(undef, length(realelems))

    for (i, e) in enumerate(realelems)
        offset = (e - 1) * Nqr
        cells[i] = MeshCell(VTKCellTypes.VTK_LAGRANGE_CURVE, offset .+ con)
    end

    if x2 == nothing
        @assert isnothing(x3)
        vtkfile =
            vtk_grid("$(base_name)", @view(x1[:]), cells; compress = false)
    elseif x3 == nothing
        vtkfile = vtk_grid(
            "$(base_name)",
            @view(x1[:]),
            @view(x2[:]),
            cells;
            compress = false,
        )
    else
        vtkfile = vtk_grid(
            "$(base_name)",
            @view(x1[:]),
            @view(x2[:]),
            @view(x3[:]),
            cells;
            compress = false,
        )
    end
    for (name, v) in fields
        vtk_point_data(vtkfile, v, name)
    end
    outfiles = vtk_save(vtkfile)
end

#=
This is the 2D WriteMesh routine
=#
function writemesh(
    base_name,
    x1,
    x2;
    x3 = nothing,
    fields = (),
    realelems = 1:size(x1)[end],
)
    @assert size(x1) == size(x2)
    (Nqr, Nqs, _) = size(x1)
    @assert Nqr == Nqs
    con = vtk_connectivity_map(Nqr, Nqs)

    M = MeshCell{VTKCellTypes.VTKCellType, Array{Int, 1}}
    cells = Array{M, 1}(undef, length(realelems))
    for (i, e) in enumerate(realelems)
        offset = (e - 1) * Nqr * Nqs
        cells[i] =
            MeshCell(VTKCellTypes.VTK_LAGRANGE_QUADRILATERAL, offset .+ con[:])
    end

    if x3 == nothing
        vtkfile = vtk_grid(
            "$(base_name)",
            @view(x1[:]),
            @view(x2[:]),
            cells;
            compress = false,
        )
    else
        vtkfile = vtk_grid(
            "$(base_name)",
            @view(x1[:]),
            @view(x2[:]),
            @view(x3[:]),
            cells;
            compress = false,
        )
    end
    for (name, v) in fields
        vtk_point_data(vtkfile, v, name)
    end
    outfiles = vtk_save(vtkfile)
end

#=
This is the 3D WriteMesh routine
=#
function writemesh(
    base_name,
    x1,
    x2,
    x3;
    fields = (),
    realelems = 1:size(x1)[end],
)
    (Nqr, Nqs, Nqt, _) = size(x1)
    @assert Nqr == Nqs == Nqt
    con = vtk_connectivity_map(Nqr, Nqs, Nqt)
    M = MeshCell{VTKCellTypes.VTKCellType, Array{Int, 1}}
    cells = Array{M, 1}(undef, length(realelems))
    for (i, e) in enumerate(realelems)
        offset = (e - 1) * Nqr * Nqs * Nqt
        cells[i] = MeshCell(VTKCellTypes.VTK_LAGRANGE_HEXAHEDRON, offset .+ con)
    end

    vtkfile = vtk_grid(
        "$(base_name)",
        @view(x1[:]),
        @view(x2[:]),
        @view(x3[:]),
        cells;
        compress = false,
    )
    for (name, v) in fields
        vtk_point_data(vtkfile, v, name)
    end
    outfiles = vtk_save(vtkfile)
end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
function fGn_predict_cond(X::AbstractVector{<:Real}, α::AbstractMatrix, u::Integer=1; bias::Bool=true, kwargs...)
    # estimation of bias  β = μy - α * μx
    n, k = size(α)  # n: step of prediction, k: number of samples used for prediction

    gy = (t,n) -> t .+ collect(1:n)  # grid of prediction, starting from t
    # gx = (t,k,u) -> collect((t-k*u+1):u:t)  # grid of samples, ending by t
    gx = (t,k,u) -> reverse(collect(t:-u:(t-k*u+1)))  # grid of samples, ending by t

    β::AbstractVector = if bias
        βm = hcat([X[gy(t,n)] - α*X[gx(t,k,u)] for t=(k*u):(length(X)-n)]...)  # t=(k*u):(length(X)-n)
        mean(βm, dims=2)[:]
    else
        zeros(n)
    end

    # conditional mean and covariance
    μ::AbstractVector = α * X[gx(length(X),k,u)] + β
#     Σ::AbstractMatrix = Σyy - Σyx * iΣxx * Σyx'

    return (μ=μ, β=β)
end

function fGn_predict_replict(X::AbstractVector{<:Real}, hurst::Real, u::Integer=1, k::Integer=1; kwargs...)
    N = length(X)
    μ = sign(hurst-0.5) * mean(X[end:-u:(end-k*u+1)])
#     μ = sign(hurst-0.5) * mean(X[end-k+1:end])

    return (μ=μ,)
end

"""
MLE of fGn model and prediction by replication.

# Args
- s: sub window size
- l: length of decorrelation
- k: number of samples used for prediction
"""
function fGn_MLE_estim_predict_replict(X::AbstractVector{<:Real}, d::Integer, s::Integer, l::Integer, k::Integer=1; kwargs...)
    est = fGn_MLE_estim(X, d, s, l; kwargs...)
    μ = sign(est.hurst-0.5) * mean(X[end-k+1:end])

    return (hurst=est.hurst, σ=est.σ, μ=μ)
end

"""
MLE of fGn model and prediction by conditional statistics.

# Args
- s: sub window size
- l: length of decorrelation
- k: number of samples used for prediction
- u: downsampling factor for samples
- n: step of prediction
"""
function fGn_MLE_estim_predict_cond(X::AbstractVector{<:Real}, d::Integer, s::Integer, l::Integer, k::Integer, u::Integer=1, n::Integer=d; bias::Bool=false, kwargs...)
    @assert 0 < k && 0 < n
    # @assert n+k*u <= length(X)

    # estimation of hurst and volatility
    est = fGn_MLE_estim(X, d, s, l; kwargs...)

    P = FractionalGaussianNoise(est.hurst, d)  # fGn with estimated hurst at step d
    gy = t -> t .+ collect(1:n)  # grid of prediction, starting from t
    gx = t -> collect(reverse(t:-u:(t-k*u+1)))  # grid of samples, ending by t

    # covariance matrices
    Σyy = covmat(P, gy(0))
    Σyx = covmat(P, gy(0), gx(0))
    Σxx = covmat(P, gx(0))
    iΣxx = pinv(Matrix(Σxx))
    Σ = Σyy - Σyx * iΣxx * Σyx'

    α::AbstractMatrix = Σyx * iΣxx  # kernels of conditional mean
    # estimation of bias  β = μy - α * μx
    β::AbstractVector = if bias
        # βm = hcat([X[gy(t)] - α*X[gx(t)] for t=1:length(X) if gx(t)[1]>0 && gy(t)[end]<=length(X)]...)  # t=(k*u):(length(X)-n)
        βm = hcat([X[gy(t)] - α*X[gx(t)] for t=k*u:length(X)-n]...)  # t=(k*u):(length(X)-n)
        mean(βm, dims=2)[:]
    else
        zeros(n)
    end

    # conditional mean and covariance
    μ = α * X[gx(length(X))] + β

    return (hurst=est.hurst, σ=est.σ, μ=μ, Σ=Σ, α=α, β=β)
end

"""
# Args
- X: each row corresponds to a symbol
"""
function fGn_MLE_estim_predict_cond(X::AbstractMatrix{<:Real}, d::Integer, s::Integer, l::Integer, k::Integer, u::Integer=1, n::Integer=d; bias::Bool=false, kwargs...)
    @assert 0 < k && 0 < n
    # @assert n+k*u <= length(X)

    # estimation of hurst and volatility
    res = []  # results of estimation
    for r=1:size(X,1)
        push!(res, fGn_MLE_estim(X[r,:], d, s, l; kwargs...))
    end

    Hs = [x.hurst for x in res]  # vector of Hurst
    # Hμ = mean(Hs)
    # Hσ = std(Hs)
    Vs = [x.σ for x in res]  # vector of volatility
    # Vμ = mean(Vs)
    # Vσ = std(Vs)

    # hurst = mean(Hs)
    hurst = median(Hs)

    P = FractionalGaussianNoise(hurst, d)  # fGn with estimated hurst at step d
    gy = t -> t .+ collect(1:n)  # grid of prediction, starting from t
    gx = t -> collect(reverse(t:-u:(t-k*u+1)))

    # covariance matrices
    Σyy = covmat(P, gy(0))
    Σyx = covmat(P, gy(0), gx(0))
    Σxx = covmat(P, gx(0))
    iΣxx = pinv(Matrix(Σxx))
    Σ = Σyy - Σyx * iΣxx * Σyx'

    α::AbstractMatrix = Σyx * iΣxx  # kernels of conditional mean
    # estimation of bias  β = μy - α * μx
    β::AbstractMatrix = if bias
        βm = [X[:,gy(t)]' - α*X[:,gx(t)]' for t=k*u:size(X,2)-n]  # t=(k*u):(length(X)-n)
        mean(βm)
    else
        zeros(size(X,1),n)
    end

    # conditional mean and covariance
    μs = X[:,gx(size(X,2))] * α' + β
    # for r=1:size(X,1)
    #     push!(μs, α*X[r,gx(size(X,2))] + β)
    # end

    return (hurst=Hs, σ=Vs, μ=μs, Σ=Σ, α=α, β=β)
end

function fGn_MLE_estim_predict_old(X::AbstractVector{<:Real}, d::Integer, s::Integer, l::Integer, k::Integer, u::Integer, n::Integer; dmode::Symbol=:regular, kwargs...)
    N = length(X)  # number of samples

    # when sub window equals to rolling window then `rolling_vectorize()` has no effect
    est = fGn_MLE_estim(X, d, s, l; kwargs...)

    # convention of time arrow: from left to right
    sgrid::AbstractVector{<:Integer} = Int[]  # grid of historical samples
    S::AbstractVector{<:Real} = Real[]  # value of historical samples

    if dmode == :logscale  # logscale downsampling
        tidx = findall(reverse(logtrain(N, N÷u)))
        # select the last k samples
        sgrid = (k==0 || k>=length(tidx)) ? tidx : tidx[end-k+1:end]
        S = X[sgrid]
    else # regular downsampling
        Sm = vec2mat(X, u, keep=:tail)
        S0 = mean(Sm, dims=1)[:]
        # S0 = Sm[1,:]
        S = (k==0 || k>=length(S0)) ? S0 : S0[end-k+1:end]
        sgrid = u*(1:length(S))
    end

    # conditional mean and covariance
    pgrid = sgrid[end] .+ (1:n)  # grid of prediction
    prd = cond_mean_cov(FractionalGaussianNoise(est.hurst, d), pgrid, sgrid, S)
    C = prd.C[end,:]

    return (hurst=est.hurst, σ=est.σ, μ=prd.μ, C=C, Σ=prd.Σ)
    # return hurst, σ, μ, σ^2*Σ
end

"""
Power-law estimator and predictor.

# Args
- X: input matrix. The first row is the observation of fBm and the others are fGn
- lags: integer time lags (increment step) used to compute each component of `X` starting from the second row. Values in `lags` must be all distinct.
- s: time lag for prediction
- d: (average) downsampling factor

# Returns
- H, σ, (μr, Vr), (μy, Vy): estimation of hurst, volatility, conditional mean and covariance of fGn and fBm.

# Notes
- This function is intended to be used with `rolling_apply` where `X` is some observation on a time window. - The time arrow is on the second dimension (i.e. horizontal) from left to right.
- μr and

"""
function powlaw_estim_predict(X::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}, s::Integer, d::Integer, k::Integer=0; kwargs...)
    # @assert s in lags
    # @assert "all values in lags must be distinct"

    N = size(X, 2)  # number of rows
    Y = view(X, 1, :)  # first row is an observation of fBm
    R = view(X, 2:size(X,1), :)  # the other rows are observations of fGn, corresponding to `lags`
    sidx = findall(lags.==s)[1]  # index of `s` in `lags`

    H, σ, opm = powlaw_estim(R, lags; kwargs...)

    # convention of time arrow: from left to right
    # method 1: logscale downsampling
    tidx0 = findall(reverse(logtrain(N, N÷d)))
    # @assert k <= length(tidx0)
    tidx = k>0 ? tidx0[end-k+1:end] : tidx0  # select the last k samples
    S = R[sidx,tidx]

    # # method 2: regular downsampling with averaging
    # S0 = mean(vec2mat(R[sidx,:], d, keep=:tail), dims=1)[:]
    # # S0 = transpose(vec2mat(R[sidx,:], d, keep=:tail))[1,:]
    # S = k>0 ? S0[end-k+1:end] : S0
    # tidx = d*(1:size(S,1))

    # conditional mean and covariance
    G = tidx[end]+1:tidx[end]+s  # grid of prediction
    # fGn: `S` is the observation of fGn corresponding to `s`
    μr, Σr = cond_mean_cov(FractionalGaussianNoise(H, s), G, tidx, S)
    # fBm: `Y` is the observation of fBm
    μy, Σy = cond_mean_cov(FractionalBrownianMotion(H), G, tidx, Y[tidx])

    return H, σ, (μr, σ^2*Σr), (μy, σ^2*Σy)
end

# """
# EXPERIMENTAL: prediction by recursive conditional mean
# No visible improvements over classical conditional mean for long range prediction.
# """
# function powlaw_estim_predict(X::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}, s::Integer, k::Integer; kwargs...)
#     # @assert k>0 && d>0
#     H, σ, opm = FracFin.powlaw_estim(X, lags; kwargs...)
#     l = size(X,2)÷4
#     Cv = cond_mean_coeff(FractionalGaussianNoise(H, lags[s]), k, l)
#     μc = Cv * X[s, end-l+1:end]
#     return H, σ, μc
# end

max_stars_repo_name: yanncalec/FracFin

max_issues_repo_name: yanncalec/FracFin

max_forks_repo_name: yanncalec/FracFin


text:
import Term.tree: Tree

tree_dict = Dict(
    "nestedasdasdsadasdasdsadsadasdasdasdsadasasdasdassfsdfdsfdsfdsfsdfsdfdsfsdfd" => Dict(
        "n1"=>1,
        "n2"=>2,
    ),
)

tree_dict_1 = Dict(
    "nested" => Dict(
        "n1"=>1,
        "n2"=>2,
    ),
    "nested2" => Dict(
        "n1"=>"a",
        "n2"=>2,
    ),
)

tree_dict_2 = Dict(
    "nested" => Dict(
        "n1"=>1,
        "n2"=>2,
    ),
    "leaf2" => 2,
    "leaf" => 2,
    "leafme" => "v",
    "canopy" => "test",
)

tree_dict_3 = Dict(
    "nested" => Dict(
        "deeper"=>Dict(
            "aleaf"=>"unbeliefable",
            "leaflet"=>"level 3"
        ),
        "n2"=>Int,
        "n3"=>1 + 2,
    ),
    "nested2" => Dict(
        "n1"=>"a",
        "n2"=>2,
    ),
)

tree_dict_4 = Dict(
    "nested" => Dict(
        "deeper"=>Dict(
            "aleaf"=>"unbeliefable",
            "leaflet"=>"level 3",
            "sodeep" => Dict(
                "a"=>4
            )
        ),
        "n2"=>Int,
        "n3"=>1 + 2,
        "adict"=>Dict(
            "x"=>2
        )
    ),
    "nested2" => Dict(
        "n1"=>"a",
        "n2"=>2,
    ),
)

@testset "\e[34mTree" begin
    # creation
    testpanel(
        Tree(tree_dict), 28, 6
    )

    testpanel(
        Tree(tree_dict_1), 15, 9
    )

    testpanel(
        Tree(tree_dict_2), 18, 11
    )

    testpanel(
        Tree(tree_dict_3), 33, 13
    )

    testpanel(
        Tree(tree_dict_4), 33, 19
    )

    # styling
    for guides_type in (:standardtree, :boldtree, :asciitree)
         @test_nothrow Tree(tree_dict; title=string(guides_type), guides_type=guides_type)
    end

    testpanel(
        Tree(
            tree_dict,
            title="my colors",
            title_style="bold red",
            node_style="blue underline",
            leaf_style="green",
            guides_style="red dim"
        ), 30, 6
    )
end


text:
const density = Property("D")
const pressuretemperature = InputPair("PT_INPUTS")

function isothermalwork(fluid::WorkingFluid,
                        temperature,
                        pressure_initial,
                        pressure_final;
                        rtol=1e-8)
    function ρ(fluid,pressure,temperature)
        update(fluid, pressuretemperature, pressure, temperature)
        keyed_output(fluid, density)
    end
    integeral, err = quadgk(x->1/ρ(fluid,x,temperature),pressure_initial,pressure_final,rtol=rtol)
    integeral
end

function adiabaticwork(fluid::WorkingFluid,
                       pressure_initial,
                       temperature_initial,
                       temperature_final;
                       rtol=1e-8)

end


text:
#module BibTex

#
# Example of a BibTex entry:
#
# @article{name,
#   author = {Lastname1, Firstname1 and Lastname2, Firstname2},
#   title = {Title of the article},
#   journal = {Fancy journal},
#   year = {1982},
#   volume = {36},
#   pages = {1000-1010},
#   url = {http://arxiv.org/abs/1801.00315},
#   doi = {10.1103/PhysRevB.90.155136}
#   }
#

mutable struct BTEntry
  ref_num::Int
  name::String
  authors ::Array{String,1}
  journal ::String
  title ::String
  doi::String
  url::String
  eprint::String
  year::String
  pages::String
  volume::String

  function BTEntry(name_::String)
    return new(0,
               name_,
               [], #authors
               "", #journal
               "", #title
               "", #doi
               "", #url
               "", #eprint
               "",  #year
               "", #pages
               "")  #volume
  end
end

function convertToMD(bt::BTEntry)::String
  md = ""
  if bt.url != ""
    md *= "[_$(bt.title)_]($(bt.url))"
  else
    md *= "_$(bt.title)_"
  end
  for a in bt.authors
    md *= ", $a"
  end
  if bt.journal != ""
    md *= ", "*bt.journal*""
    if bt.volume != ""
      md *= " $(bt.volume)"
    end
    if bt.pages != ""
      md *= ", $(bt.pages)"
    end
  end
  if bt.year != ""
    md *= " ($(bt.year))"
  end
  if bt.eprint != ""
    md *= ", "*bt.eprint
  end
  return md
end

function getEntry(contents::String,start::Int)::String
  lev = 1
  n = start
  while (lev > 0 && n < length(contents))
    if (contents[n] == '{') lev += 1 end
    if (contents[n] == '}') lev -= 1 end
    #@show n,contents[n],lev
    n += 1
  end
  return contents[start:n-1]
end

function parseBibTex(fname::String)
  contents = open(fname) do file read(file,String) end

  entries = Dict{String,BTEntry}()

  btre = r"@article{(.+?),"s

  function getField(re,source::String)::String
    if occursin(re,source)
      return strip(match(re,source).captures[1])
    end
    return ""
  end

  for m in eachmatch(btre,contents)
    name = convert(String,m.captures[1])
    bt = BTEntry(name)
    start = m.offset+length(m.match)
    entry = getEntry(contents,start)

    bt.title = getField(r"title\W*=\W*{(.+?)}"is,entry)
    bt.journal = getField(r"journal\W*=\W*{(.+?)}"is,entry)
    bt.volume = getField(r"volume\W*=\W*{(.+?)}"is,entry)
    bt.pages = getField(r"pages\W*=\W*{(.+?)}"is,entry)
    bt.year = getField(r"year\W*=\W*{(.+?)}"is,entry)
    bt.doi = getField(r"doi\W*=\W*{(.+?)}"is,entry)
    bt.eprint = getField(r"eprint\W*=\W*{(.+?)}"is,entry)
    bt.url = getField(r"url\W*=\W*{(.+?)}"is,entry)

    all_authors = getField(r"author[s]*\W*=\W*{(.+?)}"is,entry)
    if all_authors != ""
      authors = split(all_authors,"and")
      for a in authors
        a = strip(a)
        rev = r"(\w+?), (.+)"
        if occursin(rev,a)
          m = match(rev,a)
          push!(bt.authors,"$(m.captures[2]) $(m.captures[1])")
        else
          push!(bt.authors,a)
        end
      end
    end

    entries[name] = bt
  end
  return entries
end

#export parseBibTex
#end


text:
module CustomUnitRanges

filename_for_urange = joinpath(dirname(@__FILE__), "URange.jl")
filename_for_zerorange = joinpath(dirname(@__FILE__), "ZeroRange.jl")

end # module


text:
module Dec7Part1
child_to_parent = Dict()

function extract_chunks(x::String, pattern::Regex)::Array{String,1}
    matches = split(x, pattern)
    return strip.(String.(matches))
end

function capture_groups(x::String, pattern::Regex)::Tuple{String,String}
    matching = match(pattern, x)
    return strip(String(matching[1])), strip(String(matching[2]))
end

function extract_colors_number(x::String)::Tuple{Union{Nothing,String},Union{Nothing,String}}
    if x == "no other bags."
        return nothing, nothing
    end
    count, color = capture_groups(x, r"(\d)\s?([a-zA-Z ]+)\s?bags?")
    return count, color
end

function add_child_parent(parent::String, child::Tuple{Union{Nothing,String},Union{Nothing,String}})
    global child_to_parent
    child_name = child[2]
    count = child[1]
    value = get(child_to_parent, child_name, nothing)
    if value === nothing
        child_to_parent[child_name] = [(count, parent)]
    else
        push!(child_to_parent[child_name], (count, parent))
    end
end

function get_child_to_parent(parent::String, child_count)
    add_child_parent.(parent, child_count)
end
function get_child_parents(x::String)::Nothing
    strings = extract_chunks(x, r"bags contain|,")
    child_count = extract_colors_number.(strings[2:end])
    parent = strings[1]
    get_child_to_parent(parent, child_count)
    return nothing
end

function get_nested_bags(x::Dict, key::String)::Union{Nothing,Array{String,1}}
    if get(x, key, nothing) === nothing
        return nothing
    end
    parents = map(y -> y[2], x[key])
    actual_parents = []
    for parent in parents
        nested_parents = get_nested_bags(x, parent)
        if nested_parents ≠ nothing
            append!(actual_parents, nested_parents)
        end
        push!(actual_parents, parent)
    end
    return actual_parents
end

    function main(file_path::String)::Int64
    baggage_conditions = readlines(file_path)
    child_parents = get_child_parents.(baggage_conditions)
    bags = Set(get_nested_bags(child_to_parent, "shiny gold"))
    return length(bags)
end
end


text:
module JLLDateTime

import Dates;

"""
    unix2datetime(unix_timestamp)

Convert a Unix timestamp (seconds since midnight, Jan 1 1970) into a DateTime.
"""
function unix2datetime(unix_timestamp)
    return Dates.DateTime(1970, 1, 1) + Dates.Second(unix_timestamp);
end

end # JLLDateTime

max_stars_repo_name: joshua-laughner/Julia-JLLUtils

max_issues_repo_name: joshua-laughner/Julia-JLLUtils

max_forks_repo_name: joshua-laughner/Julia-JLLUtils


text:
import DataFrames
import StatsModels

function MutableDataFrame2DecisionTreeTransformer(
        feature_names::AbstractVector,
        single_label_name::Symbol;
        levels::AbstractVector = [],
        )
    dffeaturecontrasts = FeatureContrastsNotYetGenerated()
    result = MutableDataFrame2DecisionTreeTransformer(
        feature_names,
        single_label_name,
        levels,
        dffeaturecontrasts,
        )
    return result
end

"""
"""
function set_feature_contrasts!(
        x::MutableDataFrame2DecisionTreeTransformer,
        feature_contrasts::AbstractFeatureContrasts,
        )
    x.dffeaturecontrasts = feature_contrasts
    return nothing
end

"""
"""
function get_underlying(
        x::MutableDataFrame2DecisionTreeTransformer;
        saving::Bool = false,
        loading::Bool = false,
        )
    result = x.dffeaturecontrasts
    return result
end

"""
"""
function transform(
        transformer::MutableDataFrame2DecisionTreeTransformer,
        features_df::DataFrames.AbstractDataFrame,
        labels_df::DataFrames.AbstractDataFrame;
        kwargs...
        )
    single_label_name = transformer.single_label_name
    labelsarray = convert(Array, labels_df[single_label_name])
    my_formula = transformer.dffeaturecontrasts.formula_without_intercept
    my_schema = transformer.dffeaturecontrasts.schema_without_intercept
    my_formula = StatsModels.apply_schema(my_formula, my_schema)
    response, featuresarray = StatsModels.modelcols(my_formula,
                                                    features_df)
    return featuresarray, labelsarray
end

"""
"""
function transform(
        transformer::MutableDataFrame2DecisionTreeTransformer,
        features_df::DataFrames.AbstractDataFrame;
        kwargs...
        )
    my_formula = transformer.dffeaturecontrasts.formula_without_intercept
    my_schema = transformer.dffeaturecontrasts.schema_without_intercept
    my_formula = StatsModels.apply_schema(my_formula, my_schema)
    response, featuresarray = StatsModels.modelcols(my_formula,
                                                    features_df)
    return featuresarray
end

"""
"""
function fit!(
        transformer::MutableDataFrame2DecisionTreeTransformer,
        features_df::DataFrames.AbstractDataFrame,
        labels_df::DataFrames.AbstractDataFrame;
        kwargs...
        )
    return transform(transformer, features_df, labels_df)
end

"""
"""
function predict(
        transformer::MutableDataFrame2DecisionTreeTransformer,
        features_df::DataFrames.AbstractDataFrame,
        varargs...;
        kwargs...
        )
    return (transform(transformer, features_df), varargs...)
end

"""
"""
function predict_proba(
        transformer::MutableDataFrame2DecisionTreeTransformer,
        features_df::DataFrames.AbstractDataFrame;
        kwargs...
        )
    return transform(transformer, features_df)
end


text:
# This file is a part of DensityInterface.jl, licensed under the MIT License (MIT).

"""
    hasdensity(d)::Bool

Return `true` if `d` is compatible with the `DensityInterface` interface.

`hasdensity(d) == true` implies that `d` is either a density itself or has an
associated density, e.g. a probability density function or a Radon–Nikodym
derivative with respect to an implicit base measure. It also implies that the
value of that density at given points can be calculated via
[`logdensityof`](@ref) and [`densityof`](@ref).
```
"""
function hasdensity end
export hasdensity

@inline hasdensity(::Any) = false

function check_hasdensity(d)
    hasdensity(d) || throw("Object of type $(typeof(d)) is not compatible with DensityInterface")
end

"""
    logdensityof(d, x)::Real

Compute the logarithmic value of density `d` or it's associated density
at a given point `x`.

```jldoctest a
julia> hasdensity(d)
true

julia> logy = logdensityof(d, x); logy isa Real
true
```

See also [`hasdensity`](@ref) and [`densityof`](@ref).
"""
function logdensityof end
export logdensityof

"""
    logdensityof(d)

Return a function that computes the logarithmic value of density `d`
or its associated density at a given point.

```jldoctest a
julia> log_f = logdensityof(d); log_f isa Function
true

julia> log_f(x) == logdensityof(d, x)
true
```

`logdensityof(d)` defaults to `Base.Fix1(logdensityof, d)`, but may be
specialized. If so, [`logfuncdensity`](@ref) will typically have to be
specialized for the return type of `logdensityof` as well.

[`logfuncdensity`](@ref) is the inverse of `logdensityof`, so
`logfuncdensity(log_f)` must be equivalent to `d`.
"""
function logdensityof(d)
    check_hasdensity(d)
    Base.Fix1(logdensityof, d)
end

"""
    logfuncdensity(log_f)

Return a `DensityInterface`-compatible density that is defined by a given
log-density function `log_f`:

```jldoctest
julia> d = logfuncdensity(log_f);

julia> hasdensity(d) == true
true

julia> logdensityof(d, x) == log_f(x)
true
```

`logfuncdensity(log_f)` returns an instance of [`DensityInterface.LogFuncDensity`](@ref)
by default, but may be specialized to return something else depending on the
type of `log_f`). If so, [`logdensityof`](@ref) will typically have to be
specialized for the return type of `logfuncdensity` as well.

`logfuncdensity` is the inverse of `logdensityof`, so the following must
hold true:

* `logfuncdensity(logdensityof(d))` is equivalent to `d`
* `logdensityof(logfuncdensity(log_f))` is equivalent to `log_f`.

See also [`hasdensity`](@ref).
"""
function logfuncdensity end
export logfuncdensity

logfuncdensity(log_f) = LogFuncDensity(log_f)

logfuncdensity(log_f::Base.Fix1{typeof(logdensityof)}) = log_f.x

InverseFunctions.inverse(::typeof(logfuncdensity)) = logdensityof
InverseFunctions.inverse(::typeof(logdensityof)) = logfuncdensity

"""
    struct DensityInterface.LogFuncDensity{F}

Wraps a log-density function `log_f` to make it compatible with `DensityInterface`
interface. Typically, `LogFuncDensity(log_f)` should not be called
directly, [`logfuncdensity`](@ref) should be used instead.
"""
struct LogFuncDensity{F}
    _log_f::F
end
LogFuncDensity

@inline hasdensity(::LogFuncDensity) = true

@inline logdensityof(d::LogFuncDensity, x) = d._log_f(x)
@inline logdensityof(d::LogFuncDensity) = d._log_f

function Base.show(io::IO, d::LogFuncDensity)
    print(io, nameof(typeof(d)), "(")
    show(io, d._log_f)
    print(io, ")")
end

"""
    densityof(d, x)::Real

Compute the value of density `d` or its associated density at a given point
`x`.

```jldoctest a
julia> hasdensity(d)
true

julia> densityof(d, x) == exp(logdensityof(d, x))
true
```

`densityof(d, x)` defaults to `exp(logdensityof(d, x))`, but
may be specialized.

See also [`hasdensity`](@ref) and [`logdensityof`](@ref).
"""
densityof(d, x) = exp(logdensityof(d, x))
export densityof

"""
    densityof(d)

Return a function that computes the value of density `d` or its associated
density at a given point.

```jldoctest a
julia> f = densityof(d);

julia> f(x) == densityof(d, x)
true
```

`densityof(d)` defaults to `Base.Fix1(densityof, d)`, but may be specialized.
"""
function densityof(d)
    check_hasdensity(d)
    Base.Fix1(densityof, d)
end


text:
using RobustDMD, Printf, Random

# Generate DMD Synthetic Data
#------------------------------------------------------------------------------
# dimensions
m = 500;    # temporal dimension
n = 1000;    # spatial dimension
k = 3;      # number of modes
T = Float64;
# generate data
X, t, at, Bt = simDMD(m, n, k, T; seed=123456);

# additive noise
sigma = 1.0e-6 # background noise level
lds = 0.1 # larger deviation size
ldf = 0.05 # frequency of larger deviations
Xdat = (X .+ randn(real(eltype(X)),size(X))*sigma
      .+ (rand(real(eltype(X)),size(X)) .< ldf)*lds)
# specify the loss
kappa = 10.0*sigma
lossFunc1(r) = huber_func(r,kappa)
lossGrad1(r) = huber_grad!(r,kappa)
lossFunc2(r) = l2_func(r)
lossGrad2(r) = l2_grad!(r)

a0, B0 = dmdexactestimate(m, n, k, Xdat, t, dmdtype="trap");

#
# Create object
#------------------------------------------------------------------------------
params1 = DMDParams(k, Xdat, t, lossFunc1, lossGrad1);
params2 = DMDParams(k, Xdat, t, lossFunc2, lossGrad2, inner_directl2=true);
copyto!(params1.a,a0)
copyto!(params2.a,a0)

#
# Apply Solver
#------------------------------------------------------------------------------
tau = 10;
eta = 2e1;
itm = 1000;
tol = 1e-7;
ptf = 100;

function prox_stab(ar)
    # project for numerical stability
    k = length(ar) >> 1;
    @inbounds for i = 1:k
        ar[i<<1-1] = min(2.0, ar[i<<1-1]);
    end
end

opts = DMD_SVRG_Options(tau, eta, itm=itm, tol=tol, ptf=ptf,
                        prox=prox_stab);

@time obj_his1, err_his1 = solveDMD_withSVRG(params1, opts);
@time obj_his2, err_his2 = solveDMD_withSVRG(params2, opts);

@printf "%7.4e abs err evals w/ Huber\n" besterrperm(params1.a,at)
@printf "%7.4e abs err evals w/ L2\n" besterrperm(params2.a,at)


text:
module BinaryBuilder
using Libdl, LibGit2, Random, JSON
using BinaryBuilderBase
using ObjectFile
using GitHub
using Pkg, Base.BinaryPlatforms, Pkg.Artifacts
using ghr_jll

# Re-export useful stuff from Base.BinaryPlatforms:
export HostPlatform, platform_dlext, valid_dl_path, arch, libc, nbits,
       libgfortran_version, libstdcxx_version, cxxstring_abi, detect_libgfortran_version,
       detect_libstdcxx_version, detect_cxxstring_abi, call_abi, wordsize, triplet,
       select_platform, platforms_match, AbstractPlatform, Platform, os
# BinaryBuilderBase/src/Prefix.jl
export Prefix, bindir, libdirs, includedir, logdir, activate, deactivate,
       isinstalled, install, uninstall, list_tarball_files, verify, temp_prefix
# BinaryBuilderBase/src/Rootfs.jl
export supported_platforms, expand_gfortran_versions, expand_cxxstring_abis, expand_microarchitectures
# BinaryBuilderBase/src/Platforms.jl
export AnyPlatform
# BinaryBuilderBase/src/Products.jl
export Product, LibraryProduct, FileProduct, ExecutableProduct, FrameworkProduct, satisfied,
       locate, write_deps_file, variable_name
# BinaryBuilderBase/src/Dependency.jl
export Dependency, BuildDependency, HostBuildDependency
# BinaryBuilderBase/src/Sources.jl
export ArchiveSource, FileSource, GitSource, DirectorySource
# Auditor.jl
export audit, collect_files, collapse_symlinks

# Autocomplete BinaryBuilder.runshell
const runshell = BinaryBuilderBase.runshell

include("Auditor.jl")
include("Wizard.jl")

using OutputCollectors, BinaryBuilderBase, .Auditor, .Wizard

# Autocomplete BinaryBuilder.run_wizard
const run_wizard = Wizard.run_wizard

include("AutoBuild.jl")
include("Declarative.jl")
include("Logging.jl")

function __init__()
    # If we're running on Azure, enable azure logging:
    if !isempty(get(ENV, "AZP_TOKEN", ""))
        enable_azure_logging()
    end
end

get_bb_version() =
    BinaryBuilderBase.get_bbb_version(@__DIR__, "12aac903-9f7c-5d81-afc2-d9565ea332ae")
versioninfo() = BinaryBuilderBase.versioninfo(; name=@__MODULE__, version=get_bb_version())

# Precompilation ahoy!
# include("../deps/SnoopCompile/precompile/precompile_BinaryBuilder.jl")
# _precompile_()

end # module


text:
# This file contains the implementation of the general convex interface

export
  putnlcallbacks,
  clearnlcallbacks

type MSKnlinfo
  numvar :: Int
  numcon :: Int

  nlgetva :: Ptr{Void}
  nlgetsp :: Ptr{Void}

  # these are stored with 0-base indexes
  grdobjsub  :: Array{Int32,1}
  grdconsub  :: Array{Int32,1}
  grdconptr  :: Array{Int32,1}

  hessubi    :: Array{Int32,1}
  hessubj    :: Array{Int32,1}

  evalobj  :: Function
  evalconi :: Function
  grdlag   :: Function
  heslag   :: Function
  grdobj   :: Function
  grdconi  :: Function
end

# This wraps the callback method returning the sparsity of the non-linear terms.
# Depending on which arguments are NULL it should return gradient and hessian sparsity.
# NOTE: This is not documented (at all!), but in reality it is only called in a few ways.
#       While the arguments allow fetching only a portion of the hessian pattern, we ignore this
#       and always return the whole thing. It doesn't really matter since the sizes are only
#       used for pre-allocating space in MOSEK - only a reasonable upper bound is really needed.
# NOTE: The 'hessian' here is the second derivative of the lagrangian
function msk_nl_getsp_wrapper(nlhandle::    Ptr{Void},
                              numgrdobjnz:: Ptr{Int32}, # number of nonzeros in gradient of objective
                              grdobjsub::   Ptr{Int32}, # subscripts of nonzeros in gradient of objective
                              i_::          Int32,      # constraint index
                              convali::     Ptr{Bool},  # 0/1 indicating whether constraint i is non-linear
                              grdconinz::   Ptr{Int32}, # number of nonzeros in gradient of constraint i
                              grdconisub::  Ptr{Int32}, # subscripts of nonzeros in gradient of constraint i
                              yo::          Int32,      # 0/1 include objective in computation of hessian pattern
                              numycnz::     Int32,      # number of constraints to include in computation of the hessian pattern
                              ycsub::       Ptr{Int32}, # indexes of constraints to include in computation of the hessian pattern
                              maxnumhesnz:: Int32,      # lengths of hessubi and hessubj
                              numhesnz_::   Ptr{Int32}, # number of hessian nonzeros
                              hessubi::     Ptr{Int32}, # column subscrips of hessian non-zeros
                              hessubj::     Ptr{Int32}) # row subscripts of hessian non-zeros
  nlinfo = unsafe_pointer_to_objref(nlhandle) :: MSKnlinfo

  i = i_+1

  grdobjlen = length(nlinfo.grdobjsub)

  if numgrdobjnz != C_NULL
    unsafe_store!(numgrdobjnz, convert(Int32,grdobjlen))
  end

  if grdobjsub != C_NULL
    if grdobjlen > 0
      grdobjsub_a = unsafe_wrap(Array{Int32,1},grdobjsub,(grdobjlen,))
      grdobjsub_a[1:grdobjlen] = nlinfo.grdobjsub
    end
  end

  if i <= nlinfo.numcon
    if convali != C_NULL
      if nlinfo.grdconptr[i+1] > nlinfo.grdconptr[i]
        unsafe_store!(convali, convert(Int32,1))
      else
        unsafe_store!(convali, convert(Int32,0))
      end
    end

    if grdconinz != C_NULL
      unsafe_store!(grdconinz, convert(Int32, nlinfo.grdconptr[i+1] - nlinfo.grdconptr[i]))
    end

    if grdconisub != C_NULL
      num = nlinfo.grdconptr[i+1] - nlinfo.grdconptr[i]
      if num > 0
        grdconisub_a = unsafe_wrap(Array{Int32,1},grdconisub,num)
        grdconisub_a[1:num] = nlinfo.grdconsub[nlinfo.grdconptr[i]+1:nlinfo.grdconptr[i+1]]
      end
    end
  end

  numhesnz = length(nlinfo.hessubi)

  if numhesnz_ != C_NULL
    unsafe_store!(numhesnz_, convert(Int32, numhesnz))
  end

  if hessubi != C_NULL && hessubj != C_NULL && maxnumhesnz >= numhesnz
    hessubi_a = unsafe_wrap(Array{Int32,1},hessubi,(numhesnz,))
    hessubj_a = unsafe_wrap(Array{Int32,1},hessubj,(numhesnz,))

    hessubi_a[1:numhesnz] = nlinfo.hessubi
    hessubj_a[1:numhesnz] = nlinfo.hessubj
  end

  return Int32(0) :: Int32
end

function msk_nl_getva_wrapper(nlhandle    :: Ptr{Void},
                              xx_         :: Ptr{Float64}, # input
                              yo          :: Float64,
                              yc_         :: Ptr{Float64}, # input, length = numcon
                              objval      :: Ptr{Float64},
                              numgrdobjnz :: Ptr{Int32},
                              grdobjsub   :: Ptr{Int32},
                              grdobjval   :: Ptr{Float64},
                              numi_       :: Int32,
                              subi_       :: Ptr{Int32},   # input
                              conval      :: Ptr{Float64},
                              grdconptrb_ :: Ptr{Int32},   # input
                              grdconptre_ :: Ptr{Int32},   # input
                              grdconsub_  :: Ptr{Int32},   # input
                              grdconval   :: Ptr{Float64},
                              grdlag      :: Ptr{Float64},
                              maxnumhesnz :: Int32,
                              numhesnz    :: Ptr{Int32},
                              hessubi     :: Ptr{Int32},
                              hessubj     :: Ptr{Int32},
                              hesval      :: Ptr{Float64})
  nlinfo = unsafe_pointer_to_objref(nlhandle) :: MSKnlinfo

  numi = convert(Int,numi_)
  xx = unsafe_wrap(Array{Float64,1},xx_,(nlinfo.numvar,))
  yc = unsafe_wrap(Array{Float64,1},yc_,(nlinfo.numcon,))
  subi = unsafe_wrap(Array{Int32,1},subi_,(numi,)) .+ 1

  if objval != C_NULL
    unsafe_store!(objval, nlinfo.evalobj(xx))
  end

  ngrdobjnz = length(nlinfo.grdobjsub)

  if numgrdobjnz != C_NULL
    unsafe_store!(numgrdobjnz, convert(Int32,ngrdobjnz))
  end

  if grdobjsub != C_NULL && grdobjval != C_NULL
    grdobjval_a = unsafe_wrap(Array{Float64,1},grdobjval,(ngrdobjnz,))
    grdobjsub_a = unsafe_wrap(Array{Int32,1},grdobjsub,(ngrdobjnz,))

    grdobjsub = nlinfo.grdobjsub .+ 1
    nlinfo.grdobj(xx,
                  grdobjsub,
                  grdobjval_a)
    for i in 1:length(grdobjsub_a)
      grdobjsub_a[i] = nlinfo.grdobjsub[i]
    end
  end

  if numi > 0 && conval != C_NULL
    conv   = unsafe_wrap(Array{Float64,1},conval,(numi,))
    for i=1:numi
      conv[i] = nlinfo.evalconi(xx,subi[i])
    end
  end

  if grdconval != C_NULL
    grdconptrb = unsafe_wrap(Array{Int32,1},grdconptrb_,(numi,))
    grdconptre = unsafe_wrap(Array{Int32,1},grdconptre_,(numi,))
    for i=1:numi
      ptrb = grdconptrb[i]
      n    = grdconptre[i] - grdconptrb[i]
      nlinfo.grdconi(xx,subi[i],
                     unsafe_wrap(Array{Int32,1},grdconsub_+ptrb*4,(n,)) .+ 1,
                     unsafe_wrap(Array{Float64,1},grdconval +ptrb*4,(n,)))
    end
  end

  if grdlag != C_NULL
    nlinfo.grdlag(xx,yo,yc,subi,unsafe_wrap(Array{Float64,1},grdlag,(nlinfo.numvar,)))
  end

  nhesnz = length(nlinfo.hessubi)

  if numhesnz != C_NULL
    unsafe_store!(numhesnz,convert(Int32,nhesnz))
  end

  if maxnumhesnz > 0 && hessubi != C_NULL && hessubj != C_NULL && hesval != C_NULL
    hessubi_a = unsafe_wrap(Array{Int32,1},hessubi,(nhesnz,))
    hessubj_a = unsafe_wrap(Array{Int32,1},hessubj,(nhesnz,))
    nlinfo.heslag(xx,yo,yc,subi,
                  hessubi_a,
                  hessubj_a,
                  unsafe_wrap(Array{Float64,1},hesval, (nhesnz,)))

    # Hum... Very strange! If I use "hessubi_a -= 1" it appears that the
    # underlying data of the array is not the native pointer anymore, but
    # if I do it this way it is.
    for i=1:length(hessubi_a)
      hessubi_a[i] = hessubi_a[i]-convert(Int32,1)
      hessubj_a[i] = hessubj_a[i]-convert(Int32,1)
    end

    hesval_a  = unsafe_wrap(Array{Float64,1},hesval,(nhesnz,))
  end
  return convert(Int32,0) :: Int32
end

function putnlcallbacks(task::MSKtask,
                        grdobjsub :: Array{Int,1},
                        grdconsub :: Array{Int,1},
                        grdconptr :: Array{Int,1},
                        hessubi   :: Array{Int,1},
                        hessubj   :: Array{Int,1},
                        evalobj   :: Function,
                        evalconi  :: Function,
                        grdlag    :: Function,
                        grdobj    :: Function,
                        grdconi   :: Function,
                        heslag    :: Function)
  nvar = convert(Int,getnumvar(task))
  ncon = convert(Int,getnumcon(task))

  #print( ncon,length(grdconptr))
  if ( length(hessubi) !=  length(hessubj))
    error("Arrays hessubi and hessubj have mismatching lengths")
  end
  if ( length(grdconptr) != ncon+1 )
    error("Length of grdconptr should match number of constraints")
  end

  nlgetsp = cfunction(msk_nl_getsp_wrapper,
                      Int32,
                      (Ptr{Void},Ptr{Int32},Ptr{Int32},Int32,Ptr{Bool},Ptr{Int32},Ptr{Int32},Int32,Int32,Ptr{Int32},Int32,Ptr{Int32},Ptr{Int32},Ptr{Int32}))
  nlgetva = cfunction(msk_nl_getva_wrapper,
                      Int32,
                      ( Ptr{Void}, # nlhandle
                        Ptr{Float64},Float64,Ptr{Float64}, # xx,yo,yc
                        Ptr{Float64},Ptr{Int32},Ptr{Int32},Ptr{Float64}, # objval,numgrdobjnz,grdobjsub,grdobjval
                        Int32,Ptr{Int32},Ptr{Float64}, # numi,subi,conval
                        Ptr{Int32},Ptr{Int32},Ptr{Int32},Ptr{Float64},Ptr{Float64},# grdconptrb,grdconptre,grdconsub,grdconval,grdlag
                        Int32,Ptr{Int32},Ptr{Int32},Ptr{Int32},Ptr{Float64}))# maxnumhesnz, numhesnz.hessubi,hessubj,hesval

  nlinfo = MSKnlinfo(nvar,ncon,
                     nlgetsp,nlgetva,
                     round(Int32,grdobjsub .- 1),
                     round(Int32,grdconsub .- 1),
                     round(Int32,grdconptr .- 1),
                     round(Int32,hessubi .- 1),
                     round(Int32,hessubj .- 1),
                     evalobj,evalconi, grdlag,heslag,grdobj,grdconi)

  @msk_ccall("putnlfunc",
             Int32, (Ptr{Void},Any,Ptr{Void},Ptr{Void}),
             task.task, nlinfo, nlgetsp, nlgetva)
  task.nlinfo = nlinfo
end

function clearnlcallbacks(task::MSKtask)
  @msk_ccall("putnlfunc",
             Int32, (Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void}),
             task.task, C_NULL,C_NULL,C_NULL,C_NULL)
  task.nlinfo = nothing
end


text:
"""implementation of Simultaneous Sampling-and-Search Planning (SSSP)"""
module LibSSSP
export SSSP

import Printf: @sprintf, @printf
import Base: @kwdef
import ...MRMP: AbsState, Node, now, elapsed_sec, gen_uniform_sampling, get_mid_status, dist
import DataStructures: PriorityQueue, enqueue!, dequeue!
import ..Solvers: gen_g_func, get_distance_tables, get_distance_table

@kwdef mutable struct SuperNode{State<:AbsState}
    Q::Vector{Node{State}}  # set of search nodes
    next::Int64  # next agent, 0 -> fixed agents
    id::String = get_Q_id(Q, next)
    parent_id::Union{Nothing,String} = nothing  # parent node
    g::Float64 = 0.0  # g-value
    h::Float64 = 0.0  # h-value
    f::Float64 = g + h  # f-value
    depth::Int64 = 1  # depth
end

"""
    SSSP(
        config_init::Vector{State},
        config_goal::Vector{State},
        connect::Function,
        collide::Function,
        check_goal::Function;
        g_func::Function = gen_g_func(greedy = true),
        steering_depth::Int64 = 2,
        num_vertex_expansion::Int64 = 10,
        init_min_dist_thread::Float64 = 0.1,
        decreasing_rate_min_dist_thread::Float64 = 0.99,
        epsilon::Union{Float64,Nothing} = nothing,
        TIME_LIMIT::Union{Nothing,Real} = 30,
        VERBOSE::Int64 = 0,

        use_random_h_func::Bool = false,         # for ablation study
        no_roadmap_at_beginning::Bool = false,   # for ablation study

        no_fast_collision_check::Bool = false,   # use 'slow' collision checker
    )::Tuple{
        Union{Nothing,Vector{Vector{Node{State}}}},  # solution
        Vector{Vector{Node{State}}},  # roadmap
    } where {State<:AbsState}

implementation of SSSP
"""
function SSSP(
    config_init::Vector{State},
    config_goal::Vector{State},
    connect::Function,
    collide::Function,
    check_goal::Function;
    g_func::Function = gen_g_func(greedy = true),
    steering_depth::Int64 = 2,
    num_vertex_expansion::Int64 = 10,
    init_min_dist_thread::Float64 = 0.1,
    decreasing_rate_min_dist_thread::Float64 = 0.99,
    epsilon::Union{Float64,Nothing} = nothing,
    TIME_LIMIT::Union{Nothing,Real} = 30,
    VERBOSE::Int64 = 0,

    # for ablation study
    use_random_h_func::Bool = false,
    no_roadmap_at_beginning::Bool = false,

    # just for fairness of experiments
    no_fast_collision_check::Bool = false,
)::Tuple{
    Union{Nothing,Vector{Vector{Node{State}}}},  # solution
    Vector{Vector{Node{State}}},  # roadmap
} where {State<:AbsState}
    """get initial solution / refinement"""

    t_s = now()
    elapsed() = elapsed_sec(t_s)
    timeover() = TIME_LIMIT != nothing && elapsed() > TIME_LIMIT

    # number of agents
    N = length(config_init)

    # define sampler
    sampler = gen_uniform_sampling(config_init[1])

    # used in steering
    conn(q_from::State, q_to::State, i::Int64) = begin
        (isnothing(epsilon) || dist(q_from, q_to) <= epsilon) && connect(q_from, q_to, i)
    end

    # get initial roadmap by RRT-connect
    roadmaps = (
        no_roadmap_at_beginning ?
        map(
            i -> begin
                v_init = Node{State}(config_init[i], 1, [])
                v_goal = Node{State}(config_goal[i], 2, [])
                conn(v_init.q, v_goal.q, i) && push!(v_init.neighbors, v_goal.id)
                conn(v_goal.q, v_init.q, i) && push!(v_goal.neighbors, v_init.id)
                [v_init, v_goal]
            end,
            1:N,
        ) :
        gen_RRT_connect_roadmaps(
            config_init,
            config_goal,
            connect;
            steering_depth = steering_depth,
            epsilon = epsilon,
            TIME_LIMIT = (isnothing(TIME_LIMIT) ? nothing : TIME_LIMIT - elapsed()),
        )
    )
    if isnothing(roadmaps)
        VERBOSE > 0 &&
            @info @sprintf("\t%6.4f sec: failed to construct initial roadmaps\n", elapsed())
        return (nothing, map(i -> Vector{Node{State}}(), 1:N))
    end

    VERBOSE > 0 && !no_roadmap_at_beginning && @info ("\tdone, setup initial roadmaps")

    # setup distance tables
    distance_tables = get_distance_tables(roadmaps)

    # verbose
    print_progress! =
        (S::SuperNode{State}, loop_cnt::Int64; force::Bool = false) -> begin
            (VERBOSE < 1 || (!force && (loop_cnt % 100 != 0))) && return
            @printf(
                "\r\t\t%6.4f sec, explored node: %08d, f: %.4f, depth: %04d",
                elapsed(),
                loop_cnt,
                S.f,
                S.depth
            )
        end

    # setup heuristic function
    h_func(Q::Vector{Node{State}}) =
        (use_random_h_func ? rand() : sum(map(i -> distance_tables[i][Q[i].id], 1:N)) / N)

    # initial configuration
    Q_init = [roadmaps[i][1] for i = 1:N]

    # initial search node
    S_init = SuperNode(Q = Q_init, next = 1, id = get_Q_id(Q_init, 0), h = h_func(Q_init))

    k = 0
    while !timeover()
        k += 1

        # open list
        OPEN = PriorityQueue{SuperNode{State},Float64}()

        # discovered list to avoid duplication
        VISITED = Dict{String,SuperNode{State}}()

        # threshold of space-filling metric
        min_dist_thread = init_min_dist_thread * (decreasing_rate_min_dist_thread^(k - 1))

        # setup initail node
        enqueue!(OPEN, S_init, S_init.f)
        VISITED[S_init.id] = S_init

        loop_cnt = 0
        while !isempty(OPEN) && !timeover()
            loop_cnt += 1

            # pop
            S = dequeue!(OPEN)

            # check goal
            if check_goal(S.Q)
                print_progress!(S, loop_cnt, force = true)
                VERBOSE > 0 && @info @sprintf("\n\t%6.4f sec: found solution\n", elapsed())
                return (backtrack(S, VISITED), roadmaps)
            end

            # initial search or update for refine agents
            i = S.next
            j = mod1(S.next + 1, N)

            v = S.Q[i]

            # explore new states
            expand!(
                (q_from::State, q_to::State) -> conn(q_from, q_to, i),
                sampler,
                v,
                roadmaps[i],
                min_dist_thread,
                num_vertex_expansion,
                steering_depth,
            ) && (distance_tables[i] = get_distance_table(roadmaps[i]))

            # expand search node
            for p_id in vcat(v.neighbors, v.id)

                # create new configuration
                p = roadmaps[i][p_id]
                Q = copy(S.Q)
                Q[i] = p

                # check duplication and collision
                Q_id = get_Q_id(Q, j)
                haskey(VISITED, Q_id) && continue
                !no_fast_collision_check && collide(S.Q, p.q, i) && continue
                no_fast_collision_check && collide(S.Q, Q) && continue

                # create new search node
                S_new = SuperNode(
                    Q = Q,
                    next = j,
                    id = Q_id,
                    parent_id = S.id,
                    h = h_func(Q),
                    g = S.g + g_func(S.Q, Q),
                    depth = S.depth + 1,
                )

                # insert
                enqueue!(OPEN, S_new, S_new.f)
                VISITED[S_new.id] = S_new
            end
            print_progress!(S, loop_cnt, force = isempty(OPEN))
        end

        VERBOSE > 1 && println()
    end

    VERBOSE > 0 && @info @sprintf("\t%6.4f sec: failed to find solution\n", elapsed())
    return (nothing, roadmaps)
end

"""implementation of vertex expansion"""
function expand!(
    connect::Function,
    sampler::Function,
    v_from::Node{State},
    roadmap::Vector{Node{State}},
    min_dist_thread::Float64,
    num_vertex_expansion::Int64,
    steering_depth::Int64,
)::Bool where {State<:AbsState}

    updated = false
    for _ = 1:num_vertex_expansion
        # steering
        q_new = steering(connect, sampler(), v_from.q, steering_depth)
        # check space-filling metric
        if minimum(v -> dist(v.q, q_new), roadmap) > min_dist_thread
            # add vertex and edges
            u = Node(
                q_new,
                length(roadmap) + 1,
                map(v -> v.id, filter(v -> connect(q_new, v.q), roadmap)),
            )
            push!(roadmap, u)
            # update neighbors
            foreach(
                v -> push!(v.neighbors, u.id),
                filter(v -> connect(v.q, q_new), roadmap),
            )
            updated = true
        end
    end
    return updated
end

function steering(
    connect::Function,
    q_rand::State,
    q_near::State,
    steering_depth::Int64,
)::State where {State<:AbsState}

    q_h = q_rand  # probably unsafe -> eventually safe
    if !connect(q_near, q_h)
        # binary search
        q_l = q_near  # safe
        for _ = 1:steering_depth
            q = get_mid_status(q_l, q_h)
            if connect(q_near, q)
                q_l = q
            else
                q_h = q
            end
        end
        q_h = q_l
    end
    return q_h
end

"""generate initial roadmaps via RRT-Connect"""
function gen_RRT_connect_roadmaps(
    config_init::Vector{State},
    config_goal::Vector{State},
    connect::Function;
    steering_depth::Int64 = 2,
    epsilon::Union{Float64,Nothing} = 0.2,
    TIME_LIMIT::Union{Nothing,Real} = 30,
)::Union{Nothing,Vector{Vector{Node{State}}}} where {State<:AbsState}

    # utilities for timeout
    t_s = now()
    elapsed() = elapsed_sec(t_s)
    timeover() = TIME_LIMIT != nothing && elapsed() > TIME_LIMIT

    N = length(config_init)
    roadmaps = Vector{Vector{Node{State}}}()
    for i = 1:N
        rmp = gen_RRT_connect_roadmap(
            i,
            config_init[i],
            config_goal[i],
            connect;
            steering_depth = steering_depth,
            epsilon = epsilon,
            TIME_LIMIT = (isnothing(TIME_LIMIT) ? nothing : TIME_LIMIT - elapsed()),
        )
        isnothing(rmp) && return nothing
        push!(roadmaps, rmp)
    end
    return roadmaps
end

function gen_RRT_connect_roadmap(
    i::Int64,
    q_init::State,
    q_goal::State,
    connect::Function;
    steering_depth::Int64 = 2,
    epsilon::Union{Float64,Nothing} = 0.2,
    TIME_LIMIT::Union{Nothing,Real} = 30,
)::Union{Nothing,Vector{Node{State}}} where {State<:AbsState}

    # utilities for timeout
    t_s = now()
    elapsed() = elapsed_sec(t_s)
    timeover() = TIME_LIMIT != nothing && elapsed() > TIME_LIMIT

    # define sampler
    sampler = gen_uniform_sampling(q_init)

    # re-define connect function
    conn(q_from::State, q_to::State)::Bool = begin
        # check distance
        !(isnothing(epsilon) || dist(q_from, q_to) <= epsilon) && return false

        !connect(q_from, q_to, i) && return false

        return true
    end

    conn(q::State)::Bool = connect(q, i)

    # special case
    if conn(q_init, q_goal)
        v_init = Node(q_init, 1, [2])
        v_goal = Node(q_goal, 2, Vector{Int64}())
        conn(q_goal, q_init) && push!(v_goal.neighbors, v_init.id)
        return [v_init, v_goal]
    end

    # store all samples
    V1 = [q_init]
    V2 = [q_goal]

    V = V1

    # main loop
    iter = 0
    while !timeover()
        iter += 1
        from_start = (iter % 2 == 1)

        (flg, q_new) = extend!(conn, sampler(), V1, from_start, steering_depth)
        if flg != :trapped &&
           extend!(conn, q_new, V2, !from_start, steering_depth)[1] == :reached
            # fix tree
            if from_start
                V1, V2 = V2, V1
            end
            # create roadmap, insert vertices
            roadmap = [Node{State}(q_init, 1, []), Node{State}(q_goal, 2, [])]
            foreach(q -> push!(roadmap, Node{State}(q, length(roadmap) + 1, [])), V1[2:end])
            foreach(q -> push!(roadmap, Node{State}(q, length(roadmap) + 1, [])), V2[2:end])
            # create roadmap, insert edges
            for (k, v_from) in enumerate(roadmap)
                for (l, v_to) in enumerate(roadmap)
                    if k != l && conn(v_from.q, v_to.q)
                        push!(v_from.neighbors, v_to.id)
                    end
                end
                timeover() && return nothing
            end
            return roadmap
        end

        # swap tree
        V1, V2 = V2, V1
    end

    return nothing
end

"""used in RRT-Connect"""
function extend!(
    connect::Function,
    q_rand::State,  # sampled point
    V::Vector{State},  # vertices
    from_start::Bool,  # true -> tree rooted C_init, false -> tree rooted C_goal
    steering_depth::Int64,
)::Tuple{Symbol,State} where {State<:AbsState}

    # find nearest existing vertex
    q_near =
        (from_start ? argmin(q -> dist(q, q_rand), V) : argmin(q -> dist(q_rand, q), V))
    flg = :reached

    q_l = q_near  # safe sample
    q_h = q_rand  # probably unsafe sample -> eventually safe

    # steering, binary search
    if from_start ? !connect(q_near, q_h) : !(connect(q_h) && connect(q_h, q_near))
        flg = :trapped
        for _ = 1:steering_depth
            q = from_start ? get_mid_status(q_l, q_h) : get_mid_status(q_h, q_l)
            if from_start ? connect(q_near, q) : (connect(q) && connect(q, q_near))
                flg = :advanced
                q_l = q
            else
                q_h = q
            end
        end
        q_h = q_l
    end

    flg != :trapped && push!(V, q_h)
    return (flg, q_h)
end

"""generate id of search nodes"""
function get_Q_id(Q::Vector{Node{State}}, next::Int64)::String where {State<:AbsState}
    return @sprintf("%s_%d", join([v.id for v in Q], "-"), next)
end

"""obtain solution from search nodes by backtracking"""
function backtrack(
    S_fin::SuperNode{State},
    VISITED::Dict{String,SuperNode{State}},
)::Vector{Vector{Node{State}}} where {State<:AbsState}

    S = S_fin
    solution = Vector{Vector{Node{State}}}()
    while S.parent_id != nothing
        pushfirst!(solution, S.Q)
        S = VISITED[S.parent_id]
    end
    pushfirst!(solution, S.Q)
    return solution
end

end


text:
function secret_handshake(code)
    messages = ["wink", "double blink", "close your eyes", "jump"]
    answer = String[]
    for i in 1:4
        q = div(code, 2)
        if 2q != code
            push!(answer, messages[i])
        end
        code = q
    end
    q = div(code, 2)
    return 2q != code ? reverse(answer) : answer
end

max_stars_repo_name: twolodzko/Learning

max_issues_repo_name: twolodzko/Learning

max_forks_repo_name: twolodzko/Learning


text:
using Test
using LubyTransformCode

@testset "Soliton ideal distribution" begin
    N = 10
    soliton_ideal = SolitonIdeal(N)
    @test length(soliton_ideal) == N
    @test sum(soliton_ideal.probabilities) ≈ 1.0
end

@testset "Soliton robust distribution" begin
    N = 20
    soliton_robust = SolitonIdeal(N)
    @test length(soliton_robust) == N
    @test sum(soliton_robust.probabilities) ≈ 1.0
end

max_stars_repo_name: gyk/TrivialSolutions

max_issues_repo_name: gyk/TrivialSolutions

max_forks_repo_name: gyk/TrivialSolutions


text:
# TODO Move this to GPLikelihoods.jl
@doc raw"""
    StudentTLikelihood(ν::Real, σ::Real)

Likelihood with a Student-T likelihood:
```math
    p(y|f,\sigma, \nu) = \frac{\Gamma\left(\frac{\nu+1}{2}\right)}{\Gamma\left(\frac{\nu}{2}\right)\sqrt{\pi\nu}\sigma}\left(1 + \frac{1}{\nu}\left(\frac{x-\nu}{\sigma}\right)^2\right)^{-\frac{\nu+1}{2}}.
```

## Arguments
- `ν::Real`, number of degrees of freedom, should be positive and larger than 0.5 to be able to compute moments
- `σ::Real`, scaling of the inputs.
"""
struct StudentTLikelihood{Tν<:Real,Tσ<:Real,Thalfν<:Real} <: AbstractLikelihood
    ν::Tν
    σ::Tσ
    σ²::Tσ
    halfν::Thalfν
end

StudentTLikelihood(ν::Real, σ::Real) = StudentTLikelihood(ν, σ, abs2(σ), ν / 2)

function (lik::StudentTLikelihood)(f::Real)
    return Distributions.AffineDistribution(f, lik.σ, TDist(lik.ν))
end

_α(lik::StudentTLikelihood) = (lik.ν + 1) / 2

function (lik::StudentTLikelihood)(f::AbstractVector{<:Real})
    return Product(lik.(f))
end

aux_field(::StudentTLikelihood, Ω) = getproperty(Ω, :ω)

function init_aux_variables(rng::AbstractRNG, ::StudentTLikelihood, n::Int)
    return TupleVector((; ω=rand(rng, Gamma(), n)))
end

function init_aux_posterior(T::DataType, lik::StudentTLikelihood, n::Int)
    α = _α(lik)
    return For(TupleVector(; β=zeros(T, n))) do φ
        NTDist(Gamma(α, inv(φ.β))) # Distributions uses a different parametrization
    end
end

function aux_full_conditional(lik::StudentTLikelihood, y::Real, f::Real)
    return NTDist(Gamma(_α(lik), 2 / (lik.ν / abs2(lik.σ) + abs2(y - f))))
end

function aux_posterior!(
    qΩ, lik::StudentTLikelihood, y::AbstractVector, qf::AbstractVector{<:Normal}
)
    φ = only(qΩ.inds)
    map!(φ.β, y, qf) do yᵢ, fᵢ
        (lik.ν / abs2(lik.σ) + second_moment(fᵢ, yᵢ)) / 2
    end
    return qΩ
end

function auglik_potential(::StudentTLikelihood, Ω, y::AbstractVector)
    return (y .* Ω.ω,)
end

function auglik_precision(::StudentTLikelihood, Ω, ::AbstractVector)
    return (Ω.ω,)
end

function expected_auglik_potential(::StudentTLikelihood, qΩ, y::AbstractVector)
    return (tvmean(qΩ).ω .* y,)
end

function expected_auglik_precision(::StudentTLikelihood, qΩ, ::AbstractVector)
    return (tvmean(qΩ).ω,)
end

function logtilt(::StudentTLikelihood, ω::Real, y::Real, f::Real)
    return logpdf(Normal(f, sqrt(inv(ω))), y)
end

function expected_logtilt(::StudentTLikelihood, qωᵢ::NTDist{<:Gamma}, yᵢ::Real, qfᵢ::Normal)
    θ = ntmean(qωᵢ)
    return logpdf(Normal(yᵢ, sqrt(inv(θ.ω))), mean(qfᵢ)) - var(qfᵢ) * θ.ω / 2
end

function aux_prior(lik::StudentTLikelihood, y::AbstractVector{<:Real})
    return For(length(y)) do _
        NTDist(aux_prior(lik))
    end
end

aux_prior(lik::StudentTLikelihood) = Gamma(lik.halfν, lik.σ² / lik.halfν)


text:
function bool2int(vec::BitVector)

    num::Int = 0
    reversed_vec = reverse(vec)

    for n = 1 : length(vec)
        num += reversed_vec[n] * 2 ^ (n - 1)
    end

    num

end

function main()

    # open file and store content into a string
    input = open("inputs//day3.txt", "r") do io
        read(io, String)
    end

    # split and collect input into a data Vector{Vector{Char}} of '0's and '1's
    data = collect.(split(input))
    n_data = length(data)
    n_bits = length(data[1])

    # build a Matrix{Bool} from data
    matrix = Matrix{Bool}(undef, n_data, n_bits)
    for m = 1 : n_data
        for n = 1 : n_bits
            matrix[m, n] = parse(Bool, data[m][n])
        end
    end

    ########## part 1 ##########

    gamma_rate_bits = BitVector(undef, n_bits)
    for n = 1 : n_bits
        gamma_rate_bits[n] = sum(matrix[:, n]) >= n_data / 2
    end
    gamma_rate = bool2int(gamma_rate_bits)
    epsilon_rate = bool2int(gamma_rate_bits .== false)

    println("part 1: ", gamma_rate * epsilon_rate)

    ########## part 2 ##########

    matrix_O2 = copy(matrix)
    for n = 1 : n_bits
        n_data = size(matrix_O2, 1)
        if n_data == 1
            break
        end
        bit_criteria = sum(matrix_O2[:, n]) >= n_data / 2
        new_matrix_O2 = Matrix{Bool}(undef, 0, n_bits)
        for m = 1 : n_data
            if matrix_O2[m, n] == bit_criteria
                new_matrix_O2 = vcat(new_matrix_O2, matrix_O2[m, :]')
            end
        end
        matrix_O2 = copy(new_matrix_O2)
    end
    O2_rate_bits = BitVector(vec(matrix_O2))
    O2_rate = bool2int(O2_rate_bits)

    matrix_CO2 = copy(matrix)
    for n = 1 : n_bits
        n_data = size(matrix_CO2, 1)
        if n_data == 1
            break
        end
        bit_criteria = sum(matrix_CO2[:, n]) < n_data / 2
        new_matrix_CO2 = Matrix{Bool}(undef, 0, n_bits)
        for m = 1 : n_data
            if matrix_CO2[m, n] == bit_criteria
                new_matrix_CO2 = vcat(new_matrix_CO2, matrix_CO2[m, :]')
            end
        end
        matrix_CO2 = copy(new_matrix_CO2)
    end
    CO2_rate_bits = BitVector(vec(matrix_CO2))
    CO2_rate = bool2int(CO2_rate_bits)

    println("part 2: ", O2_rate * CO2_rate)

end

main()


text:
type BasicDiscUnvParameterNState{NI<:Integer, NR<:Real} <: ParameterNState{Discrete, Univariate}
  value::Vector{NI}
  loglikelihood::Vector{NR}
  logprior::Vector{NR}
  logtarget::Vector{NR}
  diagnosticvalues::Matrix
  n::Integer
  diagnostickeys::Vector{Symbol}
  copy::Function

  function BasicDiscUnvParameterNState(
    n::Integer,
    monitor::Vector{Bool}=[true; fill(false, 3)],
    diagnostickeys::Vector{Symbol}=Symbol[],
    ::Type{NI}=Int64,
    ::Type{NR}=Float64,
    diagnosticvalues::Matrix=Array(Any, length(diagnostickeys), isempty(diagnostickeys) ? 0 : n)
  )
    instance = new()

    t = [NI; fill(NR, 3)]

    l = Array(Integer, 4)
    for i in 1:4
      l[i] = (monitor[i] == false ? zero(Integer) : n)
    end

    fnames = fieldnames(BasicDiscUnvParameterNState)
    for i in 1:4
      setfield!(instance, fnames[i], Array(t[i], l[i]))
    end

    instance.diagnosticvalues = diagnosticvalues

    instance.n = n
    instance.diagnostickeys = diagnostickeys

    instance.copy = eval(codegen(:copy, instance, monitor))

    instance
  end
end

BasicDiscUnvParameterNState{NI<:Integer, NR<:Real}(
  n::Integer,
  monitor::Vector{Bool}=[true; fill(false, 3)],
  diagnostickeys::Vector{Symbol}=Symbol[],
  ::Type{NI}=Int64,
  ::Type{NR}=Float64,
  diagnosticvalues::Matrix=Array(Any, length(diagnostickeys), isempty(diagnostickeys) ? 0 : n)
) =
  BasicDiscUnvParameterNState{NI, NR}(n, monitor, diagnostickeys, NI, NR, diagnosticvalues)

function BasicDiscUnvParameterNState{NI<:Integer, NR<:Real}(
  n::Integer,
  monitor::Vector{Symbol},
  diagnostickeys::Vector{Symbol}=Symbol[],
  ::Type{NI}=Int64,
  ::Type{NR}=Float64,
  diagnosticvalues::Matrix=Array(Any, length(diagnostickeys), isempty(diagnostickeys) ? 0 : n)
)
  fnames = fieldnames(BasicDiscUnvParameterNState)
  BasicDiscUnvParameterNState(
    n, [fnames[i] in monitor ? true : false for i in 1:4], diagnostickeys, NI, NR, diagnosticvalues
  )
end

typealias DiscUnvMarkovChain BasicDiscUnvParameterNState

codegen(f::Symbol, nstate::BasicDiscUnvParameterNState, monitor::Vector{Bool}) = codegen(Val{f}, nstate, monitor)

function codegen(::Type{Val{:copy}}, nstate::BasicDiscUnvParameterNState, monitor::Vector{Bool})
  body = []
  fnames = fieldnames(BasicDiscUnvParameterNState)
  local f::Symbol # f must be local to avoid compiler errors. Alternatively, this variable declaration can be omitted

  for j in 1:4
    if monitor[j]
      f = fnames[j]
      push!(body, :(getfield(_nstate, $(QuoteNode(f)))[_i] = getfield(_state, $(QuoteNode(f)))))
    end
  end

  if !isempty(nstate.diagnosticvalues)
    push!(body, :(_nstate.diagnosticvalues[:, _i] = _state.diagnosticvalues))
  end

  @gensym _copy

  quote
    function $_copy(_nstate::BasicDiscUnvParameterNState, _state::BasicDiscUnvParameterState, _i::Integer)
      $(body...)
    end
  end
end

Base.eltype{NI<:Integer, NR<:Real}(::Type{BasicDiscUnvParameterNState{NI, NR}}) = (NI, NR)
Base.eltype{NI<:Integer, NR<:Real}(::BasicDiscUnvParameterNState{NI, NR}) = (NI, NR)

Base.(:(==)){S<:BasicDiscUnvParameterNState}(z::S, w::S) =
  reduce(&, [getfield(z, n) == getfield(w, n) for n in fieldnames(S)[1:7]])

Base.isequal{S<:BasicDiscUnvParameterNState}(z::S, w::S) =
  reduce(&, [isequal(getfield(z, n), getfield(w, n)) for n in fieldnames(S)[1:7]])

function Base.show{NI<:Integer, NR<:Real}(io::IO, nstate::BasicDiscUnvParameterNState{NI, NR})
  fnames = fieldnames(BasicDiscUnvParameterNState)
  fbool = map(n -> !isempty(getfield(nstate, n)), fnames[1:4])
  indentation = "  "

  println(io, "BasicDiscUnvParameterNState:")

  println(io, indentation*"eltype: $(eltype(nstate))")
  println(io, indentation*"number of states = $(nstate.n)")

  print(io, indentation*"monitored components:")
  if !any(fbool)
    println(io, " none")
  else
    print(io, "\n")
    for i in 1:4
      if fbool[i]
        println(io, string(indentation^2, fnames[i]))
      end
    end
  end

  print(io, indentation*"diagnostics:")
  if isempty(nstate.diagnostickeys)
    print(io, " none")
  else
    for k in nstate.diagnostickeys
      print(io, "\n")
      print(io, string(indentation^2, k))
    end
  end
end


text:
# developed with Julia 1.4.2
#
# EMS simulation with a MPC controller
# the online computing of optimal controls is written in LP form
#
# this example can be run with any JuMP compatible LP solver
# e.g replace Clp with CPLEX

using EMSx
using JuMP, Clp

include("arguments.jl")

args = parse_commandline()

mutable struct Mpc <: EMSx.AbstractController
    model::Model
    horizon::Int64
    Mpc() = new()
end

const controller = Mpc()

function EMSx.initialize_site_controller(controller::Mpc, site::EMSx.Site, prices::EMSx.Prices)

    controller = Mpc()

    model = Model(Clp.Optimizer)
    set_optimizer_attribute(model, "LogLevel", 0)

    horizon = 96
    battery = site.battery

    @variable(model, 0 <= u_c[1:horizon])
    @variable(model, 0 <= u_d[1:horizon])
    @variable(model, 0 <= x[1:horizon+1])
    @variable(model, 0 <= z[1:horizon])
    @variable(model, w[1:horizon])
    @variable(model, x0)

    @expression(model, u,  u_c - u_d)

    @constraint(model, u_c .<= battery.power*0.25)
    @constraint(model, u_d .<= battery.power*0.25)
    @constraint(model, x .<= battery.capacity)
    @constraint(model, u.+w .<= z)
    @constraint(model, x[1] == x0)
    @constraint(model, dynamics, diff(x) .== u_c*battery.charge_efficiency .-
        u_d/battery.discharge_efficiency)

    controller.model = model
    controller.horizon = horizon

    return controller

end

function EMSx.compute_control(controller::Mpc, information::EMSx.Information)

    fix(controller.model[:x0], information.soc*information.battery.capacity)
    fix.(controller.model[:w], information.forecast_load - information.forecast_pv)

    # set prices, padding out of test period prices with zero values
    prices = information.prices
    price_window = information.t:min(information.t+controller.horizon-1, size(prices.buy, 1))
    if length(price_window) != controller.horizon
        padding = controller.horizon - length(price_window)
        prices = EMSx.Prices(prices.name, vcat(prices.buy[price_window], zeros(padding)),
            vcat(prices.sell[price_window], zeros(padding)))
    else
        prices = EMSx.Prices(prices.name, prices.buy[price_window], prices.sell[price_window])
    end

    @objective(controller.model, Min,
        sum(prices.buy.*controller.model[:z]-
            prices.sell.*(controller.model[:z]-controller.model[:u]-controller.model[:w])))

    optimize!(controller.model)

    return value(controller.model[:u][1]) / (information.battery.power*0.25)

end

EMSx.simulate_sites(controller,
    joinpath(args["save"], "mpc"),
    args["price"],
    args["metadata"],
    args["test"])


text:
export Deque

mutable struct Deque{T}
    data::Vector{T}
    Deque{T}(d::T) where T = new(T[d])
end

function Base.show(io::IO, deque::Deque)
    print(io, "Deque($(deque.data))")
end

function Base.push!(deque::Deque, d)
    push!(deque.data, d)
end

function Base.pushfirst!(deque::Deque, d)
    pushfirst!(deque.data, d)
end

function Base.pop!(deque::Deque)
    data = pop!(deque.data)
    return data
end

function Base.popfirst!(deque::Deque)
    data = popfirst!(deque.data)
    return data
end


text:
export build_SOSrelaxation

function build_SOSrelaxation(relaxctx::RelaxationContext, mmtrelax_pb::SDPDual{T}; debug=false) where T<:Number
    sdpblocks = Dict{Tuple{Moment, String, Exponent, Exponent}, T}()
    sdplinsym = Dict{Tuple{Moment, String, Exponent}, T}()
    sdplin = Dict{Tuple{Moment, Exponent}, T}()
    sdpcst = Dict{Moment, T}()
    block_to_vartype = Dict{String, Symbol}()

    ## Build blocks dict
    for ((cstrname, cliquename), mmt) in mmtrelax_pb.constraints
        block_name = get_blockname(cstrname, cliquename, mmtrelax_pb)
        block_to_vartype[block_name] = mmt.matrixkind

        for ((γ, δ), poly) in mmt.mm
            for (moment, λ) in poly
                debug && (expo = product(moment.conj_part, moment.expl_part))
                debug && (@assert expo.degree.explvar == moment.expl_part.degree.explvar)
                debug && (@assert expo.degree.conjvar == moment.conj_part.degree.conjvar)
                # Check the current monomial has correct degree
                if (relaxctx.hierarchykind==:Complex) && ((moment.expl_part.degree.explvar > relaxctx.di[cstrname]) || (moment.conj_part.degree.conjvar > relaxctx.di[cstrname]))
                    warn(LOGGER, "build_SOSrelaxation(): Found exponent pair of degree $(expo.degree) > $(relaxctx.di[cstrname]) for Complex hierarchy.\n($(expo), at $((γ, δ)) of MM matrix)")
                elseif (relaxctx.hierarchykind==:Real) && ((moment.expl_part.degree.explvar > 2*relaxctx.di[cstrname]) || (moment.conj_part.degree.conjvar != 0))
                    warn(LOGGER, "build_SOSrelaxation(): Found exponent pair of degree $(expo.degree) > 2*$(relaxctx.di[cstrname]) for Real hierarchy.\n($(expo), at $((γ, δ)) of MM matrix)")
                end
                !isnan(λ) || warn(LOGGER, "build_SOSrelaxation(): isNaN ! constraint $cstrname - clique $blocname - mm entry $((γ, δ)) - moment $(moment)")

                # Add the current coeff to the SDP problem
                # Constraints are fα - ∑ Bi.Zi = 0
                if mmt.matrixkind == :SDP || mmt.matrixkind == :SDPC
                    key = (moment, block_name, γ, δ)
                    debug && (@assert !haskey(sdpblocks, key))

                    sdpblocks[key] = -λ
                elseif mmt.matrixkind == :Sym || mmt.matrixkind == :SymC
                    ## TODO: look into ht_keyindex2!, ht_keyindex for avoiding two dict table lookup
                    # Maybe implement this operation (if haskey add, else set) using 'setindex!(h::Dict{K,V}, v0, key::K) where V where K' as inspiration
                    key = (moment, block_name, product(γ, δ))
                    val = -λ * (γ!=δ ? 2 : 1)

                    addindex!(sdplinsym, val, key)

                    # haskey(sdplinsym, key) || (sdplinsym[key] = 0)
                    # sdplinsym[key] += val
                else
                    error(LOGGER, "build_SOSrelaxation(): Unhandled matrix kind $(mmt.matrixkind) for ($cstrname, $cliquename)")
                end

            end
        end
    end

    ## Build linear dict
    # Enforce clique coupling constraints on moments
    for (expo, cliques) in mmtrelax_pb.moments_overlap
        expo == Exponent() && continue
        # print_with_color(:light_cyan, "-> $expo - $(collect(cliques))\n")
        @assert length(cliques)>1
        cliqueref = first(cliques)

        refmoment = Moment(expo, cliqueref)
        for clique in setdiff(cliques, [cliqueref])
            curmoment = Moment(expo, clique)
            var = Exponent(get_ccmultvar(relaxctx, expo, cliqueref, clique))

            @assert !haskey(sdplin, (refmoment, var))
            @assert !haskey(sdplin, (curmoment, var))
            sdplin[(refmoment, var)] =  1
            sdplin[(curmoment, var)] = -1
        end
    end

    ## Build constants dict
    for (moment, fαβ) in mmtrelax_pb.objective
        # Determine which moment to affect the current coefficient.

        # Constraints are fα - ∑ Bi.Zi = 0
        addindex!(sdpcst, fαβ, moment)

        # if !haskey(sdpcst, moment)
        #     sdpcst[moment] = 0.0
        # end

        # # Constraints are fα - ∑ Bi.Zi = 0
        # sdpcst[moment] += fαβ
    end

    sosrelaxation = SDPPrimal{T}(block_to_vartype, sdpblocks, sdplinsym, sdplin, sdpcst)

    print_build_SOSrelax(relaxctx, sosrelaxation)
    return sosrelaxation
end

# """
#     α, β = split_expo(expo::Exponent)

#     Split the exponent into two exponents of conjugated and explicit variables in the complex case.
#     Real case is not supported yet.
# """
# function split_expo(relaxctx::RelaxationContext, expo::Exponent)
#     α, β = Exponent(), Exponent()

#     for (var, deg) in expo
#         product!(α, Exponent(Dict(var=>Degree(0, deg.conjvar))))
#         product!(β, Exponent(Dict(var=>Degree(deg.explvar, 0))))
#     end

#     if (relaxctx.hierarchykind == :Real) && (α.degree != Degree(0,0))
#         error(LOGGER, "split_expo(): Inconsistent degree $α, $β found for $(relaxctx.hierarchykind) hierarchy.")
#     end
#     return α, β
# end


text:
using Cassette, ReplMaker

Cassette.@context Immutabilitizer

maybefreeze(x) = x
maybefreeze(A::Array) = freeze(A)

# zeros was considered but removed after considering usage in dict's rehash!
# also removed ones, fill as well due to the above
for f in (rand, randn, *, +, -, /, Broadcast.materialize)
    @eval Cassette.overdub(::Immutabilitizer, ::typeof($f), args...) = maybefreeze($f(args...))
end

# Treat .= as an optional mutation
Cassette.overdub(::Immutabilitizer, ::typeof(Broadcast.materialize!), dest::ImmutableArray, x) = freeze(Broadcast.materialize!(thaw(dest), x))

# Create a REPL mode
function parse_to_expr(s)
    p = Meta.parse(s)
    isa(p, Expr) && p.head in (:using, :import) && return p
    quote $Cassette.@overdub($(Immutabilitizer()), $p) end
end
initrepl(parse_to_expr,
                prompt_text="Immutable> ",
                prompt_color = :blue,
                start_key=')',
                mode_name="Immutable_mode")

# Interesting examples:
# Dict rehash! fails because it uses zeros as a calloc
# Plots.plot(1:10, rand(10)) - fails with a Cassette bug
# Flux.train!()
# - fails because TrackedArray requires ::Arrays
# - Moved to the sf/Zygote branch, and then...
# - fails because it tries to materialize! into a reshaped immutable array. Fixed by making reshape a copy (that's potentially elided)
# - uses both zeros and randn to initialize its weights — making an asymmetry between immutable and not between W (from rand) and b (from zeros)
# svd(rand(10,10)) \ rand(10) -- returns an immutable array
#
#
# A common idiom:
# function f(A)
#     result = similar(A)
#     for i in eachindex(result)
#         result[i] = op(A[i])
#     end
#     return freeze(result)
# end


text:
using IterTools, ProgressBars
using JLD
using LaTeXStrings
using PyPlot

include("synthetic.jl")

SEED    = 42
Random.seed!(SEED)

RUNS    = 1
EPOCHS  = 1000
T       = n * EPOCHS

# synthetic d-dimensional least squares problem
wgen = rand(d);

# Quasi-random distributions
qrngs = Dict([("IID Uniform"     , UniformSample()),
              ("Sobol"           , SobolSample())])
qrng_names = collect(keys(qrngs))

function gen_syn_data_all()
    gaussian_samples = Dict()
    # for each distribution, pre-generate the stream of examples
    # to be used in T iterations of online SGD
    # for a specified number of runs
    iter_qrngs = ProgressBar(qrngs)
    for (qrng_name, ~) in iter_qrngs
        println(iter_qrngs, "generating data for $qrng_name...")

        data = []
        iter_runs = ProgressBar(1:RUNS)
        for run in iter_runs
            println(iter_runs, "working on run $run")
            data_run = gen_syn_data(T, qrng_name, d+1, SEED+run)
            push!(data, data_run)
        end
        gaussian_samples[qrng_name] = data
    end
    return gaussian_samples
end

function grad_errors(T::Int64, wtau::Vector{Float64}, fullgrad::Vector{Float64},
             samples::Vector{Any}, ordering::Vector{Int64})
    errs = zeros(T)
    acc = zeros(d)
    for m in 1:T
        (xi, yi) = samples[ordering[m]]
        acc += (xi'wtau - yi) * xi
        errs[m] = norm(acc/m - fullgrad)^2
    end
    return errs
end

println("generating data...")
samples_online = gen_syn_data_all()

println("begin computing gradient errors...")
result = Dict()

w0 = zeros(d)
for qrng_name in qrng_names
    err_list = zeros(RUNS, T)

    for run in 1:RUNS
        samples = samples_online[qrng_name][run]
        full_grad = w0 - wgen
        err_list[run, :] = grad_errors(T, w0, full_grad, samples, collect(1:T))
    end
    result[qrng_name] = err_list
end

println("completed")
println("plotting...")

figure(figsize = (6,6))
for qrng_name in ["IID Uniform", "Sobol"]
    err_list = result[qrng_name]
    err_mean = mean(err_list, dims=1)'[1000:1000:T]

    loglog(collect(1000:1000:T), err_mean, label=qrng_name)
end

legend(loc="lower left", fontsize=18)

xlabel(L"m", fontsize=20)
ylabel(L"\Vert \frac{1}{m}\Sigma_{t=0}^{m-1}\nabla R(w_0; x_t) - \nabla R(w_0)\Vert^2", fontsize=20)
title("Online", fontsize=24)
xticks(fontsize=20)
yticks(fontsize=20)

ax = gca()
locmin = matplotlib.ticker.LogLocator(base=10.0, subs=(0.25,0.5,0.75), numticks=50)
ax.yaxis.set_minor_locator(locmin)
ax.yaxis.set_minor_formatter(matplotlib.ticker.NullFormatter())

tight_layout()
plot_fn = joinpath(pwd(),"figs", "qmc-graderr-online.pdf")
savefig(plot_fn)
println("plot saved at $plot_fn")

max_stars_repo_name: EugeneLYC/qmc-ordering

max_issues_repo_name: EugeneLYC/qmc-ordering

max_forks_repo_name: EugeneLYC/qmc-ordering


text:
# this user supplied file creates a dictionary of arguments
# if a key is repeated, the last use of the key is used
# it is a little bit dangerous letting the user run arbitrary code
# as part of the solver
# now that this file is read inside a function, it is better encapsulated

arg_dict = Dict{String, Any}(
"physics" => "Euler",
"use_DG" => true,
"Flux_name" => "RoeFlux",
"run_type" => 5,
"jac_method" => 1,
"jac_type" => 1,
"order" => 1,
"IC_name" => "ICIsentropicVortex",
"numBC" => 1,
"BC1" => [0,1,2,3],
"BC1_name" => "isentropicVortexBC",
#"BC2_name" => "isentropicVortexBC",
"delta_t" => 0.005,
"t_max" => 10.000,
"smb_name" => "SRCMESHES/squarevortex_small.smb",
"res_abstol" => 1e-9,
"res_reltol" => -1.0,
"step_tol" => 1e-10,
"itermax" => 20,
"write_sparsity" => true,
"write_jac" => true,
"write_edge_vertnums" => true,
)