text:
### Tools for dealing with multiple algorithms
"""
Algorithm information for CUDNN algorithms that offer implementation selections.

$(FIELDS)

Methods: [`enums`](@ref), [`times`](@ref), [`timeat`](@ref), [`bytesat`](@ref)
"""
struct CUDNNAlgorithm
    """
    The enum value used to select this algorithm.

    This is a detail of cuDNN
    """
    enum::UInt32

    """
    The execution time of the algorithm - type match with the return value from the C++
    code.
    """
    time::Float32

    "The number of bytes required for working space"
    bytes::UInt64
end

"""
$(SIGNATURES)

Return the CUDNN algorithm enums from `algs`.
"""
enums(algs::Vector{CUDNNAlgorithm}) = map(x -> x.enum, algs)

"""
$(SIGNATURES)

Return the CUDNN algorithm times from `algs`.
"""
times(algs::Vector{CUDNNAlgorithm}) = map(x -> x.time, algs)

"""
$(SIGNATURES)

Return the time from `algs` for the algorithm witn enum encoding `enum`.
"""
function timeat(algs::Vector{CUDNNAlgorithm}, enum::Integer)
    ind = something(findfirst(x -> x.enum == enum, algs))
    return algs[ind].time
end

"""
$(SIGNATURES)

Return the working space requirement in bytes  from `algs` for the algorithm witn enum
encoding `enum`.
"""
function bytesat(algs::Vector{CUDNNAlgorithm}, e::Integer)
    ind = something(findfirst(x -> x.enum == e, algs))
    return algs[ind].bytes
end

#####
##### XTensor
#####

# We're going to abstract the notion of a Tensor so we can add a bunch of metadata to it
@enum TensorRole Arg Constant Intermediate

const TENSOR_GROUP_NUM = Ref(0)

"""
Wrapper around `TensorDescriptors` that includes extra metadata useful for generating
the ILP.

$(FIELDS)
"""
mutable struct XTensor{T}
    # Can possible map a single XTensor to multiple tensor descriptors if two args are
    # merged.
    "The `TensorDescriptor` that this `XTensor` represents."
    tensor::TensorDescriptor

    """
    Ordered collection of `XNode`s that reference this XTensor. The first entry is the
    producer of this tensor. All the rest read this tensor.
    """
    users::Vector{T}

    """
    The role of the tensor.
    """
    role::TensorRole

    """
    Flag for it this tensor is fixed at some specific location or not. Helpful for reducing
    the number of variables generated for the ILP.

    If the tensor is not fixed, this should be set to `nothing`.
    """
    fixed_at::Union{Nothing,TensorLocation}

    "Size of the tensor in bytes."
    size::Int

    "Collection of `TensorLocations` that this tensor can occupy."
    locations::Vector{TensorLocation}

    # Metadata for parameters and results
    "Indicate that all tensors in this group are actually the same tensor."
    inplace::Bool

    """
    Indicate the group number of this tensor. All tensors in the same group need to have the
    same PMM/DRAM assignment.

    For now, this should only be tweaked for IO tensors.
    """
    group::Int
end

isfixed(x::XTensor) = !isnothing(x.fixed_at)
fixed_location(x::XTensor) = something(x.fixed_at)

function XTensor(tensor::TensorDescriptor, ::Type{T}; fixed_at = nothing) where {T <: nGraph.AbstractBackendType}
    r = role(tensor)

    # For now - treat constants and inputs as zero size for the PMM case.
    #
    # CPU models are big enough that the extra memory consumption doesn't throw off
    # calculations that much.
    #
    # Accounting for IO size is done in later analysis steps.
    #
    # On the other hand, it's necessary to include this for the GPU code.
    # if isa(T, nGraph.CPU)
    #     if in(r, (Constant, Arg))
    #         sz = 0
    #     else
    #         sz = sizeof(tensor)
    #     end
    # else
    #     sz = sizeof(tensor)
    # end

    if in(r, (Constant, Arg))
        sz = 0
    else
        sz = sizeof(tensor)
    end

    # Add in tensors that may also live in PMEM
    xtensor = XTensor(
        tensor,
        XNode[],
        r,
        fixed_at,
        sz,
        TensorLocation[],
        false,
        TENSOR_GROUP_NUM[],
    )

    # Update the group number.
    TENSOR_GROUP_NUM[] += 1

    # Check if this tensor is fixed at a location. If so - just append that location
    if isfixed(xtensor)
        push!(xtensor.locations, fixed_at)
    else
        push!(xtensor.locations, DRAM)
        if !isconstant(xtensor)
            push!(xtensor.locations, PMEM)
        end
    end
    return xtensor
end

isinplace(t::XTensor) = t.inplace
makeinplace(t::XTensor) = (t.inplace = true)

"""
$(SIGNATURES)

Merge `XTensor`s `a` and `b` into the same tensor group, which will be the smallest group
of either `a` or `b`.

CAUTION: Currently doesn't merge all tensors in `a`'s or `b`'s group.
"""
function Base.merge!(a::XTensor, b::XTensor)
    # Sanitize inputs.
    if size(unx(a)) != size(unx(b))
        throw(DimensionMismatch("""
            Trying to merge two XTensors with different tensor sizes: $(size(a)), $(size(b))
            """)
        )
    end
    if eltype(unx(a)) != eltype(unx(b))
        throw(ArgumentError("""
            Trying to merge two XTensors with different eltypes: $(eltype(a)), $(eltype(b))
            """)
        )
    end

    # Set the group of the tensors to whatever the minimum is.
    # This makes merging two groups somewhat difficult since we don't have a way to collect
    # all tensors that belong to a group - but I don't think that will generally be a use
    # case.
    #
    # If that does become a ues-case, we'll have to keep a global tracker to all the live
    # allocated XTensors and use finalizers to clean up the list as XTensors are GC'd.
    n = min(a.group, b.group)
    a.group = n
    b.group = n
end

"""
$(SIGNATURES)

Merge all `XTensor`s in `A` into the same group.
"""
function Base.merge!(A::Vector{<:XTensor})
    # Find the lowest group number and merge everything with that.
    _, indmin = findmin(map(x -> x.group, A))
    map(x -> merge!(x, A[indmin]), A)
    return nothing
end

"""
$(SIGNATURES)

Return the `TensorRole` of a `TensorDescriptor` from it's name.
"""
function role(tensor::TensorDescriptor)
    tensorname = nGraph.name(tensor)
    if isparam(tensorname) || isresult(tensorname)
        role = Arg
    elseif isconstant(tensorname)
        role = Constant
    else
        role = Intermediate
    end
    return role
end

Utils.isconstant(t::XTensor) = (t.role == Constant)
isarg(t::XTensor) = t.role == Arg
is_persistent(t::XTensor) = nGraph.is_persistent(unx(t))
adduser!(t::XTensor, n) = !in(n, users(t)) && push!(t.users, n)

"""
$(SIGNATURES)

Return all users of `xtensor`.
"""
users(xtensor::XTensor) = xtensor.users

"""
$(SIGNATURES)

Return the XNode the produces `xtensor`.
"""
producer(xtensor::XTensor) = first(xtensor.users)

"""
$(SIGNATURES)

Return the last `XNode` that uses `xtensor`.
"""
consumer(xtensor::XTensor) = last(xtensor.users)

"""
$(SIGNATURES)

Return the allocation size of `xtensor`.
"""
Base.sizeof(xtensor::XTensor) = xtensor.size

"""
$(SIGNATURES)

Get the underlying `nGraph.TensorDescriptor` from `xtensor`.
"""
unx(xtensor::XTensor) = xtensor.tensor
locations(t::XTensor) = t.locations
nGraph.name(t::XTensor) = nGraph.name(t.tensor)
Base.show(io::IO, x::XTensor) = println(io, "XTensor: ", nGraph.name(x))
getoffset(t::XTensor) = nGraph.get_pool_offset(t.tensor)

#####
##### XNode
#####

"""
Wrapper for `nGraph.NodeDescriptor` with extra metadata.

$(FIELDS)
"""
mutable struct XNode
    "The `NodeDescriptor` this `XNode` represents"
    node::NodeDescriptor

    "The execution number of this node in the topological sort of the graph"
    index::Int

    """
    Timings for each possible `IOConfig` of the node. The collection of `IOConfigs` is
    decided at runtime based off the possible locations of the input and output tensors
    of this node.
    """
    timings::Dict{IOConfig, Union{Float64, Vector{CUDNNAlgorithm}}}

    "Ordered collection of output tensors for this node."
    outputs::Vector{XTensor}

    "Ordered collection of input tensors for this node."
    inputs::Vector{XTensor}

    "Tensors who begin life at this node."
    newlist::Vector{XTensor}

    "Tensors that may be freed after the execution of this node."
    freelist::Vector{XTensor}
end

function XNode(node::NodeDescriptor, index)
    # Don't determine the inputs or outputs yet.
    #
    # Need to make sure all XTensors and XNodes are consistent.
    # We leave this to the FunctionData constructor.
    return XNode(
        node,
        index,
        Dict{IOConfig, Union{Float64, Vector{CUDNNAlgorithm}}}(),
        XTensor[],
        XTensor[],
        XTensor[],
        XTensor[],
    )
end

_length_one(x) = (@assert(isone(length(x))); x)

can_select_algo(n::XNode, c::IOConfig) = isa(n.timings[c], Vector{CUDNNAlgorithm})
can_select_algo(n::XNode) = any(c -> can_select_algo(n, c), configs_for(n))

settime!(n::XNode, c::IOConfig, time) = n.timings[c] = time
configs_for(n::XNode) = keys(n.timings)
config_for(n::XNode) = configs_for(n) |> collect |> _length_one |> first

gettime(n::XNode, c::IOConfig) = n.timings[c]
gettime(n::XNode) = gettime(n, config_for(n))
unx(n::XNode) = n.node

Utils.hasprofile(n::XNode) = hasprofile(n.node)

nGraph.outputs(n::XNode) = n.outputs
nGraph.inputs(n::XNode) = n.inputs
nGraph.name(n::XNode) = nGraph.name(n.node)
Base.show(io::IO, n::XNode) = println(io, "XNode: ", nGraph.name(n))

Base.isless(a::XNode, b::XNode) = isless(a.index, b.index)

#####
##### Profile Data
#####

struct FunctionData{T}
    tensors::Set{XTensor}
    nodes::Vector{XNode}
end

nodes(f::FunctionData) = f.nodes
tensors(f::FunctionData) = f.tensors

FunctionData(fn::nGraph.NFunction, ::nGraph.Backend{T}) where {T} = FunctionData(fn, T)
function FunctionData(fn::nGraph.NFunction, ::Type{T}) where {T}
    tensors = Set{XTensor{XNode}}()
    nodes = XNode[]

    for (index, unwrapped_op) in enumerate(fn)
        op = NodeDescriptor(unwrapped_op)
        xnode = XNode(op, index)
        push!(nodes, xnode)

        for tensor in outputs(op)
            if isresult(op) || isparam(op)
                xtensor = XTensor(tensor, T; fixed_at = DRAM)
            else
                xtensor = XTensor(tensor, T)
            end
            # Create an XTensor from this output.
            #xtensor = XTensor(tensor, T)
            push!(tensors, xtensor)

            # Register the producing node as the first user
            adduser!(xtensor, xnode)

            # Register the xtensor as an output of the xnode
            push!(xnode.outputs, xtensor)
        end

        # Do the same thing for the node inputs
        for tensor in inputs(op)
            # Get the xtensor we made previously from the set of all xtensors
            xtensor = _get(tensors, tensor)

            # Register users and inputs
            adduser!(xtensor, xnode)
            push!(xnode.inputs, xtensor)
        end
    end

    # Run liveness analysis on the set of nodes.
    liveness!(nodes, tensors)

    return FunctionData{T}(tensors, nodes)
end

iosize(F::FunctionData) = sum(sizeof, Iterators.filter(x -> isarg(x) || isconstant(x), F.tensors))

function _get(a::Set{XTensor{XNode}}, t::TensorDescriptor)
    for i in a
        i.tensor == t && return i
    end
    throw(KeyError(t))
end

#####
##### Liveness Analysis
#####

function liveness!(nodes::Vector{XNode}, tensors::Set{XTensor{XNode}})
    # forward pass
    for op in nodes
        empty!(op.newlist)
        append!(op.newlist, filter(x -> !isarg(x) && !isconstant(x), outputs(op)))
    end

    # add all of the argument tensors to the liveness list for the first op
    ind = findfirst(hasprofile, nodes)
    append!(nodes[ind].newlist, filter(x -> isarg(x), tensors))

    # backward pass
    freed_tensors = Set{XTensor{XNode}}()
    for op in reverse(nodes)
        empty!(op.freelist)
        for tensor in inputs(op)
            if !in(tensor, freed_tensors) && !isarg(tensor) && !isconstant(tensor)
                push!(op.freelist, tensor)
                push!(freed_tensors, tensor)
            end
        end
    end
    ind = findlast(hasprofile, nodes)
    append!(nodes[ind].freelist, filter(x -> isarg(x), tensors))

    # Now, we need to do a cleanup phase.
    #
    # The BatchNorm pass can mess up some aspects of outputs from liveness analysis.
    # Here, we check to see if any tensor shows up in the new list but not the free list.
    # If so, we set its free point to the place where it was created.
    #
    # I hate batchnorm
    #
    # TODO: Fix this because it's harming liveness analysis for arguments.
    tensor_start = Dict{XTensor, Int}()
    for (index, op) in enumerate(nodes), tensor in op.newlist
        tensor_start[tensor] = index
    end

    for op in nodes, tensor in op.freelist
        @assert haskey(tensor_start, tensor)
        delete!(tensor_start, tensor)
    end

    for (tensor, index) in tensor_start
        isarg(tensor) || push!(nodes[index].freelist, tensor)
    end

    return nothing
end

# Convenience for iterating over live tensors
struct LiveTensorIterator
    data::FunctionData
    live_tensors::Set{XTensor{XNode}}
end

live_tensors(data::FunctionData) = LiveTensorIterator(data, Set{XTensor{XNode}}())
Base.length(L::LiveTensorIterator) = length(L.data.nodes)
function Base.iterate(L::LiveTensorIterator, s = 1)
    s > length(L) && return nothing

    # Free tensors from the previous iteration
    if !isone(s)
        for tensor in L.data.nodes[s-1].freelist
            delete!(L.live_tensors, tensor)
        end
    end

    # Add new tensors for this iteration
    for tensor in L.data.nodes[s].newlist
        push!(L.live_tensors, tensor)
    end

    return L.live_tensors, s+1
end

"""
    allocation_bounds(data::FunctionData)

Return upper and lower bounds on the amount of DRAM required for input, output,
constant, and intermediate tensors.

Upper bound is determined by the maximum tensors concurrently live.

Lower bound is determined by the total size of input, output, and constant tensors.
"""
function allocation_bounds(data::FunctionData)
    lower_bound = 0

    # Compute Upper Bound
    upper_bound = 0
    for tensors in live_tensors(data)
        if !isempty(tensors)
            upper_bound = max(upper_bound, sum(sizeof(n) for n in tensors))
        end
    end

    return (upper_bound = upper_bound, lower_bound = lower_bound)
end

#####
##### Valid locations that a tensor can live
#####

function possible_configs(backend::nGraph.Backend, data::FunctionData)
    configs = Set{Tuple{XNode, IOConfig}}()
    for node in nodes(data)
        hasprofile(node) || continue

        # Get possible configs based on backend type
        pc = possible_configs(backend, node)
        for config in pc
            push!(configs, (node, config))
        end
    end
    return configs
end

function possible_configs(::nGraph.Backend{nGraph.CPU}, node::XNode)
    config_inputs = [locations(t) for t in inputs(node)]
    config_outputs = [locations(t) for t in outputs(node)]

    # First - generate all of the configs generated by the Cartesian product of input and
    # output locations.
    configs = vec(map(x -> IOConfig(x...),
        Iterators.product(
            Iterators.product(config_inputs...),
            Iterators.product(config_outputs...),
        )
    ))

    # Filter out configs where some tensors belong to the same group and hence whose
    # locations should vary together>
    tensors = vcat(inputs(node), outputs(node))
    unique_groups = unique(map(x -> x.group, tensors))
    group_indices = [findall(x -> x.group == g, tensors) for g in unique_groups]
    for indices in group_indices
        filter!(x -> isone(length(unique(x[indices]))), configs)
    end

    return configs
end

function possible_configs(::nGraph.Backend{nGraph.GPU}, node::XNode)
    config = IOConfig(
        ntuple(_ -> DRAM, length(inputs(node))),
        ntuple(_ -> DRAM, length(outputs(node))),
    )
    return (config,)
end

function getconfig(n::nGraph.Node)
    f = x -> nGraph.is_persistent(x) ? PMEM : DRAM
    input = map(f, inputs(n)) |> Tuple
    output = map(f, outputs(n)) |> Tuple

    return IOConfig(input, output)
end

#####
##### Setup and cleanup code
#####

function _setup!(node::nGraph.Node, config::IOConfig)
    # Outputs
    for (i, location) in enumerate(config.outputs)
        if location == PMEM
            nGraph.make_persistent(nGraph.output(node, i))
        end
    end

    # Inputs
    for (i, location) in enumerate(config.inputs)
        if location == PMEM
            nGraph.make_persistent(nGraph.input(node, i))
        end
    end
    return nothing
end

max_stars_repo_name: darchr/AutoTM

max_issues_repo_name: darchr/AutoTM

max_forks_repo_name: darchr/AutoTM


text:
module TestBase

using GeoIP
using Sockets: IPv4, @ip_str
using Test

TEST_DATADIR = joinpath(dirname(@__FILE__), "data")

ENV["GEOIP_DATADIR"] = TEST_DATADIR
@testset "ZipFile loading" begin
    ENV["GEOIP_ZIPFILE"] = "GeoLite2-City-CSV.zip"

    geodata = load()

    @testset "Known result" begin
        ip1 = IPv4("1.0.8.1")
        geoip1 = geolocate(geodata, ip1)
        @test geoip1["country_iso_code"] == "CN"
        @test geoip1["time_zone"] == "Asia/Shanghai"
        @test ceil(Int, geoip1["location"].lat) == 24

        # String indexing
        geoip1 = geolocate(geodata, "1.0.8.1")
        @test geoip1["country_iso_code"] == "CN"
        @test geoip1["time_zone"] == "Asia/Shanghai"
        @test ceil(Int, geoip1["location"].lat) == 24
    end

    @testset "Null results" begin
        @test isempty(geolocate(geodata, ip"0.0.0.0"))
        @test isempty(geolocate(geodata, ip"127.0.0.1"))
    end

    @testset "Broadcasing" begin
        result = geolocate.(geodata, [ip"1.0.16.1", ip"1.0.8.1"])
        @test length(Set(result)) == 2
        @test !isempty(result[1])
        @test !isempty(result[2])
    end

    @testset "Dict indexing" begin
        geoip1 = geodata[ip"1.0.8.1"]
        @test geoip1["country_iso_code"] == "CN"
        @test geoip1["time_zone"] == "Asia/Shanghai"
        @test ceil(Int, geoip1["location"].lat) == 24
    end
end

end # module

ext: jl


text:
module SArrayMod

export SArray, SArrayTyped, nummissing, clip, clip_lower, clip_upper, fillna, count_ngrams, count_words

using Cxx
using Lazy
import Base: getindex, setindex!, show, +, -, *, /, .+, .-, ./, .*, .>, .<, .==, .>=, .<=

import SFrames: FlexibleTypeMod
import SFrames.FlexibleTypeMod: FlexibleType, Cflexible_type
import SFrames: head, tail, materialize, ismaterialized, sample, dropna, save, load, unpack
import SFrames.Util: cstring

abstract SAbstractArray

immutable SArray <: SAbstractArray
    val
end

Base.eltype(::SArray) = FlexibleType

immutable SArrayTyped{T} <: SAbstractArray
    array::SArray

    SArrayTyped(array::SArray) = new(array)
    SArrayTyped(val::Cxx.CppValue) = new(SArray(val))
end

cval(s::SArray) = s.val
cval(s::SArrayTyped) = s.array.val

Base.eltype{T}(::SArrayTyped{T}) = T

function SArray{T<:Union{Int,Float64}}(x::AbstractVector{T})
    if T==Int
        writer = icxx"new gl_sarray_writer(flex_type_enum::INTEGER, 1);"
    elseif T==Float64
        writer = icxx"new gl_sarray_writer(flex_type_enum::FLOAT, 1);"
    end
    for _ in x
        icxx"$writer->write($_, 0);"
    end
    array = icxx"$writer->close();"
    icxx"delete $writer;"
    SArray(array)
end

function SArray{T<:AbstractString}(x::AbstractVector{T})
    writer = icxx"return new gl_sarray_writer(flex_type_enum::STRING, 1);"
    for _ in x
        icxx"$writer->write($(cstring(_)), 0);"
    end
    array = icxx"$writer->close();"
    icxx"delete $writer;"
    SArray(array)
end

function SArray{T<:Associative}(x::AbstractVector{T})
    writer = icxx"return new gl_sarray_writer(flex_type_enum::DICT, 1);"
    for _ in x
        icxx"$writer->write($(FlexibleType(_).val), 0);"
    end
    array = icxx"$writer->close();"
    icxx"delete $writer;"
    SArray(array)
end

function getindex{T}(s::SArray, ::Type{T})
    SArrayTyped{T}(s)
end

function getindex{T<:Union{ASCIIString, UTF8String, AbstractString}}(s::SArray, ::Type{T})
    SArrayTyped{UTF8String}(s)
end

function getindex(x::SArray, idx::Int)
    FlexibleType(icxx"$(x.val)[$idx-1];")
end

function getindex(s::SArrayTyped, idx::Int)
    s.array[idx] |> eltype(s)
end

function getindex(s::SAbstractArray, idx::UnitRange)
    icxx"$(cval(s))[{$(idx.start-1), $(idx.stop)}];" |> typeof(s)
end

function getindex(s::SAbstractArray, idx::StepRange)
    icxx"$(cval(s))[{$(idx.start-1), $(idx.step), $(idx.stop)}];" |> typeof(s)
end

function Base.start(s::SAbstractArray)
    ra = icxx"$(cval(s)).range_iterator();"
    (icxx"$ra.begin();", icxx"$ra.end();")
end

function Base.next(s::SArray, state)
    iter, enditer = state
    value = icxx"auto val=*$iter; return val;" |> FlexibleType
    (value, (icxx"$iter+1;", enditer))
end

function Base.next(s::SArrayTyped, state)
    value, next_state = next(s.array, state)
    (eltype(s)(value), next_state)
end

function Base.done(s::SAbstractArray, state)
    iter, enditer = state
    icxx"$iter == $enditer;"
end

function Base.show(io::IO, s::SAbstractArray)
    sstream = icxx"stringstream();"
    icxx"$sstream<<$(cval(s));"
    c = icxx"$sstream.str().c_str();" |> bytestring
    write(io, c)
end

Base.size(x::SAbstractArray) = (Int(icxx"$(cval(x)).size();"),)
Base.ndims(::SAbstractArray) = 1
Base.length(x::SAbstractArray) = size(x)[1]

for (julia_op, c_op) in [
    (:+, "+"),
    (:-, "-"),
    (:/, "/"),
    (:*, "*"),
    (:.>, ">"),
    (:.>=, ">="),
    (:.==, "=="),
    (:.<, "<"),
    (:.<=, "<=")
    ]
    cstr = "\$(cval(s1)) $c_op \$(cval(s2));"
    cstr_scalar = "\$(cval(s1)) $c_op \$s2;"
    for op in [julia_op, Symbol(".$julia_op")]
        @eval  function $op(s1::SAbstractArray, s2::SAbstractArray)
            typeof(s1)($(Expr(:macrocall, Symbol("@icxx_str"), cstr)))
        end
    end

    @eval function $julia_op(s1::SAbstractArray, s2::Int)
        typeof(s1)($(Expr(:macrocall, Symbol("@icxx_str"), cstr_scalar)))
    end

    @eval $julia_op(s1::Int, s2::SAbstractArray) = $julia_op(s2, s1)
end

_get(f::FlexibleType) = get(f)
_get(x) = x

function Base.convert(::Type{Vector}, s::SAbstractArray)
    res = Int[]
    for i=1:length(s)
        push!(res, s[i]|>_get)
    end
    res
end

function save(s::SAbstractArray, dir, format="binary")
    icxx"$(cval(s)).save($dir, $format);"
end

function head(s::SAbstractArray, n)
    icxx"$(cval(s)).head($n);" |> typeof(s)
end

function tail(s::SAbstractArray, n)
    icxx"$(cval(s)).tail($n);" |> typeof(s)
end

function ismaterialized(s::SAbstractArray)
    icxx"$(cval(s)).is_materialized();"
end

function materialize(s::SAbstractArray)
    icxx"$(cval(s)).materialize();"
end

function sample(s::SAbstractArray, fraction)
    icxx"$(cval(s)).sample($fraction);" |> typeof(s)
end

function Base.all(s::SAbstractArray)
    icxx"$(cval(s)).all();"
end

function Base.any(s::SAbstractArray)
    icxx"$(cval(s)).any();"
end

function nummissing(s::SAbstractArray)
    icxx"$(cval(s)).num_missing();" |> Int
end

function clip(s::SAbstractArray, lower, upper)
    icxx"$(cval(s)).clip($lower, $upper);" |> typeof(s)
end

clip_lower(s::SAbstractArray, threshold) =
    icxx"$(cval(s)).clip_lower($threshold);" |> typeof(s)

clip_upper(s::SAbstractArray, threshold) =
    icxx"$(cval(s)).clip_upper($threshold);" |> typeof(s)

function dropna(s::SAbstractArray)
    icxx"$(cval(s)).dropna();" |> typeof(s)
end

function append(s1::SAbstractArray, s2::SAbstractArray)
    icxx"$(s1.val).append($(s2.val));" |> typeof(s)
end

function fillna(s::SAbstractArray, value)
    icxx"$(cval(s)).fillna(value);" |> typeof(s)
end

function Base.unique(s::SAbstractArray)
    icxx"$(cval(s)).unique();"
end

function Base.nnz(s::SAbstractArray)
    icxx"$(cval(s)).nnz();" |> Int
end

function generic_apply_int(f)
    x = Int(FlexibleType(f))
    r = (global_f::Function)(x)
    icxx"return flexible_type($r);"::Cflexible_type
end

const generic_apply_int_ptr = cfunction(generic_apply_int, Cflexible_type, (Cflexible_type,))

function generic_apply_float(f)
    x = Float64(FlexibleType(f))
    r = (global_f::Function)(x)
    icxx"return flexible_type($r);"::Cflexible_type
end

const generic_apply_float_ptr = cfunction(generic_apply_float, Cflexible_type, (Cflexible_type,))

function Base.apply(s::SArrayTyped{Int}, f, skip_undefined=true)
    global global_f = f
    res = icxx"$(cval(s)).apply((long (*)(long))$(generic_apply_int_ptr), flex_type_enum::INTEGER);"
    SArrayTyped{Int}(res)
end

function Base.apply(s::SArrayTyped{Float64}, f, skip_undefined=true)
    global global_f = f
    res = icxx"$(cval(s)).apply((double (*)(double))$(generic_apply_float_ptr), flex_type_enum::FLOAT);"
    SArrayTyped{Float64}(res)
end

# Base.map(f, s::SAbstractArray) = apply(s, f)

Base.maximum(s::SArray) = icxx"$(cval(s)).max();" |> FlexibleType
Base.minimum(s::SArray) = icxx"$(cval(s)).min();" |> FlexibleType
Base.mean(s::SArray) = icxx"$(cval(s)).mean();" |> FlexibleType
Base.std(s::SArray) = icxx"$(cval(s)).std();" |> FlexibleType

Base.maximum(s::SArrayTyped) = icxx"$(cval(s)).max();" |> FlexibleType |> eltype(s)
Base.minimum(s::SArrayTyped) = icxx"$(cval(s)).min();" |> FlexibleType |> eltype(s)
Base.mean(s::SArrayTyped) = icxx"$(cval(s)).mean();" |> FlexibleType |> Float64
Base.std(s::SArrayTyped) = icxx"$(cval(s)).std();" |> FlexibleType |> Float64

Base.endof(s::SAbstractArray) = length(s)

dict_keys(s::SArray) = icxx"$(s.val).dict_keys();" |> SArray
dict_values(s::SArray) = icxx"$(s.val).dict_values();" |> SArray

function count_ngrams(s::SArray, n=2, method="word", to_lower=true, ignore_space=true)
    icxx"$(s.val).count_ngrams($n, $(cstring(method)), $to_lower, $ignore_space);" |>
    SArray
end

function count_words(s::SArray, to_lower=true)
    icxx"$(s.val).count_words($to_lower);" |> SArray
end

function Base.get(s::SArray)
    map(get, s)
end

end


text:
version = v"1.5.0"
api_version = v"1.2.43"

# Collection of sources required to build this package
#
# The level zero repository contains both the API headers and a loader,
# which are versioned independently.
sources = [
    GitSource("https://github.com/oneapi-src/level-zero.git",
              "556cbac1a2adce87ff28c32813f23543700d95f2"),
]

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# Conjugate models for multivariate normal

using Distributions
using Base.Test

# MvNormal -- Normal (known covariance)

n = 3
p = 4
X = reshape(Float64[1:12;], p, n)
w = rand(n)
Xw = X * diagm(w)

Sigma = 0.75 * eye(p) + fill(0.25, 4, 4)
ss = suffstats(MvNormalKnownCov(Sigma), X)
ssw = suffstats(MvNormalKnownCov(Sigma), X, w)

s_t = sum(X, 2)
ws_t = sum(Xw, 2)
tw_t = length(w)
wtw_t = sum(w)

@test_approx_eq ss.sx s_t
@test_approx_eq ss.tw tw_t

@test_approx_eq ssw.sx ws_t
@test_approx_eq ssw.tw wtw_t

# Posterior
n = 100
mu_true = [2., 3.]
Sig_true = eye(2)
Sig_true[1,2] = Sig_true[2,1] = 0.25
mu0 = [2.5, 2.5]
Sig0 = eye(2)
Sig0[1,2] = Sig0[2,1] = 0.5
X = rand(MultivariateNormal(mu_true, Sig_true), n)

pri = MultivariateNormal(mu0, Sig0)

post = posterior((pri, Sig_true), MvNormal, X)
@test isa(post, FullNormal)

@test_approx_eq post.μ inv(inv(Sig0) + n*inv(Sig_true))*(n*inv(Sig_true)*mean(X,2) + inv(Sig0)*mu0)
@test_approx_eq post.Σ.mat inv(inv(Sig0) + n*inv(Sig_true))

# posterior_sample

ps = posterior_randmodel((pri, Sig_true), MvNormal, X)

@test isa(ps, FullNormal)
@test insupport(ps, ps.μ)
@test insupport(InverseWishart, ps.Σ.mat)

# NormalInverseWishart - MvNormal

mu_true = [2., 2.]
Sig_true = eye(2)
Sig_true[1,2] = Sig_true[2,1] = 0.25

X = rand(MultivariateNormal(mu_true, Sig_true), n)
Xbar = mean(X,2)
Xm = X .- mean(X,2)

mu0 = [2., 3.]
kappa0 = 3.
nu0 = 4.
T0 = eye(2)
T0[1,2] = T0[2,1] = .5
pri = NormalInverseWishart(mu0, kappa0, T0, nu0)

post = posterior(pri, MvNormal, X)

@test_approx_eq post.mu (kappa0*mu0 + n*Xbar)./(kappa0 + n)
@test_approx_eq post.kappa kappa0 + n
@test_approx_eq post.nu nu0 + n
@test_approx_eq (post.Lamchol[:U]'*post.Lamchol[:U]) T0 + A_mul_Bt(Xm, Xm) + kappa0*n/(kappa0+n)*(Xbar-mu0)*(Xbar-mu0)'

ps = posterior_randmodel(pri, MultivariateNormal, X)

@test isa(ps, MultivariateNormal)
@test insupport(ps, ps.μ) && insupport(InverseWishart, ps.Σ.mat)

# NormalWishart - MvNormal

mu_true = [2., 2.]
Lam_true = eye(2)
Lam_true[1,2] = Lam_true[2,1] = 0.25

X = rand(MvNormal(mu_true, inv(Lam_true)), n)
Xbar = mean(X,2)
Xm = X .- Xbar

mu0 = [2., 3.]
kappa0 = 3.
nu0 = 4.
T0 = eye(2)
T0[1,2] = T0[2,1] = .5
pri = NormalWishart(mu0, kappa0, T0, nu0)

post = posterior(pri, MvNormal, X)

@test_approx_eq post.mu (kappa0*mu0 + n*Xbar)./(kappa0 + n)
@test_approx_eq post.kappa kappa0 + n
@test_approx_eq post.nu nu0 + n
@test_approx_eq (post.Tchol[:U]'*post.Tchol[:U]) T0 + A_mul_Bt(Xm, Xm) + kappa0*n/(kappa0+n)*(Xbar-mu0)*(Xbar-mu0)'

ps = posterior_randmodel(pri, MvNormal, X)

@test isa(ps, MultivariateNormal)
@test insupport(ps, ps.μ)
@test insupport(InverseWishart, ps.Σ.mat)  # InverseWishart on purpose


text:
if lowercase(get(ENV, "CI", "false")) == "true"
    include("install_dependencies.jl")
end

using Test
using CausalityTools
using Distributions
using UncertainData

# @testset "Plot recipes" begin
# 	include("plot_recipes.jl")
# end

@testset "Discrete example systems" begin
    include("systems/discrete/test_discrete_systems.jl")
end

@testset "Continuous example systems" begin
    include("systems/continuous/test_continuous_systems.jl")
end

@testset "High level wrappers" begin
    include("test_wrappers_te.jl")
end

@testset "Causality tests on scalar time series" begin
    include("causality_tests/test_scalarseries_CrossMappingTest.jl")
    include("causality_tests/test_scalarseries_ConvergentCrossMappingTest.jl")
    include("causality_tests/test_scalarseries_JointDistanceDistributionTest.jl")
    include("causality_tests/test_scalarseries_JointDistanceDistributionTTest.jl")
    include("causality_tests/test_scalarseries_TransferOperatorGridTest.jl")
    include("causality_tests/test_scalarseries_VisitationFrequencyTest.jl")
    include("causality_tests/test_scalarseries_ApproximateSimplexIntersectionTest.jl")
    include("causality_tests/test_scalarseries_ExactSimplexIntersectionTest.jl")
    include("causality_tests/test_scalarseries_PredictiveAsymmetryTest.jl")
    include("causality_tests/test_scalarseries_SMeasure.jl")
end

@testset "DynamicalSystems.jl integration" begin
    include("causality_tests/integration_dynamicalsystems/test_integration_discrete_system.jl")
    include("causality_tests/integration_dynamicalsystems/test_integration_continuous_system.jl")
end

@testset "UncertainData.jl integration" begin
    include("causality_tests/integration_uncertaindata/test_uncertaindata_integration_JointDistanceDistributionTest.jl")
    include("causality_tests/integration_uncertaindata/test_uncertaindata_integration_JointDistanceDistributionTTest.jl")
    include("causality_tests/integration_uncertaindata/test_uncertaindata_integration_CrossMappingTest.jl")
    include("causality_tests/integration_uncertaindata/test_uncertaindata_integration_ConvergentCrossMappingTest.jl")
    include("causality_tests/integration_uncertaindata/test_uncertaindata_integration_TransferOperatorGridTest.jl")
    include("causality_tests/integration_uncertaindata/test_uncertaindata_integration_VisitationFrequencyTest.jl")
    include("causality_tests/integration_uncertaindata/test_uncertaindata_integration_ApproximateSimplexIntersectionTest.jl")
    include("causality_tests/integration_uncertaindata/test_uncertaindata_integration_ExactSimplexIntersectionTest.jl")
    include("causality_tests/integration_uncertaindata/test_uncertaindata_integration_PredictiveAsymmetryTest.jl")

    include("causality_tests/integration_uncertaindata/test_UncertainIndexValueDataset.jl")

    # High-level tests
    include("causality_tests/highlevel_tests/test_ConstrainedTest.jl")
    include("causality_tests/highlevel_tests/test_InterpolateBinTest.jl")
    include("causality_tests/highlevel_tests/test_RandomSequencesTest.jl")

    include("causality_tests/integration_uncertaindata/test_uncertain_indexvalue_dataset_with_schemes.jl")

end


text:
function todo_input_ui(store)
    todos = get_state(store).todos
    flag = CImGui.ImGuiInputTextFlags_AutoSelectAll |
           CImGui.ImGuiInputTextFlags_EnterReturnsTrue
    buffer_len = 128
    spacing = CImGui.GetStyle().ItemInnerSpacing.x

    all_complete = Ref(all(x->x.completed, todos))
    CImGui.Checkbox("###todo_all", all_complete) && dispatch!(store, Todo.CompleteAllTodos())
    CImGui.SameLine(0.0, spacing)
    str = Vector{Cuchar}("What needs to be done?")
    append!(str, ['\0' for i = 1:buffer_len-length(str)])
    CImGui.PushItemWidth(-1)
    if CImGui.InputText("###todo_input", str, buffer_len, flag)
        new_str = GC.@preserve str unsafe_string(pointer(str))
        dispatch!(store, Todo.AddTodo(new_str))
    end
    CImGui.PopItemWidth()
end

function todo_list_ui(store)
    todos = get_state(store).todos
    visi = get_state(store).visibility
    flag = CImGui.ImGuiInputTextFlags_EnterReturnsTrue
    spacing = CImGui.GetStyle().ItemInnerSpacing.x
    buffer_len = 128

    CImGui.PushItemWidth(-1)
    if CImGui.ListBoxHeader("List", -2)
        for i = 1:length(todos)
            id = todos[i].id
            status = todos[i].completed
            visi.filter == TodoFilter.SHOW_ACTIVE && status && continue
            visi.filter == TodoFilter.SHOW_COMPLETED && !status && continue
            # todo marker
            check = Ref(status)
            CImGui.Checkbox("###todo_$id", check) && dispatch!(store, Todo.CompleteTodo(id))
            CImGui.SameLine(0.0, spacing)
            # todo text
            old_text = todos[i].text
            if check[]
                CImGui.TextDisabled(old_text)
                CImGui.SameLine(CImGui.GetWindowContentRegionMax().x-18, spacing)
            else
                CImGui.PushItemWidth(-18)
                str = Vector{Cuchar}(old_text)
                append!(str, ['\0' for i = 1:buffer_len-length(str)])
                if CImGui.InputText("###todolist_$i", str, buffer_len, flag)
                    new_str = GC.@preserve str unsafe_string(pointer(str))
                    dispatch!(store, Todo.EditTodo(id, new_str))
                end
                CImGui.PopItemWidth()
                CImGui.SameLine(0.0, spacing)
            end
            # delete todo
            CImGui.PushID(i)
            CImGui.PushStyleColor(CImGui.ImGuiCol_Button, CImGui.HSV(0.0, 0.3, 0.3))
            CImGui.PushStyleColor(CImGui.ImGuiCol_ButtonHovered, CImGui.HSV(0.0, 0.6, 0.6))
            CImGui.PushStyleColor(CImGui.ImGuiCol_ButtonActive, CImGui.HSV(0.0, 0.9, 0.9))
            CImGui.PushItemWidth(-1)
            CImGui.Button("X") && dispatch!(store, Todo.DeleteTodo(id))
            CImGui.PopItemWidth()
            CImGui.PopStyleColor(3)
            CImGui.PopID()
        end
        CImGui.ListBoxFooter()
    end
    CImGui.PopItemWidth()
end

function todo_mvc_ui(store)
    # flag = 0
    flag = CImGui.ImGuiWindowFlags_NoTitleBar |
           CImGui.ImGuiWindowFlags_NoResize |
           CImGui.ImGuiWindowFlags_NoSavedSettings |
           CImGui.ImGuiWindowFlags_NoFocusOnAppearing |
           CImGui.ImGuiWindowFlags_NoNav
    todos = get_state(store).todos
    visi = get_state(store).visibility
    CImGui.Begin("TodoMVC", Ref(true), flag)
        todo_input_ui(store)
        CImGui.Separator()
        todo_list_ui(store)
        CImGui.Separator()
        # active items counter
        active_item_num = count(x -> !x.completed, todos)
        CImGui.TextColored((1.0,1.0,0.0,1.0), "$active_item_num items left.")
        CImGui.SameLine()
        # visibility filter
        CImGui.RadioButton("All", visi.filter == TodoFilter.SHOW_ALL) && dispatch!(store, TodoFilter.SetVisibilityFilter(TodoFilter.SHOW_ALL))
        CImGui.SameLine()
        CImGui.RadioButton("Active", visi.filter == TodoFilter.SHOW_ACTIVE) && dispatch!(store, TodoFilter.SetVisibilityFilter(TodoFilter.SHOW_ACTIVE))
        CImGui.SameLine()
        CImGui.RadioButton("Completed", visi.filter == TodoFilter.SHOW_COMPLETED) && dispatch!(store, TodoFilter.SetVisibilityFilter(TodoFilter.SHOW_COMPLETED))
        CImGui.SameLine()
        # clear completed
        if count(x -> x.completed, todos) != 0
            CImGui.PushItemWidth(-1)
            CImGui.PushStyleColor(CImGui.ImGuiCol_Button, CImGui.HSV(0.0, 0.5, 0.5))
            CImGui.PushStyleColor(CImGui.ImGuiCol_ButtonHovered, CImGui.HSV(0.0, 0.7, 0.7))
            CImGui.PushStyleColor(CImGui.ImGuiCol_ButtonActive, CImGui.HSV(0.0, 0.9, 0.9))
            CImGui.Button("Clear completed") && dispatch!(store, Todo.ClearCompleted())
            CImGui.PopStyleColor(3)
            CImGui.PopItemWidth()
        end
    CImGui.End()
end


text:
function jt10b()
#
#  test polynomial suggested by Jenkins and Traub
#
    a = 1e6;
    p = poly([a;1;1/a])
    z = [[a;1;1/a] ones(3)]
    p, PolyZeros(z)
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 1dd53841-569e-4501-b791-1fd07d8df47f
using Pkg, DrWatson

# ╔═╡ 3515bef0-6ac0-464a-beb0-087719922e5c
begin
	@quickactivate "NumericalMethodsforEngineers"
	using NumericalMethodsforEngineers, Test
end

# ╔═╡ ae8610bd-7f9c-4da9-8e51-ec2590cac91f
md" ## f-02-07.jl"

# ╔═╡ 12bdfd3f-3de0-4ed6-af37-c009af54360c
begin
	a = [16. 4. 8.; 4. 5. -4.; 8. -4. 22.]
	b = [4., 2., 5.]
end

# ╔═╡ 0a0eb8e6-8830-4850-84ab-a2ff12574b8f
as = sparse(a)

# ╔═╡ ab5cb4dd-8d90-44aa-84e2-214f10b6b8f2
begin
	f = cholesky(as)
	upper = f.U
end

# ╔═╡ 4093cc75-910d-420e-ac6e-2a7d19108c22
lower = f.L

# ╔═╡ 042eb9b3-0c0b-4a85-96d1-9dd80b733bf1
y = lower \ b

# ╔═╡ 26fce91d-c747-4f9b-9318-432cb33b7296
md" ##### Solution Vector:"

# ╔═╡ 4d09b4a2-4228-4128-a9b7-2a6c5e984c37
begin
	c = upper \ y
	c
end

# ╔═╡ aec7ad54-783f-4fc2-b957-5d6f98ea8deb
@test c == a \ b

# ╔═╡ fa236fb7-a681-4fba-a458-746f183756da
a1 = [16. 4. 8.; 4. 5. -4.; 8. -4. 22.]

# ╔═╡ 75dc1c4d-a0cc-4186-92f5-b241f53cf318
b1 = [4., 2., 5.]

# ╔═╡ be9d0599-52ec-4565-a8e6-f78fbe697489
a1[2, 2] += 1.0e20

# ╔═╡ ddd06285-0de7-4df0-8c41-ce2bc9ca71b9
b1[2] = a1[2, 2] * 5.0

# ╔═╡ 86ef4033-d955-450e-97bf-cb00ccea27df
F = cholesky(a1)

# ╔═╡ e1c22dbb-5115-46d1-8e1d-7de8020000fa
F.U

# ╔═╡ fa2d3458-d0c6-4c1c-af5a-4957a8abe29c
F.L

# ╔═╡ f256cf15-ca6d-40d3-9a9e-bd02a0484c4f
y1 = F.L \ b

# ╔═╡ b4887565-0f49-4782-9aae-e527bfb01a44
c1 = F.U \ y1

# ╔═╡ aea18ff7-21b2-4a8a-a953-05b9ad6549e2
md" ##### Solution Vector:"

# ╔═╡ 35fae17c-71ff-4158-a5be-4c3b8b143486
c

# ╔═╡ 0c8f15ce-1ef7-4b8a-b122-5632c0793a93
@test ≈(c, a \ b, atol=1e-5)

# ╔═╡ 9a8af48d-f8f1-434b-9640-024778bb4d0c
md" ### End of f-02-07.jl"

# ╔═╡ Cell order:
# ╟─ae8610bd-7f9c-4da9-8e51-ec2590cac91f
# ╠═1dd53841-569e-4501-b791-1fd07d8df47f
# ╠═3515bef0-6ac0-464a-beb0-087719922e5c
# ╠═12bdfd3f-3de0-4ed6-af37-c009af54360c
# ╠═0a0eb8e6-8830-4850-84ab-a2ff12574b8f
# ╠═ab5cb4dd-8d90-44aa-84e2-214f10b6b8f2
# ╠═4093cc75-910d-420e-ac6e-2a7d19108c22
# ╠═042eb9b3-0c0b-4a85-96d1-9dd80b733bf1
# ╠═26fce91d-c747-4f9b-9318-432cb33b7296
# ╠═4d09b4a2-4228-4128-a9b7-2a6c5e984c37
# ╠═aec7ad54-783f-4fc2-b957-5d6f98ea8deb
# ╠═fa236fb7-a681-4fba-a458-746f183756da
# ╠═75dc1c4d-a0cc-4186-92f5-b241f53cf318
# ╠═be9d0599-52ec-4565-a8e6-f78fbe697489
# ╠═ddd06285-0de7-4df0-8c41-ce2bc9ca71b9
# ╠═86ef4033-d955-450e-97bf-cb00ccea27df
# ╠═e1c22dbb-5115-46d1-8e1d-7de8020000fa
# ╠═fa2d3458-d0c6-4c1c-af5a-4957a8abe29c
# ╠═f256cf15-ca6d-40d3-9a9e-bd02a0484c4f
# ╠═b4887565-0f49-4782-9aae-e527bfb01a44
# ╠═aea18ff7-21b2-4a8a-a953-05b9ad6549e2
# ╠═35fae17c-71ff-4158-a5be-4c3b8b143486
# ╠═0c8f15ce-1ef7-4b8a-b122-5632c0793a93
# ╟─9a8af48d-f8f1-434b-9640-024778bb4d0c


text: Dict{String,Any}("newline" => Dict{String,Any}("value" => "lf","type" => "string"),"os" => Dict{String,Any}("value" => "unix","type" => "string"))

max_stars_repo_name: SmiteWindows/julia

max_issues_repo_name: SmiteWindows/julia

max_forks_repo_name: SmiteWindows/julia


text:
using Test
using ModelingToolkit
using Graphs
using SparseArrays
using UnPack

# Define some variables
@parameters t L g
@variables x(t) y(t) w(t) z(t) T(t)
D = Differential(t)

# Simple pendulum in cartesian coordinates
eqs = [D(x) ~ w,
       D(y) ~ z,
       D(w) ~ T*x,
       D(z) ~ T*y - g,
       0 ~ x^2 + y^2 - L^2]
pendulum = ODESystem(eqs, t, [x, y, w, z, T], [L, g], name=:pendulum)
sys = initialize_system_structure(pendulum)
StructuralTransformations.find_solvables!(sys)
sss = structure(sys)
@unpack graph, solvable_graph, fullvars, varassoc = sss
@test isequal(fullvars, [D(x), D(y), D(w), D(z), x, y, w, z, T])
@test graph.fadjlist == [[1, 7], [2, 8], [3, 5, 9], [4, 6, 9], [5, 6]]
@test graph.badjlist == 9 == length(fullvars)
@test ne(graph) == nnz(incidence_matrix(graph)) == 12
@test nv(solvable_graph) == 9 + 5
@test varassoc == [0, 0, 0, 0, 1, 2, 3, 4, 0]

se = collect(StructuralTransformations.𝑠edges(graph))
@test se == mapreduce(vcat, enumerate(graph.fadjlist)) do (s, d)
    StructuralTransformations.BipartiteEdge.(s, d)
end
@test_throws ArgumentError collect(StructuralTransformations.𝑑edges(graph))
@test_throws ArgumentError collect(StructuralTransformations.edges(graph))


text:
# this file defines iterators used for looping over a grid
struct UpdateFlags
    nodes::Bool
    coords::Bool
    celldofs::Bool
end

UpdateFlags(; nodes::Bool=true, coords::Bool=true, celldofs::Bool=true) =
    UpdateFlags(nodes, coords, celldofs)

"""
    CellIterator(grid::Grid)
    CellIterator(grid::DofHandler)

Return a `CellIterator` to conveniently loop over all the cells in a grid.

# Examples
```julia
for cell in CellIterator(grid)
    coords = getcoordinates(cell) # get the coordinates
    dofs = celldofs(cell)         # get the dofs for this cell
    reinit!(cv, cell)             # reinit! the FE-base with a CellIterator
end
```
"""
struct CellIterator{dim,C,T}
    flags::UpdateFlags
    grid::Grid{dim,C,T}
    current_cellid::ScalarWrapper{Int}
    nodes::Vector{Int}
    coords::Vector{Vec{dim,T}}
    cellset::Vector{Int}
    dh::Union{DofHandler{dim,C,T}, MixedDofHandler{dim,C,T}} #Future: remove DofHandler and rename MixedDofHandler->DofHandler
    celldofs::Vector{Int}

    function CellIterator{dim,C,T}(dh::Union{DofHandler{dim,C,T}, MixedDofHandler{dim,C,T}}, cellset::AbstractVector{Int}, flags::UpdateFlags) where {dim,C,T}
        isconcretetype(C) || _check_same_celltype(grid, cellset)
        N = nnodes_per_cell(dh.grid, first(cellset))
        cell = ScalarWrapper(0)
        nodes = zeros(Int, N)
        coords = zeros(Vec{dim,T}, N)
        n = ndofs_per_cell(dh)
        celldofs = zeros(Int, n)
        return new{dim,C,T}(flags, dh.grid, cell, nodes, coords, cellset, dh, celldofs)
    end

    function CellIterator{dim,C,T}(grid::Grid{dim,C,T}, cellset::AbstractVector{Int}, flags::UpdateFlags) where {dim,C,T}
        isconcretetype(C) || _check_same_celltype(grid, cellset)
        N = nnodes_per_cell(grid, first(cellset))
        cell = ScalarWrapper(0)
        nodes = zeros(Int, N)
        coords = zeros(Vec{dim,T}, N)
        return new{dim,C,T}(flags, grid, cell, nodes, coords)
    end
end

CellIterator(grid::Grid{dim,C,T}, cellset::AbstractVector{Int}=1:getncells(grid), flags::UpdateFlags=UpdateFlags()) where {dim,C,T} =
    CellIterator{dim,C,T}(grid, cellset, flags)
CellIterator(dh::Union{DofHandler{dim,C,T}, MixedDofHandler{dim,C,T}}, cellset::AbstractVector{Int}=1:getncells(dh.grid), flags::UpdateFlags=UpdateFlags()) where {dim,C,T} =
    CellIterator{dim,C,T}(dh, cellset, flags)

# iterator interface
function Base.iterate(ci::CellIterator, state = 1)
    if state > length(ci.cellset)
        return nothing
    else
        return (reinit!(ci, state), state+1)
    end
end
Base.length(ci::CellIterator)  = length(ci.cellset)

Base.IteratorSize(::Type{T})   where {T<:CellIterator} = Base.HasLength() # this is default in Base
Base.IteratorEltype(::Type{T}) where {T<:CellIterator} = Base.HasEltype() # this is default in Base
Base.eltype(::Type{T})         where {T<:CellIterator} = T

# utility
@inline getnodes(ci::CellIterator) = ci.nodes
@inline getcoordinates(ci::CellIterator) = ci.coords
@inline nfaces(ci::CellIterator) = nfaces(eltype(ci.grid.cells))
@inline onboundary(ci::CellIterator, face::Int) = ci.grid.boundary_matrix[face, ci.current_cellid[]]
@inline cellid(ci::CellIterator) = ci.current_cellid[]
@inline celldofs!(v::Vector, ci::CellIterator) = celldofs!(v, ci.dh, ci.current_cellid[])
@inline celldofs(ci::CellIterator) = ci.celldofs

function reinit!(ci::CellIterator{dim,C}, i::Int) where {dim,C}
    ci.current_cellid[] = ci.cellset[i]

    ci.flags.nodes  && cellnodes!(ci.nodes, ci.dh, ci.current_cellid[])
    ci.flags.coords && cellcoords!(ci.coords, ci.dh, ci.current_cellid[])

    if isdefined(ci, :dh) && ci.flags.celldofs # update celldofs
        celldofs!(ci.celldofs, ci.dh, ci.current_cellid[])
    end
    return ci
end

function check_compatible_geointerpolation(cv::Union{CellValues, FaceValues}, ci::CellIterator)
    if length(getnodes(ci)) != getngeobasefunctions(cv)
        msg = """The given CellValues and CellIterator are incompatiblet.
        Likely an appropriate geometry interpolate must be passed when constructing the CellValues.
        See also issue #265: https://github.com/KristofferC/JuAFEM.jl/issues/265"""
        throw(ArgumentError(msg))
    end
end

@inline function reinit!(cv::CellValues{dim,T}, ci::CellIterator{dim,N,T}) where {dim,N,T}
    check_compatible_geointerpolation(cv, ci)
    reinit!(cv, ci.coords)
end

@inline function reinit!(fv::FaceValues{dim,T}, ci::CellIterator{dim,N,T}, face::Int) where {dim,N,T}
    check_compatible_geointerpolation(fv, ci)
    reinit!(fv, ci.coords, face)
end

function _check_same_celltype(grid::AbstractGrid, cellset::AbstractVector{Int})
    celltype = typeof(grid.cells[first(cellset)])
    for cellid in cellset
        if celltype != typeof(grid.cells[cellid])
            error("You are trying to use a CellIterator to loop over a cellset with different celltypes.")
        end
    end
end


text:
using Base.Threads
using Base.Threads: @spawn, fetch
const VERBOSE = false

struct Partition{k}
    nt::Int
    n::Int
end
chunksize(pt::Partition{k}) where {k} = k
chunksize(pt::Type{Partition{k}}) where {k} = k
function index_slow(pt, i)
    j, i′ = divrem(i-1, chunksize(pt))
    j + 1, i′ + 1
end
@generated function index_fast(pt, i)
  pow = round(Int, log2(chunksize(pt)))
  :((i-1) >> $pow + 1, (i-1) & $(2^pow-1) + 1)
end

@generated function (pt::Partition{k})(i) where {k}
  ispow2(k) ?
    :(index_fast(pt, i)) :
    :(index_slow(pt, i))
end

Partition(nt,n) = Partition{div(n, nt)}(nt, n)
Base.length(pt::Partition) = pt.nt
#(pt::Partition)(i) = div((i-1), pt.k) + 1, i - pt.k*div((i-1), pt.k)
(pt::Partition{k})(q1, q2) where{k} = (q1-1)*k + q2
each(pt::Partition) = 1:pt.nt

function parallel_innermost!(rng, partition, G, G1, G2, ∇ϕ, i, t, x, θ, t′, Q, c, b, t_old,
    F::Union{ZigZag,FactBoomerang}, (factor, adapt), args...)
    t, x, θ = smove_forward!(G, i, t, x, θ, t′, F)
    ∇ϕi = ∇ϕ(x, i, args...)
    l, lb = sλ(∇ϕi, i, x, θ, F), sλ̄(b[i], t[i] - t_old[i])
    num = 1
    if rand(rng)*lb < l
        if l >= lb
            !adapt && error("Tuning parameter `c` too small.")
            adapt!(c, i, factor)
        end
        t, x, θ = smove_forward!(G2, i, t, x, θ, t′, F)
        θ = reflect!(i, ∇ϕi, x, θ, F)
        for j in neighbours(G1, i)
            b[j] = ab(G1, j, x, θ, c, F)
            t_old[j] = t[j]
            q1, q2 = partition(j)
            Q[q1][q2] = t[j] + poisson_time(b[j], rand(rng))
        end
        return true
    else
        b[i] = ab(G1, i, x, θ, c, F)
        t_old[i] = t[i]
        q1, q2 = partition(i)
        Q[q1][q2] = t[i] + poisson_time(b[i], rand(rng))
        return false
    end
end

function parallel_spdmp_inner!(latch, wakeup, ret, events, partition, ti, (t0, Δ), inner, G, G1, G2, ∇ϕ, t, x, θ, Q, c, b, t_old,
        F::Union{ZigZag,FactBoomerang}, (factor, adapt), args...)
    acc = num = 0
    VERBOSE && println("$ti starts.")
    tnext = t0 + Δ
    rng = Rng()
    while true
        num += 1
        ii, t′ = peek(Q[ti])
        i = partition(ti, ii)
        if !inner[i] ||  t′ > tnext # need neighbours at t′, or just wait
            tnext = t′ + Δ
            ret[] = i, t′, acc, num
            u = UInt(1) << (ti - 1)
            ac = Threads.atomic_and!(latch.active, ~u)
            done = false
            lock(wakeup) do
                if ac == u # last one turns the light off
                    VERBOSE && println("Notify $ti")
                    lock(latch.condition) do
                        notify(latch.condition, nothing; all = true, error = false)
                    end
                end

                VERBOSE && println("Sleep: $ti at $(t′) ($ac, $u).")
                done = wait(wakeup)
            end
            if done
                println("Done: $ti.")
                return
            end
            acc = num = 0
        else
            success = parallel_innermost!(rng, partition, G, G1, G2, ∇ϕ, i, t, x, θ, t′, Q, c, b, t_old, F, (factor, adapt), args...)
            success || continue
            acc += 1
            push!(events, event(i, t, x, θ, F))
        end
    end
end

function parallel_spdmp(partition, ∇ϕ, t0, x0, θ0, T, c, G, F::Union{ZigZag,FactBoomerang}, args...;
    factor=1.8, adapt=false, Δ = 0.1, progress=true, progress_stops = 20)
    nthr = length(partition)
    n = length(x0)
    t′ = fill(t0, nthr)
    t = fill(t0, size(θ0)...)
    t_old = copy(t)
    if G === nothing
        G = [i => rowvals(F.Γ)[nzrange(F.Γ, i)] for i in eachindex(θ0)]
    end
    inner = [all(j -> partition(j)[1] == partition(i)[1], js) for  (i,js) in G]

    G1 = [i => rowvals(F.Γ)[nzrange(F.Γ, i)] for i in eachindex(θ0)]

    @assert all(a.second ⊇ b.second for (a,b) in zip(G, G1))

    G2 = [i => setdiff(union((G1[j].second for j in G1[i].second)...), G[i].second) for i in eachindex(G1)]

    if !all(all( first(partition(g[1])) .== first.(partition.(g[2]))) for g in G2)
        error("Upper bounds may not depend across chunks.")

    end
    x, θ = copy(x0), copy(θ0)

    Q = [SPriorityQueue{Int,Float64}() for thr in 1:nthr]
    b = [ab(G1, i, x, θ, c, F) for i in eachindex(θ)]
    for i in eachindex(θ)
        q1, q2 = partition(i)
        enqueue!(Q[q1], q2 => poisson_time(b[i], rand()))
    end
    Ξ = Trace(t0, x0, θ0, F)
    task = Vector{Task}(undef, length(partition))
    evtime = zeros(nthr)
    perm = collect(1:nthr)
    waitfor = zeros(Int, nthr)
    events = [resize!([event(1, 0., x, θ, F)], 0) for ts in each(partition)]
    res = [Ref((1, 1.0, 1, 1)) for ts in each(partition)]
    latch = (;active = Threads.Atomic{UInt}(1), condition=Threads.Condition())
    wakeup = [Threads.Condition() for _ in each(partition)]
    if progress
        prg = Progress(progress_stops, 1)
    else
        prg = missing
    end
    parallel_spdmp_loop(t′, T, task, waitfor, latch, wakeup, evtime, perm, res, Ξ, events, partition, inner, G, G1, G2, ∇ϕ, t, x, θ, Q,
    c, b, t_old, F, (Δ, factor, adapt), prg, args...)
end
function parallel_spdmp_loop(t′, T, task, waitfor, latch, wakeup, evtime, perm, res, Ξ, events, partition, inner, G, G1, G2, ∇ϕ, t, x, θ, Q,
    c, b, t_old, F, (Δ, factor, adapt), prg, args...)

    tmin = minimum(t′)
    for ti in each(partition)
        Threads.atomic_or!(latch.active, UInt(1) << (ti - 1))
        task[ti] = Threads.@spawn parallel_spdmp_inner!(latch, wakeup[ti], res[ti], events[ti], partition, ti, (tmin, Δ), inner, G, G1, G2, ∇ϕ, t, x, θ, Q,
            c, b, t_old, F, (factor, adapt), args...)
    end
    task_outer = Threads.@spawn parallel_spdmp_outer!(tmin, t′, T, task, waitfor, latch, wakeup, evtime, perm, res, Ξ, events, partition, inner, G, G1, G2, ∇ϕ, t, x, θ, Q,
    c, b, t_old, F, (Δ, factor, adapt), prg, args...)

    for ti in each(partition)
        wait(task[ti])
    end

    (acc, num), (run, runs), (run2, runs2) = fetch(task_outer)
    sort!(Ξ.events, by=ev->ev[1])
    ismissing(prg) || ProgressMeter.finish!(prg)
    println("events per spawn: ", run/runs, " spawns: ", runs)
    println("wakeups per round: ", run2/runs2, " rounds: ", runs2)

    Ξ, (t, x, θ), (acc, num)
end
function parallel_spdmp_outer!(tmin, t′, T, task, waitfor, latch, wakeup, evtime, perm, res, Ξ, events, partition, inner, G, G1, G2, ∇ϕ, t, x, θ, Q,
    c, b, t_old, F, (Δ, factor, adapt), prg, args...)
    acc = num = 0
    run = runs = 0
    run2 = runs2 = 0
    stops = ismissing(prg) ? 0 : max(prg.n - 1, 0) # allow one stop for cleanup
    tstop = T/stops
    rng = Rng()
    while tmin < T
        if latch.active[] !== 0
            VERBOSE && println("Waiting $tmin")
            lock(latch.condition) do
                while latch.active[] != 0
                    wait(latch.condition)
                end
            end
        end
        for ti in each(partition)
            if waitfor[ti] == 0
                run += res[ti][][end]
                runs += 1
                #println(res[ti][end])
                append!(Ξ.events, events[ti])
                resize!(events[ti], 0)
                evtime[ti] = res[ti][][2]
            end
        end
        sortperm!(perm, evtime, alg=InsertionSort)
        alldone = true
        for ti in perm
            i, t′_, acc_, num_ = res[ti][]
            if waitfor[ti] == 0
                num += num_
                acc += acc_
                t′[ti] = t′_
            end

            waitfor[ti] = 0
            for j in G[i][2]
                i == j && continue
                q1, q2 = partition(j)
                if t′[q1] < t′_
                    waitfor[ti] = i
                end
            end

            waitfor[ti] != 0 && continue

            success = parallel_innermost!(rng, partition, G, G1, G2, ∇ϕ, i, t, x, θ, t′_, Q, c, b, t_old, F, (factor, adapt), args...)

            if success
                acc += 1
                push!(Ξ, event(i, t, x, θ, F))
            end
        end

        tmin = minimum(t′)
        runs2 += 1
        if tmin > tstop
            tstop += T/stops
            next!(prg)
        end
        for ti in each(partition)
            if waitfor[ti] .== 0 ||  tmin >= T
                run2 += 1
                Threads.atomic_or!(latch.active, UInt(1) << (ti - 1))
                lock(wakeup[ti]) do
                    notify(wakeup[ti], tmin >= T)
                end
            end
        end

        tmin >= T && return (acc, num), (run, runs),  (run2, runs2)

    end
end