text:
module JournaledJets

using Distributed, Jets, LinearAlgebra, DistributedOperations, Serialization, Schedulers

const registry = Dict{Tuple, Any}()
let JID::Int = 1
    global next_jid
    next_jid() = (id = JID; JID += 1; (myid(),id))
end

"""
    struct JArray

Journal Array Struct for collecting Jet blocks in a 'Journal' for fault tolerance
"""
mutable struct JArray{T,N,A<:AbstractArray{T}} <: AbstractArray{T,N}
    id::Tuple{Int,Int}
    pids::Vector{Int}
    blockmap::Array{Int,N}
    localblocks::Array{Union{A,Nothing}, N}
    indices::Array{NTuple{N,UnitRange{Int}},N}
    journal::Vector{Expr}
end

getblockmap_from_id(id) = registry[id].blockmap

"""
    getblockmap_from_id!(id, _blockmap)

Pull out a block map from the registry based on id
"""
function getblockmap_from_id!(id, _blockmap)
    blockmap = localpart(_blockmap)
    blockmap .= registry[id].blockmap
end

"""
    setblockmap_from_id!(id, _blockmap)

Add a block map to the registry based using id as the 'key'
"""
setblockmap_from_id!(id, _blockmap) = registry[id].blockmap .= localpart(_blockmap)

"""
    indices(A::JArray)

return indices from JArray passed in
"""
indices(A::JArray) = A.indices

"""
    indices_from_id!(id, _indices)

return indices from registry based on passed in id
"""
indices_from_id!(id, _indices) = registry[id].indices = _indices

"""
    blocklength_from_id(id, iblock)

compute local block length based on blockmap id and local block id
"""
blocklength_from_id(id, iblock) = length(registry[id].localblocks[iblock])

"""
    blocksize_from_id(id, iblock, idim)

compute local block size based on blockmap id and local block id and a given dimension
"""
blocksize_from_id(id, iblock, idim) = size(registry[id].localblocks[iblock], idim)

"""
    initialize_local_part(id, pids, nblocks::NTuple{N,Int}, ::Type{A}) where {N,A}

Create an empty block map for a given id in the registry
"""
function initialize_local_part(id, pids, nblocks::NTuple{N,Int}, ::Type{A}) where {N,A}
    blockmap = zeros(Int, nblocks)
    localblocks = Union{Nothing,A}[nothing for idx in CartesianIndices(nblocks)]
    indices = Array{NTuple{N,UnitRange{Int}},N}(undef, nblocks)
    journal = Expr[]
    x = JArray{eltype(A),N,A}(id, pids, blockmap, localblocks, indices, journal)
    registry[id] = x
    nothing
end

"""
    initialize_local_part(id, pids, nblocks::NTuple{N,Int}, ::Type{A}) where {N,A}

Create an empty block map for a given id in the registry
"""
function fill_local_part(iblock, id, f)
    x = registry[id]
    x.blockmap[iblock] = myid()
    x.localblocks[iblock] = f(iblock)
    nothing
end

"""
    finish_local_part(id, indices)
"""
function finish_local_part(id, indices)
    x = registry[id]
    x.indices .= indices
    nothing
end

"""
    JArray(f::Function, nblocks::NTuple{N,Int}) where {N}

"""
function JArray(f::Function, nblocks::NTuple{N,Int}) where {N}
    pids = workers()
    id = next_jid()

    A = Base.return_types(f, (CartesianIndex{N},))[1]
    T = eltype(A)
    @sync for pid in pids
        @async remotecall_fetch(initialize_local_part, pid, id, pids, nblocks, A)
    end

    epmap(fill_local_part, CartesianIndices(nblocks), id, f)

    blockmap = zeros(Int, nblocks)
    x = ArrayFutures(blockmap)
    @sync for pid in pids
        @async remotecall_fetch(getblockmap_from_id!, pid, id, x)
    end
    reduce!(x)
    x = bcast(blockmap)
    @sync for pid in pids
        @async remotecall_fetch(setblockmap_from_id!, pid, id, x)
    end

    indices = getindices(id, nblocks, blockmap)

    @sync for pid in pids
        @async remotecall_fetch(finish_local_part, pid, id, indices)
    end

    local x
    if myid() ∈ pids
        x = registry[id]
    else
        localblocks = Union{Nothing,A}[nothing for idx in CartesianIndices(nblocks)]
        x = JArray{T,N,A}(id, pids, blockmap, localblocks, indices, Expr[])
        registry[id] = x
    end

    finalizer(close, x)

    x
end

"""
    close_by_id(id)

Delete the work id from the registry if it hasnt already been deleted. Used when work finishes
"""
function close_by_id(id)
    #ccall(:printf, Cvoid, (Cstring,), "close_by_id\n")
    if haskey(registry, id)
        delete!(registry, id)
    end
    nothing
end

"""
    Base.close(x::JArray)
"""
function Base.close(x::JArray)
    #ccall(:printf, Cvoid, (Cstring,), "close\n")
    @sync for pid in workers()
        @async remotecall_fetch(close_by_id, pid, x.id)
    end
    if haskey(registry, x.id)
        delete!(registry, x.id)
    end
    nothing
end

"""
    getindices(id, nblocks::NTuple{1}, blockmap)
"""
function getindices(id, nblocks::NTuple{1}, blockmap)
    indices = Vector{NTuple{1,UnitRange{Int}}}(undef, nblocks)
    i1 = 1
    for iblock in CartesianIndices(nblocks)
        i2 = i1 + remotecall_fetch(blocklength_from_id, blockmap[iblock], id, iblock) - 1
        indices[iblock] = (i1:i2,)
        i1 = i2 + 1
    end
    indices
end

"""
    getindices(id, nblocks::NTuple{N}, blockmap) where {N}
"""
function getindices(id, nblocks::NTuple{N}, blockmap) where {N}
    i1 = [Int[] for i=1:N]
    i2 = [Int[] for i=1:N]
    for idim = 1:N
        _i1 = 1
        for iblock = 1:nblocks[idim]
            _i2 = _i1 + remotecall_fetch(blocksize_from_id, blockmap[iblock], id, iblock, idim) - 1
            push!(i1[idim], _i1)
            push!(i2[idim], _i2)
            _i1 = _i2 + 1
        end
    end

    indices = Array{NTuple{N,UnitRange{Int}},N}(undef, nblocks)
    for (i,idx) in enumerate(CartesianIndices(nblocks))
        indices[idx] = ntuple(idim->i1[idim][idx[idim]]:i2[idim][idx[idim]], N)
    end

    indices
end

"""
    Serialization.serialize(S::AbstractSerializer, x::JArray{T,N}) where {T,N}

Serilization implementation for JArray struct
"""
function Serialization.serialize(S::AbstractSerializer, x::JArray{T,N}) where {T,N}
    _where = worker_id_from_socket(S.io)

    Serialization.serialize_type(S, typeof(x))
    if _where ∈ x.pids || _where == x.id[1]
        serialize(S, (true, x.id))
    else
        serialize(S, (false, x.id))
        for n in [:pids, :blockmap, :indices, :journal]
            serialize(S, getfield(x, n))
        end
    end
end

"""
    Serialization.deserialize(S::AbstractSerializer, ::Type{J}) where {T,N,A,J<:JArray{T,N,A}}

Deserilization implementation for JArray struct
"""
function Serialization.deserialize(S::AbstractSerializer, ::Type{J}) where {T,N,A,J<:JArray{T,N,A}}
    what = deserialize(S)
    id_only = what[1]
    id = what[2]

    local x
    if id_only
        x = registry[id]
    else
        pids = deserialize(S)
        blockmap = deserialize(S)
        indices = deserialize(S)
        journal = deserialize(S)
        localblocks = Union{A,Nothing}[nothing for idx in CartesianIndices(blockmap)]
        x = J(id, pids, blockmap, localblocks, indices, journal)
        registry[id] = x
    end

    x
end

"""
    getblock_and_delete_fromid!(_id, _whence, iblock)

Pop a block off the registry based on the block map id and block key
"""
function getblock_and_delete_fromid!(_id, _whence, iblock)
    x = registry[_id]
    blocks = x.localblocks
    block = copy(blocks[iblock]) # TODO... is this copy needed?
    blocks[iblock] = nothing

    blockmap = x.blockmap
    blockmap[iblock] = _whence

    if _whence ∉ x.pids
        sort!(push!(x.pids, _whence))
    end

    block
end

"""
    block_delete_fromid!(_id, _whence, iblock)

Remove a block from the registry based on the block map id and block key
"""
function block_delete_fromid!(_id, _whence, iblock)
    x = registry[_id]
    blocks = x.localblocks
    blocks[iblock] = nothing

    blockmap = x.blockmap
    blockmap[iblock] = _whence

    if _whence ∉ x.pids
        sort!(push!(x.pids, _whence))
    end
    nothing
end

"""
    update_blockmap_and_pids_fromid!(_id, _whence, iblock)

Update a block in the registry based on the block map id and block key with _whence
"""
function update_blockmap_and_pids_fromid!(_id, _whence, iblock)
    if haskey(registry, _id)
        x = registry[_id]
        x.blockmap[iblock] = _whence
        if _whence ∉ x.pids
            sort!(push!(x.pids, _whence))
        end
    end
    nothing
end

"""
    Jets.getblock(x::JArray{T,N,A}, iblock::CartesianIndex) where {T,N,A}

Get a block from a remote worked based on the iblock index
"""
function Jets.getblock(x::JArray{T,N,A}, iblock::CartesianIndex) where {T,N,A}
    _where = x.blockmap[iblock]
    _myid = myid()

    if _myid == _where
        return x.localblocks[iblock]::A
    end

    _id = x.id

    block = remotecall_fetch(getblock_and_delete_fromid!, _where, _id, _myid, iblock)::A
    x.blockmap[iblock] = _myid
    x.localblocks[iblock] = block

    if _myid ∉ x.pids
        sort!(push!(x.pids, _myid))
    end

    for pid in procs()
        if pid != _myid && pid != _where
            remotecall_fetch(update_blockmap_and_pids_fromid!, pid, _id, _myid, iblock)
        end
    end

    block
end
Jets.getblock(x::JArray, iblock::Int...) = getblock(x, CartesianIndex(iblock))

"""
    setblock_from_id!(id, xblock, iblock)

Set a block in the registry based on the blockmap id, and block index
"""
function setblock_from_id!(id, xblock, iblock)
    x = registry[id]
    _xblock = x.localblocks[iblock]
    _xblock .= xblock
    nothing
end

"""
    setblock_from_id!(x::JArray, xblock, iblock::CartesianIndex)

Set a block from remote worker based on the blockmap id, and block index
"""
function Jets.setblock!(x::JArray, xblock, iblock::CartesianIndex)
    if myid() == x.blockmap[iblock]
        _xblock = x.localblocks[iblock]
        _xblock .= xblock
        return nothing
    end

    remotecall_fetch(setblock_from_id!, x.blockmap[iblock], x.id, xblock, iblock)
    nothing
end
Jets.setblock!(x::JArray, xblock, iblock::Int...) = setblock!(x, xblock, CartesianIndex(iblock))

size_dim(x::JArray{T,N}, idim) where {T,N} = x.indices[ntuple(_idim->_idim==idim ? size(x.indices)[idim] : 1, N)...][idim][end]
Base.size(x::JArray{T,N}) where {T,N} = ntuple(idim->size_dim(x, idim), N)

"""
    getindex_from_id(id, iblock, i, N)

"""
function getindex_from_id(id, iblock, i, N)
    x = registry[id]
    iₒ = ntuple(idim->indices(x)[iblock][idim][1], N)
    δ = CartesianIndex(ntuple(idim->i[idim]-iₒ[idim]+1, N))
    a = x.localblocks[iblock][δ]
    a
end

function Base.getindex(x::JArray{T,N}, i::CartesianIndex) where {T,N}
    iblock = findfirst(rng->mapreduce(idim->i[idim]∈rng[idim], &, 1:N), x.indices)
    remotecall_fetch(getindex_from_id, x.blockmap[iblock], x.id, iblock, i, N)
end
Base.getindex(x::JArray, i::Int...) = getindex(x, CartesianIndex(i))

Base.similar(x::Nothing, ::Type{S}) where {S} = nothing
function similar_localpart_from_id(id, pids, similar_id, ::Type{A}, N, ::Type{S}) where {A,S}
    if haskey(registry, id)
        x = registry[id]
        blockmap = copy(x.blockmap)
        localblocks = Union{Nothing,A}[similar(x.localblocks[idx], S) for idx in CartesianIndices(size(x.blockmap))]
        indices = copy(x.indices)
        _x = JArray{S,N,A}(similar_id, pids, blockmap, localblocks, indices, Expr[])
        registry[similar_id] = _x
    end
    nothing
end

function Base.similar(x::JArray{T,N,A}, ::Type{S}) where {T,N,A,S}
    id = x.id
    pids = x.pids
    similar_id = next_jid()
    _A = typeof(similar(A(undef,ntuple(_->0,N)), S))
    for pid in workers()
        remotecall_fetch(similar_localpart_from_id, pid, id, pids, similar_id, _A, N, S)
    end

    if haskey(registry, similar_id)
        _x = registry[similar_id]
    else
        blockmap = copy(x.blockmap)
        localblocks = Union{Nothing,_A}[similar(x.localblocks[idx], S) for idx in CartesianIndices(size(x.blockmap))]
        indices = copy(x.indices)
        _x = JArray{S,N,_A}(similar_id, pids, blockmap, localblocks, indices, Expr[])
        registry[similar_id] = _x
    end

    finalizer(close, _x)

    _x
end

Distributed.procs(x::JArray) = x.pids
Jets.nblocks(x::JArray) = length(x.localblocks)

function Base.collect(x::JArray{T,1,A}) where {T,A}
    _x = Vector{A}(undef, length(x.localblocks))
    n = 0
    @sync for iblock = 1:length(_x)
        @async begin
            _x[iblock] = remotecall_fetch(getblock, x.blockmap[iblock], x, iblock)
        end
    end
    Jets.BlockArray(_x, [x.indices[i][1] for i=1:length(x.indices)])
end

Base.convert(::Type{S}, x::JArray{T,1,A}) where {S<:Jets.BlockArray,T,A} = collect(x)
Base.convert(::Array, x::JArray{T,1,A}) where {T,A} = convert(Array, collect(x))

JArray_local_norm(iblock, x::JArray, p) = norm(getblock(x, iblock), p)
function LinearAlgebra.norm(x::JArray{T}, p::Real=2) where {T}
    _T = float(real(T))
    z = epmap(JArray_local_norm, 1:length(x.localblocks), x, p)
    if p == Inf
        maximum(z)
    elseif p == -Inf
        minimum(z)
    elseif p == 0 || p == 1
        sum(z)
    else
        _p = _T(p)
        mapreduce(_z->_z^p, +, z)^(one(_T)/_p)
    end
end

JArray_local_dot(iblock, x, y) = dot(getblock(x, iblock), getblock(y, iblock))
function LinearAlgebra.dot(x::JArray, y::JArray)
    z = epmap(JArray_local_dot, 1:length(x.localblocks), x, y)
    sum(z)
end

JArray_local_extrema(iblock, x) = extrema(getblock(x, iblock))
function Base.extrema(x::JArray)
    mnmx = epmap(JArray_local_extrema, 1:length(x.localblocks), x)
    mn, mx = mnmx[1]
    for i = 2:length(mnmx)
        _mn, _mx = mnmx[i]
        _mn < mn && (mn = _mn)
        _mx > mx && (mx = _mx)
    end
    mn,mx
end

JArray_local_fill!(iblock, x, a) = fill!(getblock(x, iblock), a)
function Base.fill!(x::JArray, a)
    epmap(JArray_local_fill!, 1:length(x.localblocks), x, a)
    x
end

function Base.copy(x::JArray{T,N,A}) where {T,N,A}
    _x = JArray(i->copy(getblock(x,i))::A, size(x.localblocks))
    _x.journal = deepcopy(x.journal)
    _x
end
Base.deepcopy(x::JArray) = copy(x)
# -->

# JArray broadcasting implementation --<
struct JArrayStyle <: Broadcast.AbstractArrayStyle{1} end
Base.BroadcastStyle(::Type{<:JArray}) = JArrayStyle()
JArrayStyle(::Val{1}) = JArrayStyle()

Base.similar(bc::Broadcast.Broadcasted{JArrayStyle}, ::Type{T}) where {T} = similar(find_jarray(bc), T)
find_jarray(bc::Broadcast.Broadcasted) = find_jarray(bc.args)
find_jarray(args::Tuple) = find_jarray(find_jarray(args[1]), Base.tail(args))
find_jarray(x) = x
find_jarray(a::JArray, rest) = a
find_jarray(::Any, rest) = find_jarray(rest)

Jets.getblock(A::JArray, ::Type{<:Any}, iblock, indices) = getblock(A, iblock)

function JArray_bcast_local_copyto!(iblock, dest, bc, S)
    _x = getblock(bc, S, iblock, dest.indices[iblock][1])
    setblock!(dest, _x, iblock)
    nothing
end
function Base.copyto!(dest::JArray{T,1,<:AbstractArray{T,N}}, bc::Broadcast.Broadcasted{JArrayStyle}) where {T,N}
    S = Broadcast.DefaultArrayStyle{N}
    epmap(JArray_bcast_local_copyto!, 1:length(dest.localblocks), dest, bc, S)
    dest
end
# -->

#
# JetJSpace
#
struct JetJSpace{T,S<:Jets.JetAbstractSpace{T}} <: Jets.JetAbstractSpace{T,1}
    spaces::JArray{S,1,Array{S,1}}
    indices::Vector{UnitRange{Int}}
end

JetJSpace_length(iblock, spaces) = length(getblock(spaces, iblock)[1])

function JetJSpace(spaces::JArray{S,1,A}) where {S<:JetAbstractSpace, A<:AbstractArray{S}}
    n = epmap(JetJSpace_length, 1:length(spaces.localblocks), spaces)
    iₒ = 1
    indices = Vector{UnitRange{Int}}(undef, length(spaces))
    for iblock = 1:length(indices)
        indices[iblock] = iₒ:(iₒ+n[iblock]-1)
        iₒ = indices[iblock][end] + 1
    end
    JetJSpace(spaces, indices)
end

Jets.indices(R::JetJSpace) = R.indices

Base.size(R::JetJSpace) = (R.indices[end][end],)
Base.eltype(R::JetSpace{T}) where {T} = T
Base.eltype(R::Type{JetJSpace{T}}) where {T} = T
Base.eltype(R::Type{JetJSpace{T,S}}) where {T,S} = T

Jets.space(R::JetJSpace, iblock::Integer) = getblock(R.spaces, iblock)[1]
Jets.nblocks(R::JetJSpace) = length(R.spaces)
Distributed.procs(R::JetJSpace) = procs(R.spaces)

for f in (:Array, :ones, :rand, :zeros)
    @eval (Base.$f)(R::JetJSpace) = JArray(iblock->($f)(space(R, iblock[1])), (nblocks(R),))
end

#
# Block operator
#
JetJBlock_dom(iblock, ops) = domain(getblock(ops, 1, iblock[1])[1])
JetJBlock_rng(iblock, ops) = range(getblock(ops, iblock[1], 1)[1])
JetJBlock_spc_length(iblock, blkspaces) = length(getblock(blkspaces, iblock)[1])
function Jets.JetBlock(ops::JArray{T,2}) where {T<:Jop}
    n1,n2 = size(ops)

    local dom
    if n2 == 1
        dom = remotecall_fetch(JetJBlock_dom, ops.blockmap[1,1], 1, ops)
    else
        indices_dom = Vector{UnitRange{Int}}(undef, n2)
        blkspaces_dom = JArray(iblock->[JetJBlock_dom(iblock, ops)], (n2,))
        n = epmap(Int, JetJBlock_spc_length, 1:n2, blkspaces_dom)
        i1 = 1
        for i = 1:length(n)
            i2 = i1 + n[i] - 1
            indices_dom[i] = i1:i2
            i1 = i2 + 1
        end
        dom = JetJSpace(blkspaces_dom, indices_dom)
    end

    local rng
    if n1 == 1
        rng = remotecall_fetch(JetJBlock_rng, ops.blockmap[1,1], 1, ops)
    else
        indices_rng = Vector{UnitRange{Int}}(undef, n1)
        blkspaces_rng = JArray(iblock->[JetJBlock_rng(iblock, ops)], (n1,))
        n = epmap(Int, JetJBlock_spc_length, 1:n1, blkspaces_rng)
        i1 = 1
        for i = 1:length(n)
            i2 = i1 + n[i] - 1
            indices_rng[i] = i1:i2
            i1 = i2 + 1
        end
        rng = JetJSpace(blkspaces_rng, indices_rng)
    end
    Jet(dom = dom, rng = rng, f! = JetJBlock_f!, df! = JetJBlock_df!, df′! = JetJBlock_df′!, s = (ops=ops,))
end

function addmasterpid(pids)
    if myid() ∉ pids
        return [myid();pids]
    end
    pids
end

function JetJBlock_local_f!(iblock, d, _m, ops)
    dᵢ = getblock(d, iblock)
    opsᵢ = getblock(ops, iblock)[1]
    m = localpart(_m)
    mul!(dᵢ, opsᵢ, m)
    nothing
end

function JetJBlock_f!(d::JArray, m::AbstractArray; ops, kwargs...)
    _m = bcast(m, procs())
    epmap(JetJBlock_local_f!, 1:nblocks(d), d, _m, ops)
    d
end

function JetJBlock_local_df!(iblock, d, _m, ops)
    dᵢ = getblock(d, iblock)
    opsᵢ = getblock(ops, iblock)[1]
    m = localpart(_m)
    mul!(dᵢ, JopLn(opsᵢ), m)
    nothing
end

function JetJBlock_df!(d::JArray, m::AbstractArray; ops, kwargs...)
    _m = bcast(m, procs())
    epmap(JetJBlock_local_df!, 1:nblocks(d), d, _m, ops)
    d
end

function JetJBlock_local_df′!(iblock, _m, d, ops)
    m = localpart(_m)
    dᵢ = getblock(d, iblock)
    opsᵢ = JopLn(getblock(ops, iblock)[1])
    m .+= opsᵢ' * dᵢ
    nothing
end

function JetJBlock_df′!(m::AbstractArray, d::JArray; ops, kwargs...)
    m .= 0
    epmapreduce!(m, JetJBlock_local_df′!, 1:nblocks(d), d, ops)
end

function JetJBlock_local_point!(iblock, j, mₒ)
    op = getblock(state(j).ops, iblock)[1]
    Jets.point!(jet(op), getblock(mₒ, iblock))
    nothing
end

function Jets.point!(jet::Jet{D,R,typeof(JetJBlock_f!)}, mₒ::AbstractArray) where {D<:Jets.JetAbstractSpace, R<:Jets.JetAbstractSpace}
    epmap(JetJBlock_local_point!, CartesianIndices(size(state(jet).ops)), jet, mₒ)
    jet
end

Jets.getblock(jet::Jet{D,R,typeof(JetJBlock_f!)}, i, j) where {D,R} = state(jet).ops[i,j]
Jets.getblock(A::JopLn{T}, i, j) where {T<:Jet{<:JetAbstractSpace,<:JetAbstractSpace,typeof(JetJBlock_f!)}} = JopLn(getblock(jet(A), i, j))
Jets.getblock(A::JopNl{T}, i, j) where {T<:Jet{<:JetAbstractSpace,<:JetAbstractSpace,typeof(JetJBlock_f!)}} = getblock(jet(A), i, j)
Jets.getblock(A::T, i, j) where {J<:Jet{<:JetAbstractSpace,<:JetAbstractSpace,typeof(JetJBlock_f!)},T<:JopAdjoint{J}} = getblock(A.op, j, i)'
Jets.getblock(::Type{JopNl}, A::Jop{T}, i, j) where {T<:Jet{<:JetAbstractSpace,<:JetAbstractSpace,typeof(JetJBlock_f!)}} = getblock(jet(A), i, j)::JopNl
Jets.getblock(::Type{JopLn}, A::Jop{T}, i, j) where {T<:Jet{<:JetAbstractSpace,<:JetAbstractSpace,typeof(JetJBlock_f!)}} = JopLn(getblock(jet(A), i, j))

#
# Journal
#

# partial operator for journal replay
function JetJSubBlock_f!(d::AbstractArray, m::AbstractArray; ops, iblock, kwargs...)
    opsᵢ = getblock(ops, iblock, 1)[1]
    mul!(d, opsᵢ, m)
    d
end

function JetJSubBlock_df!(d::AbstractArray, m::AbstractArray; ops, iblock, kwargs...)
    opsᵢ = getblock(ops, iblock, 1)[1]
    mul!(d, JopLn(jet(opsᵢ)), m)
    d
end

journal_getblock(x::JArray{T}, iblock::CartesianIndex) where {T} = getblock(x, iblock)
journal_getblock(x::JopNl, iblock::CartesianIndex) = JopNl(journal_getblock(jet(x), iblock))
journal_getblock(x::JopLn, iblock::CartesianIndex) = JopLn(journal_getblock(jet(x), iblock))
#journal_getblock(x, iblock::CartesianIndex) = x
journal_getblock(x, iblock::Int) = journal_getblock(x, CartesianIndex(iblock))

function journal_getblock(x::Jet, iblock::CartesianIndex)
    dom = domain(x)
    rng = range(getblock(state(x).ops, iblock.I[1], 1)[1])
    Jet(;dom = dom, rng = rng, f! = JetJSubBlock_f!, df! = JetJSubBlock_df!, df′! = _->@info("not implemented"), s=(ops=state(x).ops, iblock=iblock.I[1]))
end

journal!(x::JArray, expr) = push!(x.journal, expr)
journal!(x::Jop, expr) = journal!(state(x).ops, expr)

emptyjournal!(x::JArray) = x.journal = x.journal[1:1]
emptyjournal!(x::Jop) = emptyjournal!(state(x).ops)

journal(x::JArray) = x.journal
journal(x::Jop) = state(x).ops.journal

macro journal(expr::Expr)
    sexpr = string(expr)

    quote
        x = $(esc(expr))
        journal!(x, Meta.parse($sexpr))
        x
    end
end

macro journal(x::Symbol, expr::Expr)
    sexpr = string(expr)
    quote
        journal!($(esc(x)), Meta.parse($sexpr))
        $(esc(expr))
    end
end

macro journal(r::QuoteNode, x::Symbol, expr::Expr)
    sexpr = string(expr)
    quote
        $(esc(r)) == :reset && emptyjournal!($(esc(x)))
        journal!($(esc(x)), Meta.parse($sexpr))
        $(esc(expr))
    end
end

function replay_construct!(mod, x::JArray, iblock::CartesianIndex)
    expr = x.journal[1]

    # using the JArray constructor is a special case
    for iexpr = 1:length(expr.args)
        if expr.args[iexpr] == :JArray
            _where = x.blockmap[iblock]
            _myid = myid()
            _id = x.id

            if _myid != _where
                remotecall_fetch(block_delete_fromid!, _where, _id, _myid, iblock)
                x.blockmap[iblock] = _myid
            end

            if _myid ∉ x.pids
                sort!(push!(x.pids, _myid))
            end

            for pid in procs()
                if pid != _myid && pid != _where
                    remotecall_fetch(update_blockmap_and_pids_fromid!, pid, _id, _myid, iblock)
                end
            end

            _expr = expr.args[iexpr+1]
            x.localblocks[iblock] = @eval mod $(_expr)($iblock)
            return nothing
        end
    end

    _expr = metablocks(mod, expr, iblock)
    x.localblocks[iblock] = @eval mod $(_expr)

    nothing
end
replay_construct!(mod, x::JArray, iblock) = replay_construct!(mod, x, CartesianIndex(iblock))

function replay!(mod, symx::String, x::JArray, iblock)
    replay_construct!(mod, x, iblock)
    for expr in x.journal[2:end]
        if block_in_play(Symbol(symx), expr, iblock)
            _expr = metablocks(mod, expr, iblock)
            @eval mod $(_expr)
        end
    end
end
replay!(mod, symx, x::Jop, iblock::Int) = replay!(mod, symx, state(x).ops, iblock)

macro replay(x, iblock)
    _x = string(x)
    :(replay!(@__MODULE__, $_x, $(esc(x)), $(esc(iblock))))
end

function block_in_play(symx::Symbol, expr::Expr, iblock)
    inplay = true
    args = expr.args
    for (iarg, arg) in enumerate(args)
        if arg == :getblock
            isx = args[iarg+1] == symx
            isi = CartesianIndex(iblock) == CartesianIndex(eval(args[iarg+2]))
            (isx && isi) && (return true)
            (isx && !isi) && (return false)
        end
        inplay = block_in_play(symx, arg, iblock)
        inplay == false && (return inplay)
    end
    true
end
block_in_play(symx::Symbol, notexpr, iblock) = true

function metablocks!(mod::Module, expr::Expr, iblock)
    if expr.head == :call
        _args = @view expr.args[2:end]
        metablocks!(mod, _args, iblock)
    else
        metablocks!(mod, expr.args, iblock)
    end
    expr
end

function metablocks!(mod::Module, arg::Symbol, iblock)
    e = :(@isdefined $arg)
    isdefined = @eval mod $e
    if isdefined
        _arg = @eval mod $arg
        T = Union{
            JArray,
            Jet{<:JetAbstractSpace,<:JetAbstractSpace,typeof(JetJBlock_f!)},
            JopNl{<:Jet{<:JetAbstractSpace,<:JetAbstractSpace,typeof(JetJBlock_f!)}},
            JopLn{<:Jet{<:JetAbstractSpace,<:JetAbstractSpace,typeof(JetJBlock_f!)}}
            }
        isblocked = isa(_arg, T)
        if isblocked
            return :(JournaledJets.journal_getblock($arg, $iblock))
        end
    end
    arg
end

function metablocks!(mod::Module, args::AbstractArray, iblock)
    for iarg in 1:length(args)
        args[iarg] = metablocks!(mod, args[iarg], iblock)
    end
    args
end

metablocks!(mod, arg, iblock) = arg

function metablocks(mod::Module, expr::Expr, iblock)
    _expr = deepcopy(expr)
    metablocks!(mod, _expr, iblock)
    _expr
end

export JArray, JetJSpace, @journal, @replay

end


text:
module curved_cantilever_examples
##  Curved (45 deg circular arc) cantilever, Transverse force at the tip
#
# Large-deflection problem. Solved many times:
#
# Bathe, Boulourchi 1979
# Simo, Vu-Quoc 1986
# Cardona, Geradin 1988
# Krysl 1993
# Krysl, FAESOR script curved_cantilever
#
# Present calculation refers to the data of
# Reference: EULERIAN FORMULATION FOR LARGE-DISPLACEMENT ANALYSIS OF SPACE
# FRAMES, by B.  A.  lzzuddin I and A.  S.  Elnashai,
# Journal of Engineering Mechanics, Vol.  119, No.  3,  March, 1993.

using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.CrossSectionModule: CrossSectionRectangle
using FinEtoolsFlexStructures.MeshFrameMemberModule: frame_member, merge_members
using FinEtoolsFlexStructures.FEMMCorotBeamModule
using FinEtoolsFlexStructures.FEMMCorotBeamModule: FEMMCorotBeam
stiffness = FEMMCorotBeamModule.stiffness
geostiffness = FEMMCorotBeamModule.geostiffness
mass = FEMMCorotBeamModule.mass
distribloads_global = FEMMCorotBeamModule.distribloads_global
restoringforce = FEMMCorotBeamModule.restoringforce
gyroscopic = FEMMCorotBeamModule.gyroscopic
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using LinearAlgebra: dot
using Arpack
using LinearAlgebra
using SparseArrays
using VisualStructures: plot_points, plot_nodes, plot_midline, render, plot_space_box, plot_solid, space_aspectratio
using PlotlyJS
using JSON

# Comparison with reference solution.

izzelna_ux =[    0.0317    0.6048
-0.0628    0.6048
8.7273   34.7095
16.6668   85.8664
23.7557  134.8918
32.6404  213.7587
37.8389  264.9157
41.7141  318.2041
44.9277  382.1503
47.9523  446.0965
50.4098  510.0426
51.9220  578.2519
53.9069  657.1188
53.9069  657.1188];

izzelna_uy =[-0.1574    0.6048
0.0317    2.7364
1.1659   73.0772
2.7727  132.7603
4.9466  211.6272
7.1205  284.0995
8.9164  350.1772
10.2396  414.1234
11.3738  480.2011
13.1697  548.4104
14.0203  618.7511
14.6820  657.1188
14.5875  657.1188];

izzelna_uz =[   -0.5354    2.7364
-0.2519    2.7364
-1.1970   51.7618
-2.2367   90.1295
-4.3161  124.2341
-6.8681  183.9172
-9.1366  235.0741
-11.8776  288.3626
-14.5241  354.4403
-16.9815  411.9919
-19.3445  475.9380
-21.7074  550.5419
-23.9759  618.7511
-25.3936  661.3819];

function curved_cantilever()
    # Parameters:
    E=1e7;# lb.in^-2
    nu=0.3;
    h= 1; b =1;# in
    Fmag=600;#  lb
    radius=100;# in
    ang=45;
    nel=8;# number of elements
    utol=1e-6;
    maxit = 150

    # Cross-sectional properties
    cs = CrossSectionRectangle(s -> b, s -> h, s -> [0.0, 0.0, 1.0])

    # Select the number of elements per leg.
    fens, fes = frame_member([0 0 0; nel 0 0], nel, cs)
    for i in 1:count(fens)
        k = fens.xyz[i, 1]
        fens.xyz[i, :] .= (radius*cos(k/nel*ang*2*pi/360), radius*sin(k/nel*ang*2*pi/360), 0)
    end
    clampedn = [1]
    tipn = [nel+1];

    tolerance = radius/nel/100;

    # Material properties
    material = MatDeforElastIso(DeforModelRed3D, 0.0, E, nu, 0.0)

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    for i in [1, 2, 3, 4, 5, 6]
        setebc!(dchi, clampedn, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi);

    # Additional fields
    u1 = deepcopy(u0)
    Rfield1 = deepcopy(Rfield0)
    rhs = gathersysvec(dchi);
    TMPv = deepcopy(rhs)
    utol = 1e-13*dchi.nfreedofs;

    # tbox = plot_space_box([[0.0*radius 0.0*radius 0.0*radius]; [1.15*radius 0.75*radius 2.0*radius]])
    # tshape0 = plot_solid(fens, fes; x = geom0.values, u = 0.0.*dchi.values[:, 1:3], R = Rfield0.values, facecolor = "rgb(125, 155, 125)", opacity = 0.3);
    # plots = cat(tbox,  tshape0; dims = 1)
    # pl = render(plots)
    # sleep(0.5)

    femm = FEMMCorotBeam(IntegDomain(fes, GaussRule(1, 2)), material)
    loadbdry = FESetP1(reshape(tipn, 1, 1))
    lfemm = FEMMBase(IntegDomain(loadbdry, PointRule()))

    load_parameters = 0.0:0.05:1;
    tip_displacement = fill(0.0, length(load_parameters), 3);
    step = 0
    for load_parameter in load_parameters
        applyebc!(dchi) # Apply boundary conditions
        u1.values[:] = u0.values[:]; # guess
        Rfield1.values[:] = Rfield0.values[:]; # guess
        fi = ForceIntensity(FFlt[0, 0, load_parameter*Fmag, 0, 0, 0]);

        println("Load: $load_parameter")
        iter = 1;
        while true
            F = distribloads(lfemm, geom0, dchi, fi, 3);
            Fr = restoringforce(femm, geom0, u1, Rfield1, dchi);       # Internal forces
            @. rhs = F + Fr;
            K = stiffness(femm, geom0, u1, Rfield1, dchi) + geostiffness(femm, geom0, u1, Rfield1, dchi);
            dchi = scattersysvec!(dchi, (K)\rhs); # Disp. incr
            u1.values[:] += (dchi.values[:,1:3])[:];   # increment displacement
            update_rotation_field!(Rfield1, dchi)
            print("$iter: ||du||=$(maximum(abs.(dchi.values[:])))\n")
            if maximum(abs.(dchi.values[:])) < utol # convergence check
                break;
            end
            if (iter > maxit)# bailout for failed convergence
                error("Possible failed convergence");
            end
            iter += 1;
        end
        u0.values[:] = u1.values[:];       # update the displacement
        Rfield0.values[:] = Rfield1.values[:]; # update the rotations

        step = step + 1
        tip_displacement[step, :] .= u1.values[tipn[1], :]

    end
    step = step + 1

    plots = cat(
    scatter(;x=tip_displacement[:,1], y=load_parameters.*Fmag, mode="lines", line_color = "rgb(0, 0, 0)", name="Tip-x"),
    scatter(;x=tip_displacement[:,2], y=load_parameters.*Fmag, mode="lines", line_color = "rgb(0, 0, 0)", name="Tip-y"),
    scatter(;x=tip_displacement[:,3], y=load_parameters.*Fmag, mode="lines", line_color = "rgb(0, 0, 0)", name="Tip-z"),
    scatter(;x=izzelna_ux[:,1], y=izzelna_ux[:,2], mode="markers", line_color = "rgb(255, 0, 0)", name="Ref-x"),
    scatter(;x=izzelna_uy[:,1], y=izzelna_uy[:,2], mode="markers", line_color = "rgb(0, 255, 0)", name="Ref-y"),
    scatter(;x=izzelna_uz[:,1], y=izzelna_uz[:,2], mode="markers", line_color = "rgb(0, 0, 255)", name="Ref-z");
    dims = 1)
    layout = Layout(width=700, height=500, xaxis=attr(title="Tip displacement [in]"),
    yaxis=attr(title="Load [lb]"))
    pl = plot(plots, layout)
    display(pl)

    return true
end # curved_cantilever

function curved_cantilever_thin()
    # Parameters:
    E=1e7;# lb.in^-2
    nu=0.3;
    thinness_factor = 100
    h= 1.0/thinness_factor; b = 1.0/thinness_factor;# in
    Fmag=600.0/(thinness_factor^4);#  lb
    radius=100;# in
    ang=45;
    nel=8;# number of elements
    utol=1e-6;
    maxit = 150

    # Cross-sectional properties
    cs = CrossSectionRectangle(s -> b, s -> h, s -> [0.0, 0.0, 1.0])

    # Select the number of elements per leg.
    fens, fes = frame_member([0 0 0; nel 0 0], nel, cs)
    for i in 1:count(fens)
        k = fens.xyz[i, 1]
        fens.xyz[i, :] .= (radius*cos(k/nel*ang*2*pi/360), radius*sin(k/nel*ang*2*pi/360), 0)
    end
    clampedn = [1]
    tipn = [nel+1];

    tolerance = radius/nel/100;

    # Material properties
    material = MatDeforElastIso(DeforModelRed3D, 0.0, E, nu, 0.0)

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    for i in [1, 2, 3, 4, 5, 6]
        setebc!(dchi, clampedn, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi);

    # Additional fields
    u1 = deepcopy(u0)
    Rfield1 = deepcopy(Rfield0)
    rhs = gathersysvec(dchi);
    TMPv = deepcopy(rhs)
    utol = 1e-13*dchi.nfreedofs;

    # tbox = plot_space_box([[0.0*radius 0.0*radius 0.0*radius]; [1.15*radius 0.75*radius 2.0*radius]])
    # tshape0 = plot_solid(fens, fes; x = geom0.values, u = 0.0.*dchi.values[:, 1:3], R = Rfield0.values, facecolor = "rgb(125, 155, 125)", opacity = 0.3);
    # plots = cat(tbox,  tshape0; dims = 1)
    # pl = render(plots)
    # sleep(0.5)

    femm = FEMMCorotBeam(IntegDomain(fes, GaussRule(1, 2)), material)
    loadbdry = FESetP1(reshape(tipn, 1, 1))
    lfemm = FEMMBase(IntegDomain(loadbdry, PointRule()))

    load_parameters = 0.0:0.05:1;
    tip_displacement = fill(0.0, length(load_parameters), 3);
    step = 0
    for load_parameter in load_parameters
        applyebc!(dchi) # Apply boundary conditions
        u1.values[:] = u0.values[:]; # guess
        Rfield1.values[:] = Rfield0.values[:]; # guess
        fi = ForceIntensity(FFlt[0, 0, load_parameter*Fmag, 0, 0, 0]);

        println("Load: $load_parameter")
        iter = 1;
        while true
            F = distribloads(lfemm, geom0, dchi, fi, 3);
            Fr = restoringforce(femm, geom0, u1, Rfield1, dchi);       # Internal forces
            @. rhs = F + Fr;
            K = stiffness(femm, geom0, u1, Rfield1, dchi) + geostiffness(femm, geom0, u1, Rfield1, dchi);
            dchi = scattersysvec!(dchi, (K)\rhs); # Disp. incr
            u1.values[:] += (dchi.values[:,1:3])[:];   # increment displacement
            update_rotation_field!(Rfield1, dchi)
            print("$iter: ||du||=$(maximum(abs.(dchi.values[:])))\n")
            if maximum(abs.(dchi.values[:])) < utol # convergence check
                break;
            end
            if (iter > maxit)# bailout for failed convergence
                error("Possible failed convergence");
            end
            iter += 1;
        end
        u0.values[:] = u1.values[:];       # update the displacement
        Rfield0.values[:] = Rfield1.values[:]; # update the rotations

        step = step + 1
        tip_displacement[step, :] .= u1.values[tipn[1], :]

    end
    step = step + 1

    plots = cat(
    scatter(;x=tip_displacement[:,1], y=load_parameters.*Fmag, mode="lines", line_color = "rgb(0, 0, 0)", name="Tip-x"),
    scatter(;x=tip_displacement[:,2], y=load_parameters.*Fmag, mode="lines", line_color = "rgb(0, 0, 0)", name="Tip-y"),
    scatter(;x=tip_displacement[:,3], y=load_parameters.*Fmag, mode="lines", line_color = "rgb(0, 0, 0)", name="Tip-z"),
    scatter(;x=izzelna_ux[:,1], y=izzelna_ux[:,2]/(thinness_factor^4), mode="markers", line_color = "rgb(255, 0, 0)", name="Ref-x"),
    scatter(;x=izzelna_uy[:,1], y=izzelna_uy[:,2]/(thinness_factor^4), mode="markers", line_color = "rgb(0, 255, 0)", name="Ref-y"),
    scatter(;x=izzelna_uz[:,1], y=izzelna_uz[:,2]/(thinness_factor^4), mode="markers", line_color = "rgb(0, 0, 255)", name="Ref-z");
    dims = 1)
    layout = Layout(width=700, height=500, xaxis=attr(title="Tip displacement [in]"),
    yaxis=attr(title="Load [lb]"))
    pl = plot(plots, layout)
    display(pl)

    return true
end # curved_cantilever_thin

function allrun()
    println("#####################################################")
    println("# curved_cantilever ")
    curved_cantilever()
    println("#####################################################")
    println("# curved_cantilever_thin ")
    curved_cantilever_thin()
    return true
end # function allrun

@info "All examples may be executed with "
println("using .$(@__MODULE__); $(@__MODULE__).allrun()")

end # module
nothing

ext: jl


text:
using ApproxInferenceProblems
using Test


text:
module day08

include("inputs.jl")

ex = inputs.ex1_d8()
input = inputs.input_d8()

function parse_tape(input)
    map(input) do line
        op, arg = split(line)
        (Symbol(op), parse(Int, arg))
    end
end

function run(tape)
    accum = 0
    lines_visited = Set{Int}()
    ptr = 1

    while ptr ∉ lines_visited
        ptr > length(tape) && return (true, accum)

        op, arg = tape[ptr]
        push!(lines_visited, ptr)

        if op == :nop
            ptr += 1
        elseif op == :acc
            accum += arg
            ptr += 1
        elseif op == :jmp
            ptr += arg
        end
    end
    return (false, accum)
end

# Part 1
part1(input) = run(parse_tape(input))[2]

@assert part1(ex) == 5
@assert part1(input) == 1930
@info "Day 8, Part 1 answer: $(part1(input))"

# Part 2
function part2(input)
    tape = parse_tape(input)
    for i in eachindex(tape)
        op, arg = tape[i]
        if op == :jmp || op == :nop
            tape[i] = op == :jmp ? (:nop, arg) : (:jmp, arg)
            (was_terminated, accum) = run(tape)
            if was_terminated
                return accum
            else
                tape[i] = (op, arg)
            end
        end
    end
end

@assert part2(ex) == 8
@assert part2(input) == 1688
@info "Day 8, Part 2 answer: $(part2(input))"
end


text:
using Plots, LaTeXStrings, Statistics

function Deposition(;len, tot_time)
    surf = [0 for i=1:len]
    AvgList = [0.0 for i=1:tot_time]
    for n in 2:tot_time+1
        randsurf = rand(1:len)
        for i in randsurf
            index = FindLeast(surf, i, len)
            surf[index] += 1
        end
        AvgList[n-1] = mean(surf)
    end
    return AvgList
end

function sides(n, L)
    if n == L
        return n-1 , 1
    elseif n == 1
        return L , n+1
    else
        return n-1, n+1
    end
end

function FindLeast(surface, index_,L_surf)
    i1 , i2 = sides(index_, L_surf)
    lens = Dict(surface[i1]=>i1,
        surface[index_]=>index_,
        surface[i2]=>i2)
    minlen = min(surface[i1],surface[index_],surface[i2])
    return lens[minlen]
end

function Linear_fit(;Time, AvgList, tot_time)
    A = [hcat(Time[1:tot_time]) reshape(ones(tot_time), tot_time, 1)]
    b = reshape(AvgList[1:tot_time], tot_time, 1)
    line = (A \ b)
    return line
end

iternum = 1000
Parameters = Dict(:len => 300,
                    :tot_time => 50,)
allAvg = [ [0.0 for i in 1:Parameters[:tot_time]] for j = 1:iternum]
meanAvg = [0.0 for i in 1:Parameters[:tot_time]]
vars = [0.0 for i in 1:Parameters[:tot_time]]
for i in 1:iternum
    AvgList = Deposition(;Parameters...)
    allAvg[i] = AvgList
    meanAvg += AvgList
    print("\r$i")
end
meanAvg /= iternum
for i in 1:Parameters[:tot_time]
    vars[i] = std(hcat(allAvg...)[i,:])
end

Time = 0:Parameters[:tot_time]-1

Paraline = Dict(
                :Time => Time[1:end],
                :AvgList => meanAvg[1:end],
                :tot_time => 30)
Line = Linear_fit(;Paraline...)
X = Time[1]:Time[50]
Y = X .* Line[1] .+ Line[2]
Line[2]
plot(X,Y,c=:black,label = L"y = %$(round(Line[1],digits= 3))x + %$(round(Line[2],digits= 3))")
scatter!(Time[1:end], meanAvg[1:end],
    # xlims = (1, Parameters[:tot_time]),
    c = :steelblue,
    xlabel= L"Time",
    ylabel= L"H_{(t)}",
    title= L"~H_{(t)}-Time~\ (L = %$(Parameters[:len]))",
    label = L"Data\ point",
    yerror = vars)
    # legend = nothing)

savefig("../../Figs/Q1/H-t(L=$(Parameters[:len])).pdf")


text:
include("../src/SFFM.jl")
using LinearAlgebra, Plots, JLD2, SparseArrays

T = [-1.0 1.0 0; 0.5 -1.0 0.5; 0.5 0.5 -1]
C = [2.0; -2.0; -1.0]

r = (
    r = function (x)
        [-2*(x.>-1) 2*(x.>-1) -1*(x.>-1)] # [ones(size(x)) ones(size(x))]#
    end,
    R = function (x)
        [-2*(x.>-1).*x 2*(x.>-1).*x -1*(x.>-1).*x] # [x x]
    end,
)

Bounds = [0 6; -Inf Inf]
model = SFFM.Model(T = T, C = C, r = r, Bounds = Bounds)

Δ = 2
Nodes = collect(Bounds[1,1]:Δ:Bounds[1,2])
NBases = 3
Basis = "lagrange"
mesh = SFFM.MakeMesh(model = model, Nodes = Nodes, NBases = NBases, Basis=Basis)

V = SFFM.vandermonde(NBases=NBases)
Matrices2 = SFFM.MakeMatrices2(model=model, mesh=mesh)
MatricesR = SFFM.MakeMatricesR(model=model, mesh=mesh)
MR = SparseArrays.spzeros(Float64, model.NPhases * mesh.TotalNBases, model.NPhases * mesh.TotalNBases)
Minv =
    SparseArrays.spzeros(Float64, model.NPhases * mesh.TotalNBases, model.NPhases * mesh.TotalNBases)
FGR = SparseArrays.spzeros(Float64, model.NPhases * mesh.TotalNBases, model.NPhases * mesh.TotalNBases)
for i in 1:model.NPhases
    idx = (1:mesh.TotalNBases) .+ (i-1)*mesh.TotalNBases
    MR[idx,idx] = MatricesR.Global.M[i]
    Minv[idx,idx] = Matrices2.Global.MInv
    FGR[idx,idx] = model.C[i]*(MatricesR.Global.F[i]+MatricesR.Global.G[i])
end

T = SparseArrays.kron(model.T, SparseArrays.sparse(LinearAlgebra.I,mesh.TotalNBases,mesh.TotalNBases))
DR = MR * T * Minv + FGR * Minv

All = SFFM.MakeAll(model = model, mesh = mesh, approxType = "interpolation")
Matrices = SFFM.MakeMatrices(model=model,mesh=mesh,probTransform=false)
MatricesR = SFFM.MakeMatricesR(model=model,mesh=mesh)
B = SFFM.MakeB(model=model,mesh=mesh,Matrices=Matrices,probTransform=false)
Dr = SFFM.MakeDR(
    Matrices=Matrices,
    MatricesR=MatricesR,
    model=model,
    mesh=mesh,
    B=B,
)

Db = All.R.R[3:end-1,3:end-1]*B.B[3:end-1,3:end-1]

max_stars_repo_name: angus-lewis/SFFM

max_issues_repo_name: angus-lewis/SFFM

max_forks_repo_name: angus-lewis/SFFM


text:
export free_cells, free_string
export add_check, add_card, rm_card

"""
`free_cells` is a set of (up to) four cards
"""
free_cells = Set{Card}

function free_string(FC::free_cells)
    cards = sort(collect(FC))
    nc = length(cards)
    result = ""

    for k = 1:nc
        result *= string(cards[k]) * " "
    end
    for k = nc+1:4
        result *= "-- "
    end
    return result
end

"""
`add_card(FC::free_cells, PC::Card)::Bool` adds the card `PC`
to the free cells, returning `true` is ok.
"""
function add_card(FC::free_cells, PC::Card)::Bool
    if length(FC) > 3
        return false
    end
    push!(FC, PC)
    return true
end

"""
`rm_card(FC::free_cells, PC::Card)::Bool` removes the card `PC`
from the free cells `FC`, returning `true` if successful.
"""
function rm_card(FC::free_cells, PC::Card)::Bool
    if PC ∉ FC
        return false
    end
    delete!(FC, PC)
    return true
end

"""
`add_check(FC::free_cells,PC::Card)::Bool` checks if there's room
in the freecells to add the card.
"""
function add_check(c::free_cells, ::Card)::Bool
    return length(c) < 4
end

"""
`score(c::free_cells)` measures gives a penalty for too many
occupied free cells.
"""
function score(c::free_cells)
    svec = [10 10 0 -20 -50 -500]
    return svec[length(c)+1]
end


text:
#
# @description Implements the CBC (Cipher-Block Chaining) mode of operation algorithm.
#
# @author Alejandro Claro (alejandro.claro@gmail.com)
#
# Copyright 2015 All rights reserved.
# Use of this source code is governed by a MIT-style license that can be found in the LICENSE file.
#
export CbcModeCipher
export block_size, key_size, iv_size, iv, reset!, encrypt_next_blocks!, encrypt_last_block!, decrypt_next_blocks!, decrypt_last_block!

# @description Defines the CBC cipher mode data struture.
struct CbcModeCipher <: MultiBlockCipher
  block_cipher::BlockCipher
  padder::PaddingMethod
  iv::Vector{UInt8}

  # Operation state.
  previous_block::Vector{UInt8}
  partial_block::Vector{UInt8}

  # @description Initialize the CBC cipher with a random initialization vector.
  #
  # @param {BlockCipher}   block_cipher The block cipher.
  # @param {PaddingMethod} padder       The padding method.
  function CbcModeCipher(block_cipher::BlockCipher, padder::PaddingMethod)
    # TODO: Generate random IV using cryptographic secure RNG.
    iv = rand(UInt8, block_size(block_cipher))
    return reset!(new(block_cipher, padder, iv, UInt8[], UInt8[]))
  end

  # @description Initialize the CBC cipher.
  #
  # @param {BlockCipher}   block_cipher The block cipher.
  # @param {PaddingMethod} padder       The padding method.
  # @param {Vector{UInt8}} iv           The initialization vector.
  function CbcModeCipher(block_cipher::BlockCipher, padder::PaddingMethod, iv::Vector{UInt8})
    if length(iv) != block_size(block_cipher)
      throw(ArgumentError("Invalid initialization vector length. $(length(iv) * 8) bits recived, but $(block_size(block_cipher) * 8) bits expected."))
    end

    return reset!(new(block_cipher, padder, iv, UInt8[], UInt8[]))
  end
end

# @description Gets the cipher block size.
#
# @param {CbcModeCipher} self The cipher data struture.
#
# @return {Integer} The block size in bytes.
function block_size(self::CbcModeCipher)
  return block_size(self.block_cipher)
end

# @description Gets the cipher key size.
#
# @param {CbcModeCipher} self The cipher data struture.
#
# @return {Integer} The key size in bytes.
function key_size(self::CbcModeCipher)
  return key_size(self.block_cipher)
end

# @description Gets the cipher initialization vector size.
#
# @param {CbcModeCipher} self The cipher data struture.
#
# @return {Integer} The initialization vector size in bytes.
function iv_size(self::CbcModeCipher)
  return block_size(self.block_cipher)
end

# @description Gets the cipher initialization vector.
#
# @param {CbcModeCipher} self The cipher data struture.
#
# @return {Vector{UInt8}} The initialization vector.
function iv(self::CbcModeCipher)
  return self.iv
end

# @description Resets the state of the CBC data structure.
#
# @param {CbcModeCipher} self The hash algorithm.
#
# @return The reference to the resetted cipher data structure.
function reset!(self::CbcModeCipher)
  empty!(self.partial_block)

  resize!(self.previous_block, block_size(self))
  map!(x -> x, self.previous_block, self.iv)

  return self
end

# @description Continues a multiple-part encryption operation processing another data part.
#
# @param {CbcModeCipher} self The cipher data struture.
# @param {Vector{UInt8}} data The message part to encrypt.
#
# @return {Vector{UInt8}} The resulting ciphertext blocks. It can be a empty vector if the input data is too short to
# result in a new block.
function encrypt_next_blocks!(self::CbcModeCipher, data::Vector{UInt8})
  data   = vcat(self.partial_block, data)
  step   = block_size(self)
  len    = length(data)
  result = UInt8[]

  for i in 1 : step : (len - 1)
    xor        = map((x,y) -> x ⊻ y, self.previous_block, data[i : (i + step - 1)])
    ciphertext = encrypt(self.block_cipher, xor)

    map!(x -> x, self.previous_block, ciphertext)
    append!(result, ciphertext)
  end

  empty!(self.partial_block)
  append!(self.partial_block, data[(end - (len % step) + 1) : end])

  return result
end

# @description Encrypts the given last plaintext block.
#
# @param {CbcModeCipher} self  The cipher data struture.
# @param {Vector{UInt8}} block The last block from the message to encrypt.
#
# @return {Vector{UInt8}} The resulting last ciphertext blocks. This could be more than the size of a block.
function encrypt_last_block!(self::CbcModeCipher, block::Vector{UInt8})
  result = encrypt_next_blocks!(self, block)

  padded_block = pad(self.padder, self.partial_block, block_size(self))
  xor          = map((x,y) -> x ⊻ y, self.previous_block, padded_block)
  ciphertext   = encrypt(self.block_cipher, xor)

  append!(result, ciphertext)
  reset!(self)

  return result
end

# @description Continues a multiple-part decryption operation processing another data part.
#
# @param {CbcModeCipher} self  The cipher data struture.
# @param {Vector{UInt8}} data  The message part to decrypt.
#
# @return {Vector{UInt8}} The resulting plaintext blocks. It can be a empty vector if the input data is too short to
# result in a new block.
function decrypt_next_blocks!(self::CbcModeCipher, data::Vector{UInt8})
  data   = vcat(self.partial_block, data)
  step   = block_size(self)
  len    = length(data)
  result = UInt8[]

  for i in 1 : step : (len - 1)
    ciphertext = data[i : (i + step - 1)]
    plaintext  = decrypt(self.block_cipher, ciphertext)
    xor        = map((x,y) -> x ⊻ y, self.previous_block, plaintext)

    map!(x -> x, self.previous_block, ciphertext)
    append!(result, xor)
  end

  empty!(self.partial_block)
  append!(self.partial_block, data[(end - (len % step) + 1) : end])

  return result
end

# @description Decrypts the given last ciphertext block.
#
# @param {CbcModeCipher} self  The cipher data struture.
# @param {Vector{UInt8}} block The last block from the message to decrypt.
#
# @return {Vector{UInt8}} The resulting last plaintext block. This could be less than the size of a block.
function decrypt_last_block!(self::CbcModeCipher, block::Vector{UInt8})
  result = decrypt_next_blocks!(self, block)

  if isempty(self.partial_block) == false
    throw(BlockSizeError())
  end

  result = unpad(self.padder, result)
  reset!(self)

  return result
end


text:
a1 = [1,2,3,4,5,6]
println(a1)
res = a1 .* 2.143
println(res)

f(x) = 3x
res = f.(a1)

max_stars_repo_name: rpillar/JuliaExamples

max_issues_repo_name: rpillar/JuliaExamples

max_forks_repo_name: rpillar/JuliaExamples


text:
using JuMP, Ipopt
# using KNITRO
using SparseArrays # ybus
# using DataFrames # for comparison

include(string(pwd(), "/src/struct_network.jl"))
include(string(pwd(), "/src/readnetwork.jl"))
include(string(pwd(), "/src/ybus.jl"))
#----------------------------------------------------------------------------------#
solver = "Ipopt"
# solver = "KNITRO"

testsystem = "case_ieee30"
buses, lines, generators, datamat = readnetwork(testsystem)
#----------------------------------------------------------------------------------#

lineset = 1:length(lines)
busset = 1:length(buses)
genset = 1:length(generators)

B_g = []
for g in genset
    push!(B_g, generators[g].location)
end

B_gn = Array{Array{Int64,1}}(undef, length(buses))
for ii in 1:length(buses)
    B_gn[ii] = Int64[]
end
for g in genset
    push!(B_gn[generators[g].location], g)
end

Ybus, yff, yft, ytf, ytt = ybus(buses, lines)

yff_r = real.(yff)
yff_i = imag.(yff)
ytt_r = real.(ytt)
ytt_i = imag.(ytt)
yft_r = real.(yft)
yft_i = imag.(yft)
ytf_r = real.(ytf)
ytf_i = imag.(ytf)

#---------------------------------------------------------------------------------#
if solver == "Ipopt"
    m = Model(Ipopt.Optimizer)
    set_optimizer_attribute(m, "max_iter", 1_000_000)
    # set_optimizer_attribute(m, "linear_solver", "ma86")
    # set_optimizer_attribute(m, "print_level", 5) # default = 5 (0~12)
elseif solver == "KNITRO"
    m = Model(KNITRO.Optimizer)
    # set_optimizer_attribute(m, "option_file", "knitro.opt")
    #    set_optimizer_attribute(m, "ms_enable", 0)
    #    set_optimizer_attribute(m, "ms_maxsolves", 300)
    #    set_optimizer_attribute(m, "ms_terminate", 1) # 0 = maxsolve | 1 = first local optimum
    #    set_optimizer_attribute(m, "par_blasnumthreads", 10)
    #    set_optimizer_attribute(m, "maxit", 100000)
    # set_optimizer_attribute(m, "algorithm", 5)
    # set_optimizer_attribute(m, "convex", 1)
end

@variables(m, begin
    buses[b].Vmin <= v[b in busset] <= buses[b].Vmax, (start = 1.0)
    - π <= θ[b in busset] <= π
    generators[g].Pmin <= pg[g in genset] <= generators[g].Pmax  # real power output of generator
    generators[g].Qmin <= qg[g in genset] <= generators[g].Qmax # reactive power output of generator
end)

for b in busset
    if buses[b].btype == 3
        fix(θ[b], 0; force = true) # θ[1] = 0 at root node
        # fix(v[b], 1; force=true)
    end
end

@NLexpression(m, p_ft[l in lineset],
    +v[lines[l].fbus]^2 * yff_r[l]
    +
    v[lines[l].fbus] * v[lines[l].tbus] * (
        +yft_r[l] * cos(θ[lines[l].fbus] - θ[lines[l].tbus])
        +
        yft_i[l] * sin(θ[lines[l].fbus] - θ[lines[l].tbus])
    )
)

@NLexpression(m, q_ft[l in lineset],
    -v[lines[l].fbus]^2 * yff_i[l]
    +
    v[lines[l].fbus] * v[lines[l].tbus] * (
        -yft_i[l] * cos(θ[lines[l].fbus] - θ[lines[l].tbus])
        +
        yft_r[l] * sin(θ[lines[l].fbus] - θ[lines[l].tbus])
    )
)

@NLexpression(m, p_tf[l in lineset],
    +v[lines[l].tbus]^2 * ytt_r[l]
    +
    v[lines[l].tbus] * v[lines[l].fbus] * (
        +ytf_r[l] * cos(θ[lines[l].tbus] - θ[lines[l].fbus])
        +
        ytf_i[l] * sin(θ[lines[l].tbus] - θ[lines[l].fbus])
    )
)

@NLexpression(m, q_tf[l in lineset],
    -v[lines[l].tbus]^2 * ytt_i[l]
    +
    v[lines[l].tbus] * v[lines[l].fbus] * (
        -ytf_i[l] * cos(θ[lines[l].tbus] - θ[lines[l].fbus])
        +
        ytf_r[l] * sin(θ[lines[l].tbus] - θ[lines[l].fbus])
    )
)

@NLconstraint(m, PBalance[b in busset],
    sum(p_ft[l] for l in buses[b].outline)
    + sum(p_tf[l] for l in buses[b].inline)
    +
 - sum(pg[g] for g in B_gn[b])
    + buses[b].Pd
    + buses[b].Gs * v[b]^2
    ==
    0)

@NLconstraint(m, QBalance[b in busset],
    sum(q_ft[l] for l in buses[b].outline)
    + sum(q_tf[l] for l in buses[b].inline)
    +
 - sum(qg[g] for g in B_gn[b])
    +
    buses[b].Qd
    -
    buses[b].Bs * v[b]^2
    ==
    0)

@NLconstraint(m, LineCapFW[l in lineset; lines[l].u != 0], p_ft[l]^2 + q_ft[l]^2 <= (lines[l].u)^2)
@NLconstraint(m, LineCapBW[l in lineset; lines[l].u != 0], p_tf[l]^2 + q_tf[l]^2 <= (lines[l].u)^2)

@constraint(m, AngDiffmax[l in lineset], θ[lines[l].fbus] - θ[lines[l].tbus] <= min(lines[l].angmax, deg2rad(60)))
@constraint(m, AngDiffmin[l in lineset], θ[lines[l].fbus] - θ[lines[l].tbus] >= max(lines[l].angmin, -deg2rad(60)))

@expression(m, OCgen[g in genset],
    generators[g].cost[1] * (pg[g] * datamat["baseMVA"])^2
    + generators[g].cost[2] * (pg[g] * datamat["baseMVA"])
    + generators[g].cost[3])
@expression(m, OC, sum(OCgen[g] for g in genset))

@objective(m, Min, OC)

optimize!(m)
@show termination_status(m)
@show objective_value(m)

max_stars_repo_name: JipKim/ACOPF

max_issues_repo_name: JipKim/ACOPF

max_forks_repo_name: JipKim/ACOPF


text:
@testset "denman_beavers" begin
    A = [0.1 0.3; 0.2 0.7]
    (graph, cref) = graph_denman_beavers(10)
    @test sqrt(A) ≈ eval_graph(graph, A)
end


text:
module Reflected

using Test
using ExproniconLite
using Configurations
using OrderedCollections
using Configurations: @option, Reflect, has_duplicated_reflect_type

@option struct OptionA
    type::Reflect
    name::String
    age::Int
end

@option struct OptionB
    type::Reflect
    name::String
end

@option struct OptionC{T}
    type::Reflect
    typevar::T
end

# NOTE: this is not supported yet
# @option "aaa" struct OptionD{T}
#     type::Reflect
#     typevar::T
# end

@option "bbb" struct OptionE
    type::Reflect
    name::String
end

@option struct Composite
    person::Union{OptionA,OptionB,OptionC,OptionE}
end

@option struct VectorOfUnion
    data::Vector{Union{OptionA,OptionB}}
end

@testset "Reflect Type" begin
    @testset "OptionA" begin
        opt = OptionA(Reflect(), "Sam", 2)
        @test opt == OptionA(; name="Sam", age=2)
        d = to_dict(opt)
        @test from_dict(OptionA, d) == opt
    end

    @testset "OptionC{T}" begin
        opt = OptionC(Reflect(), 2)
        @test opt == OptionC(; typevar=2)
        d = to_dict(opt)

        @test_throws ArgumentError from_dict(OptionC, d)
        @test_throws ArgumentError from_dict(OptionC{Float32}, d)
        @test from_dict(OptionC{Int}, d) == OptionC(Reflect(), 2)
    end

    @testset "Composite" begin
        opt = Composite(OptionA(Reflect(), "Sam", 2))

        d = to_dict(opt)
        @test from_dict(Composite, d) == opt

        opt = Composite(OptionB(Reflect(), "Sam"))
        d = to_dict(opt)
        @test from_dict(Composite, d) == opt

        opt = Composite(OptionC(Reflect(), "Sam"))
        d = to_dict(opt)
        @test from_dict(Composite, d) == opt

        d = Dict(
            "person" => Dict{String,Any}(
                "type" => "Main.Reflected.OptionC{Float32}", "typevar" => 2
            ),
        )
        opt = from_dict(Reflected.Composite, d)
        @test typeof(opt.person) <: Reflected.OptionC{Float32}
    end

    @testset "duplicated Reflect" begin
        ex = @expr struct OptionError
            type_a::Reflect
            type_b::Reflect
            age::Int
        end

        @test_throws ArgumentError Configurations.option_m(Reflected, ex)

        ex = @expr JLKwStruct struct OptionError
            type_a::Reflect
            type_b::Reflect
            age::Int
        end
        @test has_duplicated_reflect_type(Reflected, ex)

        ex = @expr JLKwStruct struct OptionB end
        @test !has_duplicated_reflect_type(Reflected, ex)
    end

    @testset "type_alias as reflection" begin
        opt = Reflected.Composite(Reflected.OptionE(; name="aaa"))
        d = to_dict(opt)
        @test d == OrderedDict{String,Any}(
            "person" => OrderedDict{String,Any}("type" => "bbb", "name" => "aaa")
        )
        @test from_dict(Reflected.Composite, d) == opt
    end

    @testset "vector of union" begin
        d = Dict{String,Any}(
            "data" => [
                Dict{String,Any}(
                    "type" => "Main.Reflected.OptionA", "name" => "A", "age" => 1
                ),
                Dict{String,Any}(
                    "type" => "Main.Reflected.OptionA", "name" => "B", "age" => 2
                ),
                Dict{String,Any}("type" => "Main.Reflected.OptionB", "name" => "C"),
                Dict{String,Any}("type" => "Main.Reflected.OptionB", "name" => "D"),
            ],
        )

        @test from_dict(Reflected.VectorOfUnion, d) == Reflected.VectorOfUnion(
            Union{Reflected.OptionA,Reflected.OptionB}[
                Reflected.OptionA(Reflect(), "A", 1),
                Reflected.OptionA(Reflect(), "B", 2),
                Reflected.OptionB(Reflect(), "C"),
                Reflected.OptionB(Reflect(), "D"),
            ],
        )

        opt = Reflected.VectorOfUnion(;
            data=[
                Reflected.OptionA(; name="A", age=1),
                Reflected.OptionA(; name="B", age=2),
                Reflected.OptionB(; name="C"),
                Reflected.OptionB(; name="D"),
            ],
        )
        @test from_dict(Reflected.VectorOfUnion, to_dict(opt)) == opt
    end
end

end # Reflected


text:
export Gaussian, prod!, convert

# Convert parameterizations
function convert(::Type{ProbabilityDistribution{V, GaussianMeanPrecision}}, dist::ProbabilityDistribution{V, GaussianMeanVariance}) where V<:VariateType
    w = cholinv(dist.params[:v])
    m = deepcopy(dist.params[:m])

    return ProbabilityDistribution(V, GaussianMeanPrecision, m=m, w=w)
end

function convert(::Type{ProbabilityDistribution{V, GaussianMeanPrecision}}, dist::ProbabilityDistribution{V, GaussianWeightedMeanPrecision}) where V<:VariateType
    w = deepcopy(dist.params[:w])
    m = cholinv(w)*dist.params[:xi]

    return ProbabilityDistribution(V, GaussianMeanPrecision, m=m, w=w)
end

function convert(::Type{ProbabilityDistribution{V, GaussianMeanVariance}}, dist::ProbabilityDistribution{V, GaussianMeanPrecision}) where V<:VariateType
    v = cholinv(dist.params[:w])
    m = deepcopy(dist.params[:m])

    return ProbabilityDistribution(V, GaussianMeanVariance, m=m, v=v)
end

function convert(::Type{ProbabilityDistribution{V, GaussianMeanVariance}}, dist::ProbabilityDistribution{V, GaussianWeightedMeanPrecision}) where V<:VariateType
    v = cholinv(dist.params[:w])
    m = v*dist.params[:xi]

    return ProbabilityDistribution(V, GaussianMeanVariance, m=m, v=v)
end

function convert(::Type{ProbabilityDistribution{V, GaussianWeightedMeanPrecision}}, dist::ProbabilityDistribution{V, GaussianMeanPrecision}) where V<:VariateType
    w = deepcopy(dist.params[:w])
    xi = w*dist.params[:m]

    return ProbabilityDistribution(V, GaussianWeightedMeanPrecision, xi=xi, w=w)
end

function convert(::Type{ProbabilityDistribution{V, GaussianWeightedMeanPrecision}}, dist::ProbabilityDistribution{V, GaussianMeanVariance}) where V<:VariateType
    w = cholinv(dist.params[:v])
    xi = w*dist.params[:m]

    return ProbabilityDistribution(V, GaussianWeightedMeanPrecision, xi=xi, w=w)
end

# Convert VariateTypes
convert(::Type{ProbabilityDistribution{Multivariate, GaussianMeanVariance}}, dist::ProbabilityDistribution{Univariate, GaussianMeanVariance}) =
    ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=[dist.params[:m]], v=mat(dist.params[:v]))
convert(::Type{ProbabilityDistribution{Multivariate, GaussianMeanPrecision}}, dist::ProbabilityDistribution{Univariate, GaussianMeanPrecision}) =
    ProbabilityDistribution(Multivariate, GaussianMeanPrecision, m=[dist.params[:m]], w=mat(dist.params[:w]))
convert(::Type{ProbabilityDistribution{Multivariate, GaussianWeightedMeanPrecision}}, dist::ProbabilityDistribution{Univariate, GaussianWeightedMeanPrecision}) =
    ProbabilityDistribution(Multivariate, GaussianWeightedMeanPrecision, xi=[dist.params[:xi]], w=mat(dist.params[:w]))

function prod!(
    x::ProbabilityDistribution{Univariate, <:Gaussian},
    y::ProbabilityDistribution{Univariate, <:Gaussian},
    z::ProbabilityDistribution{Univariate, GaussianWeightedMeanPrecision}=ProbabilityDistribution(Univariate, GaussianWeightedMeanPrecision, xi=0.0, w=1.0))

    z.params[:xi] = unsafeWeightedMean(x) + unsafeWeightedMean(y)
    z.params[:w] = unsafePrecision(x) + unsafePrecision(y)

    return z
end

@symmetrical function prod!(
    x::ProbabilityDistribution{Univariate, <:Gaussian},
    y::ProbabilityDistribution{Univariate, PointMass},
    z::ProbabilityDistribution{Univariate, PointMass}=ProbabilityDistribution(Univariate, PointMass, m=0.0))

    z.params[:m] = y.params[:m]
    return z
end

function prod!(
    x::ProbabilityDistribution{Multivariate, <:Gaussian},
    y::ProbabilityDistribution{Multivariate, <:Gaussian},
    z::ProbabilityDistribution{Multivariate, GaussianWeightedMeanPrecision}=ProbabilityDistribution(Multivariate, GaussianWeightedMeanPrecision, xi=[NaN], w=transpose([NaN])))

    z.params[:xi] = unsafeWeightedMean(x) + unsafeWeightedMean(y)
    z.params[:w] = unsafePrecision(x) + unsafePrecision(y)

    return z
end

@symmetrical function prod!(
    x::ProbabilityDistribution{Multivariate, <:Gaussian},
    y::ProbabilityDistribution{Multivariate, PointMass},
    z::ProbabilityDistribution{Multivariate, PointMass}=ProbabilityDistribution(Multivariate, PointMass, m=[NaN]))

    z.params[:m] = deepcopy(y.params[:m])

    return z
end

function sample(dist::ProbabilityDistribution{Univariate, <:Gaussian})
    isProper(dist) || error("Cannot sample from improper distribution")
    (m,v) = unsafeMeanCov(dist)
    return sqrt(v)*randn() + m
end

function sample(dist::ProbabilityDistribution{Univariate, <:Gaussian}, n_samples::Int64)
    isProper(dist) || error("Cannot sample from improper distribution")
    (m,v) = unsafeMeanCov(dist)

    return sqrt(v).*randn(n_samples) .+ m
end

function sample(dist::ProbabilityDistribution{Multivariate, <:Gaussian})
    isProper(dist) || error("Cannot sample from improper distribution")
    (m,V) = unsafeMeanCov(dist)
    return (cholesky(default_cholesky_mode, V)).U' *randn(dims(dist)) + m
end

function sample(dist::ProbabilityDistribution{Multivariate, <:Gaussian}, n_samples::Int64)
    isProper(dist) || error("Cannot sample from improper distribution")
    (m,V) = unsafeMeanCov(dist)
    U = (cholesky(default_cholesky_mode, V)).U
    d = dims(dist)

    return [U' *randn(d) + m for i in 1:n_samples]
end

function naturalParams(dist::ProbabilityDistribution{<:VariateType, <:Gaussian})
    (xi, w) = unsafeWeightedMeanPrecision(dist)
    return vcat(xi, -0.5*vec(w))
end

standardDistribution(::Type{Univariate}, ::Type{<:Gaussian}; η::Vector) = ProbabilityDistribution(Univariate, GaussianWeightedMeanPrecision, xi=η[1], w=-2*η[2])
function standardDistribution(::Type{Multivariate}, ::Type{<:Gaussian}; η::Vector)
    d = Int(-0.5 + 0.5*sqrt(1 + 4*length(η))) # Extract dimensionality
    η_1 = η[1:d]
    η_2 = reshape(η[d+1:end], d, d)
    return ProbabilityDistribution(Multivariate, GaussianWeightedMeanPrecision, xi=η_1, w=-2.0*η_2)
end

logNormalizer(::Type{Univariate}, ::Type{<:Gaussian}; η::Vector) = -η[1]^2/(4*η[2]) - 0.5*log(-2*η[2])
function logNormalizer(::Type{Multivariate}, ::Type{<:Gaussian}; η::Vector)
    d = Int(-0.5 + 0.5*sqrt(1 + 4*length(η))) # Extract dimensionality
    η_1 = η[1:d]
    η_2 = reshape(η[d+1:end], d, d)
    return η_1'*cholinv(-4*η_2)*η_1 - 0.5*logdet(-2*η_2)
end

logPdf(V::Type{Univariate}, ::Type{F}, x::Number; η::Vector) where F<:Gaussian = -0.5*log(2pi) + vcat(x, x^2)'*η - logNormalizer(V, F; η=η)
function logPdf(V::Type{Multivariate}, ::Type{F}, x::Vector; η::Vector) where F<:Gaussian
    d = Int(-0.5 + 0.5*sqrt(1 + 4*length(η))) # Extract dimensionality
    return -0.5*d*log(2pi) + vcat(x, vec(x*x'))'*η - logNormalizer(V, F; η=η)
end

# Entropy functional
function differentialEntropy(dist::ProbabilityDistribution{Univariate, <:Gaussian})
    return  0.5*log(unsafeCov(dist)) +
            0.5*log(2*pi) +
            0.5
end

function differentialEntropy(dist::ProbabilityDistribution{Multivariate, <:Gaussian})
    d = dims(dist)[1]
    return  0.5*logdet(unsafeCov(dist)) +
            (d/2)*log(2*pi) +
            d/2
end


text:
#https://discourse.julialang.org/t/encapsulating-enum-access-via-dot-syntax/11785/3
# @scopedenum Fruit APPLE=1 PEAR=2 BANANA=3
# Fruit.APPLE
# Fruit.PEAR
# # access enum value
# Fruit.APPLE.value
# # make an APPLE from string
# Fruit.Enum("APPLE")
# # restricting type signatures
# f(x::Fruit.Enum) = # do stuff w/ x
macro scopedenum(T, args...)
    defs = Expr(:block)
    append!(defs.args, collect(:(const $(x.args[1]) = Enum($(x.args[2]))) for x in args))
    names = Dict(x.args[2]=>String(x.args[1]) for x in args)
    str2val = Dict(String(x.args[1])=>x.args[2] for x in args)
    push!(defs.args, quote
        function name(e::Enum)
            nms = $names
            return nms[e.value]
        end
        Enum(str::String) = Enum($(str2val)[str])
        Base.show(io::IO, e::E) where {E <: Enum} = print(io, "$(Base.datatype_module(E)).$(name(e)) = $(e.value)")
    end)
    blk = esc(:(module $T; struct Enum{T}; value::T; end; Enum{T}(e::Enum{T}) where {T} = e; $defs; end))
    return Expr(:toplevel, blk)
end

#translated from auto-generated Thrift.jl to more ergonomic enum usage via macro
@scopedenum TDatumType SMALLINT=Int32(0) INT=Int32(1) BIGINT=Int32(2) FLOAT=Int32(3) DECIMAL=Int32(4) DOUBLE=Int32(5) STR=Int32(6) TIME=Int32(7) TIMESTAMP=Int32(8) DATE=Int32(9) BOOL=Int32(10) INTERVAL_DAY_TIME=Int32(11) INTERVAL_YEAR_MONTH=Int32(12) POINT=Int32(13) LINESTRING=Int32(14) POLYGON=Int32(15) MULTIPOLYGON=Int32(16) TINYINT=Int32(17) GEOMETRY=Int32(18) GEOGRAPHY=Int32(19)

#this function used to translate enum to Julia types
#if a Julia type isn't defined, mark as String
function getcolumntype(x::Int32)
    d = Dict(0 => Int16,  #SMALLINT
             1 => Int32,  #INT
             2 => Int64,  #BIGINT
             3 => Float32, #FLOAT
             4 => Dec64,  #DECIMAL
             5 => Float64, #DOUBLE
             6 => String,  #STR
             7 => Time,  #TIME
             8 => DateTime,  #TIMESTAMP
             9 => Date,  #DATE
             10 => Bool,  #BOOL
             11 => String,  #INTERVAL_DAY_TIME
             12 => String,  #INTERVAL_YEAR_MONTH
             13 => Point,  #POINT
             14 => LineString,  #LINESTRING
             15 => Polygon,  #POLYGON
             16 => MultiPolygon,  #MULTIPOLYGON
             17 => Int8,  #TINYINT
             18 => String,  #GEOMETRY
             19 => String  #GEOGRAPHY
    )

    return d[x]
end

@scopedenum TEncodingType NONE=Int32(0) FIXED=Int32(1) RL=Int32(2) DIFF=Int32(3) DICT=Int32(4) SPARSE=Int32(5) GEOINT=Int32(6) DATE_IN_DAYS=Int32(7)

function getencodingtype(x::Int32)
    d = Dict(0 => "None",
             1 => "Fixed",
             2 => "RL",
             3 => "Diff",
             4 => "Dict",
             5 => "Sparse",
             6 => "GeoInt"
    )

    d[x]
end

@scopedenum TExecuteMode GPU=Int32(1) CPU=Int32(2)

@scopedenum TDeviceType CPU=Int32(0) GPU=Int32(1)

@scopedenum TTableType DELIMITED=Int32(0) POLYGON=Int32(1)


text:
@testset "FastBSplineSpace" begin
    P1 = FastBSplineSpace(2, Knots([1, 3, 5, 6, 8, 9]))
    @test bsplinesupport(P1, 2) == 3..8
    @test dim(P1) == 3
    @test properdim(P1) == 3
    @test isproper(P1) == true

    P2 = FastBSplineSpace(1, Knots([1, 3, 3, 3, 8, 9]))
    @test isproper(P2) == false
    @test properdim(P2) == 3

    i = 2
    k = Knots([5, 12, 13, 13, 14])
    p = 2
    P = FastBSplineSpace(p, k)
    @test bsplinesupport(P) == [5..13, 12..14]
    @test bsplinesupport(P, i) == 12..14

    @test isproper(FastBSplineSpace(2, Knots([1, 3, 5, 6, 8, 9])))
    @test !isproper(FastBSplineSpace(1, Knots([1, 3, 3, 3, 8, 9])))

    @test dim(FastBSplineSpace(2, Knots([1, 3, 5, 6, 8, 9]))) == 3

    P1 = FastBSplineSpace(1, Knots([1, 3, 5, 8]))
    P2 = FastBSplineSpace(1, Knots([1, 3, 5, 6, 8, 9]))
    P3 = FastBSplineSpace(2, Knots([1, 1, 3, 3, 5, 5, 8, 8]))
    @test P1 ⊆ P2
    @test P1 ⊆ P3
    @test P2 ⊈ P3
end


text:
using Revise, BenchmarkTools, DataFrames
using Distributions, StatsPlots, FFTDists, Plots.Measures
include("Benchmarks.jl")
############################################################
#####################  Define Models  ######################
############################################################
model = Gamma(5, .05) + LogNormal(-1, .5)
N = 10 .^(0:7)
############################################################
#####################  Benchmark  Model  ###################
############################################################
results = benchmarker(model, N)
############################################################
#####################  Plot Model  #########################
############################################################
pyplot()
options = (norm=true, legend=false, ylabel="Mean Time (ms)", xlabel="N")
@df results plot(:N, :time, xaxis=:log; options...)

max_stars_repo_path: Benchmarks/Run Benchmarks.jl

max_issues_repo_path: Benchmarks/Run Benchmarks.jl

max_forks_repo_path: Benchmarks/Run Benchmarks.jl


text:
struct SSubHandler end

gethandler(broadcasted, ::Type{typeof(-)}, ::Type{<:SBitstreamLike}, ::Type{<:SBitstreamLike}) = SSubHandler()
init_state(::SSubHandler) = (id = 0,)

function (handler::SSubHandler)(buffer, netlist, state, inputs, outputs)
    # update netlist with inputs
    setsigned!(netlist, inputs[1], true)
    setsigned!(netlist, inputs[2], true)

    # compute output size
    lname, rname = handle_broadcast_name(name(inputs[1]), name(inputs[2]),
                                         netsize(inputs[1]), netsize(inputs[2]))
    outsize = netsize(outputs[1])

    # add output net to netlist
    setsigned!(netlist, outputs[1], true)

    # add internal nets to netlist
    push!(netlist, Net(name = "sub$(state.id)_out_pp", size = outsize))
    push!(netlist, Net(name = "sub$(state.id)_out_pm", size = outsize))
    push!(netlist, Net(name = "sub$(state.id)_out_mp", size = outsize))
    push!(netlist, Net(name = "sub$(state.id)_out_mm", size = outsize))

    write(buffer, """
        $stdcomment
        // BEGIN sub$(state.id)
        stoch_signed_sub_mat #(
                .NUM_ROWS($(outsize[1])),
                .NUM_COLS($(outsize[2]))
            ) sub$(state.id) (
                .CLK(CLK),
                .nRST(nRST),
                .A_p($(lname("_p"))),
                .A_m($(lname("_m"))),
                .B_p($(rname("_p"))),
                .B_m($(rname("_m"))),
                .Y_p($(name(outputs[1]))_p),
                .Y_m($(name(outputs[1]))_m)
            );
        // END sub$(state.id)
        \n""")

    return buffer, (id = state.id + 1,)
end


text:
using DifferentialEquations
using ParameterizedFunctions
using PyPlot

f = @ode_def_bare BallBounce begin
  dy = v
  dv = -g
end g=9.81

abstract type DiscreteTimeSystem end

mutable struct ControlSystem{T<:Number} <: DiscreteTimeSystem
  Ts::T
  Tprev::T
end

det = ControlSystem(0.15,0.0)

function condition(det::ControlSystem,t,u,integrator)
  t - (det.Tprev + det.Ts)
end

function affect!(integrator,det::ControlSystem)
  integrator.u[2] += 1
  det.Tprev = integrator.t
end

cb = ContinuousCallback((t, u, int) -> condition(det,t,u,int),
                          (int) -> affect!(int,det))

u0 = [50.0,-10.0]
tspan = (0.0,1.5)
prob = ODEProblem(f,u0,tspan)
sol = solve(prob,Tsit5(),callback=cb,force_dtmin=true)

idx = 1
v = Vector{Float64}(length(sol.u))
for x in sol.u
  v[idx] = x[2]
  idx += 1
end

#plot(sol.t, v)


text:
"takes a json biomd metadata file and converts to DF"
function jsonfn_to_df(fn)
    json = read(fn, String);
    json = JSON3.read(json)
    haskey(json, :models) ? DataFrame(jsontable(json.models)) : missing
end

# "https://www.ebi.ac.uk/biomodels/search?query=&offset=0&numResults=10"
"curls the json files"
function curl_biomd_metadata(query, meta_dir="$(datadir)/biomd_ode_meta")
    offsets = 0:100:2200 # ideally uses a while loop
    urls = .*("https://www.ebi.ac.uk/biomodels/search?query=", query,"&offset=", string.(offsets), "&numResults=100&format=json")
    @sync for i in 1:length(offsets)
        @async run(`curl $(urls[i]) -o "$(meta_dir)/sbml_$(i).json"`)
    end
    read_biomd_metadf(meta_dir)
end

"returns a dataframe with metadata about the biomodels"
function read_biomd_metadf(meta_dir)
    fns = readdir(meta_dir; join=true)
    dfs = filter(!ismissing, jsonfn_to_df.(fns))
    vcat(dfs...)
end

"this downloads the xmls directly, without needing zip "
function curl_biomd_xmls(ids; verbose=false)
    base = "https://www.ebi.ac.uk/biomodels/model/download/"
    @sync Threads.@threads for id in ids
        verbose && @info("downloading $id")
        url = "$(base)$(id)?filename=$(id)_url.xml"
        fn = "$(datadir)/biomd/$(id).xml"
        run(`curl $(url) -o "$(fn)"`)
    end
end

"creates a metadata dataframe from queries to the biomodels api and then downloads the actual models"
function biomodels(
    meta_dir="$(datadir)/biomd_meta",
    biomd_dir="$(datadir)/biomd/"; # the xml files are put here
    odes_only=false,
    curl_meta=false,
    limit=nothing, # still gets all the metadata, just limits for curling zips
    verbose=true
    )

    verbose && @info("in biomodels()")
    mkpath.([meta_dir, biomd_dir])
    query = odes_only ? "*%3A*%20AND%20modellingapproach%3A%22Ordinary%20differential%20equation%20model%22&domain=biomodels" : "sbml"
    df = curl_meta ? curl_biomd_metadata(query, meta_dir) : read_biomd_metadf(meta_dir)
    verbose && display(df)
    CSV.write("$(datadir)/sbml_biomodels.csv", df)
    limit === nothing ? curl_biomd_xmls(df.id; verbose=verbose) : curl_biomd_xmls(df.id[1:limit]; verbose=verbose)
    df
end

function get_biomd_fns()
    readdir(joinpath(datadir, "biomd/"); join=true)
end


text:
using NonUniformRandomVariateGeneration
using Distributions
using BenchmarkTools

function benchNative(n::Int64, ps::Vector{Float64}, N::Int64)
  m::Vector{Int64} = zeros(Int64, length(ps))
  v::Vector{Int64} = Vector{Int64}(undef, length(ps))
  for i in 1:N
    sampleMultinomial!(n, ps, v)
    m .+= v
  end
  return m / N / n
end

function benchDistributions(n::Int64, ps::Vector{Float64}, N::Int64)
  multinom = Distributions.Multinomial(n, ps)
  m::Vector{Int64} = zeros(Int64, length(ps))
  v::Vector{Int64} = Vector{Int64}(undef, length(ps))
  for i in 1:N
    Distributions._rand!(multinom, v)
    m .+= v
  end
  return m / N / n
end

n = 100000
ps = [0.5; 0.1; 0.25; 0.15]

@btime benchNative(n, ps, 1024*1024*8)
@btime benchDistributions(n, ps, 1024*1024*8)


text:
using SIGPROC
using Documenter

ENV["GKSwstype"] = "100"

DocMeta.setdocmeta!(SIGPROC, :DocTestSetup, :(using SIGPROC); recursive=true)

makedocs(;
         modules=[SIGPROC],
         authors="Kiran Shila <me@kiranshila.com> and contributors",
         repo="https://github.com/kiranshila/SIGPROC.jl/blob/{commit}{path}#{line}",
         sitename="SIGPROC.jl",
         format=Documenter.HTML(;
                                prettyurls=get(ENV, "CI", "false") == "true",
                                canonical="https://kiranshila.github.io/SIGPROC.jl",
                                assets=String[]),
         pages=["Home" => "index.md",
                "Examples" => "examples.md",
                "API" => "api.md"])

deploydocs(;
           repo="github.com/kiranshila/SIGPROC.jl",
           devbranch="main")


text:
"""
`SolutionODE`: Solution of an ordinary differential equation

Contains all fields necessary to store the solution of an ODE.

### Fields

* `nd`: dimension of the dynamical variable ``q``
* `nt`: number of time steps to store
* `ni`: number of initial conditions
* `t`:  time steps
* `q`:  solution `q[nd, nt+1, ni]` with `q[:,0,:]` the initial conditions
* `ntime`: number of time steps to compute
* `nsave`: store every nsave'th time step (default: 1)
* `nwrite`: save data to disk after every nwrite'th time step (default: ntime)
* `counter`: counter for copied solution entries
* `woffset`: counter for file offset
* `h5`: HDF5 file for storage
"""
mutable struct SolutionODE{dType, tType, N} <: DeterministicSolution{dType, tType, N}
    nd::Int
    nt::Int
    ni::Int
    t::TimeSeries{tType}
    q::SDataSeries{dType,N}
    ntime::Int
    nsave::Int
    nwrite::Int
    counter::Vector{Int}
    woffset::Int
    h5::HDF5File

    function SolutionODE{dType, tType, N}(nd, nt, ni, t, q, ntime, nsave, nwrite) where {dType <: Number, tType <: Real, N}
        new(nd, nt, ni, t, q, ntime, nsave, nwrite, zeros(Int, ni), 0)
    end
end

function SolutionODE(equation::Union{ODE{DT,TT,FT},SODE{DT,TT,FT}}, Δt::TT, ntime::Int, nsave::Int=DEFAULT_NSAVE, nwrite::Int=DEFAULT_NWRITE; filename=nothing) where {DT,TT,FT}
    @assert nsave > 0
    @assert ntime == 0 || ntime ≥ nsave
    @assert nwrite == 0 || nwrite ≥ nsave
    @assert mod(ntime, nsave) == 0

    if nwrite > 0
        @assert mod(nwrite, nsave) == 0
        @assert mod(ntime, nwrite) == 0
    end

    N  = (equation.n > 1 ? 3 : 2)
    nd = equation.d
    ni = equation.n
    nt = div(ntime, nsave)
    nt = (nwrite == 0 ? nt : div(nwrite, nsave))
    nw = (nwrite == 0 ? ntime : nwrite)

    @assert nd > 0
    @assert ni > 0
    @assert nt > 0
    @assert nw > 0

    t = TimeSeries{TT}(nt, Δt, nsave)
    q = SDataSeries(DT, nd, nt, ni)
    s = SolutionODE{DT,TT,N}(nd, nt, ni, t, q, ntime, nsave, nw)
    set_initial_conditions!(s, equation)

    if !isnothing(filename)
        isfile(filename) ? @warn("Overwriting existing HDF5 file.") : nothing
        create_hdf5(s, filename)
    end

    return s
end

function SolutionODE(t::TimeSeries{TT}, q::SDataSeries{DT,N}, ntime::Int) where {DT,TT,N}
    # extract parameters
    nd = q.nd
    ni = q.ni
    nt = t.n
    ns = div(ntime, nt)

    @assert mod(ntime, nt) == 0

    # create solution
    SolutionODE{DT,TT,N}(nd, nt, ni, t, q, ntime, ns, 0)
end

function SolutionODE(file::String)
    # open HDF5 file
    get_config(:verbosity) > 1 ? @info("Reading HDF5 file ", file) : nothing
    h5 = h5open(file, "r")

    # read attributes
    ntime = read(attrs(h5)["ntime"])
    nsave = read(attrs(h5)["nsave"])

    # reading data arrays
    t = TimeSeries(read(h5["t"]), nsave)
    q = SDataSeries(read(h5["q"]))

    # need to close the file
    close(h5)

    # create solution
    SolutionODE(t, q, ntime)
end

Base.:(==)(sol1::SolutionODE, sol2::SolutionODE) = (
                                sol1.nd == sol2.nd
                             && sol1.nt == sol2.nt
                             && sol1.ni == sol2.ni
                             && sol1.t  == sol2.t
                             && sol1.q  == sol2.q
                             && sol1.ntime == sol2.ntime
                             && sol1.nsave == sol2.nsave
                             && sol1.counter == sol2.counter)

hdf5(sol::SolutionODE)  = sol.h5
time(sol::SolutionODE)  = sol.t.t
ntime(sol::SolutionODE) = sol.ntime
nsave(sol::SolutionODE) = sol.nsave
offset(sol::SolutionODE) = sol.woffset

function set_initial_conditions!(sol::SolutionODE{DT,TT}, equ::Union{ODE{DT,TT},SODE{DT,TT}}) where {DT,TT}
    set_initial_conditions!(sol, equ.t₀, equ.q₀)
end

function set_initial_conditions!(sol::SolutionODE{DT,TT}, t₀::TT, q₀::Union{Array{DT}, Array{TwicePrecision{DT}}}) where {DT,TT}
    set_data!(sol.q, q₀, 0)
    compute_timeseries!(sol.t, t₀)
    sol.counter .= 1
end

function get_initial_conditions!(sol::SolutionODE{DT,TT}, q::SolutionVector{DT}, k, n=1) where {DT,TT}
    get_solution!(sol, q, n-1, k)
end

function get_initial_conditions(sol::SolutionODE, k, n=1)
    get_solution(sol, n-1, k)
end

function CommonFunctions.get_solution!(sol::SolutionODE{DT,TT}, q::SolutionVector{DT}, n, k) where {DT,TT}
    q .= sol.q[:, n, k]
end

function CommonFunctions.get_solution(sol::SolutionODE, n, k)
    (sol.t[n], sol.q[:, n, k])
end

function CommonFunctions.set_solution!(sol::SolutionODE, t, q, n, k)
    set_solution!(sol, q, n, k)
end

function CommonFunctions.set_solution!(sol::SolutionODE{DT,TT}, q::SolutionVector{DT}, n, k) where {DT,TT}
    @assert n <= sol.ntime
    @assert k <= sol.ni
    if mod(n, sol.nsave) == 0
        if sol.counter[k] > sol.nt
            @error("Solution overflow. Call write_to_hdf5() and reset!() before continuing the simulation.")
        end
        set_data!(sol.q, q, sol.counter[k], k)
        sol.counter[k] += 1
    end
end

function CommonFunctions.reset!(sol::SolutionODE)
    reset!(sol.q)
    compute_timeseries!(sol.t, sol.t[end])
    sol.counter .= 1
    sol.woffset += sol.nt
end

function Base.close(solution::SolutionODE)
    # close(solution.t)
    # close(solution.q)
    close(solution.h5)
end

"Creates HDF5 file and initialises datasets for ODE solution object."
function create_hdf5(solution::SolutionODE{DT,TT,2}, file::AbstractString) where {DT,TT}
    # create HDF5 file and save attributes and common parameters
    solution.h5 = createHDF5(solution, file)

    # save attributes
    save_attributes(solution)

    # create dataset
    nt = div(solution.ntime, solution.nsave)

    # nt can be used to set the expected total number of timesteps
    # so that the size of the array does not need to be adapted dynamically.
    # Right now, the size has to be set to nt as dynamical size adaptation is
    # not yet working.
    t = d_create(solution.h5, "t", datatype(TT), dataspace((nt+1,)), "chunk", (1,))
    q = d_create(solution.h5, "q", datatype(DT), dataspace(solution.nd, solution.nt+1), "chunk", (solution.nd,1))

    # copy initial conditions
    t[1] = solution.t[0]
    q[:, 1] = solution.q[:, 0]

    return solution.h5
end

"Creates HDF5 file and initialises datasets for ODE solution object."
function create_hdf5(solution::SolutionODE{DT,TT,3}, file::AbstractString) where {DT,TT}
    # create HDF5 file and save attributes and common parameters
    solution.h5 = createHDF5(solution, file)

    # create dataset
    nt = div(solution.ntime, solution.nsave)

    # nt can be used to set the expected total number of timesteps
    # so that the size of the array does not need to be adapted dynamically.
    # Right now, the size has to be set to nt as dynamical size adaptation is
    # not yet working.
    t = d_create(solution.h5, "t", datatype(TT), dataspace((nt+1,)), "chunk", (1,))
    q = d_create(solution.h5, "q", datatype(DT), dataspace(solution.nd, nt+1, solution.ni), "chunk", (solution.nd,1,1))

    # copy initial conditions
    t[1] = solution.t[0]
    q[:, 1, :] = solution.q[:, 0, :]

    return solution.h5
end

"Append solution to HDF5 file."
function CommonFunctions.write_to_hdf5(solution::SolutionODE{DT,TT,2}, h5::HDF5File, offset=0) where {DT,TT}
    # set convenience variables and compute ranges
    n  = solution.nt
    j1 = offset+2
    j2 = offset+1+n

    # TODO # extend dataset if necessary
    # if size(x, 2) < j2
    #     set_dims!(x, (d, j2))
    # end

    # copy data from solution to HDF5 dataset
    h5["t"][j1:j2] = solution.t[1:n]
    h5["q"][:, j1:j2] = solution.q[:, 1:n]

    return nothing
end

"Append solution to HDF5 file."
function CommonFunctions.write_to_hdf5(solution::SolutionODE{DT,TT,3}, h5::HDF5File, offset=0) where {DT,TT}
    # set convenience variables and compute ranges
    n  = solution.nt
    j1 = offset+2
    j2 = offset+1+n

    # TODO # extend dataset if necessary
    # if size(x, 2) < j2
    #     set_dims!(x, (d, j2))
    # end

    # copy data from solution to HDF5 dataset
    h5["t"][j1:j2] = solution.t[1:n]
    h5["q"][:, j1:j2, :] = solution.q[:, 1:n, :]

    return nothing
end


text:
using Weiqi

import Weiqi: empty, white, black, cb
import Weiqi: NewPosition, bp, wp, nextplayer, pass, neighbors, liberties

@testset "assign var `np` to NewPosition(player, coords, stone)" begin
newpos = NewPosition(bp, (1,2), black)
newpos = NewPosition(wp, (1,2), white)
# board state should decide when overwrite can happen

newpos = NewPosition(bp, (1,2), empty) # cannot happen
newpos = NewPosition(wp, (1,2), empty) # cannot happen

newpos.coords
newpos.stone
end

@testset "call `function nextplayer(np)`" begin
newpos.player
nextplayer(newpos)
end

@testset "set a pass, then unset it" begin
newpos = NewPosition(bp, (0,0), black)
pass(newpos)

newpos = NewPosition(bp, (1,1), black)
pass(newpos)
end

@testset "test `function neighbors(cb, row::Int64, col::Int64)`" begin
# assume cb = createboard(19)
neighbors(cb, 1, 1)
neighbors(cb, 19, 19)
neighbors(cb, 7, 13)
end

@testset "test `function liberties(cb, row::Int64, col::Int64)``" begin
cb.array[3:5,3:5] = black; # create black block
cb
liberties(cb, 4, 5) # black block
liberties(cb, 3, 3) # black block
liberties(cb, 5, 4) # black block
end


text:
(function (x::Int64)
  x = 10
  println(10)
  x = "Did something print?"
  println(x)
end)(10)

max_stars_repo_name: idrougge/ple

max_issues_repo_name: idrougge/ple

max_forks_repo_name: idrougge/ple


text:
using VectorBackedLists

include("test_lists.jl")


text:
"""
Computes the Bayesian score component for the given target variable index and
    Dirichlet prior counts given in alpha

INPUT:
    i       - index of the target variable
    parents - list of indeces of parent variables (should not contain self)
    r       - list of instantiation counts accessed by variable index
              r[1] gives number of discrete states variable 1 can take on
    data - matrix of sufficient statistics / counts
              d[j,k] gives the number of times the target variable took on its kth instantiation
              given the jth parental instantiation

OUTPUT:
    the Bayesian score, Float64
"""
function bayesian_score_component(
    i::Int,
    parents::AbstractVector{I},
    ncategories::AbstractVector{Int},
    data::AbstractMatrix{Int},
    alpha::AbstractMatrix{Float64}, # ncategories[i]×prod(ncategories[parents])
) where {I<:Integer}

    (n, m) = size(data)
    if !isempty(parents)
        Np = length(parents)
        stridevec = fill(1, Np)
        for k in 2:Np
            stridevec[k] = stridevec[k-1] * ncategories[parents[k-1]]
        end
        js = (data[parents,:] - 1)' * stridevec + 1
    else
        js = fill(1, m)
    end

    N = sparse(vec(data[i,:]), vec(js), 1, size(alpha)...) # note: duplicates are added together
    sum(SpecialFunctions.loggamma.(alpha + N)) - sum(SpecialFunctions.loggamma.(alpha)) + sum(SpecialFunctions.loggamma.(sum(alpha,dims=1))) - sum(SpecialFunctions.loggamma.(sum(alpha,dims=1) + sum(N,dims=1)))::Float64
end
function bayesian_score_component_uniform(
    i::Int,
    parents::AbstractVector{I},
    ncategories::AbstractVector{Int},
    data::AbstractMatrix{Int},
    prior::DirichletPrior,
) where {I<:Integer}

    (n, m) = size(data)
    if !isempty(parents)
        Np = length(parents)
        stridevec = fill(1, Np)
        for k in 2:Np
            stridevec[k] = stridevec[k-1] * ncategories[parents[k-1]]
        end
        js = (data[parents,:] .- 1)' * stridevec .+ 1
    else
        js = fill(1, m)
    end

    n = prod(ncategories[parents])
    N = sparse(data[i,:], js, 1, ncategories[i], n...) # note: duplicates are added together

    u = prior.α
    p = SpecialFunctions.loggamma(u)

    # Given a sparse N, we can be clever in our calculation and not waste time
    # computing the same SpecialFunctions.loggamma values by exploiting the sparse structure.
    sum0 = sum(SpecialFunctions.loggamma.(nonzeros(N) .+ u)) .+ p * (ncategories[i] * n - nnz(N))
    sum1 = n * ncategories[i] * p
    sum2 = n * SpecialFunctions.loggamma(ncategories[i] * u)
    cc = ncategories[i] * u

    @static if Base.VERSION.major == 0 && Base.VERSION.minor < 5
        sN = sparsevec(sum(N[1:ncategories[i],:], 1)) # Slower, but should be supported by Julia 0.4
    else
        sN = sum(N[i,:] for i=1:ncategories[i])
    end
    sum3 = sum(SpecialFunctions.loggamma.(nonzeros(sN) .+ cc)) .+ (size(N, 2) - nnz(sN)) * SpecialFunctions.loggamma(cc)
    sum0 - sum1 + sum2 - sum3::Float64
end
function bayesian_score_component(
    i::Int,
    parents::AbstractVector{I},
    ncategories::AbstractVector{Int},
    data::AbstractMatrix{Int},
    prior::DirichletPrior,
) where {I<:Integer}

    if typeof(prior) == UniformPrior
        return bayesian_score_component_uniform(i, parents, ncategories, data, prior)
    end
    alpha = get(prior, i, ncategories, parents)
    bayesian_score_component(i, parents, ncategories, data, alpha)
end

function bayesian_score(
    parent_list::Vector{Vector{Int}},
    ncategories::AbstractVector{Int},
    data::AbstractMatrix{Int},
    prior::DirichletPrior,
    )

    tot = 0.0
    for (i, p) in enumerate(parent_list)
        tot += bayesian_score_component(i, p, ncategories, data, prior)
    end
    tot
end
function bayesian_score(bn::DiscreteBayesNet, data::DataFrame, prior::DirichletPrior=UniformPrior())

    n = length(bn)
    parent_list = Array{Vector{Int}}(undef, n)
    ncategories = Array{Int}(undef, n)
    datamat = Matrix{Int}(data)'

    for (i,cpd) in enumerate(bn.cpds)
        parent_list[i] = inneighbors(bn.dag, i)
        ncategories[i] = infer_number_of_instantiations(convert(Vector{Int}, data[!,i]))
    end

    bayesian_score(parent_list, ncategories, datamat, prior)
end

function bayesian_score_component(
    i::Int,
    parents::AbstractVector{Int},
    ncategories::AbstractVector{Int},
    data::AbstractMatrix{Int},
    prior::DirichletPrior,
    cache::ScoreComponentCache,
    )

    if !haskey(cache[i], parents)
        (cache[i][parents] = bayesian_score_component(i, parents, ncategories, data, prior))
    end

    cache[i][parents]
end
function bayesian_score_components(
    parent_list::Vector{Vector{Int}},
    ncategories::AbstractVector{Int},
    data::AbstractMatrix{Int},
    prior::DirichletPrior,
    )

    score_components = Array{Float64}(undef, length(parent_list))
    for (i,p) in enumerate(parent_list)
        score_components[i] = bayesian_score_component(i, p, ncategories, data, prior)
    end
    score_components
end
function bayesian_score_components(
    parent_list::Vector{Vector{Int}},
    ncategories::AbstractVector{Int},
    data::AbstractMatrix{Int},
    prior::DirichletPrior,
    cache::ScoreComponentCache,
    )

    score_components = Array{Float64}(undef, length(parent_list))
    for (i,p) in enumerate(parent_list)
        score_components[i] = bayesian_score_component(i, p, ncategories, data, prior, cache)
    end
    score_components
end
function bayesian_score_components(bn::DiscreteBayesNet, data::DataFrame, prior::DirichletPrior=UniformPrior())

    n = length(bn)
    parent_list = Array{Vector{Int}}(undef, n)
    ncategories = Array{Int}(undef, n)
    datamat = Matrix{Int}(data)'

    for (i,cpd) in enumerate(bn.cpds)
        parent_list[i] = inneighbors(bn.dag, i)
        ncategories[i] = infer_number_of_instantiations(convert(Vector{Int}, data[!,i]))
    end

    bayesian_score_components(parent_list, ncategories, datamat, prior)
end

"""
    bayesian_score(G::DAG, names::Vector{Symbol}, data::DataFrame[, ncategories::Vector{Int}[, prior::DirichletPrior]])

Compute the bayesian score for graph structure `g`, with the data in `data`. `names` containes a symbol corresponding to each vertex in `g` that is the name of a column in `data`. `ncategories` is a vector of the number of values that each variable in the Bayesian network can take.

Note that every entry in data must be an integer greater than 0
"""
function bayesian_score(G::DAG,
                        names::Vector{Symbol},
                        data::DataFrame,
                        ncategories::Vector{Int}=Int[infer_number_of_instantiations(convert(Vector{Int}, data[!,n])) for n in names],
                        prior::DirichletPrior=UniformPrior())
    datamat = Array{Int}(undef, ncol(data), nrow(data))
    for i in 1:nv(G)
        datamat[i,:] = data[!,names[i]]
    end

    # NOTE: this is badj(G) prior to v0.6 and inneighbors(G) in v0.6
    backwards_adjacency = [inneighbors(G, i) for i in 1 : nv(G)]
    return bayesian_score(backwards_adjacency, ncategories, datamat, prior)
end


text:
using Downloads

# Link from https://mvnrepository.com/artifact/org.jacop/jacop/4.8.0
const jacop_url = "https://repo1.maven.org/maven2/org/jacop/jacop/4.8.0/jacop-4.8.0.jar"

const depsfile = joinpath(dirname(@__FILE__), "deps.jl")
if isfile(depsfile)
    rm(depsfile)
end

function write_depsfile(path)
    open(depsfile, "w") do f
        println(f, "const libjacopjava = \"$(escape_string(path))\"")
        return
    end
    return
end

const jacop_path = joinpath(dirname(@__FILE__), "jacop.jar")
Downloads.download(jacop_url, jacop_path, verbose=true)
if isfile(jacop_path)
    write_depsfile(jacop_path)
end