text:
## -----------  Access

own_name(o :: ModelObject) = own_name(get_object_id(o));

"""
	$(SIGNATURES)

Is `o` a `ModelObject`?
This can be used to apply all of the code of this package to objects that are not sub-types of `ModelObject`. But this is not fully implemented.
"""
is_model_object(o :: ModelObject) = true;
is_model_object(o) = false;

"""
	$(SIGNATURES)

Retrieve an objects ObjectId. Return `nothing` if not found.
"""
get_object_id(o :: ModelObject) = o.objId :: ObjectId
get_object_id(o) = nothing;

## ------------  Collect objects inside a ModelObject

"""
	$(SIGNATURES)

Collect all ObjectId's from an object as a `Vector{ObjectId}`.
"""
collect_object_ids(o :: ModelObject) =
    get_object_id.(collect_model_objects(o));

"""
    $(SIGNATURES)

Collect all model objects inside an object, including the object itself.
Only those that satisfy `is_model_object`.
Recursive. Also collects objects inside child objects and so on.
Returns empty `Vector` if no objects found.

# Arguments
- flatten: If `true`, return a single `Vector` with all objects. If `false`, return a nested `Vector`.

# Example
```julia
struct O1
    c1
    c2
end
v = collect_model_objects(o; flatten = true);
v == [o, o.c1, o.c1.c1, o.c2, o.c2.c1, o.c2.c2];
# Each element contains an objects and all its children (and their children).
v = collect_model_objects(o; flatten = false);
v[1] == o;
v[2] == [o.c1, o.c1.c1]
v[3] == [o.c2, o.c2.c1, o.c2.c2]
```
"""
function collect_model_objects(o :: ModelObject; flatten :: Bool = true)
    outV = Vector{Any}();
    if is_model_object(o)
        push!(outV, o);
    end

    # Objects directly contained in `o`
    childObjV = get_child_objects(o);
    if !Base.isempty(childObjV)
        for i1 = 1 : length(childObjV)
            nestedObjV = collect_model_objects(childObjV[i1]);
            if flatten
                append!(outV, nestedObjV);
            else
                push!(outV, nestedObjV);
            end
        end
    end
    return outV :: Vector
end

collect_model_objects(o; flatten :: Bool = true) = Vector{Any}();

"""
    $(SIGNATURES)

Find the child objects inside a model object.
Returns empty Vector if no objects found.
"""
function get_child_objects(o :: ModelObject)
    childV = Vector{Any}();
    for pn in propertynames(o)
        @assert isdefined(o, pn)  "$pn undefined in $o"
        obj = getproperty(o, pn);
        if isa(obj, Vector)  &&  !isempty(obj)
            # This check is not quite right. But objects should all be the same type.
            if is_model_object(obj[1])
                append!(childV, obj);
            end
        else
            if is_model_object(obj)
                push!(childV, obj);
            end
        end
    end
    return childV :: Vector
end

get_child_objects(o) = Vector{Any}();

# """
# 	$(SIGNATURES)

# Model objects as nested Vector.
# """
# function model_object_dict(o :: ModelObject)
#     # Objects directly contained in `o`
#     childObjV = get_child_objects(o);
#     if !Base.isempty(childObjV)
#         for i1 = 1 : length(childObjV)
#             nestedObjV = model_object_dict(childObjV[i1]);
#             push! = nestedObjV;
#         end
#     end
#     return d
# end

# model_object_dict(o) = Dict{Any, Any}

"""
	$(SIGNATURES)

Find child object with a given `ObjectId`. Expects that all ModelObjects have ObjectIds.
Returns `nothing` if not found.
"""
function find_object(o :: ModelObject, id :: ObjectId)
    oOut = nothing;
    objV = collect_model_objects(o);
    if !isempty(objV)
        for obj in objV
            if isequal(get_object_id(obj), id)
                oOut = obj;
                break;
            end
        end
    end
    return oOut
end

# find_object(o :: ModelObject, name :: Symbol) = find_object(o)

find_object(o, id) = nothing;

"""
	$(SIGNATURES)

Find all child objects that have a name given by a `Symbol`. Easier than having to specify an entire `ObjectId`.
"""
function find_object(o :: ModelObject, oName :: Symbol)
    outV = Vector{Any}();
    objV = collect_model_objects(o);
    if !isempty(objV)
        for obj in objV
            if own_name(obj) == oName
                push!(outV, obj);
            end
        end
    end
    return outV
end

find_object(o, oName :: Symbol) = Vector{Any}();

"""
	$(SIGNATURES)

Find the only object that matches the name `oName`. Errors if not found.
"""
function find_only_object(o :: ModelObject, oName :: Symbol)
    return only(find_object(o, oName));
end

find_only_object(o, oName :: Symbol) = nothing;

"""
	$(SIGNATURES)

Retrieve the value of a field in a `ModelObject` or its children.
Object name `oName` must be unique. This is the name in the `ObjectId`, not what the object is called as a field in the model object.
"""
function get_value(x :: ModelObject, oName :: Symbol, pName :: Symbol)
    objV = find_object(x, oName);
    @assert length(objV) == 1  "Found $(length(objV)) matches for $oName / $pName"
    return getfield(objV[1], pName)
end

## -----------  Show

"""
	$(SIGNATURES)

Show structure of a `ModelObject`.
"""
function Base.show(io :: IO,  o :: ModelObject)
    show_object_structure(o; io = io);
end

"""
	$(SIGNATURES)

Show the "tree" structure of a ModelObject. Returns a Vector of formatted strings that show objects, children, grandchildren with their types.
"""
function show_object_structure(o; io :: IO = stdout)
    idV, typeV = object_structure(o);
    if !isempty(idV)
        for j = 1 : length(idV)
            println(io,  indent_string(n_parents(idV[j])),  own_name(idV[j]),
                "\t",  typeV[j]);
        end
    end
    return nothing
end

indent_string(n) = "   " ^ n;

"""
	$(SIGNATURES)

Represent the "tree" structure of a ModelObject.

Returns ObjectIds and DataTypes.
"""
function object_structure(o)
    objV = collect_model_objects(o);
    idV = Vector{ObjectId}();
    typeV = Vector{DataType}();
    if !isempty(objV)
        for obj in objV
            # (id = objId, nparents = n_parents(objId))
            push!(idV, get_object_id(obj))
            push!(typeV, typeof(obj));
        end
    end
    return idV, typeV
end

# -------------


text:
function pagerank(
    g::AbstractGraph{U},
    α=0.85,
    n=100::Integer,
    ϵ=1.0e-6
    ) where U <: Integer

    # indegree(g, v) is estimated run-time to iterate over inneighbors(g, v)
    partitions = LightGraphs.optimal_contiguous_partition(indegree(g), nthreads(), nv(g))

    α_div_outdegree = Vector{Float64}(undef,nv(g))
    dangling_nodes = Vector{U}()
    @inbounds for v in vertices(g)
        if outdegree(g, v) == 0
            push!(dangling_nodes, v)
        end
        α_div_outdegree[v] = (α/outdegree(g, v))
    end

    nvg = Int(nv(g))
    # solution vector and temporary vector
    x = fill(1.0 / nvg, nvg)
    xlast = copy(x)
    @inbounds for _ in 1:n
        dangling_sum = 0.0
        for v in dangling_nodes
            dangling_sum += x[v]
        end
        # flow from teleprotation
        y = (1 - α + α * dangling_sum) * (1.0 / nvg)
        xlast .= y
        # flow from edges
        let x = x
            @threads for v_set in partitions
                for v in v_set
                    for u in inneighbors(g, v)
                        xlast[v] += (x[u] * α_div_outdegree[u])
                    end
                end
            end
        end

        # l1 change in solution convergence criterion
        err = 0.0
        for v in vertices(g)
            err += abs(xlast[v] - x[v])
            x[v] = xlast[v]
        end
        if (err < nvg * ϵ)
            return x
        end
    end
    error("Pagerank did not converge after $n iterations.") # TODO 0.7: change to InexactError with appropriate msg.
end

ext: jl


text:
"""
    Path{D}
Alias of `MutableLinkedList{Dims{D}}`. Used to represent the path to be
taken by an agent in a `D` dimensional [`GridSpace`](@ref).
"""
const Path{D} = MutableLinkedList{Dims{D}}

struct Pathfinder{W<:Union{BitArray,Nothing},M<:Union{CostMetric,Nothing}}
    diagonal_movement::Bool
    admissibility::Float64
    walkable::W
    cost_metric::M
end

"""
    Pathfinding.Pathfinder(; kwargs...)

Enable pathfinding using the A* algorithm by passing an instance of `Pathfinder` into
[`GridSpace`](@ref). Pathfinding works by using the functions
[`Pathfinding.set_target!`](@ref) and [`move_along_route!`](@ref) see [`Pathfinding`](@ref)
for more.

## Keywords

* `diagonal_movement = true` states that agents are allowed to move diagonally.
    Otherwise, only orthogonal directions are possible.
* `admissibility = 0` allows the algorithm to approximate paths to speed up pathfinding
    significantly. A value of `admissibility` allows paths at most `(1+admissibility)` times
    the optimal path length.
* `walkable = nothing` specifies (un)walkable regions of the space. If specified, it should
    be a `BitArray` array of the same size as the corresponding [`GridSpace`](@ref). This defaults
    to `nothing`, which allows agents to walk on any position in the space. An example usage can
    be found in [Maze Solver](@ref).
* `cost_metric` is an instance of a cost metric and specifies the method
    to use for approximating the distance between two points. This defaults
    to [`Pathfinding.DirectDistance`](@ref) with appropriate dimensionality.
"""
function Pathfinder(;
    diagonal_movement::Bool = true,
    admissibility::Float64 = 0.0,
    walkable::W = nothing,
    cost_metric::M = nothing,
) where {W<:Union{BitArray,Nothing},M<:Union{CostMetric,Nothing}}
    return Pathfinder(diagonal_movement, admissibility, walkable, cost_metric)
end

struct AStar{D,P,M}
    agent_paths::Dict{Int,Path{D}}
    grid_dims::Dims{D}
    neighborhood::Vector{CartesianIndex{D}}
    admissibility::Float64
    walkable::BitArray{D}
    cost_metric::CostMetric{D}

    function AStar{D,P,M}(
        agent_paths::Dict,
        grid_dims::Dims{D},
        neighborhood::Vector{CartesianIndex{D}},
        admissibility::Float64,
        walkable::BitArray{D},
        cost_metric::CostMetric{D},
    ) where {D,P,M}
        @assert admissibility >= 0 "Invalid value for admissibility: $admissibility ≱ 0"
        if typeof(cost_metric) == HeightMap{D}
            @assert size(cost_metric.hmap) == grid_dims "Heightmap dimensions must be same as provided space"
        elseif typeof(cost_metric) == DirectDistance{D}
            if M
                @assert length(cost_metric.direction_costs) >= D "DirectDistance direction_costs must have as many values as dimensions"
            else
                @assert length(cost_metric.direction_costs) >= 1 "DirectDistance direction_costs must have non-zero length"
            end
        end
        new(agent_paths, grid_dims, neighborhood, admissibility, walkable, cost_metric)
    end
end

"""
    AStar(dims::Dims{D}; kwargs...)
Provides pathfinding capabilities and stores agent paths. Not part of public API.

`periodic = false` specifies if the space is periodic

`diagonal_movement = true` specifies if movement can be to diagonal neighbors of a
tile, or only orthogonal neighbors.

`admissibility = 0.0` specifies how much a path can deviate from optimal, in favour of
faster pathfinding. Admissibility (`ε`) of `0.0` will always find the optimal path.
Larger values of `ε` will explore fewer nodes, returning a path length with at most
`(1+ε)` times the optimal path length.

`walkable = trues(size(space))` is used to specify (un)walkable positions of the
space. All positions are assumed to be walkable by default.

`cost_metric = DirectDistance` specifies the metric used to approximate the distance
between any two walkable points on the grid. See [`CostMetric`](@ref).

Example usage in [Maze Solver](@ref) and [Mountain Runners](@ref).
"""
function AStar(
    dims::Dims{D};
    periodic::Bool = false,
    diagonal_movement::Bool = true,
    admissibility::Float64 = 0.0,
    walkable::BitArray{D} = trues(dims),
    cost_metric::CostMetric{D} = DirectDistance{D}(),
) where {D}
    neighborhood = diagonal_movement ? moore_neighborhood(D) : vonneumann_neighborhood(D)
    return AStar{D,periodic,moore_neighbors}(
        Dict{Int,Path{D}}(),
        dims,
        neighborhood,
        admissibility,
        walkable,
        cost_metric,
    )
end

function AStar(dims::Dims{D}, periodic::Bool, pathfinder::Pathfinder) where {D}
    walkable = pathfinder.walkable === nothing ? trues(dims) : pathfinder.walkable

    metric =
        isnothing(pathfinder.cost_metric) ? DirectDistance{D}() : pathfinder.cost_metric

    neighborhood =
        pathfinder.diagonal_movement ? moore_neighborhood(D) : vonneumann_neighborhood(D)
    return AStar{D,periodic,pathfinder.diagonal_movement}(
        Dict{Int,Path{D}}(),
        dims,
        neighborhood,
        pathfinder.admissibility,
        walkable,
        metric,
    )
end

moore_neighborhood(D) = [
    CartesianIndex(a)
    for a in Iterators.product([-1:1 for φ in 1:D]...) if a != Tuple(zeros(Int, D))
]

function vonneumann_neighborhood(D)
    hypercube = CartesianIndices((repeat([-1:1], D)...,))
    [β for β ∈ hypercube if LinearAlgebra.norm(β.I) == 1]
end

function Base.show(io::IO, pathfinder::AStar{D,P,M}) where {D,P,M}
    periodic = P ? "periodic, " : ""
    moore = M ? "diagonal, " : "orthogonal, "
    s =
        "A* in $(D) dimensions, $(periodic)$(moore)ϵ=$(pathfinder.admissibility), " *
        "metric=$(pathfinder.cost_metric)"
    print(io, s)
end

"""
    position_delta(pathfinder::AStar{D}, from::NTuple{Int,D}, to::NTuple{Int,D})
Returns the absolute difference in coordinates between `from` and `to` taking into account
periodicity of `pathfinder`.
"""
position_delta(pathfinder::AStar{D,true}, from::Dims{D}, to::Dims{D}) where {D} =
    min.(abs.(to .- from), pathfinder.grid_dims .- abs.(to .- from))

position_delta(pathfinder::AStar{D,false}, from::Dims{D}, to::Dims{D}) where {D} =
    abs.(to .- from)

"""
    delta_cost(pathfinder::AStar{D}, metric::M, from, to) where {M<:CostMetric}
Calculate an approximation for the cost of travelling from `from` to `to` (both of
type `NTuple{N,Int}`. Expects a return value of `Float64`.
"""
function delta_cost(
    pathfinder::AStar{D,periodic,true},
    metric::DirectDistance{D},
    from::Dims{D},
    to::Dims{D},
) where {D,periodic}
    delta = collect(position_delta(pathfinder, from, to))

    sort!(delta)
    carry = 0
    hdist = 0
    for i in D:-1:1
        hdist += metric.direction_costs[i] * (delta[D+1-i] - carry)
        carry = delta[D+1-i]
    end
    return hdist
end

function delta_cost(
    pathfinder::AStar{D,periodic,false},
    metric::DirectDistance{D},
    from::Dims{D},
    to::Dims{D},
) where {D,periodic}
    delta = position_delta(pathfinder, from, to)

    return sum(delta) * metric.direction_costs[1]
end

delta_cost(
    pathfinder::AStar{D},
    metric::MaxDistance{D},
    from::Dims{D},
    to::Dims{D},
) where {D} = max(position_delta(pathfinder, from, to)...)

delta_cost(
    pathfinder::AStar{D},
    metric::HeightMap{D},
    from::Dims{D},
    to::Dims{D},
) where {D} =
    delta_cost(pathfinder, metric.base_metric, from, to) +
    abs(metric.hmap[from...] - metric.hmap[to...])

delta_cost(pathfinder::AStar{D}, from::Dims{D}, to::Dims{D}) where {D} =
    delta_cost(pathfinder, pathfinder.cost_metric, from, to)

struct GridCell
    f::Int
    g::Int
    h::Int
end

GridCell(g::Int, h::Int, admissibility::Float64) =
    GridCell(round(Int, g + (1 + admissibility) * h), g, h)

GridCell() = GridCell(typemax(Int), typemax(Int), typemax(Int))

"""
    find_path(pathfinder::AStar{D}, from::NTuple{D,Int}, to::NTuple{D,Int})
Calculate the shortest path from `from` to `to` using the A* algorithm.
If a path does not exist between the given positions, an empty linked list is returned.
"""
function find_path(pathfinder::AStar{D}, from::Dims{D}, to::Dims{D}) where {D}
    grid = Dict{Dims{D},GridCell}()
    parent = Dict{Dims{D},Dims{D}}()

    open_list = MutableBinaryMinHeap{Tuple{Int,Dims{D}}}()
    open_list_handles = Dict{Dims{D},Int64}()
    closed_list = Set{Dims{D}}()

    grid[from] = GridCell(0, delta_cost(pathfinder, from, to), pathfinder.admissibility)
    push!(open_list, (grid[from].f, from))

    while !isempty(open_list)
        _, cur = pop!(open_list)
        cur == to && break
        push!(closed_list, cur)

        nbors = get_neighbors(cur, pathfinder)
        for nbor in Iterators.filter(n -> inbounds(n, pathfinder, closed_list), nbors)
            nbor_cell = haskey(grid, nbor) ? grid[nbor] : GridCell()
            new_g_cost = grid[cur].g + delta_cost(pathfinder, cur, nbor)

            if new_g_cost < nbor_cell.g
                parent[nbor] = cur
                grid[nbor] = GridCell(
                    new_g_cost,
                    delta_cost(pathfinder, nbor, to),
                    pathfinder.admissibility,
                )
                if haskey(open_list_handles, nbor)
                    update!(open_list, open_list_handles[nbor], (grid[nbor].f, nbor))
                else
                    open_list_handles[nbor] = push!(open_list, (grid[nbor].f, nbor))
                end
            end
        end
    end

    agent_path = Path{D}()
    cur = to
    while true
        haskey(parent, cur) || break
        pushfirst!(agent_path, cur)
        cur = parent[cur]
    end
    return agent_path
end

@inline get_neighbors(cur, pathfinder::AStar{D,true}) where {D} =
    (mod1.(cur .+ β.I, pathfinder.grid_dims) for β in pathfinder.neighborhood)
@inline get_neighbors(cur, pathfinder::AStar{D,false}) where {D} =
    (cur .+ β.I for β in pathfinder.neighborhood)
@inline inbounds(n, pathfinder, closed) =
    all(1 .<= n .<= pathfinder.grid_dims) && pathfinder.walkable[n...] && n ∉ closed

"""
    Pathfinding.set_target!(agent, target::NTuple{D,Int}, model)
Calculate and store the shortest path to move the agent from its current position to
`target` (a grid position e.g. `(1, 5)`) for models using [`Pathfinding`](@ref).

Use this method in conjuction with [`move_along_route!`](@ref).
"""
function set_target!(
    agent::A,
    target::Dims{D},
    model::ABM{<:GridSpace{D,P,<:AStar{D}},A},
) where {D,P,A<:AbstractAgent}
    model.space.pathfinder.agent_paths[agent.id] =
        find_path(model.space.pathfinder, agent.pos, target)
end

"""
    Pathfinding.set_best_target!(agent, targets::Vector{NTuple{D,Int}}, model)

Calculate and store the best path to move the agent from its current position to
a chosen target position taken from `targets` for models using [`Pathfinding`](@ref).

The `condition = :shortest` keyword retuns the shortest path which is shortest
(allowing for the conditions of the models pathfinder) out of the possible target
positions. Alternatively, the `:longest` path may also be requested.

Returns the position of the chosen target.
"""
function set_best_target!(
    agent::A,
    targets::Vector{Dims{D}},
    model::ABM{<:GridSpace{D,P,<:AStar{D}},A};
    condition::Symbol = :shortest,
) where {D,P,A<:AbstractAgent}
    @assert condition ∈ (:shortest, :longest)
    compare = condition == :shortest ? (a, b) -> a < b : (a, b) -> a > b
    best_path = Path{D}()
    best_target = nothing
    for target in targets
        path = find_path(model.space.pathfinder, agent.pos, target)
        if isempty(best_path) || compare(length(path), length(best_path))
            best_path = path
            best_target = target
        end
    end

    model.space.pathfinder.agent_paths[agent.id] = best_path
    return best_target
end

Agents.is_stationary(
    agent::A,
    model::ABM{<:GridSpace{D,P,<:AStar{D}},A},
) where {D,P,A<:AbstractAgent} = isempty(agent.id, model.space.pathfinder)

Base.isempty(id::Int, pathfinder::AStar) =
    !haskey(pathfinder.agent_paths, id) || isempty(pathfinder.agent_paths[id])

"""
    Pathfinding.heightmap(model)
Return the heightmap of a [`Pathfinding.Pathfinder`](@ref) if the
[`Pathfinding.HeightMap`](@ref) metric is in use, `nothing` otherwise.

It is possible to mutate the map directly, for example
`Pathfinding.heightmap(model)[15, 40] = 115`
or `Pathfinding.heightmap(model) .= rand(50, 50)`. If this is mutated,
a new path needs to be planned using [`Pathfinding.set_target!`](@ref).
"""
function heightmap(model::ABM{<:GridSpace{D,P,<:AStar{D}}}) where {D,P}
    if model.space.pathfinder.cost_metric isa HeightMap
        return model.space.pathfinder.cost_metric.hmap
    else
        return nothing
    end
end

"""
    Pathfinding.walkmap(model)
Return the walkable map of a [`Pathfinding.Pathfinder`](@ref).

It is possible to mutate the map directly, for example
`Pathfinding.walkmap(model)[15, 40] = false`.
If this is mutated, a new path needs to be planned using [`Pathfinding.set_target!`](@ref).
"""
walkmap(model::ABM{<:GridSpace{D,P,<:AStar{D}}}) where {D,P} =
    model.space.pathfinder.walkable

"""
    move_along_route!(agent, model_with_pathfinding)
Move `agent` for one step along the route toward its target set by [`Pathfinding.set_target!`](@ref)
for agents on a [`GridSpace`](@ref) using a [`Pathfinding.Pathfinder`](@ref).
If the agent does not have a precalculated path or the path is empty, it remains stationary.
"""
function Agents.move_along_route!(
    agent::A,
    model::ABM{<:GridSpace{D,P,<:AStar{D}},A},
) where {D,P,A<:AbstractAgent}
    isempty(agent.id, model.space.pathfinder) && return

    move_agent!(agent, first(model.space.pathfinder.agent_paths[agent.id]), model)
    popfirst!(model.space.pathfinder.agent_paths[agent.id])
end

function Agents.kill_agent!(
    agent::A,
    model::ABM{<:GridSpace{D,P,<:AStar{D}},A},
) where {D,P,A<:AbstractAgent}
    delete!(model.space.pathfinder.agent_paths, agent.id)
    delete!(model.agents, agent.id)
    Agents.remove_agent_from_space!(agent, model)
end


text:
# # A simple FEM application
#
# ## Introduction
# The Finite Element Method is widely used to solve PDEs in Engineering applications and particularly in Structural Analysis problems [Add citation Bathe]. A specific case os structures is Truss structures, where trusses or bars are connected but not welded. Truss models are considered during the conceptual design of bridges or other structures.
#
# The stiffness matrix of a truss element in the local coordinate system is given by
# ```math
# K_L = s
# \left(
#  \begin{matrix}
#  1 & 0 & -1 & 0 \\
#  0 & 0 &  0 & 0 \\
# -1 & 0 &  1 & 0 \\
#  0 & 0 &  0 & 0
# \end{matrix}
# \right),
# ```
#
# where $s$ is the stiffness, given by $EA/L$, with $E$ being the Young modulus, $A$ the area of the cross-section and $L$ the length of that truss element.
#
# The change-of-basis matrix is given by
# ```math
# _G(Q)_L = Q =
# \left(
#   \begin{matrix}
# \cos(\alpha) & -\sin(\alpha) & 0 & 0 \\
#  sin(alpha) & cos(alpha) &  0 & 0 \\
#  0 & 0 &  cos(alpha) & sin(alpha) \\
#  0 & 0 &  sin(alpha) & cos(alpha)
# \end{matrix}
# \right),
# ```
#
# The system of equations for each element is written in local coordinates as
# ```math
# K_L d_L = f_L
# ```
# and using the change-of-basis we obtain
# ```math
# K_G d_G = f_G \qquad K_G = Q K_L Q^T
# ```
#
# The stiffness matrix can be computed using the following function, which computes the unitary stiffness matrix for an element defined by the coordinates of its first and second nodes.
function unitaryStiffnessMatrix( coordFirstNode, coordSecondNode  )
  diff      = (coordSecondNode - coordFirstNode)
  length   = sqrt( diff'*diff )
  c        = diff[1] / length ;
  s        = diff[2] / length ;
  Qloc2glo = [ c -s 0 0 ; s c 0 0 ; 0 0 c -s ; 0 0 s c ] ;
  Kloc     = [ 1 0 -1 0 ; 0 0 0 0 ; -1 0 1 0 ; 0 0 0 0 ] ;
  Kglo     = Qloc2glo * Kloc * transpose(Qloc2glo)
  return     Kglo, length
end
#
# ## Problem with fixed parameters
# In this section, a problem based on Example 4.1 from [https://github.com/JuliaIntervals/IntervalLinearAlgebra.jl/files/7271616/skalna2006.pdf] is considered. The following diagram shows the truss structure considered.
#
# \fig{../../assets/trussDiagram.png}
#
# The scalar parameters considered are given by
E = 2e11 ; # Young modulus
A = 5e-3 ; # Cross-section area
# while the coordinate matrix is given by
nodesCMatrix = [ 0. 0. ;
                 1. 1. ;
                 2. 0. ;
                 3. 1. ;
                 4. 0. ]
# and connectivity matrix is given by
connecMatrix = [ 1 2 ;
                 1 3 ;
                 2 3 ;
                 2 4 ;
                 3 4 ;
                 3 5 ;
                 4 5 ]
# and the fixed degrees of freedom (supports) are
fixedDofs     = [2 9 10 ]

# calculations
numNodes = size( nodesCMatrix )[1]
numElems = size( connecMatrix )[1]
freeDofs = zeros(Int8, 2*numNodes-length(fixedDofs))
indDof  = 1 ;
counter = 0 ;
while indDof <= (2*numNodes)
  if !(indDof in fixedDofs)
    global counter = counter + 1 ;
    freeDofs[ counter ] = indDof ;
    print(indDof)
  end
  global indDof = indDof + 1 ;
end
print(freeDofs)
KG = zeros( 2*numNodes, 2*numNodes ) ;
FG = zeros( 2*numNodes )

# assembly
for elem in 1:numElems
  print(" assembling stiffness matrix of element ", elem , "\n")
  indexFirstNode  = connecMatrix[ elem, 1 ]
  indexSecondNode = connecMatrix[ elem, 2 ]
  dofsElem = [2*indexFirstNode-1 2*indexFirstNode 2*indexSecondNode-1 2*indexSecondNode ]
  KGelem, lengthElem = unitaryStiffnessMatrix( nodesCMatrix[ indexSecondNode, : ], nodesCMatrix[ indexFirstNode, : ] )
  stiffnessParam = E * A / lengthElem ;
  for i in 1:4
    for j in 1:4
      KG[ dofsElem[i], dofsElem[j] ] = KG[ dofsElem[i], dofsElem[j] ] + stiffnessParam * KGelem[i,j]
    end
  end
end
FG[4] = -1e4 ;
KG = KG[ freeDofs, : ]
KG = KG[ :, freeDofs ]
FG = FG[ freeDofs ]

u = KG \ FG
print(u)

# ## Problem with interval parameters


text:
# we provide support for two kinds of source.  the difference is important
# mainly because of how the source interacts with the regexp matcher, which
# can be critical when trying to get best performance with large files.

# however, the parser code will interact with "any" source provided it
# provides the API here.

# utilities

function fmt_error(line, col, text, msg)
    arrow = string(repeat(" ", max(col-1, 0)), "^")
    "$(msg) at ($(line),$(col))\n$(text)\n$(arrow)\n"
end

# strings indexed by offset.  in this case, the type of the source is simply
# the string itself, so we do not need to implement the iter protocol, since
# that already exists.  but we do need some extra functions.

# given a source and iterator, return a message identifying line and
# column

nl(x) = x == '\n'

function round_up(s, i)
    i = max(i, 1)
    nl(s[i]) ? i+1 : i
end

function round_down(s, i)
    i = i == 0 ? lastindex(s) : i
    nl(s[i]) ? i-1 : i
end

function diagnostic(s::AbstractString, i, msg)
    if i < 1
        l, c, t = 0, 0, "[Before start]"
    elseif i > length(s)
        l, c, t = count(nl, s)+2, 0, "[After end]"
    else
        l = count(nl, SubString(s, 1, max(1, i-1))) + 1
        # excluding new line
        prev_index = findprev(nl, s, max(1, i-1))
        p = round_up(s, something(prev_index, 0))
        # end of line, excluding newline
        next_index = findnext(nl, s, i)
        q = round_down(s, something(next_index, 0))
        t = SubString(s, p, q)
        c = i - p + 1
    end
    fmt_error(l, c, t, msg)
end

# given a source and iterator, return following text for regexp
forwards(s::AbstractString, i) = SubString(s, i)

# originally i thought that UTF8 would require code to move forwards
# the given number of characters.  but discard() is called only from
# regex and that returns the number of *bytes*.
discard(::AbstractString, i, n) = i + n

# io instance, from which lines are read.  the lines are stored so that we
# have backtracking.  this is pretty pointless, because you might as well just
# read the whole thing into memory at the start, but it serves as the basis
# for a similar approach in try.jl where not all lines are saved.

# since it's so pointless as it is, we thrown in an extra feature - a maximum
# depth (in lines) after which data are discarded.  it might be useful
# somewhere.

# all the below is based on line_at()
abstract type LineAt end

mutable struct LineSource{S}<:LineAt
    io::IO
    zero::Int      # offset to lines (lines[x] contains line x+zero)
    limit::Int     # maximum number of lines
    lines::Vector{S}
    LineSource(io::IO, line::S; limit=-1) where {S} = new{S}(io, 0, limit, S[line])
end

function LineSource(io::IO; limit=-1)
    line = readline(io, keep=true)
    LineSource(io, line; limit=limit)
end

LineSource(s::S; limit=-1) where {S<:AbstractString} = LineSource(IOBuffer(s); limit=limit)

struct LineException<:FailureException end

@auto_hash_equals struct LineIter
    line::Int
    column::Int
end

# used in debug

isless(a::LineIter, b::LineIter) = a.line < b.line || (a.line == b.line && a.column < b.column)

# iter protocol

function line_at(s::LineSource, i::LineIter; check::Bool=true)
    if check && i.line <= s.zero || i.column < 1
        throw(LineException())
    end

    n = i.line - s.zero
    while length(s.lines) < n
        push!(s.lines, readline(s.io, keep=true))
    end
    while s.limit > 0 && length(s.lines) > s.limit
        s.zero += 1
        pop!(s.lines)
    end
    line = s.lines[i.line - s.zero]
    if check && i.column > length(line)
        throw(LineException())
    end
    return line
end

function iterate(s::LineSource, i::LineIter=LineIter(1,1))
    # i.line is always correct for further reading (eg via forwards())
    line = line_at(s, i, check=false)
    if iterate(line, i.column) == nothing && eof(s.io)
        return nothing
    end
    line = line_at(s, i)
    c, column = iterate(line, i.column)
    if iterate(line, column) == nothing && eof(s.io)
        return c, LineIter(i.line+1, 1)
    else
        return c, LineIter(i.line, column)
    end
end

firstindex(s::LineSource) = LineIter(1,1)

# other api

function diagnostic(s::LineAt, i::LineIter, msg)
    line = "[Not available]"
    try
        line = line_at(s, i)
        if nl(line[end])
            line = line[1:end-1]
        end
    catch x
        if !isa(x, FailureException)
            throw(x)
        end
    end
    fmt_error(i.line, i.column, line, msg)
end

# regexp only works within the current line
function forwards(s::LineAt, i::LineIter)
    line = line_at(s, i; check=false)
    if iterate(line, i.column) == nothing && eof(s.io)
        return ""
    end
    return SubString(line_at(s, i), i.column)
end

function discard(s::LineAt, i::LineIter, n)
    while n > 0 && iterate(s, i) != nothing
        l = line_at(s, i)
        available = length(l) - i.column + 1
        if n < available
            i = LineIter(i.line, i.column+n)
            n = 0
        else
            i = LineIter(i.line+1, 1)
            n -= available
        end
    end
    i
end


text:
# A motivating example
module Football
    # My own array-like type
    struct Scores <: AbstractVector{Float64}
        values::Vector{Float64}
    end

    # implement AbstractArray interface
    Base.size(s::Scores) = (length(s.values),)
    Base.getindex(s::Scores, i::Int) = s.values[i]
end

#--- Type 1: Redefining an existing function from a different package
#            to different logic!

# Let's crash it immediately by typing this into the REPL!
Base.:(+)(x::Int, y::Int) = "hello world"

#=
julia> Base.:(+)(x::Int, y::Int) = "hello world"

┌ Error: Error in the keymap
│   exception =
│    MethodError: no method matching Int64(::String)
│    Closest candidates are:
│      Int64(::Union{Bool, Int32, Int64, UInt32, UInt64, UInt8, Int128, Int16, Int8, UInt128, UInt16}) at boot.jl:710
│      Int64(::Ptr) at boot.jl:720
│      Int64(::Float32) at float.jl:700
│      ...
=#

#--- Type 2: Extending an existing function from a different package with types that you don't own

# In JavaScript, we can add strings to numbers! e.g. "1" + 2 = "12"
module MyModule
    import Base.+
    (+)(s::AbstractString, n::Number) = "$s$n"
end

using .MyModule
"1" + 2
#=
julia> module MyModule
           import Base.+
           (+)(s::AbstractString, n::Number) = "$s$n"
       end
Main.MyModule

julia> using .MyModule

julia> "1" + 2
"12"
=#

# Another package author thought it would be cool to handle strings
# in math operations
module AnotherModule
    import Base: +, -, *, /
    (+)(s::AbstractString, n::T) where T <: Number = parse(T, s) + n
    (-)(s::AbstractString, n::T) where T <: Number = parse(T, s) - n
    (*)(s::AbstractString, n::T) where T <: Number = parse(T, s) * n
    (/)(s::AbstractString, n::T) where T <: Number = parse(T, s) / n
end

using .AnotherModule
"1" + 2
#=
julia> module AnotherModule
           import Base: +, -, *, /
           (+)(s::AbstractString, n::T) where T <: Number = parse(T, s) + n
           (-)(s::AbstractString, n::T) where T <: Number = parse(T, s) - n
           (*)(s::AbstractString, n::T) where T <: Number = parse(T, s) * n
           (/)(s::AbstractString, n::T) where T <: Number = parse(T, s) / n
       end
Main.AnotherModule

julia> using .AnotherModule

julia> "1" + 2
3

julia> "2" * 3
6

julia> "1" + 2.3
3.3
=#

# Avoiding Type 2 piracy - create your own type!

module MyModule
    export @str_str
    import Base: +, show

    struct MyString
        value::AbstractString
    end

    macro str_str(s::AbstractString)
        MyString(s)
    end

    show(io::IO, s::MyString) = print(io, s.value)
    (+)(s::MyString, n::Number) = MyString(s.value * string(n))
    (+)(n::Number, s::MyString) = MyString(string(n) * s.value)
    (+)(s::MyString, t::MyString) = MyString(s.value * t.value)
end

using .MyModule
str"I am " + 25 + str" years old!"

#--- Type 3: Extending an existing function from a different package
#            for your own type but has a different meaning
# aka "puns"
# A Party just contains a title and guest names
struct Party
    title::String
    guests::Vector{String}
end

# constructor
Party(title) = Party(title, String[])

# idea: extending from Base to get around with name conflict issue
Base.join(name::String, party::Party) = push!(party.guests, name)

party = Party("Halloween 2019")
join("Tom", party)
join("Kevin", party)
#=
julia> join("Tom", party)
1-element Array{String,1}:
 "Tom"

julia> join("Kevin", party)
2-element Array{String,1}:
 "Tom"
 "Kevin"
=#

# what about this bad code?
# it has bad results without warning due to duck typing.
join(["Bob", "Jeff"], party)
#=
julia> join(["Bob", "Jeff"], party)
"BobParty(\"Halloween 2019\", [\"Tom\", \"Kevin\"])Jeff"
=#


text:
"""
`binary_boundary_search(fun, x_range, err)`

Find a boundary value `X` at which `fun(X)` changes its value from `0` to `1`.
That is, `fun(x) = 0` for all `x < X` and `fun(x) = 1` for all `x > X`.

**Arguments**
- `fun`: Function that returns `0` if `x < X` and `1` if `x > X`.
- `x_range`: Search range as `[x_min, x_max]`.
- `err`: Maximal deviation of the result from the real value.

**Returns**
- `X`: Boundary value, at which `fun(x)` changes from `0` to `1`.
"""
function binary_boundary_search(fun, x_range, err)

    x1 = x_range[1]
    x2 = x_range[2]

    y1 = fun(x1)
    if y1 == 1
        return x1
    end

    y2 = fun(x2)
    if y2 == 0
        return x2
    end

    while abs(x1 - x2) > 2err

        t0 = (x1 + x2) / 2
        x0 = fun(t0)

        if x0 == y1
            y1 = x0
            x1 = t0
        else
            y2 = x0
            x2 = t0
        end

    end

    X = (x1 + x2) / 2
    return X

end

max_stars_repo_name: vkumpost/stoosc

max_issues_repo_name: vkumpost/stoosc

max_forks_repo_name: vkumpost/stoosc


text:
module SQLite

using Random, Serialization
using WeakRefStrings, DBInterface

export DBInterface, SQLiteException

struct SQLiteException <: Exception
    msg::AbstractString
end

include("consts.jl")
include("api.jl")

const DBHandle = Ptr{Cvoid}   # SQLite3 DB connection handle
const StmtHandle = Ptr{Cvoid} # SQLite3 prepared statement handle

# Normal constructor from filename
sqliteexception(handle::DBHandle) = SQLiteException(unsafe_string(sqlite3_errmsg(handle)))
sqliteerror(handle::DBHandle) = throw(sqliteexception(handle))

"""
Internal wrapper that holds the handle to SQLite3 prepared statement.
It is managed by [`SQLite.DB`](@ref) and referenced by the "public" [`SQLite.Stmt`](@ref) object.

When no `SQLite.Stmt` instances reference the given `SQlite._Stmt` object,
it is closed automatically.

When `SQLite.DB` is closed or [`SQLite.finalize_statements!`](@ref) is called,
all its `SQLite._Stmt` objects are closed.
"""
mutable struct _Stmt
    handle::StmtHandle
    params::Dict{Int, Any}

    function _Stmt(handle::StmtHandle)
        stmt = new(handle, Dict{Int, Any}())
        finalizer(_close!, stmt)
        return stmt
    end
end

# close statement
function _close!(stmt::_Stmt)
    stmt.handle == C_NULL || sqlite3_finalize(stmt.handle)
    stmt.handle = C_NULL
    return
end

# _Stmt unique identifier in DB
const _StmtId = Int

"""
    `SQLite.DB()` => in-memory SQLite database
    `SQLite.DB(file)` => file-based SQLite database

Constructors for a representation of an sqlite database, either backed by an on-disk file or in-memory.

`SQLite.DB` requires the `file` string argument in the 2nd definition
as the name of either a pre-defined SQLite database to be opened,
or if the file doesn't exist, a database will be created.
Note that only sqlite 3.x version files are supported.

The `SQLite.DB` object represents a single connection to an SQLite database.
All other SQLite.jl functions take an `SQLite.DB` as the first argument as context.

To create an in-memory temporary database, call `SQLite.DB()`.

The `SQLite.DB` will be automatically closed/shutdown when it goes out of scope
(i.e. the end of the Julia session, end of a function call wherein it was created, etc.)
"""
mutable struct DB <: DBInterface.Connection
    file::String
    handle::DBHandle
    stmts::Dict{_StmtId, _Stmt} # opened prepared statements

    lastStmtId::_StmtId

    function DB(f::AbstractString)
        handle = Ref{DBHandle}()
        f = isempty(f) ? f : expanduser(f)
        if @OK sqlite3_open(f, handle)
            db = new(f, handle[], Dict{StmtHandle, _Stmt}(), 0)
            finalizer(_close, db)
            return db
        else # error
            sqliteerror(handle[])
        end
    end
end
DB() = DB(":memory:")
DBInterface.connect(::Type{DB}) = DB()
DBInterface.connect(::Type{DB}, f::AbstractString) = DB(f)
DBInterface.close!(db::DB) = _close(db)
Base.close(db::DB) = _close(db)
Base.isopen(db::DB) = db.handle != C_NULL

# close all prepared statements of db connection
function finalize_statements!(db::DB)
    for stmt in values(db.stmts)
        _close!(stmt)
    end
    empty!(db.stmts)
end

function _close(db::DB)
    finalize_statements!(db)
    # disconnect from DB
    db.handle == C_NULL || sqlite3_close_v2(db.handle)
    db.handle = C_NULL
    return
end

sqliteerror(db::DB) = sqliteerror(db.handle)
sqliteexception(db::DB) = sqliteexception(db.handle)

Base.show(io::IO, db::SQLite.DB) = print(io, string("SQLite.DB(", "\"$(db.file)\"", ")"))

# prepare given sql statement
function _Stmt(db::DB, sql::AbstractString)
    handle = Ref{StmtHandle}()
    sqliteprepare(db, sql, handle, Ref{StmtHandle}())
    return _Stmt(handle[])
end

"""
    SQLite.Stmt(db, sql) => SQL.Stmt

Prepares an optimized internal representation of SQL statement in
the context of the provided SQLite3 `db` and constructs the `SQLite.Stmt`
Julia object that holds a reference to the prepared statement.

*Note*: the `sql` statement is not actually executed, but only compiled
(mainly for usage where the same statement is executed multiple times
with different parameters bound as values).

Internally `SQLite.Stmt` constructor creates the [`SQLite._Stmt`](@ref) object that is managed by `db`.
`SQLite.Stmt` references the `SQLite._Stmt` by its unique id.

The `SQLite.Stmt` will be automatically closed/shutdown when it goes out of scope
(i.e. the end of the Julia session, end of a function call wherein it was created, etc.).
One can also call `DBInterface.close!(stmt)` to immediately close it.

All prepared statements of a given DB connection are also automatically closed when the
DB is disconnected or when [`SQLite.finalize_statements!`](@ref) is explicitly called.
"""
mutable struct Stmt <: DBInterface.Statement
    db::DB
    id::_StmtId # id of _Stmt inside db (may refer to already closed connection)

    function Stmt(db::DB, sql::AbstractString)
        _stmt = _Stmt(db, sql)
        id = (db.lastStmtId += 1)
        stmt = new(db, id)
        db.stmts[id] = _stmt # FIXME check for duplicate handle?
        finalizer(_finalize, stmt)
        return stmt
    end
end

# check if the statement is ready (not finalized due to
# _close(_Stmt) called and the statment handle removed from DB)
isready(stmt::Stmt) = haskey(stmt.db.stmts, stmt.id)

# get underlying _Stmt or nothing if not found
_stmt_safe(stmt::Stmt) = get(stmt.db.stmts, stmt.id, nothing)

# get underlying _Stmt or throw if not found
@inline function _stmt(stmt::Stmt)
    _st = _stmt_safe(stmt)
    (_st === nothing) && throw(SQLiteException("Statement $(stmt.id) not found"))
    return _st
end

# automatically finalizes prepared statement (_Stmt)
# when no Stmt objects refer to it and removes
# it from the db.stmts collection
_finalize(stmt::Stmt) = DBInterface.close!(stmt)

# explicitly close prepared statement
function DBInterface.close!(stmt::Stmt)
    _st = _stmt_safe(stmt)
    if _st !== nothing
        _close!(_st)
        delete!(stmt.db.stmts, stmt.id) # remove the _Stmt
    end
    return stmt
end

sqliteprepare(db::DB, sql::AbstractString, stmt::Ref{StmtHandle}, null::Ref{StmtHandle}) =
    @CHECK db sqlite3_prepare_v2(db.handle, sql, stmt, null)

include("UDF.jl")
export @sr_str

"""
    SQLite.clear!(stmt::SQLite.Stmt)

Clears any bound values to a prepared SQL statement
"""
function clear!(stmt::Stmt)
    _st = _stmt(stmt)
    sqlite3_clear_bindings(_st.handle)
    empty!(_st.params)
    return
end

"""
    SQLite.bind!(stmt::SQLite.Stmt, values)

bind `values` to parameters in a prepared [`SQLite.Stmt`](@ref). Values can be:

* `Vector` or `Tuple`: where each element will be bound to an SQL parameter by index order
* `Dict` or `NamedTuple`; where values will be bound to named SQL parameters by the `Dict`/`NamedTuple` key

Additional methods exist for working individual SQL parameters:

* `SQLite.bind!(stmt, name, val)`: bind a single value to a named SQL parameter
* `SQLite.bind!(stmt, index, val)`: bind a single value to a SQL parameter by index number

From the [SQLite documentation](https://www3.sqlite.org/cintro.html):

> Usually, though,
> it is not useful to evaluate exactly the same SQL statement more than once.
> More often, one wants to evaluate similar statements.
> For example, you might want to evaluate an INSERT statement
> multiple times though with different values to insert.
> To accommodate this kind of flexibility,
> SQLite allows SQL statements to contain parameters
> which are "bound" to values prior to being evaluated.
> These values can later be changed and the same prepared statement
> can be evaluated a second time using the new values.
>
> In SQLite,
> wherever it is valid to include a string literal,
> one can use a parameter in one of the following forms:
>
> - `?`
> - `?NNN`
> - `:AAA`
> - `\$AAA`
> - `@AAA`
>
> In the examples above,
> `NNN` is an integer value and `AAA` is an identifier.
> A parameter initially has a value of `NULL`.
> Prior to calling `sqlite3_step()` for the first time
> or immediately after `sqlite3_reset()``,
> the application can invoke one of the `sqlite3_bind()` interfaces
> to attach values to the parameters.
> Each call to `sqlite3_bind()` overrides prior bindings on the same parameter.

"""
function bind! end

function bind!(stmt::_Stmt, params::DBInterface.NamedStatementParams)
    nparams = sqlite3_bind_parameter_count(stmt.handle)
    (nparams == length(params)) || throw(SQLiteException("values should be provided for all query placeholders"))
    for i in 1:nparams
        name = unsafe_string(sqlite3_bind_parameter_name(stmt.handle, i))
        isempty(name) && throw(SQLiteException("nameless parameters should be passed as a Vector"))
        # name is returned with the ':', '@' or '$' at the start
        sym = Symbol(name[2:end])
        haskey(params, sym) || throw(SQLiteException("`$name` not found in values keyword arguments to bind to sql statement"))
        bind!(stmt, i, params[sym])
    end
end

function bind!(stmt::_Stmt, values::DBInterface.PositionalStatementParams)
    nparams = sqlite3_bind_parameter_count(stmt.handle)
    (nparams == length(values)) || throw(SQLiteException("values should be provided for all query placeholders"))
    for i in 1:nparams
        @inbounds bind!(stmt, i, values[i])
    end
end

bind!(stmt::Stmt, values::DBInterface.StatementParams) = bind!(_stmt(stmt), values)

bind!(stmt::Union{_Stmt, Stmt}; kwargs...) = bind!(stmt, kwargs.data)

# Binding parameters to SQL statements
function bind!(stmt::_Stmt, name::AbstractString, val::Any)
    i::Int = sqlite3_bind_parameter_index(stmt.handle, name)
    if i == 0
        throw(SQLiteException("SQL parameter $name not found in $stmt"))
    end
    return bind!(stmt, i, val)
end

# binding method for internal _Stmt class
bind!(stmt::_Stmt, i::Integer, val::AbstractFloat)  = (stmt.params[i] = val; @CHECK stmt.db sqlite3_bind_double(stmt.handle, i, Float64(val)); return nothing)
bind!(stmt::_Stmt, i::Integer, val::Int32)          = (stmt.params[i] = val; @CHECK stmt.db sqlite3_bind_int(stmt.handle, i, val); return nothing)
bind!(stmt::_Stmt, i::Integer, val::Int64)          = (stmt.params[i] = val; @CHECK stmt.db sqlite3_bind_int64(stmt.handle, i, val); return nothing)
bind!(stmt::_Stmt, i::Integer, val::Missing)        = (stmt.params[i] = val; @CHECK stmt.db sqlite3_bind_null(stmt.handle, i); return nothing)
bind!(stmt::_Stmt, i::Integer, val::AbstractString) = (stmt.params[i] = val; @CHECK stmt.db sqlite3_bind_text(stmt.handle, i, val); return nothing)
bind!(stmt::_Stmt, i::Integer, val::WeakRefString{UInt8})   = (stmt.params[i] = val; @CHECK stmt.db sqlite3_bind_text(stmt.handle, i, val.ptr, val.len); return nothing)
bind!(stmt::_Stmt, i::Integer, val::WeakRefString{UInt16})  = (stmt.params[i] = val; @CHECK stmt.db sqlite3_bind_text16(stmt.handle, i, val.ptr, val.len*2); return nothing)
bind!(stmt::_Stmt, i::Integer, val::Bool)           = (stmt.params[i] = val; @CHECK stmt.db sqlite3_bind_int(stmt.handle, i, Int32(val)); return nothing)
bind!(stmt::_Stmt, i::Integer, val::Vector{UInt8})  = (stmt.params[i] = val; @CHECK stmt.db sqlite3_bind_blob(stmt.handle, i, val); return nothing)
# Fallback is BLOB and defaults to serializing the julia value

bind!(stmt::Stmt, param::Union{Integer, AbstractString}, val::Any) = bind!(_stmt(stmt), param, val)

# internal wrapper mutable struct to, in-effect, mark something which has been serialized
struct Serialized
    object
end

const GLOBAL_BUF = IOBuffer()
function sqlserialize(x)
    seekstart(GLOBAL_BUF)
    # deserialize will sometimes return a random object when called on an array
    # which has not been previously serialized, we can use this mutable struct to check
    # that the array has been serialized
    s = Serialized(x)
    Serialization.serialize(GLOBAL_BUF, s)
    return take!(GLOBAL_BUF)
end
# fallback method to bind arbitrary julia `val` to the parameter at index `i` (object is serialized)
bind!(stmt::_Stmt, i::Integer, val::Any) = bind!(stmt, i, sqlserialize(val))

struct SerializeError <: Exception
    msg::String
end

# magic bytes that indicate that a value is in fact a serialized julia value, instead of just a byte vector
# these bytes depend on the julia version and other things, so they are determined using an actual serialization
const SERIALIZATION = sqlserialize(0)[1:18]

function sqldeserialize(r)
    ret = ccall(:memcmp, Int32, (Ptr{UInt8}, Ptr{UInt8}, UInt),
            SERIALIZATION, r, min(sizeof(SERIALIZATION), sizeof(r)))
    if ret == 0
        try
            v = Serialization.deserialize(IOBuffer(r))
            return v.object
        catch e
            throw(SerializeError("Error deserializing non-primitive value out of database; this is probably due to using SQLite.jl with a different Julia version than was used to originally serialize the database values. The same Julia version that was used to serialize should be used to extract the database values into a different format (csv file, feather file, etc.) and then loaded back into the sqlite database with the current Julia version."))
        end
    else
        return r
    end
end
#TODO:
 #int sqlite3_bind_zeroblob(sqlite3_stmt*, int, int n);
 #int sqlite3_bind_value(sqlite3_stmt*, int, const sqlite3_value*);

function juliatype(handle, col)
    t = SQLite.sqlite3_column_decltype(handle, col)
    if t != C_NULL
        T = juliatype(unsafe_string(t))
        T !== Any && return T
    end
    x = SQLite.sqlite3_column_type(handle, col)
    if x == SQLite.SQLITE_BLOB
        val = SQLite.sqlitevalue(Any, handle, col)
        return typeof(val)
    else
        return juliatype(x)
    end
end

juliatype(x::Integer) = x == SQLITE_INTEGER ? Int : x == SQLITE_FLOAT ? Float64 : x == SQLITE_TEXT ? String : Any
juliatype(x::String) = x == "INTEGER" ? Int : x in ("NUMERIC","REAL") ? Float64 : x == "TEXT" ? String : Any

sqlitevalue(::Type{T}, handle, col) where {T <: Union{Base.BitSigned, Base.BitUnsigned}} = convert(T, sqlite3_column_int64(handle, col))
const FLOAT_TYPES = Union{Float16, Float32, Float64} # exclude BigFloat
sqlitevalue(::Type{T}, handle, col) where {T <: FLOAT_TYPES} = convert(T, sqlite3_column_double(handle, col))
#TODO: test returning a WeakRefString instead of calling `unsafe_string`
sqlitevalue(::Type{T}, handle, col) where {T <: AbstractString} = convert(T, unsafe_string(sqlite3_column_text(handle, col)))
function sqlitevalue(::Type{T}, handle, col) where {T}
    blob = convert(Ptr{UInt8}, sqlite3_column_blob(handle, col))
    b = sqlite3_column_bytes(handle, col)
    buf = zeros(UInt8, b) # global const?
    unsafe_copyto!(pointer(buf), blob, b)
    r = sqldeserialize(buf)::T
    return r
end

# conversion from Julia to SQLite3 types
sqlitetype_(::Type{<:Integer}) = "INT"
sqlitetype_(::Type{<:AbstractFloat}) = "REAL"
sqlitetype_(::Type{<:AbstractString}) = "TEXT"
sqlitetype_(::Type{Bool}) = "INT"
sqlitetype_(::Type) = "BLOB" # fallback

sqlitetype(::Type{Missing}) = "NULL"
sqlitetype(::Type{Union{T, Missing}}) where T = sqlitetype_(T)
sqlitetype(::Type{T}) where T = string(sqlitetype_(T), " NOT NULL")

"""
    SQLite.execute(db::SQLite.DB, sql::AbstractString, [params]) -> Int
    SQLite.execute(stmt::SQLite.Stmt, [params]) -> Int

An internal method that executes the SQL statement (provided either as a `db` connection and `sql` command,
or as an already prepared `stmt` (see [`SQLite.Stmt`](@ref))) with given `params` parameters
(either positional (`Vector` or `Tuple`), named (`Dict` or `NamedTuple`), or specified as keyword arguments).

Returns the SQLite status code of operation.

*Note*: this is a low-level method that just executes the SQL statement,
but does not retrieve any data from `db`.
To get the results of a SQL query, it is recommended to use [`DBInterface.execute`](@ref).
"""
function execute end

function execute(db::DB, stmt::_Stmt, params::DBInterface.StatementParams=())
    sqlite3_reset(stmt.handle)
    bind!(stmt, params)
    r = sqlite3_step(stmt.handle)
    if r == SQLITE_DONE
        sqlite3_reset(stmt.handle)
    elseif r != SQLITE_ROW
        e = sqliteexception(db)
        sqlite3_reset(stmt.handle)
        throw(e)
    end
    return r
end

execute(stmt::Stmt, params::DBInterface.StatementParams) =
    execute(stmt.db, _stmt(stmt), params)

execute(stmt::Stmt; kwargs...) = execute(stmt, kwargs.data)

function execute(db::DB, sql::AbstractString, params::DBInterface.StatementParams)
     # prepare without registering _Stmt in DB
    _stmt = _Stmt(db, sql)
    try
        return execute(db, _stmt, params)
    finally
        _close!(_stmt) # immediately close, don't wait for GC
    end
end

execute(db::DB, sql::AbstractString; kwargs...) = execute(db, sql, kwargs.data)

"""
    SQLite.esc_id(x::Union{AbstractString,Vector{AbstractString}})

Escape SQLite identifiers
(e.g. column, table or index names).
Can be either a string or a vector of strings
(note does not check for null characters).
A vector of identifiers will be separated by commas.

Example:

```julia
julia> using SQLite, DataFrames

julia> df = DataFrame(label=string.(rand("abcdefg", 10)), value=rand(10));

julia> db = SQLite.DB(mktemp()[1]);

julia> tbl |> SQLite.load!(db, "temp");

julia> DBInterface.execute(db,"SELECT * FROM temp WHERE label IN ('a','b','c')") |> DataFrame
4×2 DataFrame
│ Row │ label   │ value    │
│     │ String⍰ │ Float64⍰ │
├─────┼─────────┼──────────┤
│ 1   │ c       │ 0.603739 │
│ 2   │ c       │ 0.429831 │
│ 3   │ b       │ 0.799696 │
│ 4   │ a       │ 0.603586 │

julia> q = ['a','b','c'];

julia> DBInterface.execute(db,"SELECT * FROM temp WHERE label IN (\$(SQLite.esc_id(q)))") |> DataFrame
4×2 DataFrame
│ Row │ label   │ value    │
│     │ String⍰ │ Float64⍰ │
├─────┼─────────┼──────────┤
│ 1   │ c       │ 0.603739 │
│ 2   │ c       │ 0.429831 │
│ 3   │ b       │ 0.799696 │
│ 4   │ a       │ 0.603586 │
```
"""
function esc_id end

esc_id(x::AbstractString) = "\"" * replace(x, "\""=>"\"\"") * "\""
esc_id(X::AbstractVector{S}) where {S <: AbstractString} = join(map(esc_id, X), ',')

# Transaction-based commands
"""
    SQLite.transaction(db, mode="DEFERRED")
    SQLite.transaction(func, db)

Begin a transaction in the specified `mode`, default = "DEFERRED".

If `mode` is one of "", "DEFERRED", "IMMEDIATE" or "EXCLUSIVE" then a
transaction of that (or the default) mutable struct is started. Otherwise a savepoint
is created whose name is `mode` converted to AbstractString.

In the second method, `func` is executed within a transaction (the transaction being committed upon successful execution)
"""
function transaction end

function transaction(db::DB, mode="DEFERRED")
    execute(db, "PRAGMA temp_store=MEMORY;")
    if uppercase(mode) in ["", "DEFERRED", "IMMEDIATE", "EXCLUSIVE"]
        execute(db, "BEGIN $(mode) TRANSACTION;")
    else
        execute(db, "SAVEPOINT $(mode);")
    end
end

@inline function transaction(f::Function, db::DB)
    # generate a random name for the savepoint
    name = string("SQLITE", Random.randstring(10))
    execute(db, "PRAGMA synchronous = OFF;")
    transaction(db, name)
    try
        f()
    catch
        rollback(db, name)
        rethrow()
    finally
        # savepoints are not released on rollback
        commit(db, name)
        execute(db, "PRAGMA synchronous = ON;")
    end
end

"""
    SQLite.commit(db)
    SQLite.commit(db, name)

commit a transaction or named savepoint
"""
function commit end

commit(db::DB) = execute(db, "COMMIT TRANSACTION;")
commit(db::DB, name::AbstractString) = execute(db, "RELEASE SAVEPOINT $(name);")

"""
    SQLite.rollback(db)
    SQLite.rollback(db, name)

rollback transaction or named savepoint
"""
function rollback end

rollback(db::DB) = execute(db, "ROLLBACK TRANSACTION;")
rollback(db::DB, name::AbstractString) = execute(db, "ROLLBACK TRANSACTION TO SAVEPOINT $(name);")

"""
    SQLite.drop!(db, table; ifexists::Bool=true)

drop the SQLite table `table` from the database `db`; `ifexists=true` will prevent an error being thrown if `table` doesn't exist
"""
function drop!(db::DB, table::AbstractString; ifexists::Bool=false)
    exists = ifexists ? "IF EXISTS" : ""
    transaction(db) do
        execute(db, "DROP TABLE $exists $(esc_id(table))")
    end
    execute(db, "VACUUM")
    return
end

"""
    SQLite.dropindex!(db, index; ifexists::Bool=true)

drop the SQLite index `index` from the database `db`; `ifexists=true` will not return an error if `index` doesn't exist
"""
function dropindex!(db::DB, index::AbstractString; ifexists::Bool=false)
    exists = ifexists ? "IF EXISTS" : ""
    transaction(db) do
        execute(db, "DROP INDEX $exists $(esc_id(index))")
    end
    return
end

"""
    SQLite.createindex!(db, table, index, cols; unique=true, ifnotexists=false)

create the SQLite index `index` on the table `table` using `cols`,
which may be a single column or vector of columns.
`unique` specifies whether the index will be unique or not.
`ifnotexists=true` will not throw an error if the index already exists
"""
function createindex!(db::DB, table::AbstractString, index::AbstractString, cols::Union{S, AbstractVector{S}};
                      unique::Bool=true, ifnotexists::Bool=false) where {S <: AbstractString}
    u = unique ? "UNIQUE" : ""
    exists = ifnotexists ? "IF NOT EXISTS" : ""
    transaction(db) do
        execute(db, "CREATE $u INDEX $exists $(esc_id(index)) ON $(esc_id(table)) ($(esc_id(cols)))")
    end
    execute(db, "ANALYZE $index")
    return
end

"""
    SQLite.removeduplicates!(db, table, cols)

Removes duplicate rows from `table` based on the values in `cols`, which is an array of column names.

A convenience method for the common task of removing duplicate
rows in a dataset according to some subset of columns that make up a "primary key".
"""
function removeduplicates!(db::DB, table::AbstractString, cols::AbstractArray{T}) where {T <: AbstractString}
    colsstr = ""
    for c in cols
       colsstr = colsstr * esc_id(c) * ","
    end
    colsstr = chop(colsstr)
    transaction(db) do
        execute(db, "DELETE FROM $(esc_id(table)) WHERE _ROWID_ NOT IN (SELECT max(_ROWID_) from $(esc_id(table)) GROUP BY $(colsstr));")
    end
    execute(db, "ANALYZE $table")
    return
 end

include("tables.jl")

"""
    SQLite.tables(db, sink=columntable)

returns a list of tables in `db`
"""
tables(db::DB, sink=columntable) = DBInterface.execute(sink, db, "SELECT name FROM sqlite_master WHERE type='table';")

"""
    SQLite.indices(db, sink=columntable)

returns a list of indices in `db`
"""
indices(db::DB, sink=columntable) = DBInterface.execute(sink, db, "SELECT name FROM sqlite_master WHERE type='index';")

"""
    SQLite.columns(db, table, sink=columntable)

returns a list of columns in `table`
"""
columns(db::DB, table::AbstractString, sink=columntable) = DBInterface.execute(sink, db, "PRAGMA table_info($(esc_id(table)))")

"""
    SQLite.last_insert_rowid(db)

returns the auto increment id of the last row
"""
last_insert_rowid(db::DB) = sqlite3_last_insert_rowid(db.handle)

"""
    SQLite.enable_load_extension(db, enable::Bool=true)

Enables extension loading (off by default) on the sqlite database `db`. Pass `false` as the second argument to disable.
"""
function enable_load_extension(db::DB, enable::Bool=true)
   ccall((:sqlite3_enable_load_extension, SQLite.libsqlite), Cint, (Ptr{Cvoid}, Cint), db.handle, enable)
end

"""
    SQLite.busy_timeout(db, ms::Integer=0)

Set a busy handler that sleeps for a specified amount of milliseconds  when a table is locked. After at least ms milliseconds of sleeping, the handler will return 0, causing sqlite to return SQLITE_BUSY.
"""
function busy_timeout(db::DB, ms::Integer=0)
    sqlite3_busy_timeout(db.handle, ms)
end

end # module


text:
circle{T}(uv::Vec{2, T}) = T(0.5) - norm(uv)

type Uniforms{F}
    projection::Mat4f0
    strokecolor::Vec4f0
    glowcolor::Vec4f0
    distance_func::F
end

type TextUniforms
    projection::Mat4f0
    strokecolor::Vec4f0
    glowcolor::Vec4f0
end

vert_particles(vertex, scene, uniforms, image) = vert_particles(vertex, scene, uniforms)
function vert_particles(vertex, scene, uniforms)
    p = getposition(vertex)
    scale = getscale(vertex)
    return Vertex2Geom(
        getuv(vertex),
        getcolor(vertex),
        Vec4f0(p[1], p[2], scale[1], scale[2])
    )
end

"""
Emits a vertex with
"""
function emit_vertex(emit!, vertex, uv, arg, pos, uniforms)
    datapoint = uniforms.projection * Vec4f0(pos[1], pos[2], 0, 1)
    final_position = uniforms.projection * Vec4f0(vertex[1], vertex[2], 0, 0)
    emit!(datapoint .+ final_position, (uv, arg.color))
    return
end

function geom_particles(emit!, vertex_out, uniforms, image)
    geom_particles(emit!, vertex_out, uniforms)
    return
end
function geom_particles(emit!, vertex_out, uniforms)
    # get arguments from first face
    # (there is only one in there anywas, since primitive type is point)
    # (position, vertex_out)
    arg = vertex_out[1]
    # emit quad as triangle strip
    # v3. ____ . v4
    #    |\   |
    #    | \  |
    #    |  \ |
    #    |___\|
    # v1*      * v2
    pos_scale = arg.rect
    pos = pos_scale[Vec(1, 2)]
    scale = pos_scale[Vec(3, 4)]
    quad = Vec4f0(0f0, 0f0, scale[1], scale[2])
    uv = arg.uvrect
    emit_vertex(emit!, quad[Vec(1, 2)], uv[Vec(1, 4)], arg, pos, uniforms)
    emit_vertex(emit!, quad[Vec(1, 4)], uv[Vec(1, 2)], arg, pos, uniforms)
    emit_vertex(emit!, quad[Vec(3, 2)], uv[Vec(3, 4)], arg, pos, uniforms)
    emit_vertex(emit!, quad[Vec(3, 4)], uv[Vec(3, 2)], arg, pos, uniforms)
    return
end

function sdf2color(dist, bg_color, color)
    inside = aastep(0f0, dist)
    mix(bg_color, color, inside)
end
function frag_particles(geom_out, uniforms, image)
    uv = geom_out[1]; color = geom_out[2]
    dist = -image[uv][1]
    bg_color = Vec4f0(0f0, 0f0, 0f0, 0f0)
    (sdf2color(dist, bg_color, color), )
end
function frag_particles(geom_out, uniforms)
    uv = geom_out[1]; color = geom_out[2]
    dist = uniforms.distance_func(uv)
    bg_color = Vec4f0(0f0, 0f0, 0f0, 0f0)
    (sdf2color(dist, bg_color, color), )
end


text:
module ZeroMeanFESpacesTests

using Test
using Gridap.Arrays
using Gridap.Geometry
using Gridap.Fields
using Gridap.FESpaces
using Gridap.CellData
using Gridap.ReferenceFEs

domain = (0,1,0,1)
partition = (4,4)
model = CartesianDiscreteModel(domain,partition)

order = 2

trian = get_triangulation(model)
degree = order
dΩ = Measure(trian,degree)

_V = FESpace(model,ReferenceFE(lagrangian,Float64,order);conformity=:L2)

V = ZeroMeanFESpace(_V,dΩ)

cellmat = [rand(4,4) for cell in 1:num_cells(model)]
cellvec = [rand(4) for cell in 1:num_cells(model)]
cellmatvec = pair_arrays(cellmat,cellvec)
test_single_field_fe_space(V,cellmatvec,cellmat,cellvec,trian)

U = TrialFESpace(V)
test_single_field_fe_space(U,cellmatvec,cellmat,cellvec,trian)
@test isa(U,ZeroMeanFESpace)

fun(x) = sin(4*pi*(x[1]+x[2]^2)) + 3
uh = interpolate(fun, U)

mean1 = sum(∫(uh)*dΩ)

tol = 1.0e-10
@test abs(mean1) < tol

V = FESpace(model,ReferenceFE(lagrangian,Float64,order);conformity=:L2,constraint=:zeromean)
@test isa(V,ZeroMeanFESpace)

end # module


text:
# even_parity() computes the goal for evolution of a circuit to solve the even parity problem
# odd_parity() computes the goal for evolution of a circuit to solve the odd parity problem
function odd_parity(numinputs::Int64)
  gb = reverse(get_bits( construct_context(numinputs),numinputs))
  result = MyInt(0)
  for gbi in gb
    result <<= 1
    #result |= MyInt(1 - (count_ones(gbi) % 2))
    result |= MyInt(1 - (count_ones(gbi) & 1))
  end
  result
end

function even_parity(numinputs::Int64)
  gb = reverse(get_bits( construct_context(numinputs),numinputs))
  result = MyInt(0)
  for gbi in gb
    result <<= 1
    result |= MyInt(count_ones(gbi) & 1)
  end
  result
end

even(x::Integer) = Int64(1-(x & MyInt(1)))
odd(x::Integer) = Int64(x & MyInt(1))

function test_odd_parity( numinputs )
  blist = [  MyInt_to_binary_list(MyInt(x),numinputs) for x =  MyInt(2^numinputs-1):-1:0x0 ]
  clist = map(x->even(sum(x)),blist)
  to_bin = to_binary(even_parity(numinputs),2^numinputs)
  @assert to_bin == clist
end

function test_even_parity( numinputs )
  blist = [  MyInt_to_binary_list(MyInt(x),numinputs) for x =  MyInt(2^numinputs-1):-1:0x0 ]
  clist = map(x->odd(sum(x)),blist)
  to_bin = to_binary(odd_parity(numinputs),2^numinputs)
  @assert to_bin == clist
end

function MyInt_to_binary_list( x::MyInt, size::Int64 )
  #size = 64  # If the type of x is UInt64, must be 64.  If the type of x is UInt32, must be 32.
  if MyInt==UInt16
    @assert size <= 4
  end
  if MyInt==UInt32
    @assert size <= 5
  end
  if MyInt==UInt64
    @assert size <= 6
  end
  if MyInt==UInt128
    @assert size <= 7
  end
  #println("size: ",size)
  shift = size-1
  result = fill(MyInt(0),size)
  for i = 1:size
    result[i] = (x >> shift) & convert(MyInt,1)
    shift -= 1
  end
  result
end

function shift_extend( p::Parameters, maxsteps::Int64; parity::Bool=:true, extend::Int64=3 )
  c = random_chromosome(p)
  if parity
    goal = [even_parity(p.numinputs)]
  else
    goal= [rand_shift_extend(extend,p.numinputs)]
  end
  #println("goal: ",goal)
  (nc,steps)=neutral_evolution(c,goal,maxsteps)
  ssteps =  steps < maxsteps ? steps : 0
  failure = ssteps == 0 ? 1 : 0
  (parity,p.numinputs,p.numinteriors,p.numlevelsback,ssteps,failure)
end

function run_shift_extend( numinputs::IntRange, reps::Int64, maxsteps::Int64;
    csvfile::String="", parity::Bool=:true, extend::Int64=3 )
  number_gates_multiplier = 5
  #number_gates_multiplier = 7  # Changed from 5 to 7 on 4/29 to do experiment with 15 inputs
  df = DataFrame()
  df.parity = Bool[]
  df.numinputs = Int64[]
  df.numgates = Int64[]
  df.levelsback = Int64[]
  df.repetitions = Int64[]
  if !parity
    df.extend = Int64[]
  end
  df.average_steps = Float64[]
  df.median_steps = Float64[]
  for ni in numinputs
    if extend >= ni
      error("The keywork parameter extend=",extend,"  must be less than numinputsi=",ni)
    end
    ng = number_gates_multiplier*ni
    levelsback = Int(floor(ng/2))
    p = Parameters( ni, 1, ng, levelsback )
    funcs = default_funcs(p.numinputs)
    rows = pmap(x->shift_extend(p,maxsteps,parity=parity,extend=extend),collect(1:reps))
    sum_steps = 0
    sum_failures = 0
    steps_list = Int64[]
    for r in rows
      if r[5] > 0
        push!(steps_list,r[5])
      else
        sum_failures += r[6]
      end
    end
    r = rows[1]
    if parity
      row = (r[1],r[2],r[3],r[4],reps-sum_failures,mean(steps_list),median(steps_list))
    else
      row = (r[1],r[2],r[3],r[4],reps-sum_failures,extend,mean(steps_list),median(steps_list))
    end
    push!(df,row)
  end
  if length(csvfile) > 0
    hostname = chomp(open("/etc/hostname") do f read(f,String) end)
    open( csvfile, "w" ) do f
      println(f,"# date and time: ",Dates.now())
      println(f,"# host: ",hostname," with ",nprocs()-1,"  processes: " )
      println(f,"# funcs: ", Main.CGP.default_funcs(numinputs[end]))
      println(f,"# MyInt: ", MyInt )
      println(f,"# reps: ",reps)
      println(f,"# extend: ",extend)
      println(f,"# maxsteps: ", maxsteps )
      println(f,"# number_gates_multiplier: ",number_gates_multiplier)
      CSV.write( f, df, append=true, writeheader=true )
    end
  end
  df
end

function reverse_bits( x::MyInt, size::Int64 )
  result = MyInt(0)
  mask = MyInt(1)
  final = MyInt(1) << (size-1)
  result |= mask & x != MyInt(0) ? MyInt(1) : MyInt(0)
  while mask != final
    result <<= 1
    mask <<= 1
    result |= mask & x != MyInt(0) ? MyInt(1) : MyInt(0)
    #@printf("mask: 0x%x  result: 0x%x\n",mask,result)
  end
  result
end

function rand_shift_extend( orig_nbits::Int64, result_nbits::Int64; parity::Bool=:false )
  @assert result_nbits > orig_nbits
  if parity
    sgl = MyInt(0x96)
  else
    sgl = randgoal(orig_nbits,1)[1]
  end
  #@printf("sgl: 0x%x\n",sgl)
  rgl = sgl |= sgl | (reverse_bits(sgl,2^orig_nbits) << 2^orig_nbits)
  #@printf("rgl: 0x%x\n",rgl)
  for i = (orig_nbits+2):result_nbits
    rgl |= (rgl << 2^(i-1))
    #@printf("i: %d  rgl: 0x%x\n",i,rgl)
  end
  rgl
end


text:
"""
    assemble_matrices(model)

Assemble finite element matrices.
"""
function assemble_matrices(model)
    @unpack mesh, D, T₂, ρ = model

    # Deduce sizes
    ncompartment = length(ρ)

    # Assemble finite element matrices compartment-wise
    M_cmpts = []
    S_cmpts = []
    Mx_cmpts = [[] for _ = 1:3]
    for icmpt = 1:ncompartment
        # Finite elements
        points = mesh.points[icmpt]
        elements = mesh.elements[icmpt]
        volumes, = get_mesh_volumes(points, elements)

        # Assemble mass, stiffness and flux matrices
        push!(M_cmpts, assemble_mass_matrix(elements', volumes))
        push!(S_cmpts, assemble_stiffness_matrix(elements', points', D[icmpt]))

        # Assemble first order product moment matrices
        for dim = 1:3
            push!(Mx_cmpts[dim], assemble_mass_matrix(elements', volumes, points[dim, :]))
        end
    end

    # Assemble global finite element matrices
    M = blockdiag(M_cmpts...)
    S = blockdiag(S_cmpts...)
    R = blockdiag((M_cmpts ./ T₂)...)
    Mx = [blockdiag(Mx_cmpts[dim]...) for dim = 1:3]
    Q_blocks = assemble_flux_matrices(mesh.points, mesh.facets)
    Q = couple_flux_matrix(model, Q_blocks, false)

    (; M, S, R, Mx, Q, M_cmpts, S_cmpts, Mx_cmpts)
end


text: struct S a b c end

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
immutable Weibull <: ContinuousUnivariateDistribution
    shape::Float64
    scale::Float64
    function Weibull(sh::Real, sc::Real)
    	zero(sh) < sh && zero(sc) < sc || error("Both shape and scale must be positive")
    	new(float64(sh), float64(sc))
    end
end

Weibull(sh::Real) = Weibull(sh, 1.0)

## Support
@continuous_distr_support Weibull 0.0 Inf

## Properties
mean(d::Weibull) = d.scale * gamma(1.0 + 1.0 / d.shape)
median(d::Weibull) = d.scale * log(2.0)^(1.0 / d.shape)

mode(d::Weibull) = d.shape > 1.0 ? (ik = 1.0/d.shape; d.scale * (1.0-ik)^ik) : 0.0

var(d::Weibull) = d.scale^2 * gamma(1.0 + 2.0 / d.shape) - mean(d)^2

function skewness(d::Weibull)
    tmp = gamma(1.0 + 3.0 / d.shape) * d.scale^3
    tmp -= 3.0 * mean(d) * var(d)
    tmp -= mean(d)^3
    return tmp / std(d)^3
end

function kurtosis(d::Weibull)
    λ, k = d.scale, d.shape
    μ = mean(d)
    σ = std(d)
    γ = skewness(d)
    den = λ^4 * gamma(1.0 + 4.0 / k) -
          4.0 * γ * σ^3 * μ -
          6.0 * μ^2 * σ^2 - μ^4
    num = σ^4
    return den / num - 3.0
end

function entropy(d::Weibull)
    λ, k = d.scale, d.shape
    return ((k - 1.0) / k) * -digamma(1.0) + log(λ / k) + 1.0
end

## Functions
function pdf(d::Weibull, x::Real)
    x < zero(x) && return 0.0
    a = x/d.scale
    d.shape/d.scale * a^(d.shape-1.0) * exp(-a^d.shape)
end
function logpdf(d::Weibull, x::Real)
    x < zero(x) && return -Inf
    a = x/d.scale
    log(d.shape/d.scale) + (d.shape-1.0)*log(a) - a^d.shape
end

cdf(d::Weibull, x::Real) = x <= zero(x) ? 0.0 : -expm1(-((x / d.scale)^d.shape))
ccdf(d::Weibull, x::Real) = x <= zero(x) ? 1.0 : exp(-((x / d.scale)^d.shape))
logcdf(d::Weibull, x::Real) = x <= zero(x) ? -Inf : log1mexp(-((x / d.scale)^d.shape))
logccdf(d::Weibull, x::Real) = x <= zero(x) ? 0.0 : -(x / d.scale)^d.shape

quantile(d::Weibull, p::Real) = @checkquantile p d.scale*(-log1p(-p))^(1/d.shape)
cquantile(d::Weibull, p::Real) = @checkquantile p d.scale*(-log(p))^(1/d.shape)
invlogcdf(d::Weibull, lp::Real) = lp > zero(lp) ? NaN : d.scale*(-log1mexp(lp))^(1/d.shape)
invlogccdf(d::Weibull, lp::Real) = lp > zero(lp) ? NaN : d.scale*(-lp)^(1/d.shape)

function gradlogpdf(d::Weibull, x::Real)
  insupport(Weibull, x) ? (d.shape - 1.0) / x - d.shape * x^(d.shape - 1.0) / (d.scale^d.shape) : 0.0
end

## Sampling
rand(d::Weibull) = d.scale*Base.Random.randmtzig_exprnd()^(1/d.shape)


text:
using Agents

mutable struct SchellingAgent <: AbstractAgent
    id::Int             # The identifier number of the agent
    pos::NTuple{3, Int} # The x, y, z location of the agent on a 3D grid
    mood::Bool          # whether the agent is happy in its position. (true = happy)
    group::Int          # The group of the agent, determines mood as it interacts with neighbors
end

using Random # for reproducibility
function initialize(; numagents = 320, griddims = (20, 20, 20), min_to_be_happy = 3, seed = 125)
    space = GridSpace(griddims, periodic = false)
    properties = Dict(:min_to_be_happy => min_to_be_happy)
    rng = Random.MersenneTwister(seed)
    model = ABM(
        SchellingAgent, space;
        properties, rng, scheduler = Schedulers.randomly
    )
    for n in 1:numagents
        agent = SchellingAgent(n, (1, 1, 1), false, n < numagents / 2 ? 1 : 2)
        add_agent_single!(agent, model)
    end
    return model
end

function agent_step!(agent, model)
    minhappy = model.min_to_be_happy
    count_neighbors_same_group = 0
    for neighbor in nearby_agents(agent, model)
        if agent.group == neighbor.group
            count_neighbors_same_group += 1
        end
    end
    if count_neighbors_same_group ≥ minhappy
        agent.mood = true
    else
        move_agent_single!(agent, model)
    end
    return
end

using InteractiveDynamics, GLMakie

model = initialize()
ac(agent) = (:red, :blue)[agent.group]
fig, step = abm_play(model, agent_step!, dummystep; ac, as = 0.5)


text:
import ModelSanitizer
import Test

a = Any[1, 2, 3, 4]
Test.@test eltype(a) == Any
Test.@test eltype(a) != Int
a_fixed = ModelSanitizer._fix_vector_type(a)
Test.@test eltype(a_fixed) != Any
Test.@test eltype(a_fixed) == Int

Test.@test !ModelSanitizer._compare(1.0, missing)
Test.@test !ModelSanitizer._compare(missing, 2.0)
Test.@test !ModelSanitizer._compare(missing, missing)


text:
using CatPlanarGraphs
using Test

include("planar_graph.jl")


text:
using Soss, SossGen
using Gen

const logpdf = Soss.logpdf

m = @model σ begin
    μ ~ Normal()
    y ~ Normal(μ, σ) |> iid(1000)
end;

truth = rand(m(σ=1))

g1 = SossGenerativeFunction(m)
g2 = SossGenerativeFunction(m,codegen)

args = (σ=1,)

# tr = simulate(g, (args,)) # generate a full trace
tr1, w1 = generate(g1, (args,), choicemap(:y => truth.y)) # generate a trace with y = 3

tr2, w2 = generate(g2, (args,), choicemap(:y => truth.y)) # generate a trace with y = 3

@gen function q(current, args)
    @trace(normal(current[:μ], 0.2), :μ)
end

using BenchmarkTools

q(tr1, (args,))

@btime mh(tr1, q, (args,));
@btime mh(tr2, q, (args,));

function run_mh(tr)
  for i=1:1000
    tr, = mh(tr, q, ())
  end
  tr
end

@btime get_choices(run_mh(tr1))

# for i=1:1000
#     mh(tr, Gen.select(:x))
#     mh(tr, Gen.select(:y))
# end

using Revise, SossGen
import Gen

m = @model σ begin
    μ ~ Normal()
    y ~ Normal(μ, σ) |> iid(100000)
end;

args = (σ=1.1,);

truth = rand(m(args));
pairs(truth)

# pairs(::NamedTuple) with 3 entries:
#   :σ => 1.1
#   :μ => 0.398088
#   :y => [-0.849513, 0.777053, -0.800264, 0.105648, -1.50177, 1.26158, -0.43
# 2465…

g1 = SossGenerativeFunction(m);
g2 = SossGenerativeFunction(m,codegen);

# tr = simulate(g, (args,)) # generate a full trace

tr1, w1 = generate(g1, (args,), choicemap(:y => truth.y)); # generate a trace with y = 3
tr2, w2 = generate(g2, (args,), choicemap(:y => truth.y)); # generate a trace with y = 3

@gen (static) function q(current, args)
    @trace(normal(current[:μ], 0.2), :μ)
end

Gen.load_generated_functions();

using BenchmarkTools
@btime mh($tr1, $q, $(args,));

# 891.681 μs (147 allocations: 8.64 KiB)

@btime mh($tr2, $q, $(args,));

# 24.929 μs (148 allocations: 8.67 KiB)

max_stars_repo_name: cscherrer/SossGen

max_issues_repo_name: cscherrer/SossGen

max_forks_repo_name: cscherrer/SossGen


text:
module PSF

import ForwardDiff
import Optim

"""
All the information form a psField file needed to compute a raw PSF for a point.
Attributes:
 - rrows: A matrix of flattened eigenimages.
 - rnrow: The number of rows in an eigenimage.
 - rncol: The number of columns in an eigenimage.
 - cmat: The coefficients of the weight polynomial (see get_psf_at_point()).
"""
immutable RawPSFComponents
    rrows::Array{Float64,2}
    rnrow::Int32
    rncol::Int32
    cmat::Array{Float64,3}
end

# Constants to keep the optimization stable.
const sigma_min = diagm([0.25, 0.25])
const weight_min = 0.05

"""
Convert the parameters of a mixture of 2d Gaussians
sum_k weight_vec[k] * N(mu_vec[k], sigma_vec[k])
to an unconstrained vector that can be passed to an optimizer.
"""
function wrap_parameters(
    mu_vec::Vector{Vector{Float64}},
    sigma_vec::Vector{Matrix{Float64}},
    weight_vec::Vector{Float64})

  @assert length(mu_vec) == length(sigma_vec) == length(weight_vec)
  local K = length(mu_vec)

  # Two mean parameters, three covariance parameters,
  # and one weight per component.
  local par = zeros(K * 6)
  for k = 1:K
    offset = (k - 1) * 6
    par[offset + (1:2)] = mu_vec[k]
    sigma_chol = chol(sigma_vec[k] - sigma_min)
    par[offset + 3] = log(sigma_chol[1, 1])
    par[offset + 4] = sigma_chol[1, 2]
    par[offset + 5] = log(sigma_chol[2, 2])
    par[offset + 6] = log(weight_vec[k] - weight_min)
  end

  par
end

"""
Reverse wrap_parameters().
"""
function unwrap_parameters{T <: Number}(par::Vector{T})
  local K = round(Int, length(par) / 6)
  @assert K == length(par) / 6

  local mu_vec = Array(Vector{T}, K)
  local sigma_vec = Array(Matrix{T}, K)
  local weight_vec = zeros(T, K)

  for k = 1:K
    offset = (k - 1) * 6
    mu_vec[k] = par[offset + (1:2)]
    sigma_chol_vec = par[offset + (3:5)]
    sigma_chol = T[exp(sigma_chol_vec[1]) sigma_chol_vec[2];
                   0.0                    exp(sigma_chol_vec[3])]
    sigma_vec[k] = sigma_chol' * sigma_chol + sigma_min
    weight_vec[k] = exp(par[offset + 6]) + weight_min
  end

  mu_vec, sigma_vec, weight_vec
end

"""
Using a vector of mu, sigma, and weights, get the value of the GMM at x.

Args:
  - x: A length 2 vector to evaluate the GMM at.
  - mu_vec: A vector of location vectors.
  - sigma_vec: A vector of covariance matrices.
  - weight_vec: A vector of weights.

Returns:
  - The weighted sum of the component densities at x.
"""
function evaluate_psf_at_point{T <: Number}(
    x::Vector{Float64},
    mu_vec::Vector{Vector{T}},
    sigma_vec::Vector{Matrix{T}},
    weight_vec::Vector{T})

  @assert length(mu_vec) == length(sigma_vec) == length(weight_vec)
  local K = length(mu_vec)

  @assert length(x) == 2
  local pdf = 0.0
  for k = 1:K
    z = x - mu_vec[k]
    log_pdf = -0.5 * dot(z, sigma_vec[k] \ z) -
               0.5 * logdet(sigma_vec[k]) - log(2 * pi)
    pdf += weight_vec[k] * exp(log_pdf)
  end

  pdf
end

"""
Given a psf matrix, return a matrix of the same size, where each element of
the matrix is a 2-length vector of the [x1, x2] location of that matrix cell.
The locations are chosen so that the scale is pixels, but the center of the
psf is [0, 0].
"""
function get_x_matrix_from_psf(psf::Matrix{Float64})
  psf_center = Float64[ (size(psf, i) - 1) / 2 + 1 for i=1:2 ]
  Vector{Float64}[ Float64[i, j] - psf_center for i=1:size(psf, 1), j=1:size(psf, 2) ]
end

"""
Evaluate a GMM expressed in its unconstrained form.
"""
function render_psf{T <: Number}(par::Vector{T}, x_mat::Matrix{Vector{Float64}})
  local mu_vec, sigma_vec, weight_vec
  mu_vec, sigma_vec, weight_vec = unwrap_parameters(par)
  gmm_psf =
    T[ evaluate_psf_at_point(x_mat[i, j], mu_vec, sigma_vec, weight_vec)
       for i=1:size(x_mat, 1), j=1:size(x_mat, 2) ]
  gmm_psf
end

"""
Fit a mixture of 2d Gaussians to a PSF image (evaluated at a single point).

Args:
 - psf: An matrix image of the point spread function, e.g. as returned by
        get_psf_at_point.

"""
function fit_psf_gaussians_least_squares(
    psf::Array{Float64, 2}; initial_par=Float64[],
    grtol = 1e-9, ftol = 1e-9, iterations = 5000, verbose=false,
    K=2, optim_method=:nelder_mead)

  if (any(psf .< 0))
      if verbose
          warn("Some psf values are negative.")
      end
      psf[ psf .< 0 ] = 0
  end

  x_mat = get_x_matrix_from_psf(psf);

  if length(initial_par) == 0
    if verbose
      println("Using default initialization.")
    end
    psf_starting_mean  =
      sum([ psf[i, j]  * x_mat[i, j] for
          i in 1:size(psf, 1), j=1:size(psf, 2) ]) / sum(psf)

    mu_vec = Array(Vector{Float64}, K)
    sigma_vec = Array(Matrix{Float64}, K)
    weight_vec = zeros(Float64, K)

    for k=1:K
      mu_vec[k] = psf_starting_mean
      sigma_vec[k] = Float64[ sqrt(2 * k) 0; 0 sqrt(2 * k)]
      weight_vec[k] = 1 / K
    end
    initial_par = wrap_parameters(mu_vec, sigma_vec, weight_vec)
  end

  @assert length(initial_par) == K * 6

  function evaluate_fit{T <: Number}(par::Vector{T})
    gmm_psf = render_psf(par, x_mat)
    local fit = sum((psf .- gmm_psf) .^ 2)
    if verbose && T == Float64
      println("-------------------")
      println("Fit: $fit")
      mu_vec, sigma_vec, weight_vec = unwrap_parameters(par)
      println(mu_vec)
      println(sigma_vec)
      println(weight_vec)
    end
    fit
  end

  optim_result =
    Optim.optimize(evaluate_fit, initial_par, method=optim_method,
                   iterations=iterations, grtol=grtol, ftol=ftol)

  optim_result, unwrap_parameters(optim_result.minimum)...
end

"""
Using data from a psField file, evaluate the PSF for a source at given point.

Args:
 - row: The row of the point source (may be a float).
 - col: The column of the point source (may be a float).
 - raw_psf_comp: Raw PSF components (e.g. as read from a psField file)

Returns:
 - An rnrow x rncol image of the PSF at (row, col)

The PSF is represented as a weighted combination of "eigenimages" (stored
in rrows), where the weights vary smoothly across points (row, col) in the image
as a polynomial of the form
weight[k](row, col) = sum_{i,j} cmat[i, j, k] * (rcs * row) ^ i (rcs * col) ^ j
...where row and col are zero-indexed.

This function is based on the function sdss_psf_at_points in astrometry.net:
https://github.com/dstndstn/astrometry.net/blob/master/util/sdss_psf.py
"""
function get_psf_at_point(row::Float64, col::Float64,
                          raw_psf_comp::RawPSFComponents)

    # This is a coordinate transform to keep the polynomial coefficients
    # to a reasonable size.
    const rcs = 0.001

    # rrows' image data is in the first column a flattened form.
    # The second dimension is the number of eigen images, which should
    # match the number of coefficient arrays.
    cmat = raw_psf_comp.cmat

    k_tot = size(raw_psf_comp.rrows)[2]
    @assert k_tot == size(cmat)[3]

    nrow_b = size(cmat)[1]
    ncol_b = size(cmat)[2]

    # Get the weights.  The row and column are intended to be
    # zero-indexed.
    coeffs_mat = [ ((row - 1) * rcs) ^ i * ((col - 1) * rcs) ^ j for
                    i=0:(nrow_b - 1), j=0:(ncol_b - 1)]
    weight_mat = zeros(k_tot)
    for k = 1:k_tot, i = 1:nrow_b, j = 1:ncol_b
        weight_mat[k] += cmat[i, j, k] * coeffs_mat[i, j]
    end

    # Weight the images in rrows and reshape them into matrix form.
    # It seems I need to convert for reshape to work.  :(
    psf = reshape(sum([ raw_psf_comp.rrows[:, i] * weight_mat[i] for i=1:k_tot]),
                  (convert(Int64, raw_psf_comp.rnrow),
                   convert(Int64, raw_psf_comp.rncol)))

    psf
end

end


text:
using Test
import ParallelStencil
using ParallelStencil.ParallelKernel
import ParallelStencil.ParallelKernel: @reset_parallel_kernel, @is_initialized, SUPPORTED_PACKAGES, PKG_CUDA, PKG_THREADS
import ParallelStencil.ParallelKernel: @require, longnameof, @prettyexpand, prettystring, @gorgeousexpand, gorgeousstring
import ParallelStencil.ParallelKernel: checkargs_hide_communication, hide_communication, hide_communication_cuda
using ParallelStencil.ParallelKernel.Exceptions
TEST_PACKAGES = SUPPORTED_PACKAGES
@static if PKG_CUDA in TEST_PACKAGES
    import CUDA
    if !CUDA.functional() TEST_PACKAGES = filter!(x->x≠PKG_CUDA, TEST_PACKAGES) end
end

@static for package in TEST_PACKAGES  eval(:(
    @testset "$(basename(@__FILE__)) (package: $(nameof($package)))" begin
        @testset "1. hide_communication macro" begin
            @init_parallel_kernel($package, Float64)
            @require @is_initialized()
            @testset "@hide_communication boundary_width block (macro expansion)" begin
                @static if $package == $PKG_CUDA
                    block = string(@gorgeousexpand(1, @hide_communication(boundary_width, begin @parallel f(A, B); communication(); end)))
                    @test occursin("ranges_outer = ParallelStencil.ParallelKernel.get_ranges_outer(boundary_width, ParallelStencil.ParallelKernel.get_ranges(A, B))", block)
                    @test occursin("ranges_inner = ParallelStencil.ParallelKernel.get_ranges_inner(boundary_width, ParallelStencil.ParallelKernel.get_ranges(A, B))", block)
                    @test occursin("@hide_communication ranges_outer ranges_inner begin\n            @parallel f(A, B)\n            communication()\n        end\nend", block)
                end;
            end;
            @testset "@hide_communication" begin
                @require @is_initialized()
                @parallel_indices (ix,iy,iz) function add_indices!(A)
                    A[ix,iy,iz] = A[ix,iy,iz] + ix + (iy-1)*size(A,1) + (iz-1)*size(A,1)*size(A,2); # NOTE: $ix, $iy, $iz come from ParallelStencil.INDICES.
                    return
                end
                function communication!(A)
                    A[  1,  :,  :] .= A[end-1,    :,    :]
                    A[end,  :,  :] .= A[    2,    :,    :]
                    A[  :,  1,  :] .= A[    :,end-1,    :]
                    A[  :,end,  :] .= A[    :,    2,    :]
                    A[  :,  :,  1] .= A[    :,    :,end-1]
                    A[  :,  :,end] .= A[    :,    :,    2]
                    return A
                end
                function communication_y!(A)
                    A[  :,  1,  :] .= A[    :,end-1,    :]
                    A[  :,end,  :] .= A[    :,    2,    :]
                    return A
                end
                function communication_z!(A)
                    A[  :,  :,  1] .= A[    :,    :,end-1]
                    A[  :,  :,end] .= A[    :,    :,    2]
                    return A
                end
                @testset "@hide_communication boundary_width block" begin
                    A  = @zeros(6, 7, 8)
                    @hide_communication (2,2,3) begin
                        @parallel add_indices!(A);
                        communication!(A);
                    end
                    @test all(Array(A) .== communication!([ix + (iy-1)*size(A,1) + (iz-1)*size(A,1)*size(A,2) for ix=1:size(A,1), iy=1:size(A,2), iz=1:size(A,3)]))
                end;
                @testset "@hide_communication boundary_width block" begin  # This test verifies that the results are correct even though in the current version (using CUDA.jl < v2.0), it cannot overlap
                    A  = @zeros(6, 7, 8)
                    @hide_communication (2,2,3) begin
                        @parallel add_indices!(A);
                        A[  1,  :,  :] .= A[end-1,    :,    :]
                        if true
                            A[end,  :,  :] .= A[    2,    :,    :]
                        end
                        communication_y!(A);
                        communication_z!(A);
                    end
                    @test all(Array(A) .== communication!([ix + (iy-1)*size(A,1) + (iz-1)*size(A,1)*size(A,2) for ix=1:size(A,1), iy=1:size(A,2), iz=1:size(A,3)]))
                end;
                @testset "@hide_communication boundary_width block" begin
                    A  = @zeros(6, 7, 8)
                    @parallel_indices (iy,iz) bc_xl(A) = (A[  1,iy,iz]=A[end-1,iy,iz]; return)
                    @parallel_indices (iy,iz) bc_xr(A) = (A[end,iy,iz]=A[    2,iy,iz]; return)
                    @hide_communication (2,2,3) begin
                        @parallel add_indices!(A);
                        @parallel (1:size(A,2), 1:size(A,3)) bc_xl(A);
                        if true
                            @parallel (1:size(A,2), 1:size(A,3)) bc_xr(A);
                        end
                        communication_y!(A);
                        communication_z!(A);
                    end
                    @test all(Array(A) .== communication!([ix + (iy-1)*size(A,1) + (iz-1)*size(A,1)*size(A,2) for ix=1:size(A,1), iy=1:size(A,2), iz=1:size(A,3)]))
                end;
                @testset "@hide_communication ranges_outer ranges_inner block" begin
                    A  = @zeros(6, 7, 8)
                    ranges_outer = ParallelStencil.ParallelKernel.get_ranges_outer((1, 1, 2), ParallelStencil.ParallelKernel.get_ranges(A))
                    ranges_inner = ParallelStencil.ParallelKernel.get_ranges_inner((1, 1, 2), ParallelStencil.ParallelKernel.get_ranges(A))
                    @hide_communication ranges_outer ranges_inner begin
                        @parallel add_indices!(A);
                        communication!(A);
                    end
                    @test all(Array(A) .== communication!([ix + (iy-1)*size(A,1) + (iz-1)*size(A,1)*size(A,2) for ix=1:size(A,1), iy=1:size(A,2), iz=1:size(A,3)]))
                end;
            end;
            @reset_parallel_kernel()
        end;
        @testset "2. Exceptions" begin
            @init_parallel_kernel($package, Float64)
            @require @is_initialized
            @testset "arguments @hide_communication" begin
                @test_throws ArgumentError checkargs_hide_communication(:boundary_width, :block)               # Error: the last argument must be a code block.
                @test_throws ArgumentError checkargs_hide_communication(:ranges_outer, :ranges_inner, :block)  # Error: the last argument must be a code block.
                @static if $package == $PKG_CUDA
                    @test_throws ArgumentError hide_communication_cuda(:boundary_width, :(@parallel f()))                                                                      # Error: missing (bc and) communication code.
                    @test_throws ArgumentError hide_communication_cuda(:boundary_width, :(communication()))                                                                    # Error: missing @parallel call.
                    @test_throws ArgumentError hide_communication_cuda(:boundary_width, :(begin @parallel f(); end))                                                           # Error: missing (bc and) communication code.
                    @test_throws ArgumentError hide_communication_cuda(:boundary_width, :(begin communication(); end))                                                         # Error: missing @parallel call.
                    block = :(begin do_something(); @parallel f(); communication(); end);          @test_throws ArgumentError hide_communication_cuda(:boundary_width, block)  # Error: block does not start with @parallel call
                    block = :(begin @parallel stream=mystream f(); communication(); end);          @test_throws ArgumentError hide_communication_cuda(:boundary_width, block)  # Error: invalid keyword argument 'stream'.
                    block = :(begin @parallel ranges f(); communication(); end);                   @test_throws ArgumentError hide_communication_cuda(:boundary_width, block)  # Error: invalid optional argument 'ranges'.
                    block = :(begin @parallel nblocks nthreads f(); communication(); end);         @test_throws ArgumentError hide_communication_cuda(:boundary_width, block)  # Error: invalid optional arguments 'nblocks' and 'nthreads'.
                    block = :(begin @parallel ranges nblocks nthreads f(); communication(); end);  @test_throws ArgumentError hide_communication_cuda(:boundary_width, block)  # Error: invalid optional arguments 'ranges', 'nblocks' and 'nthreads'.

                    @test_throws ArgumentError hide_communication_cuda(:ranges_outer, :ranges_inner, :(@parallel f()))                                                                      # Error: missing (bc and) communication code.
                    @test_throws ArgumentError hide_communication_cuda(:ranges_outer, :ranges_inner, :(communication()))                                                                    # Error: missing @parallel call.
                    @test_throws ArgumentError hide_communication_cuda(:ranges_outer, :ranges_inner, :(begin @parallel f(); end))                                                           # Error: missing (bc and) communication code.
                    @test_throws ArgumentError hide_communication_cuda(:ranges_outer, :ranges_inner, :(begin communication(); end))                                                         # Error: missing @parallel call.
                    block = :(begin do_something(); @parallel f(); communication(); end);          @test_throws ArgumentError hide_communication_cuda(:ranges_outer, :ranges_inner, block)  # Error: block does not start with @parallel call
                    block = :(begin @parallel stream=mystream f(); communication(); end);          @test_throws ArgumentError hide_communication_cuda(:ranges_outer, :ranges_inner, block)  # Error: invalid keyword argument 'stream'.
                    block = :(begin @parallel ranges f(); communication(); end);                   @test_throws ArgumentError hide_communication_cuda(:ranges_outer, :ranges_inner, block)  # Error: invalid optional arguments 'ranges'.
                    block = :(begin @parallel nblocks nthreads f(); communication(); end);         @test_throws ArgumentError hide_communication_cuda(:ranges_outer, :ranges_inner, block)  # Error: invalid optional arguments 'nblocks' and 'nthreads'.
                    block = :(begin @parallel ranges nblocks nthreads f(); communication(); end);  @test_throws ArgumentError hide_communication_cuda(:ranges_outer, :ranges_inner, block)  # Error: invalid optional arguments 'ranges', 'nblocks' and 'nthreads'.

                    block = :(begin @parallel f(); @parallel g(); communication(); end);                          @test_throws ArgumentError hide_communication_cuda(:ranges_outer, :ranges_inner, block)  # Error: missing arguments 'ranges' in @parallel call for bc computations.
                    block = :(begin @parallel f(); @parallel stream=mystream g(); communication(); end);          @test_throws ArgumentError hide_communication_cuda(:ranges_outer, :ranges_inner, block)  # Error: invalid keyword argument in @parallel call for bc computations.
                    block = :(begin @parallel f(); @parallel nblocks nthreads g(); communication(); end);         @test_throws ArgumentError hide_communication_cuda(:ranges_outer, :ranges_inner, block)  # Error: invalid optional arguments 'nblocks' and 'nthreads' in @parallel call for bc computations.
                    block = :(begin @parallel f(); @parallel ranges nblocks nthreads g(); communication(); end);  @test_throws ArgumentError hide_communication_cuda(:ranges_outer, :ranges_inner, block)  # Error: invalid optional arguments 'ranges', 'nblocks' and 'nthreads' in @parallel call for bc computations.
                end
            end;
            @reset_parallel_kernel()
        end;
    end;
)) end == nothing || true;