text:
using DataFrames, SQLite
export ContinuousSpace, index!, update_space!

#######################################################################################
# Continuous space structure
#######################################################################################
struct ContinuousSpace{F, E} <: AbstractSpace
  D::Int
  update_vel!::F
  periodic::Bool
  extend::E
  metric::Symbol
  db::SQLite.DB
  insertq::SQLite.Stmt
  searchq::SQLite.Stmt
  deleteq::SQLite.Stmt
  updateq::SQLite.Stmt
  searchqNoId::SQLite.Stmt
end

const COORDS = 'a':'z' # letters representing coordinates in database

"""
    ContinuousSpace(D::Int [, update_vel!]; kwargs...)
Create a `ContinuousSpace` of dimensionality `D`.
In this case, your agent positions (field `pos`) should be of type `NTuple{D, F}`
where `F <: AbstractFloat`.
In addition, the agent type should have a third field `vel::NTuple{D, F}` representing
the agent's velocity to use [`move_agent!`](@ref).

The optional argument `update_vel!` is a **function**, `update_vel!(agent, model)` that updates
the agent's velocity **before** the agent has been moved, see [`move_agent!`](@ref).
You can of course change the agents' velocities
during the agent interaction, the `update_vel!` functionality targets arbitrary force
fields acting on the agents (e.g. some magnetic field).
By default no update is done this way.

Notice that if you need to write your own custom `move_agent` function, call
[`update_space!`](@ref) at the end, like in e.g. the [Bacterial Growth](@ref) example.

## Keywords
* `periodic = true` : whether continuous space is periodic or not
* `extend::NTuple{D} = ones` : Extent of space. The `d` dimension starts at 0
  and ends at `extend[d]`. If `periodic = true`, this is also when
  periodicity occurs. If `periodic ≠ true`, `extend` is only used at plotting.
* `metric = :cityblock` : metric that configures distances for finding nearest neighbors
  in the space. The other option is `:euclidean` but cityblock is faster (due to internals).

**Note:** if your model requires linear algebra operations for which tuples are not supported,
a performant solution is to convert between Tuple and SVector using
[StaticArrays.jl](https://github.com/JuliaArrays/StaticArrays.jl)
as follows: `s = SVector(t)` and back with `t = Tuple(s)`.
"""
function ContinuousSpace(D::Int, update_vel! = defvel;
  periodic = true, extend = Tuple(1.0 for i in 1:D), metric = :cityblock)

  @assert metric ∈ (:cityblock, :euclidean)
  @assert typeof(extend) <: NTuple{D} "extend must be a tuple of length(D)"
  @assert all(typeof(e) <: Real && e > 0 for e in extend) "extend must contain positive integers or floating point numbers"
  db, q, q2, q3, q4, q5 = prepare_database(D)
  ContinuousSpace(D, update_vel!, periodic, extend, metric, db, q, q2, q3, q4, q5)
end

# Deprecate Space constructor
@deprecate Space(D::Int, update_vel!::Function) ContinuousSpace(D::Int, update_vel!::Function)

function prepare_database(D)
  db = SQLite.DB()
  dimexpression = join("$x REAL, " for x in COORDS[1:D])
  stmt = "CREATE TABLE tab ("*dimexpression*"id INTEGER PRIMARY KEY)"
  DBInterface.execute(db, stmt)
  insertedxpression = join("$x, " for x in COORDS[1:D])
  qmarks = join("?, " for _ in 1:D)
  insertstmt = "INSERT INTO tab ($(insertedxpression)id) VALUES ($(qmarks)?)"
  q = DBInterface.prepare(db, insertstmt)
  searchexpr = join("$x BETWEEN ? AND ? AND " for x in COORDS[1:D])
  searchq = "SELECT id FROM tab WHERE $(searchexpr)id != ?"
  q2 = DBInterface.prepare(db, searchq)
  deleteq = "DELETE FROM tab WHERE id = ?"
  q3 = DBInterface.prepare(db, deleteq)
  updateexpr = strip(join("$x = ?," for x in COORDS[1:D]), ',')
  updateq = "UPDATE tab SET $updateexpr WHERE id = ?"
  q4 = DBInterface.prepare(db, updateq)
  searchexpr2 = join("$x BETWEEN ? AND ? AND " for x in COORDS[1:D])
  searchq2 = "SELECT id FROM tab WHERE $(searchexpr2)"[1:end-4]
  q5 = DBInterface.prepare(db, searchq2)
  return db, q, q2, q3, q4, q5
end

defvel(a, m) = nothing

function Base.show(io::IO, abm::ContinuousSpace)
    s = "$(abm.D)-dimensional $(abm.periodic ? "periodic " : "")ContinuousSpace"
    abm.update_vel! ≠ defvel && (s *= " with velocity updates")
    print(io, s)
end

#######################################################################################
# SQLite database functions
#######################################################################################
"Collect IDs from an SQLite.Query where IDs are stored in `colname`"
function collect_ids(q::SQLite.Query, colname=:id)
  output = Vector{Int}(undef, 0)
  for row in skipmissing(q)
    push!(output, row[colname])
  end
  return output
end

"""
    index!(model)
Index the database underlying the `ContinuousSpace` of the model.

This can drastically improve performance for finding neighboring agents, but adding new
data can become slower because after each addition, index needs to be called again.

Lack of index won't be noticed for small databases. Only use it when you have
many agents and not many additions of agents.
"""
function index!(model)
  D = model.space.D
  expr = join("$x," for x in COORDS[1:D])
  DBInterface.execute(model.space.db, "CREATE INDEX pos ON tab ($(expr)id)")
end

#######################################################################################
# add, move and kill
#######################################################################################
# central, low level function that is always called by all others!
function add_agent_pos!(agent::A, model::ABM{A, <: ContinuousSpace}) where {A<:AbstractAgent}
  DBInterface.execute(model.space.insertq, (agent.pos..., agent.id))
  model[agent.id] = agent
  return agent
end

function randompos(space::ContinuousSpace)
  pos = Tuple(rand(space.D))
  space.extend ≠ nothing && (pos = pos .* space.extend)
  return pos
end

function add_agent!(agent::A, model::ABM{A, <: ContinuousSpace}) where {A<:AbstractAgent}
  agent.pos = randompos(model.space)
  add_agent_pos!(agent, model)
end

function add_agent!(agent::A, pos, model::ABM{A, <: ContinuousSpace}) where {A<:AbstractAgent}
  agent.pos = pos
  add_agent_pos!(agent, model)
end

# versions that create the agents
function add_agent!(model::ABM{A, <: ContinuousSpace}, args...; kwargs...) where {A<:AbstractAgent}
  add_agent!(randompos(model.space), model, args...; kwargs...)
end

function add_agent!(pos::Tuple, model::ABM{A, <: ContinuousSpace}, args...; kwargs...) where {A<:AbstractAgent}
  id = nextid(model)
  agent = A(id, pos, args...; kwargs...)
  add_agent_pos!(agent, model)
end

"""
    move_agent!(agent::A, model::ABM{A, ContinuousSpace}, dt = 1.0)
Propagate the agent forwards one step according to its velocity,
_after_ updating the agent's velocity (see [`ContinuousSpace`](@ref)).
Also take care of periodic boundary conditions.

For this continuous space version of `move_agent!`, the "evolution algorithm"
is a trivial Euler scheme with `dt` the step size, i.e. the agent position is updated
as `agent.pos += agent.vel * dt`.

Notice that if you want the agent to instantly move to a specified position, do
`agent.pos = pos` and then `update_space!(agent, model)`.
"""
function move_agent!(agent::A, model::ABM{A, S, F, P}, dt = 1.0) where {A<:AbstractAgent, S <: ContinuousSpace, F, P}
  model.space.update_vel!(agent, model)
  agent.pos = agent.pos .+ dt .* agent.vel
  if model.space.periodic
    agent.pos = mod.(agent.pos, model.space.extend)
  end
  update_space!(model, agent)
  return agent.pos
end

"""
    move_agent!(agent::A, model::ABM{A, ContinuousSpace}, vel::NTuple{D, N}, dt = 1.0)
Propagate the agent forwards one step according to `vel` and the model's space, with `dt` as the time step. (`update_vel!` is not used)
"""
function move_agent!(agent::A, model::ABM{A,S,F,P}, vel::NTuple{D, N}, dt = 1.0) where {A <: AbstractAgent, S <: ContinuousSpace, F, P, D, N <: AbstractFloat}
      agent.pos = agent.pos .+ dt .* vel
      if model.space.periodic
        agent.pos = mod.(agent.pos, model.space.extend)
      end
      update_space!(model, agent)
      return agent.pos
end

"""
    update_space!(model::ABM{A, ContinuousSpace}, agent)
Update the internal representation of continuous space to match the new position of
the agent (useful in custom `move_agent` functions).
"""
function update_space!(model::ABM{A,<:ContinuousSpace}, agent) where {A}
  DBInterface.execute(model.space.updateq, (agent.pos..., agent.id))
end

function kill_agent!(agent::AbstractAgent, model::ABM{A,<:ContinuousSpace}) where {A}
  DBInterface.execute(model.space.deleteq, (agent.id,))
  delete!(model.agents, agent.id)
  return model
end

function genocide!(model::ABM{A,<:ContinuousSpace}, n::Int) where {A}
  ids = strip(join("$id," for id in keys(model.agents) if id > n), ',')
  DBInterface.execute(model.space.db, "DELETE FROM tab WHERE id IN ($ids)")
  for id in keys(model.agents)
    id > n && (delete!(model.agents, id))
  end
  return model
end

function genocide!(model::ABM{A,<:ContinuousSpace}, f::Function) where {A}
  ids = strip(join("$(agent.id)," for agent in values(model.agents) if f(agent)), ',')
  DBInterface.execute(model.space.db, "DELETE FROM tab WHERE id IN ($ids)")
  for agent in values(model.agents)
    f(agent) && (delete!(model.agents, agent.id))
  end
  return model
end

function genocide!(model::ABM{A,<:ContinuousSpace}) where {A}
  DBInterface.execute(model.space.db, "DELETE FROM tab")
  for agent in values(model.agents)
    delete!(model.agents, agent.id)
  end
end

#######################################################################################
# neighboring agents
#######################################################################################
function space_neighbors(pos::Tuple, model::ABM{A, <:ContinuousSpace}, r::Real) where {A}
  left = pos .- r
  right = pos .+ r
  res = interlace(left, right)
  ids = collect_ids(DBInterface.execute(model.space.searchqNoId, res))
  if model.space.metric == :euclidean
    filter!(i -> sqrt(sum(abs2.(model[i].pos .- pos))) ≤ r, ids)
  end
  ids
end

function space_neighbors(agent::A, model::ABM{A, <:ContinuousSpace}, r::Real) where {A<:AbstractAgent}
  left = agent.pos .- r
  right = agent.pos .+ r
  res = interlace(left, right)
  ids = collect_ids(DBInterface.execute(model.space.searchq, (res...,agent.id)))
  if model.space.metric == :euclidean
    filter!(i -> sqrt(sum(abs2.(model[i].pos .- agent.pos))) ≤ r, ids)
  end
  ids
end

@generated function interlace(left::NTuple{D}, right::NTuple{D}) where {D}
  a = [[:(left[$i]), :(right[$i])] for i=1:D]
  b = vcat(a...)
  quote
    tuple($(b...))
  end
end

#######################################################################################
# Continuous space exclusive
#######################################################################################
export nearest_neighbor, elastic_collision!, interacting_pairs

"""
    nearest_neighbor(agent, model, r) → nearest
Return the agent that has the closest distance to given `agent`, according to the
space's metric. Valid only in continuous space.
Return `nothing` if no agent is within distance `r`.
"""
function nearest_neighbor(agent, model, r)
  n = space_neighbors(agent, model, r)
  length(n) == 0 && return nothing
  d, j = Inf, 0
  for i in 1:length(n)
    if model.space.metric == :euclidean
      @inbounds dnew = sqrt(sum(abs2.(agent.pos .- model[n[i]].pos)))
    elseif model.space.metric == :cityblock
      @inbounds dnew = sum(abs.(agent.pos .- model[n[i]].pos))
    end
    _, j = findmin(d)
    if dnew < d
      d, j = dnew, i
    end
  end
  return @inbounds model[n[j]]
end

using LinearAlgebra

"""
    elastic_collision!(a, b, f = nothing)
Resolve a (hypothetical) elastic collision between the two agents `a, b`.
They are assumed to be disks of equal size touching tangentially.
Their velocities (field `vel`) are adjusted for an elastic collision happening between them.
This function works only for two dimensions.
Notice that collision only happens if both disks face each other, to avoid
collision-after-collision.

If `f` is a `Symbol`, then the agent property `f`, e.g. `:mass`, is taken as a mass
to weight the two agents for the collision. By default no weighting happens.

One of the two agents can have infinite "mass", and then acts as an immovable object
that specularly reflects the other agent. In this case of course momentum is not
conserved, but kinetic energy is still conserved.
"""
function elastic_collision!(a, b, f = nothing)
  # Do elastic collision according to
  # https://en.wikipedia.org/wiki/Elastic_collision#Two-dimensional_collision_with_two_moving_objects
  v1, v2, x1, x2 = a.vel, b.vel, a.pos, b.pos
  length(v1) != 2 && error("This function works only for two dimensions.")
  r1 = x1 .- x2; r2 = x2 .- x1
  m1, m2 = f == nothing ? (1.0, 1.0) : (getfield(a, f), getfield(b, f))
  # mass weights
  m1 == m2 == Inf && return false
  if m1 == Inf
    @assert v1 == (0, 0) "An agent with ∞ mass cannot have nonzero velocity"
    dot(r1, v2) ≤ 0 && return false
    v1 = ntuple(x -> zero(eltype(v1)), length(v1))
    f1, f2 = 0.0, 2.0
  elseif m2 == Inf
    @assert v2 == (0, 0) "An agent with ∞ mass cannot have nonzero velocity"
    dot(r2, v1) ≤ 0 && return false
    v2 = ntuple(x -> zero(eltype(v1)), length(v1))
    f1, f2 = 2.0, 0.0
  else
    # Check if disks face each other, to avoid double collisions
    !(dot(r2, v1) > 0 && dot(r2, v1) > 0) && return false
    f1 = (2m2/(m1+m2))
    f2 = (2m1/(m1+m2))
  end
  ken = norm(v1)^2 + norm(v2)^2
  dx = a.pos .- b.pos
  dv = a.vel .- b.vel
  n = norm(dx)^2
  n == 0 && return false # do nothing if they are at the same position
  a.vel = v1 .- f1 .* ( dot(v1 .- v2, r1) / n ) .* (r1)
  b.vel = v2 .- f2 .* ( dot(v2 .- v1, r2) / n ) .* (r2)
  return true
end

"""
    interacting_pairs(model, r, method; scheduler = model.scheduler)
Return an iterator that yields unique pairs of agents `(a1, a2)` that are close
neighbors to each other, within some interaction radius `r`.

This function is usefully combined with `model_step!`, when one wants to perform
some pairwise interaction across all pairs of close agents once
(and does not want to trigger the event twice, both with `a1` and with `a2`, which
is unavoidable when using `agent_step!`).

The argument `method` provides three pairing scenarios
- `:all`: return every pair of agents that are within radius `r` of each other,
  not only the nearest ones.
- `:nearest`: agents are only paired with their true nearest neighbor
  (existing within radius `r`).
  Each agent can only belong to one pair, therefore if two agents share the same nearest
  neighbor only one of them (sorted by id) will be paired.
- `:scheduler`: agents are scanned according to the given keyword `scheduler`
  (by default the model's scheduler), and each scanned
  agent is paired to its nearest neighbor. Similar to `:nearest`, each agent can belong
  to only one pair. This functionality is useful e.g. when you want some agents to be
  paired "guaranteed", even if some other agents might be nearest to each other.
- `:types`: For mixed agent models only. Return every pair of agents within radius `r`
  (similar to `:all`), only capturing pairs of differing types. For example, a model of
  `Union{Sheep,Wolf}` will only return pairs of `(Sheep, Wolf)`. In the case of multiple
  agent types, *e.g.* `Union{Sheep, Wolf, Grass}`, skipping pairings that involve
  `Grass`, can be achived by a [`scheduler`](@ref Schedulers) that doesn't schedule `Grass`
  types, *i.e.*: `scheduler = [a.id for a in allagents(model) of !(a isa Grass)]`.
"""
function interacting_pairs(model::ABM, r::Real, method; scheduler = model.scheduler)
    @assert method ∈ (:scheduler, :nearest, :all, :types)
    pairs = Tuple{Int,Int}[]
    if method == :nearest
        true_pairs!(pairs, model, r)
    elseif method == :scheduler
        scheduler_pairs!(pairs, model, r, scheduler)
    elseif method == :all
        all_pairs!(pairs, model, r)
    elseif method == :types
        type_pairs!(pairs, model, r, scheduler)
    end
    return PairIterator(pairs, model.agents)
end

function scheduler_pairs!(pairs::Vector{Tuple{Int,Int}}, model::ABM, r::Real, scheduler)
    #TODO: This can be optimized further I assume
    for id in scheduler(model)
        # Skip already checked agents
        any(isequal(id), p[2] for p in pairs) && continue
        a1 = model[id]
        a2 = nearest_neighbor(a1, model, r)
        # This line ensures each neighbor exists in only one pair:
        if a2 ≠ nothing && !any(isequal(a2.id), p[2] for p in pairs)
            push!(pairs, (id, a2.id))
        end
    end
end

function all_pairs!(pairs::Vector{Tuple{Int,Int}}, model::ABM, r::Real)
    for a in allagents(model)
        for nid in space_neighbors(a, model, r)
            # Sort the pair to overcome any uniqueness issues
            new_pair = isless(a.id, nid) ? (a.id, nid) : (nid, a.id)
            new_pair ∉ pairs && push!(pairs, new_pair)
        end
    end
end

function true_pairs!(pairs::Vector{Tuple{Int,Int}}, model::ABM, r::Real)
    distances = Vector{Float64}(undef, 0)
    for a in allagents(model)
        nn = nearest_neighbor(a, model, r)
        nn == nothing && continue
        # Sort the pair to overcome any uniqueness issues
        new_pair = isless(a.id, nn.id) ? (a.id, nn.id) : (nn.id, a.id)
        if new_pair ∉ pairs
            # We also need to check if our current pair is closer to each
            # other than any pair using our first id already in the list,
            # so we keep track of nn distances.
            dist = pair_distance(a.pos, model[nn.id].pos, model.space.metric)

            idx = findfirst(x -> first(new_pair) == x, first.(pairs))
            if idx == nothing
                push!(pairs, new_pair)
                push!(distances, dist)
            elseif idx != nothing && distances[idx] > dist
                # Replace this pair, it is not the true neighbor
                pairs[idx] = new_pair
                distances[idx] = dist
            end
        end
    end
end

function type_pairs!(pairs::Vector{Tuple{Int,Int}}, model::ABM, r::Real, scheduler)
    # We don't know ahead of time what types the scheduler will provide. Get a list.
    available_types = unique(typeof(model[id]) for id in scheduler(model))
    for id in scheduler(model)
        for nid in space_neighbors(model[id], model, r)
            neigbor_type = typeof(model[nid])
            if neigbor_type ∈ available_types && neigbor_type !== typeof(model[id])
                # Sort the pair to overcome any uniqueness issues
                new_pair = isless(id, nid) ? (id, nid) : (nid, id)
                new_pair ∉ pairs && push!(pairs, new_pair)
            end
        end
    end
end

function pair_distance(pos1, pos2, metric::Symbol)
    if metric == :euclidean
        sqrt(sum(abs2.(pos1 .- pos2)))
    elseif metric == :cityblock
        sum(abs.(pos1 .- pos2))
    end
end

struct PairIterator{A}
    pairs::Vector{Tuple{Int,Int}}
    agents::Dict{Int,A}
end

Base.length(iter::PairIterator) = length(iter.pairs)
function Base.iterate(iter::PairIterator, i = 1)
    i > length(iter) && return nothing
    p = iter.pairs[i]
    id1, id2 = p
    return (iter.agents[id1], iter.agents[id2]), i + 1
end


text:
DAPR_VERSION = v"1.4.2"

using Downloads
using CodecZlib
using Tar
using gRPCClient
using ProtoBuf

download_url = "https://github.com/dapr/dapr/archive/refs/tags/v$DAPR_VERSION.tar.gz"

open(Downloads.download(download_url)) do tar_gz
    tar = GzipDecompressorStream(tar_gz)
    dir = Tar.extract(tar)
    proto_dir = joinpath(dir, "dapr-$DAPR_VERSION")

    includes = [
        proto_dir,
        joinpath(pkgdir(ProtoBuf), "gen")
    ]

    gRPCClient.generate(
        joinpath(proto_dir, "dapr", "proto", "runtime", "v1", "dapr.proto")
        ;outdir=joinpath(@__DIR__, "..", "src", "generated", "dapr"),
        includes=includes
    )

    gRPCClient.generate(
        joinpath(proto_dir, "dapr", "proto", "runtime", "v1", "appcallback.proto")
        ;outdir=joinpath(@__DIR__, "..", "src", "generated", "appcallback"),
        includes=includes
    )

    gRPCClient.generate(
        joinpath(proto_dir, "dapr", "proto", "sentry", "v1", "sentry.proto")
        ;outdir=joinpath(@__DIR__, "..", "src", "generated", "sentry"),
        includes=includes
    )

    gRPCClient.generate(
        joinpath(proto_dir, "dapr", "proto", "operator", "v1", "operator.proto")
        ;outdir=joinpath(@__DIR__, "..", "src", "generated", "operator"),
        includes=includes
    )

    gRPCClient.generate(
        joinpath(proto_dir, "dapr", "proto", "placement", "v1", "placement.proto")
        ;outdir=joinpath(@__DIR__, "..", "src", "generated", "placement"),
        includes=includes
    )

    gRPCClient.generate(
        joinpath(proto_dir, "dapr", "proto", "internals", "v1", "service_invocation.proto")
        ;outdir=joinpath(@__DIR__, "..", "src", "generated", "service_invocation"),
        includes=includes
    )

end

#=
Following up:

1. rename service related module name & filename
2. Add `using ProtoBuf:google` in each client file
3. remove unnecessary google related files
=#

ext: jl


text:
module StaticContracts

#-- Modules
# wiring diagrams
using Catlab.WiringDiagrams
import Catlab.WiringDiagrams: oapply  # Needed to overwrite oapply function

# intervals of real numbers
using IntervalSets

#-- Accessible functions
export StaticContract, oapply, StaticContractDict

#-- Contracts are defined via intervals:
struct StaticContract{T<:Real}
    input::Vector{Interval}
    output::Vector{Interval}

    # Constructor
    function StaticContract{T}(input::Vector, output::Vector) where T<:Real
        # cannot have empty contract
        for contract in [input; output]
            if isempty(contract)
                error("the interval $contract is empty or backwards")
            end
        end
        new{T}(input, output)
    end
end

#-- Display struct as product of intervals:
function formatInterval(contract::Interval)
    # Switch symbol for infinity
    if -contract.left == contract.right == Inf
        output = "ℝ"
    else
        # check left bound
        if( contract.left == -Inf )
            left = "-∞"
            lbracket = "("
        else
            left = contract.left
            lbracket = isleftopen(contract) ? "(" : "["
        end

        # check right bound
        if( contract.right == Inf )
            right = "∞"
            rbracket = ")"
        else
            right = contract.right
            rbracket = isrightopen(contract) ? ")" : "]"
        end

        # Concatenate strings
        output = "$lbracket" * "$left,$right" * "$rbracket"
    end
    return output
end

function Base.show(io::IO, vf::StaticContract)
    # Check all contracts
    list = [vf.input; vf.output]
    output = ""

    # iterate across all contracts
    len = length(list)

    for i in 1:len
        contract = list[i]
        output *= formatInterval(contract)
        # add product operator [spaces make it easier to read]
        if i != len
            output *= " × "
        end
    end
    # display combined string
    print("StaticContract( $output )")
end

#-- Compose contracts with diagram:
function oapply(d::WiringDiagram, ms::Vector{StaticContract{T}}) where T
    # boxes in diagram
    box = boxes(d)

    # Ensure machines fill diagram
    n_contract = length(ms)

    if nboxes(d) != n_contract
        error("there are $nboxes(d) boxes but $length(ms) machines")
    end

    # Boxes must have unique names to be identified
    if n_contract > 1
        name = map( m -> m.value, box )
        for i in 1:(n_contract-1)
            curr_name = name[i]
            if curr_name in name[ (i+1):n_contract ]
                error("Two or more boxes share the name: $curr_name")
            end
        end
    end

    # Each wire must have a contract
    for id in 1:nboxes(d)
        curr_box = box[id]
        if length(curr_box.input_ports) != length(ms[id].input) || length(curr_box.output_ports) != length(ms[id].output)
            name = curr_box.value
            error("number of ports do not match number of contracts at $name (id=$id)")
        end
    end

    # Check all contracts inside the diagram (excluding wires that enter or exit the diagram)
    for w in wires(d, :Wire)
        # id number of boxes
        source_id = w.source.box
        target_id = w.target.box
        # box names
        source_name = box[source_id].value
        target_name = box[target_id].value

        # 1. Variable names, see section 3.3.3, P.16
        source_var = box[source_id].output_ports[w.source.port]
        target_var = box[target_id].input_ports[w.target.port]

            # ensure target and source name match.
        if target_var != source_var
            error("wire \"$target_var\" of $source_name (id=$source_id) " *
                  "does not match wire \"$source_var\" of $target_name (id=$target_id)")
        end

        # 2. Induced contract, see section 3.3.1, P.13
        contract_source = ms[source_id].output[w.source.port]
        contract_target = ms[target_id].input[w.target.port]
        overlap = intersect(contract_source, contract_target)

        format_source = formatInterval(contract_source)
        format_target = formatInterval(contract_target)

            # check if source contract is compatible with target contract
        if isempty(overlap) == true
            error("Incompatible contract between $source_name (id=$source_id) " *
                  "and $target_name (id=$target_id) at wire \"$target_var\": " *
                  "$format_source ∩ $format_target = ∅")

            # Check for undefined behavior
        elseif overlap != contract_source
            println("Intervals do not overlap between $source_name (id=$source_id) " *
                    "and $target_name (id=$target_id) at wire \"$target_var\": " *
                    "$format_source ∩ $format_target ≠ $format_source")
        end
    end

    # Assign external contracts
    input = map( w -> ms[w.target.box].input[w.target.port], wires(d, :InWire))
    output = map( w -> ms[w.source.box].output[w.source.port], wires(d, :OutWire))

    # Create new machine
    return StaticContract{T}(input, output)
end

# Compose using a library
const StaticContractDict{T} = Dict{ Symbol, StaticContract{T} }

function oapply(d::WiringDiagram, ms::StaticContractDict{T}) where T
    box = boxes(d)
    name_contract = keys(ms)

    # add unassigned contracts
    if length(box) != length(name_contract)
        for b in box
            name = b.value
            if !(name in name_contract)
                # create default box
                n_input = length(b.input_ports)
                n_output = length(b.output_ports)
                default_contract = StaticContract{T}(repeat([-Inf..Inf], n_input),
                                                     repeat([-Inf..Inf], n_output) )
                # add to dictionary
                ms[name] = default_contract
            end
        end
    end
    # compose with array
    oapply(d, map(box -> ms[box.value], box) )
end

end # module


text:
using BlackBoxOptim

rosenbrock(x) = sum( 100*( x[2:end] .- x[1:end-1].^2 ).^2 .+ ( x[1:end-1] .- 1 ).^2 )

const MyFitnessGoal = 30.0

function myfitnessgoalachieved(oc)
    best_fitness(oc) < MyFitnessGoal
end

function cbearlystopping(oc)
    if myfitnessgoalachieved(oc)
        BlackBoxOptim.shutdown!(oc)
    end
end

# We give max 1000 seconds to optimize but expect this to stop much earlier
@time res = bboptimize(rosenbrock; SearchRange = (-100.0, 100.0), NumDimensions = 100,
        CallbackFunction = cbearlystopping, CallbackInterval = 0.0, MaxTime = 1000.0);

@assert best_fitness(res) < MyFitnessGoal
@assert BlackBoxOptim.stop_reason(res) == "Run explicitly stopped via shutdown method"


text:
"""
    function elo_ranking(
        players::Vector{Player},
        ranking_scheme_to_update::String,
        last_tournament_id::String
    )

Update the player's rankings using Elo Ranking system.
Computed after the completion of a tournament.
A player's ranking is not updated after the result of each match.
The ranking for each team is computed as average of individual player rankings.
K factor is calculated by dividing 800 by the number of games a player's rating is based on
plus the number of games the player completed in a tournament.

Arguments
=========
- `players`:
    an array of players in the simulation.
- `ranking_scheme_to_update`:
    the name of the ranking scheme to update for each Player.
- `last_tournament_id`:
    the last tournament that was completed.
"""

function elo_ranking(
    players::Vector{Player},
    ranking_scheme_to_update::String,
    last_tournament_id::String
)
    updatedRankings = Dict{Int64, Int64}();
    for player in players
        if haskey(player.doubles_results, last_tournament_id)
            expected_score = 0
            actual_score = 0
            teammate_id = player.doubles_results[last_tournament_id].partner_player_id
            teammate_ranking = search_elo_ranking(teammate_id, players)
            for match in player.doubles_results[last_tournament_id].matches
                RA = (player.rankings["Elo"] + teammate_ranking) / 2
                RB = search_elo_ranking(match.opponent1_player_id, players) + search_elo_ranking(match.opponent2_player_id, players)
                EA = 1 / (1 + 10^((RB-RA)/400))
                SA = match.result == "Won" ? 1 : 0
                expected_score += EA
                actual_score += SA
            end
            k_factor = 800 / (player.matches_played + length(player.doubles_results[last_tournament_id].matches))
            updatedRankings[player.id] = player.rankings["Elo"] + k_factor*(actual_score-expected_score)
        else
            # If player didn't participate in tournament, their ranking doesn't change
            updatedRankings[player.id] = player.rankings["Elo"];
        end
    end
    map(player -> player.rankings[ranking_scheme_to_update] = updatedRankings[player.id], players)
    return nothing
end

"""
    function search_elo_ranking(
        player_id::Int64,
        players::Vector{Player}
    )::Float64

Search for a player's current Elo Ranking

Arguments
=========
- `player_id`:
    the player whose ranking is under search.
- `players`:
    an array of players in the simulation.
"""

function search_elo_ranking(
    player_id::Int64,
    players::Vector{Player}
)::Float64
    player_ranking = filter(p -> p.id == player_id, players)[1].rankings["Elo"]
    return player_ranking
end

max_stars_repo_name: Danial-Hussain/Tennis-Ranking

max_issues_repo_name: Danial-Hussain/Tennis-Ranking

max_forks_repo_name: Danial-Hussain/Tennis-Ranking


text:
@deprecate(
    JLSOFile(
        version::VersionNumber,
        julia::VersionNumber,
        format::Symbol,
        compression::Symbol,
        image::String,
        pkgs::Dict{String, VersionNumber},
        objects::Dict{String, Vector{UInt8}},
    ),
    JLSOFile(
        version,
        julia,
        format,
        compression,
        image,
        Pkg.TOML.parse.(_upgrade_env(pkgs))...,
        Dict{Symbol, Vector{UInt8}}(Symbol(k) => v for (k, v) in objects),
    )
)

@deprecate JLSOFile(data; kwargs...) JLSOFile(:data => data; kwargs...)
@deprecate(
    JLSOFile(data::Dict{String, <:Any}; kwargs...),
    JLSOFile(Dict(Symbol(k) => v for (k, v) in data); kwargs...)
)

function upgrade_jlso(raw_dict::AbstractDict)
    result = copy(raw_dict)
    version = version_number(result["metadata"]["version"])

    while version < v"3"
        result = upgrade_jlso(result, Val(Int(version.major)))
        version = version_number(result["metadata"]["version"])
    end

    return result
end

@deprecate getindex(jlso::JLSOFile, name::String) getindex(jlso, Symbol(name))
@deprecate setindex!(jlso::JLSOFile, value, name::String) setindex!(jlso, value, Symbol(name))

# v1 and v2 stored version numbers but v3 stores strings to be bson parser agnostic
version_number(x::VersionNumber) = x
version_number(x::String) = VersionNumber(x)

# Metadata changes to upgrade file from v1 to v2
function upgrade_jlso(raw_dict::AbstractDict, ::Val{1})
    info(LOGGER, "Upgrading JLSO format from v1 to v2")
    metadata = copy(raw_dict["metadata"])
    version = version_number(metadata["version"])
    @assert version ∈ semver_spec("1")

    if metadata["format"] == :serialize
        metadata["format"] = :julia_serialize
    end

    metadata["compression"] = :none
    metadata["version"] = v"2"

    return Dict("metadata" => metadata, "objects" => get(raw_dict, "objects", Dict()))
end

# Metadata changes to upgrade from v2 to v3
function upgrade_jlso(raw_dict::AbstractDict, ::Val{2})
    info(LOGGER, "Upgrading JLSO format from v2 to v3")
    metadata = copy(raw_dict["metadata"])
    version = version_number(metadata["version"])
    @assert version ∈ semver_spec("2")

    # JLSO v3 expects to be loading a compressed manifest
    compression = get(metadata, "compression", :none)
    project, manifest = _upgrade_env(get(metadata, "pkgs", nothing))

    return Dict{String, Dict}(
        "metadata" => Dict{String, Union{String, Vector{UInt8}}}(
            "version" => "3",
            "julia" => string(get(metadata, "julia", VERSION)),
            "format" => string(get(metadata, "format", :julia_serialize)),
            "compression" => string(compression),
            "image" => get(metadata, "image", ""),
            "project" => project,
            "manifest" => read(compress(compression, IOBuffer(manifest))),
        ),
        "objects" => Dict{String, Vector{UInt8}}(
            k => v for (k, v) in get(raw_dict, "objects", Dict())
        ),
    )
end

# Used to convert the old "PkgName" => VersionNumber metadata to a
# Project.toml and Manifest.toml file.
function _upgrade_env(pkgs::Dict)
    src_env = Base.active_project()

    try
        mktempdir() do tmp
            Pkg.activate(tmp)

            # We construct an array of PackageSpecs to avoid ordering problems with
            # adding each package individually
            try
                Pkg.add([
                    Pkg.PackageSpec(; name=key, version=value) for (key, value) in pkgs
                ])
            catch e
                # Warn about failure and fallback to simply trying to install the pacakges
                # without version constraints.
                warn(LOGGER) do
                    "Failed to construct an environment with the provide package version " *
                    "($pkgs): $e.\n Falling back to simply adding the packages."
                end
                Pkg.add([Pkg.PackageSpec(; name=key) for (key, value) in pkgs])
            end

            return _env()
        end
    finally
        Pkg.activate(src_env)
    end
end

_upgrade_env(::Nothing) = ("", "")


text:
function masternodes(porder::Int,dim::Int,elemtype::Int)

d0 = pwd();
ii = findlast("Exasim", d0);

# filename = string(d0[1:ii[end]],"Preprocessing/masternodes.mat");
# vars = MAT.matread(filename);
#
# pelem1 = vars["pelem"][elemtype+1,porder,dim];
# telem1 = vars["telem"][elemtype+1,porder,dim];
# pface1 = vars["pface"][elemtype+1,porder,dim];
# tface1 = vars["tface"][elemtype+1,porder,dim];
# perm1 = vars["perm"][elemtype+1,porder,dim];
#
# telem1 = convert(Array{Int}, telem1);
# tface1 = convert(Array{Int}, tface1);
# perm1 = convert(Array{Int}, perm1);

fn = string(d0[1:ii[end]],"/Version0.1/Julia/Preprocessing/masternodes.bin");
tmp = reinterpret(Float64,read(fn));

ndims = Int(tmp[1]);
k1 = 2; k2 = (k1-1)+(ndims);
narrays = Int.(tmp[k1:k2]);
k1 = k2+1; k2 = (k1-1) + prod(narrays);
sz1 = reshape(Int.(tmp[k1:k2]), (narrays...,));
k1 = k2+1; k2 = (k1-1) + prod(narrays);
sz2 = reshape(Int.(tmp[k1:k2]), (narrays...,));
sz = sz1.*sz2;
lz = cumsum([0; sz[:]]);
lz1 = reshape(lz[1:end-1],(narrays...,));
lz2 = reshape(lz[2:end],(narrays...,));

e = elemtype+1;
i = 1;
m1 = k2+1 + lz1[i,e,porder,dim];
m2 = m1-1 + lz2[i,e,porder,dim]-lz1[i,e,porder,dim];
pelem = reshape(tmp[m1:m2],sz1[i,e,porder,dim],sz2[i,e,porder,dim]);

i = 2;
m1 = k2+1 + lz1[i,e,porder,dim];
m2 = m1-1 + lz2[i,e,porder,dim]-lz1[i,e,porder,dim];
telem = reshape(tmp[m1:m2],sz1[i,e,porder,dim],sz2[i,e,porder,dim]);

i = 3;
m1 = k2+1 + lz1[i,e,porder,dim];
m2 = m1-1 + lz2[i,e,porder,dim]-lz1[i,e,porder,dim];
pface = reshape(tmp[m1:m2],sz1[i,e,porder,dim],sz2[i,e,porder,dim]);

i = 4;
m1 = k2+1 + lz1[i,e,porder,dim];
m2 = m1-1 + lz2[i,e,porder,dim]-lz1[i,e,porder,dim];
tface = reshape(tmp[m1:m2],sz1[i,e,porder,dim],sz2[i,e,porder,dim]);

i = 5;
m1 = k2+1 + lz1[i,e,porder,dim];
m2 = m1-1 + lz2[i,e,porder,dim]-lz1[i,e,porder,dim];
perm = reshape(tmp[m1:m2],sz1[i,e,porder,dim],sz2[i,e,porder,dim]);

telem = convert(Array{Int}, telem);
tface = convert(Array{Int}, tface);
perm = convert(Array{Int}, perm);

# display(pelem - pelem1)
# display(telem - telem1)
# display(pface - pface1)
# display(tface - tface1)
# display(perm - perm1)

return pelem, telem, pface, tface, perm;

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
# Quadratic terms & constraints
#

function add_qpterms!(model::Model, qr::IVec, qc::IVec, qv::FVec)
    nnz = length(qr)
    (nnz == length(qc) == length(qv)) || error("Inconsistent argument dimensions.")

    if nnz > 0
        ret = @grb_ccall(addqpterms, Cint, (
            Ptr{Cvoid},    # model
            Cint,         # nnz
            Ptr{Cint},    # qrow
            Ptr{Cint},    # qcol
            Ptr{Float64}, # qval
            ),
            model, nnz, qr.-Cint(1), qc.-Cint(1), qv)

        if ret != 0
            throw(GurobiError(model.env, ret))
        end
    end
    nothing
end

function add_qpterms!(model::Model, qr::Vector, qc::Vector, qv::Vector)
    add_qpterms!(model, ivec(qr), ivec(qc), fvec(qv))
end

function add_qpterms!(model, H::SparseMatrixCSC{Float64}) # H must be symmetric
    n = num_vars(model)
    (H.m == n && H.n == n) || error("H must be an n-by-n symmetric matrix.")

    nnz_h = nnz(H)
    qr = Array{Cint}(undef, nnz_h)
    qc = Array{Cint}(undef, nnz_h)
    qv = Array{Float64}(undef, nnz_h)
    k = 0

    colptr::Vector{Int} = H.colptr
    nzval::Vector{Float64} = H.nzval

    for i = 1 : n
        qi::Cint = convert(Cint, i)
        for j = colptr[i]:(colptr[i+1]-1)
            qj = convert(Cint, H.rowval[j])

            if qi < qj
                k += 1
                qr[k] = qi
                qc[k] = qj
                qv[k] = nzval[j]
            elseif qi == qj
                k += 1
                qr[k] = qi
                qc[k] = qj
                qv[k] = nzval[j] * 0.5
            end
        end
    end

    add_qpterms!(model, qr[1:k], qc[1:k], qv[1:k])
end

function add_qpterms!(model, H::Matrix{Float64}) # H must be symmetric
    n = num_vars(model)
    size(H) == (n, n) || error("H must be an n-by-n symmetric matrix.")

    nmax = round(Int,n * (n + 1) / 2)
    qr = Array{Cint}(undef, nmax)
    qc = Array{Cint}(undef, nmax)
    qv = Array{Float64}(undef, nmax)
    k::Int = 0

    for i = 1 : n
        qi = convert(Cint, i)

        v = H[i,i]
        if v != 0.
            k += 1
            qr[k] = qi
            qc[k] = qi
            qv[k] = v * 0.5
        end

        for j = i+1 : n
            v = H[j, i]
            if v != 0.
                k += 1
                qr[k] = qi
                qc[k] = convert(Cint, j)
                qv[k] = v
            end
        end
    end

    add_qpterms!(model, qr[1:k], qc[1:k], qv[1:k])
end

function add_diag_qpterms!(model, H::Vector)  # H stores only the diagonal element
    n = num_vars(model)
    n == length(H) || error("Incompatible dimensions.")
    q = [convert(Cint,1):convert(Cint,n)]
    add_qpterms!(model, q, q, fvec(h))
end

function add_diag_qpterms!(model, hv::Real)  # all diagonal elements are H
    n = num_vars(model)
    q = [convert(Cint,1):convert(Cint,n)]
    add_qpterms!(model, q, q, fill(float64(hv), n))
end

function delq!(model::Model)
    ret = @grb_ccall(delq, Cint, (
        Ptr{Cvoid},    # model
        ),
        model)

    if ret != 0
        throw(GurobiError(model.env, ret))
    end
end

function getq(model::Model)
    nz = get_intattr(model, "NumQNZs")
    rowidx = Array{Cint}(undef, nz)
    colidx = Array{Cint}(undef, nz)
    val = Array{Float64}(undef, nz)
    nzout = Ref{Cint}()

    ret = @grb_ccall(getq, Cint, (
        Ptr{Cvoid},  # model
        Ptr{Cint},  # numqnzP
        Ptr{Cint},  # qrow
        Ptr{Cint},  # qcol
        Ptr{Float64}# qval
        ),
        model,nzout,rowidx,colidx,val)

    if ret != 0
        throw(GurobiError(model.env, ret))
    end

    return rowidx, colidx, val
end

# add_qconstr!

function add_qconstr!(model::Model, lind::IVec, lval::FVec, qr::IVec, qc::IVec, qv::FVec, rel::Cchar, rhs::Float64)
    qnnz = length(qr)
    qnnz == length(qc) == length(qv) || error("Inconsistent argument dimensions.")

    lnnz = length(lind)
    lnnz == length(lval) || error("Inconsistent argument dimensions.")

    if qnnz > 0 || lnnz > 0
        ret = @grb_ccall(addqconstr, Cint, (
            Ptr{Cvoid},    # model
            Cint,         # lnnz
            Ptr{Cint},    # lind
            Ptr{Float64}, # lval
            Cint,         # qnnz
            Ptr{Cint},    # qrow
            Ptr{Cint},    # qcol
            Ptr{Float64}, # qval
            Cchar,        # sense
            Float64,      # rhs
            Ptr{UInt8}    # name
            ),
            model, lnnz, lind.-Cint(1), lval, qnnz, qr.-Cint(1), qc.-Cint(1), qv, rel, rhs, C_NULL)

        if ret != 0
            throw(GurobiError(model.env, ret))
        end
    end
    nothing
end

function add_qconstr!(model::Model, lind::Vector, lval::Vector, qr::Vector, qc::Vector,
    qv::Vector{Float64}, rel::GChars, rhs::Real)

    add_qconstr!(model, ivec(lind), fvec(lval), ivec(qr), ivec(qc), fvec(qv), cchar(rel), Float64(rhs))
end

function delqconstrs!(model::Model, indices::Vector{Int})
    ret = @grb_ccall(delqconstrs, Cint, (Ptr{Cvoid}, Cint, Ptr{Cint}), model, length(indices), Cint.(indices.-1))
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
end

function getqconstr(model::Model, index::Int)
    # From the Gurobi documentation:
    # Typical usage is to call this routine twice. In the first call, you
    # specify the requested quadratic constraint, with NULL values for the array
    # arguments.
    affine_nnz    = Cint[0]
    quadratic_nnz = Cint[0]
    ret = @grb_ccall(getqconstr, Cint, (Ptr{Cvoid}, Cint,  # model, constraint
        Ptr{Cint}, Ptr{Cint}, Ptr{Float64},                  # affine terms
        Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Float64}        # quadratic terms
        ),
        model, Cint(index-1),
        affine_nnz, C_NULL, C_NULL,
        quadratic_nnz, C_NULL, C_NULL, C_NULL
    )
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    affine_cols     = fill(Cint(0), affine_nnz[1])
    affine_coefs    = fill(0.0, affine_nnz[1])
    quadratic_rows  = fill(Cint(0), quadratic_nnz[1])
    quadratic_cols  = fill(Cint(0), quadratic_nnz[1])
    quadratic_coefs = fill(0.0, quadratic_nnz[1])

    ret = @grb_ccall(getqconstr, Cint, (Ptr{Cvoid}, Cint,  # model, constraint
        Ptr{Cint}, Ptr{Cint}, Ptr{Float64},                  # affine terms
        Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Float64}        # quadratic terms
        ),
        model, Cint(index-1),
        affine_nnz, affine_cols, affine_coefs,
        quadratic_nnz, quadratic_rows, quadratic_cols, quadratic_coefs
    )
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    return affine_cols, affine_coefs, quadratic_rows, quadratic_cols, quadratic_coefs
end


text:
using Oceananigans.Solvers
using Oceananigans.Operators
using Oceananigans.Architectures
using Oceananigans.Grids: on_architecture
using Oceananigans.Fields: Field

using Oceananigans.Models.HydrostaticFreeSurfaceModels:
             compute_vertically_integrated_lateral_areas!,
             compute_matrix_coefficients,
             flux_div_xyᶜᶜᶜ

import Oceananigans.Models.HydrostaticFreeSurfaceModels:
             build_implicit_step_solver,
             compute_implicit_free_surface_right_hand_side!

import Oceananigans.Solvers: solve!

import Oceananigans.Architectures: architecture

struct UnifiedImplicitFreeSurfaceSolver{S, R, T}
    unified_pcg_solver :: S
    right_hand_side :: R
    storage :: T
end

architecture(solver::UnifiedImplicitFreeSurfaceSolver) =
    architecture(solver.preconditioned_conjugate_gradient_solver)

function UnifiedImplicitFreeSurfaceSolver(mrg::MultiRegionGrid, settings, gravitational_acceleration::Number; multiple_devices = false)

    # Initialize vertically integrated lateral face areas
    grid = reconstruct_global_grid(mrg)
    grid = on_architecture(CPU(), grid)

    ∫ᶻ_Axᶠᶜᶜ = Field((Face, Center, Nothing), grid)
    ∫ᶻ_Ayᶜᶠᶜ = Field((Center, Face, Nothing), grid)

    vertically_integrated_lateral_areas = (xᶠᶜᶜ = ∫ᶻ_Axᶠᶜᶜ, yᶜᶠᶜ = ∫ᶻ_Ayᶜᶠᶜ)

    compute_vertically_integrated_lateral_areas!(vertically_integrated_lateral_areas)
    fill_halo_regions!(vertically_integrated_lateral_areas)

    arch = architecture(mrg)
    right_hand_side =  unified_array(arch, zeros(eltype(grid), grid.Nx*grid.Ny))
    storage = deepcopy(right_hand_side)

    # Set maximum iterations to Nx * Ny if not set
    settings = Dict{Symbol, Any}(settings)
    maximum_iterations = get(settings, :maximum_iterations, grid.Nx * grid.Ny)
    settings[:maximum_iterations] = maximum_iterations

    coeffs = compute_matrix_coefficients(vertically_integrated_lateral_areas, grid, gravitational_acceleration)

    reduced_dim = (false, false, true)
    solver = multiple_devices ? UnifiedDiagonalIterativeSolver(coeffs; reduced_dim, grid, mrg, settings...) :
                                HeptadiagonalIterativeSolver(coeffs; reduced_dim,
                                                             template = right_hand_side,
                                                             grid = on_architecture(arch, grid),
                                                             settings...)

    return UnifiedImplicitFreeSurfaceSolver(solver, right_hand_side, storage)
end

build_implicit_step_solver(::Val{:HeptadiagonalIterativeSolver}, grid::MultiRegionGrid, settings, gravitational_acceleration) =
    UnifiedImplicitFreeSurfaceSolver(grid, settings, gravitational_acceleration)
build_implicit_step_solver(::Val{:Default}, grid::MultiRegionGrid, settings, gravitational_acceleration) =
    UnifiedImplicitFreeSurfaceSolver(grid, settings, gravitational_acceleration)
build_implicit_step_solver(::Val{:PreconditionedConjugateGradient}, grid::MultiRegionGrid, settings, gravitational_acceleration) =
    throw(ArgumentError("Cannot use PCG solver with Multi-region grids!! Select :Default or :HeptadiagonalIterativeSolver as solver_method"))

function compute_implicit_free_surface_right_hand_side!(rhs, implicit_solver::UnifiedImplicitFreeSurfaceSolver, g, Δt, ∫ᶻQ, η)
    grid = ∫ᶻQ.u.grid
    M    = length(grid.partition)
    @apply_regionally compute_regional_rhs!(rhs, grid, g, Δt, ∫ᶻQ, η, Iterate(1:M), grid.partition)
    return nothing
end

function compute_regional_rhs!(rhs, grid, g, Δt, ∫ᶻQ, η, region, partition)
    arch = architecture(grid)
    event = launch!(arch, grid, :xy,
                    implicit_linearized_unified_free_surface_right_hand_side!,
                    rhs, grid, g, Δt, ∫ᶻQ, η, region, partition,
		            dependencies = device_event(arch))

    wait(device(arch), event)
    return nothing
end

# linearized right hand side
@kernel function implicit_linearized_unified_free_surface_right_hand_side!(rhs, grid, g, Δt, ∫ᶻQ, η, region, partition)
    i, j = @index(Global, NTuple)
    Az   = Azᶜᶜᶜ(i, j, 1, grid)
    δ_Q  = flux_div_xyᶜᶜᶜ(i, j, 1, grid, ∫ᶻQ.u, ∫ᶻQ.v)
    t    = displaced_xy_index(i, j, grid, region, partition)
    @inbounds rhs[t] = (δ_Q - Az * η[i, j, 1] / Δt) / (g * Δt)
end

function solve!(η, implicit_free_surface_solver::UnifiedImplicitFreeSurfaceSolver, rhs, g, Δt)

    solver = implicit_free_surface_solver.unified_pcg_solver
    storage = implicit_free_surface_solver.storage

    sync_all_devices!(η.grid.devices)

    switch_device!(getdevice(solver.matrix_constructors[1]))
    solve!(storage, solver, rhs, Δt)

    arch = architecture(solver)
    grid = η.grid

    @apply_regionally redistribute_lhs!(η, storage, arch, grid, Iterate(1:length(grid)), grid.partition)

    fill_halo_regions!(η)

    return nothing
end

function redistribute_lhs!(η, sol, arch, grid, region, partition)

    event = launch!(arch, grid, :xy, _redistribute_lhs!, η, sol, region, grid, partition,
		            dependencies = device_event(arch))

    wait(device(arch), event)
end

# linearized right hand side
@kernel function _redistribute_lhs!(η, sol, region, grid, partition)
    i, j = @index(Global, NTuple)
    t = displaced_xy_index(i, j, grid, region, partition)
    @inbounds η[i, j, 1] = sol[t]
end


text:
module MaxHelpingHandAllBlackholesStrengthTrigger

using ..Ahorn, Maple

@mapdef Trigger "MaxHelpingHand/AllBlackholesStrengthTrigger" AllBlackholesStrengthTrigger(x::Integer, y::Integer, width::Integer=Maple.defaultTriggerWidth, height::Integer=Maple.defaultTriggerHeight; strength::String="Mild")

const placements = Ahorn.PlacementDict(
    "Black Hole Strength (All Black Holes) (max480's Helping Hand)" => Ahorn.EntityPlacement(
        AllBlackholesStrengthTrigger,
        "rectangle"
    )
)

function Ahorn.editingOptions(trigger::AllBlackholesStrengthTrigger)
    return Dict{String, Any}(
        "strength" => Maple.black_hole_trigger_strengths
    )
end

end

max_stars_repo_name: Kitty-Cats/MaxHelpingHand

max_issues_repo_name: Kitty-Cats/MaxHelpingHand

max_forks_repo_name: Kitty-Cats/MaxHelpingHand


text:
module LaplBEM

using LinearAlgebra: dot, norm, cross
using Statistics: mean
using FinEtools
using FinEtools.IntegRuleModule: TriRule, GaussRule
import FinEtools.MatrixUtilityModule: locjac!
using Base.Threads

@inline function _distance(dx, dy, dz)::Float64
    return sqrt(dx^2 + dy^2 + dz^2)
end

@inline function _G(r)
	# 3-D free-field Green's function for the Laplace equation
	# x    is the position vector of the observation point
	# x0    is the position vector of the source (same size as x0)
	#         g = 1/(4*pi)/(norm(x-x0));
	return 0.079577471545948/(r);
end

"""
    Greens(k, xsource, xfield)

Compute.
"""
@inline function Greens(xsource, xfield)
    dx, dy, dz = xfield[1]-xsource[1], xfield[2]-xsource[2], xfield[3]-xsource[3]
    return _G(_distance(dx, dy, dz))
end

function _Gn(r, nxx0)#(x,x0,n)
	#
	# Normal derivative of the 3-D free-field Green's function for the Laplace equation
	#
	# x    is the position vector of the observation point (field point)
	# x0    is the position vector of the source (same size as x0)
	#         r=(x-x0);
	#         g = 1/(4*pi)*((n*r')/(norm(r)^3));
	return 0.079577471545948*((nxx0)/(r^3));
end

"""
    dGreensdn(xsource, xfield, n)

We are assuming the normal to point INTO THE FLUID!
"""
@inline function dGreensdn(xsource, xfield, n) # Normal derivative
    dx, dy, dz = xfield[1]-xsource[1], xfield[2]-xsource[2], xfield[3]-xsource[3]
    r = _distance(dx, dy, dz);
    return  _Gn(r, (n[1]*dx + n[2]*dy + n[3]*dz))
end

"""
    trinml!(n, x)

Compute triangle normal (NOT a unit normal, it has a length of twice the area
of the triangle).
"""
function trinml!(n, x)
    a1, a2, a3 = x[2, 1] - x[1, 1], x[2, 2] - x[1, 2], x[2, 3] - x[1, 3]
    b1, b2, b3 = x[3, 1] - x[1, 1], x[3, 2] - x[1, 2], x[3, 3] - x[1, 3]
    n[1] = a2*b3-a3*b2
    n[2] = a3*b1-a1*b3
    n[3] = a1*b2-a2*b1
    return n
end

"""
    triunitnml!(n, x)

Compute the UNIT triangle normal.
"""
function triunitnml!(n, x)
    n = trinml!(n, x)
    n .= n ./ norm(n)
    return n
end

"""
    tricntr!(c, x)

Compute triangle centroid.
"""
function tricntr!(c, x)
    d = length(c)
    for i in 1:d
        c[i] = 0.0
        for j in 1:size(x, 1)
            c[i] += x[j, i]
        end
        c[i] *= (1.0/d)
    end
    return c
end

"""
    centroidareanormal!(c, n, x)

Compute triangle centroid and normal and its area.
"""
function centroidareanormal!(c, n, x)
    n = trinml!(n, x)
    twicearea = norm(n)
    n *= (1.0 / twicearea) # Unit normal
    c = tricntr!(c, x)
    return c, twicearea/2.0, n
end

"""
    surfnml!(n, J)

Compute the surface normal.
A lot of effort has been put into eliminating all allocations, because this function is called in the innermost loop.
"""
function surfnml!(n, J)
   a1, a2, a3 = J[1, 1], J[2, 1], J[3, 1]
   b1, b2, b3 = J[1, 2], J[2, 2], J[3, 2]
   n[1] = a2*b3-a3*b2
   n[2] = a3*b1-a1*b3
   n[3] = a1*b2-a2*b1
   return n
end

function _bfundata(fes, ir)
    Ns = FFltMat[];
    gradNparams = FFltMat[];
    for j in 1:ir.npts
        push!(Ns, bfun(fes,vec(ir.param_coords[j,:])));
        push!(gradNparams, bfundpar(fes,vec(ir.param_coords[j,:])));
    end
    return Ns, gradNparams
end

"""
    _eint(f::T, econn, ex, xsource, n, p, w, Ns, gradNparams, xfield, J, nfield)::Complex{Float64} where {T}

Compute the element integral.
We go to a lot of trouble here to eliminate all sources of allocations. This function is called
for each entry of the Helmholtz matrices.
"""
function _eint(f::T, econn, ex, xsource, n, p, w, Ns, gradNparams, xfield, J, nfield)::Complex{Float64} where {T}
    result = 0.0  + 0.0im
    for j in 1:length(w)
        locjac!(xfield, J, ex, Ns[j], gradNparams[j])
        Jac = norm(surfnml!(nfield, J))
        result += f(xsource, xfield, n) * (w[j] * Jac)
    end
    return result::Complex{Float64}
end

function _allelxs(xyz, conn)
    elxs = fill(Matrix{Float64}(undef, size(conn, 2), 3), size(conn, 1))
    for i in 1:size(conn, 1)
        elxs[i] = xyz[view(conn, i, :), :]
    end
    return elxs
end

"""
    allelxs(xyz, conn)

Construct array of all surface panel coordinates.
"""
function allelxs(xyz, conn)
    return _allelxs(xyz, conn)
end

function _allnormals(elxs)
    normals = fill(Vector{Float64}(undef, 3), length(elxs))
    cnormal = fill(0.0, 3)
    for i in 1:length(elxs)
        normals[i] = deepcopy(triunitnml!(cnormal, elxs[i]))
    end
    return normals
end

"""
    allnormals(elxs)

Construct array of all surface panel normals.
"""
function allnormals(elxs)
    return _allnormals(elxs)
end

function _doublelayer!(A, firstrowofset, rxyz, xyz, conn, elxs, solidangle, ir::TriRule)
	nrows = size(rxyz, 1)
    neqn = size(conn, 1)
    @assert size(A, 1) == nrows "Size of matrix needs to match number of row points"
    @assert size(A, 2) == neqn "Size of matrix needs to match number of surface panels"
    econn = (1,2,3)
    dummyfes = FESetT3(reshape(collect(econn), 1, 3))
    p, w = ir.param_coords, vec(ir.weights)
    Ns, gradNparams = _bfundata(dummyfes, ir)
    f = (xsource, xfield, n) -> (-1.0) * dGreensdn(xsource, xfield, n)
    cnormal = fill(0.0, 3)
    xfield = fill(zero(FFlt), 1, 3);
    J = fill(zero(FFlt), 3, 2);
    nfield = fill(zero(FFlt), 1, 3);
    normals = _allnormals(elxs)
    for r in 1:nrows
        rx = view(rxyz, r, :)
        for c in 1:neqn
            if r+firstrowofset != c
                A[r, c] += _eint(f, econn, elxs[c], rx, normals[c], p, w, Ns, gradNparams, xfield, J, nfield)
            else
                A[r, c] += solidangle
            end
        end
    end
    return A
end

function _singlelayer!(B, firstrowofset, rxyz, xyz, conn, elxs, iroffdiagonal::TriRule, irdiagonal::TriRule)
	nrows = size(rxyz, 1)
    neqn = size(conn, 1)
    @assert size(B, 1) == nrows "Size of matrix needs to match number of row points"
    @assert size(B, 2) == neqn "Size of matrix needs to match number of surface panels"
    econn = (1,2,3)
    dummyfes = FESetT3(reshape(collect(econn), 1, 3))
    poff, woff = iroffdiagonal.param_coords, vec(iroffdiagonal.weights)
    Nsoff, gradNparamsoff = _bfundata(dummyfes, iroffdiagonal)
    pdiagonal, wdiagonal = irdiagonal.param_coords, vec(irdiagonal.weights)
    Nsdiagonal, gradNparamsdiagonal = _bfundata(dummyfes, irdiagonal)
    f = (xsource, xfield, n) -> (-1.0) * Greens(xsource, xfield)
    ex = fill(0.0, 3, 3)
    xfield = fill(zero(FFlt), 1, 3);
    J = fill(zero(FFlt), 3, 2);
    nfield = fill(zero(FFlt), 1, 3);
    normals = _allnormals(elxs)
    ccentroid = fill(0.0, 3)
    for r in 1:nrows
        rx = view(rxyz, r, :)
        for c in 1:neqn
            if r+firstrowofset != c
                B[r, c] += _eint(f, econn, elxs[c], rx, normals[c], poff, woff, Nsoff, gradNparamsoff, xfield, J, nfield)
            else
                ex[3,:] .= tricntr!(ccentroid, elxs[c])
                ex[1,:] .= elxs[c][1,:]; ex[2,:] .= elxs[c][2,:];
                B[r, c] += _eint(f, econn, ex, rx, normals[c], pdiagonal, wdiagonal, Nsdiagonal, gradNparamsdiagonal, xfield, J, nfield)
                ex[1,:] .= elxs[c][2,:]; ex[2,:] .= elxs[c][3,:];
                B[r, c] += _eint(f, econn, ex, rx, normals[c], pdiagonal, wdiagonal, Nsdiagonal, gradNparamsdiagonal, xfield, J, nfield)
                ex[1,:] .= elxs[c][3,:]; ex[2,:] .= elxs[c][1,:];
                B[r, c] += _eint(f, econn, ex, rx, normals[c], pdiagonal, wdiagonal, Nsdiagonal, gradNparamsdiagonal, xfield, J, nfield)
            end
        end
    end
    return B
end

"""
    doublelayer!(A, xyz, conn, ir::TriRule)

Compute the double-layer matrix contribution.

Arguments
- `A` = matrix
- `xyz` = array of coordinates of the nodes
- `conn` = array of connectivities of the surface elements
- `ir` = integration rule

The contribution is _added_to the matrix `A`. This matrix is also returned for convenience.
"""
function doublelayer!(A, xyz, conn, ir::TriRule)
	neqn = size(conn, 1)
    rxyz = fill(0.0, neqn, 3)
    rcentroid = fill(0.0, 3)
    elxs = _allelxs(xyz, conn)
    for r in 1:neqn
        rxyz[r, :] .= tricntr!(rcentroid, elxs[r])
    end
    return _doublelayer!(A, 0, rxyz, xyz, conn, elxs, -0.5, ir)
end

"""
    singlelayer!(B, k, xyz, conn, iroffdiagonal::TriRule, irdiagonal::TriRule)

Compute the single-layer matrix contribution.

Arguments
- `B` = complex matrix
- `xyz` = array of coordinates of the nodes
- `conn` = array of connectivities of the surface elements
- `iroffdiagonal` = integration rule for the off-diagonal terms
- `irdiagonal` = integration rule for the diagonal terms

The contribution is _added_to the matrix `B`. This matrix is also returned for convenience.
"""
function singlelayer!(B, xyz, conn, iroffdiagonal::TriRule, irdiagonal::TriRule)
	neqn = size(conn, 1)
    rxyz = fill(0.0, neqn, 3)
    rcentroid = fill(0.0, 3)
    elxs = _allelxs(xyz, conn)
    for r in 1:neqn
        rxyz[r, :] .= tricntr!(rcentroid, elxs[r])
    end
    return _singlelayer!(B, 0, rxyz, xyz, conn, elxs, iroffdiagonal, irdiagonal)
end

# function _doublelayercolumns!(A, firstrowofset, columnrange, k, rxyz, xyz, conn, elxs, normals, solidangle, ir::TriRule)
# 	nrows = size(rxyz, 1)
#     neqn = size(conn, 1)
#     econn = (1,2,3)
#     dummyfes = FESetT3(reshape(collect(econn), 1, 3))
#     p, w = ir.param_coords, vec(ir.weights)
#     Ns, gradNparams = _bfundata(dummyfes, ir)
#     f = (xsource, xfield, n) -> dGreensdn(k, xsource, xfield, n)
#     cnormal = fill(0.0, 3)
#     xfield = fill(zero(FFlt), 1, 3);
#     J = fill(zero(FFlt), 3, 2);
#     nfield = fill(zero(FFlt), 1, 3);
#     for r = 1:nrows
#     	rx = view(rxyz, r, :)
#         for c in columnrange
#             if r+firstrowofset != c
#                 A[r, c] += (-1.0) * _eint(f, econn, elxs[c], rx, normals[c], p, w, Ns, gradNparams, xfield, J, nfield)
#             else
#                 A[r, c] += solidangle
#             end
#         end
#     end
#     return A
# end

"""
    tdoublelayer!(A, k, xyz, conn, ir::TriRule)

Compute the double-layer matrix contribution.  Threaded version.

- `A` = complex matrix
- `k` = wave number
- `xyz` = array of coordinates of the nodes
- `conn` = array of connectivities of the surface elements
- `ir` = integration rule

As `doublelayer(A, k, xyz, conn, ir::TriRule)`, but the calculation is carried
out in parallel using all available workers.
"""
# function tdoublelayer!(A, k, xyz, conn, ir::TriRule)
#     neqn = size(conn, 1)
#     rxyz = fill(0.0, neqn, 3)
#     rcentroid = fill(0.0, 3)
#     elxs = _allelxs(xyz, conn)
#     for r = 1:neqn
#     	rxyz[r, :] .= tricntr!(rcentroid, elxs[r])
#     end
#     normals = _allnormals(elxs)
#     firstrowofset = 0
#     nrows = size(rxyz, 1)
#     nth = Base.Threads.nthreads()
#     _tcolumnchunk = Int(ceil(neqn/nth))
#     c = 0
#     while c < neqn
#     	tasks = []
#     	for t in 1:nth
#     		columnrange = c+1:min(neqn, c+1+_tcolumnchunk)
#     		push!(tasks, Threads.@spawn _doublelayercolumns!(A, firstrowofset, columnrange, k, rxyz, xyz, conn, elxs, normals, 0.5, ir))
#     		c = maximum(columnrange)
#     		if c >= neqn
#     			break
#     		end
#     	end
#     	Threads.wait.(tasks);
#     end
#     return A
# end

"""
    tdoublelayer!(A, k, xyz, conn, ir::TriRule)

Compute the double-layer matrix contribution.  Threaded version.

As `doublelayer(A, k, xyz, conn, ir::TriRule)`, but the calculation is carried
out in parallel using all available threads.

Don't forget to set the number of threads with
`export JULIA_NUM_THREADS=4`.
"""
# function tdoublelayer!(A, k, chiefxyz, xyz, conn, ir::TriRule)
# 	neqn = size(conn, 1)
# 	elxs = _allelxs(xyz, conn)
# 	normals = _allnormals(elxs)
# 	firstrowofset = neqn
# 	nrows = size(chiefxyz, 1)
# 	nth = Base.Threads.nthreads()
# 	_tcolumnchunk = Int(ceil(neqn/nth))
# 	c = 0
# 	while c < neqn
# 		tasks = []
# 		for t in 1:nth
# 			columnrange = c+1:min(neqn, c+1+_tcolumnchunk)
# 			push!(tasks, Threads.@spawn _doublelayercolumns!(A, firstrowofset, columnrange, k, chiefxyz, xyz, conn, elxs, normals, 0.0, ir))
# 			c = maximum(columnrange)
# 			if c >= neqn
# 				break
# 			end
# 		end
# 		Threads.wait.(tasks);
# 	end
# 	return A
# end

# function _singlelayercolumns!(B, firstrowofset, columnrange, k, rxyz, xyz, conn, elxs, normals, iroffdiagonal::TriRule, irdiagonal::TriRule)
#     nrows = size(rxyz, 1)
#     neqn = size(conn, 1)
#     econn = (1,2,3)
#     dummyfes = FESetT3(reshape(collect(econn), 1, 3))
#     poff, woff = iroffdiagonal.param_coords, vec(iroffdiagonal.weights)
#     Nsoff, gradNparamsoff = _bfundata(dummyfes, iroffdiagonal)
#     pdiagonal, wdiagonal = irdiagonal.param_coords, vec(irdiagonal.weights)
#     Nsdiagonal, gradNparamsdiagonal = _bfundata(dummyfes, irdiagonal)
#     f = (xsource, xfield, n) -> Greens(k, xsource, xfield)
#     ex = fill(0.0, 3, 3)
#     xfield = fill(zero(FFlt), 1, 3);
#     J = fill(zero(FFlt), 3, 2);
#     nfield = fill(zero(FFlt), 1, 3);
#     ccentroid = fill(0.0, 3)
#     for r = 1:nrows
#         rx = view(rxyz, r, :)
#         for c = columnrange
#             if r+firstrowofset != c
#                 B[r, c] += _eint(f, econn, elxs[c], rx, normals[c], poff, woff, Nsoff, gradNparamsoff, xfield, J, nfield)
#             else
#                 ex[3,:] .= tricntr!(ccentroid, elxs[c])
#                 ex[1,:] .= elxs[c][1,:]; ex[2,:] .= elxs[c][2,:];
#                 B[r, c] += _eint(f, econn, ex, rx, normals[c], pdiagonal, wdiagonal, Nsdiagonal, gradNparamsdiagonal, xfield, J, nfield)
#                 ex[1,:] .= elxs[c][2,:]; ex[2,:] .= elxs[c][3,:];
#                 B[r, c] += _eint(f, econn, ex, rx, normals[c], pdiagonal, wdiagonal, Nsdiagonal, gradNparamsdiagonal, xfield, J, nfield)
#                 ex[1,:] .= elxs[c][3,:]; ex[2,:] .= elxs[c][1,:];
#                 B[r, c] += _eint(f, econn, ex, rx, normals[c], pdiagonal, wdiagonal, Nsdiagonal, gradNparamsdiagonal, xfield, J, nfield)
#             end
#         end
#     end
#     return B
# end

"""
    tsinglelayer!(B, k, xyz, conn, iroffdiagonal::TriRule, irdiagonal::TriRule)

Compute the single-layer matrix contribution.  Threaded version.

As `singlelayer!(B, k, xyz, conn, iroffdiagonal::TriRule,
irdiagonal::TriRule)`, but the calculation is carried out in parallel using
all available threads.
"""
# function tsinglelayer!(B, k, xyz, conn, iroffdiagonal::TriRule, irdiagonal::TriRule)
# 	neqn = size(conn, 1)
# 	elxs = _allelxs(xyz, conn)
# 	normals = _allnormals(elxs)
# 	rxyz = fill(0.0, neqn, 3)
# 	rcentroid = fill(0.0, 3)
# 	for r = 1:neqn
# 		rxyz[r, :] .= tricntr!(rcentroid, elxs[r])
# 	end
# 	nrows = size(rxyz, 1)
# 	firstrowofset = 0
# 	nth = Base.Threads.nthreads()
# 	_tcolumnchunk = Int(ceil(neqn/nth))
# 	c = 0
# 	while c < neqn
# 		tasks = []
# 		for t in 1:nth
# 			columnrange = c+1:min(neqn, c+1+_tcolumnchunk)
# 			push!(tasks, Threads.@spawn _singlelayercolumns!(B, firstrowofset, columnrange, k, rxyz, xyz, conn, elxs, normals, iroffdiagonal, irdiagonal))
# 			c = maximum(columnrange)
# 			if c >= neqn
# 				break
# 			end
# 		end
# 		Threads.wait.(tasks);
# 	end
# 	return B
# end

"""
    tsinglelayer!(B, k, xyz, conn, iroffdiagonal::TriRule, irdiagonal::TriRule)

Compute the single-layer matrix contribution. Threaded version.

As `singlelayer!(B, k, xyz, conn, iroffdiagonal::TriRule,
irdiagonal::TriRule)`, but the calculation is carried out in parallel using
all available threads.

Don't forget to set the number of threads with
`export JULIA_NUM_THREADS=4`.
"""
# function tsinglelayer!(B, k, chiefxyz, xyz, conn, iroffdiagonal::TriRule, irdiagonal::TriRule)
#     neqn = size(conn, 1)
#     elxs = _allelxs(xyz, conn)
#     normals = _allnormals(elxs)
#     nrows = size(chiefxyz, 1)
#     firstrowofset = neqn
#     nth = Base.Threads.nthreads()
#     _tcolumnchunk = Int(ceil(neqn/nth))
#     c = 0
#     while c < neqn
#     	tasks = []
#     	for t in 1:nth
#     		columnrange = c+1:min(neqn, c+1+_tcolumnchunk)
#     		push!(tasks, Threads.@spawn _singlelayercolumns!(B, firstrowofset, columnrange, k, chiefxyz, xyz, conn, elxs, normals, iroffdiagonal, irdiagonal))
#     		c = maximum(columnrange)
#     		if c >= neqn
#     			break
#     		end
#     	end
#     	Threads.wait.(tasks);
#     end
#     return B
# end

end

max_stars_repo_name: PetrKryslUCSD/FinEtoolsVibInFluids

max_issues_repo_name: PetrKryslUCSD/FinEtoolsVibInFluids

max_forks_repo_name: PetrKryslUCSD/FinEtoolsVibInFluids


text:
module HawkesSimulator
using StatsBase,Statistics,Distributions,LinearAlgebra,Random
using FFTW
using Colors # to save rasters as png

#=
UnitType : specifies the kernel for Hawkes processes, or
the input type for input units.
PopulationState  : contains UnitType, n, which is the population size,
then trains contains n arrays, each is the spike of past trains
for that single unit.  trains_history are non-interactive spikes
AbstractPopulation : contains one postsynaptic (or input) population state,
which is called `state`.Contains `spike_proposals` for memory allocation
PopulationHawkes : contains one postsynaptic population, plus the connecting from
other populations into it.
PopulationInput : used for input types.

=#

abstract type UnitType end
abstract type Connection end
abstract type PlasticityRule end
struct NoPlasticity <: PlasticityRule end

struct ConnectionWeights{N,PL<:NTuple{N,PlasticityRule}} <: Connection
  weights::Matrix{Float64}
  plasticities::PL
end
function ConnectionWeights(weights::Matrix{Float64})
  plast = NTuple{0,NoPlasticity}()
  return ConnectionWeights(weights,plast)
end
function ConnectionWeights(weights::Matrix{Float64},
     (plasticity_rules::PL where PL<:PlasticityRule)...)
  return ConnectionWeights(weights,plasticity_rules)
end
function reset!(conn::ConnectionWeights)
  reset!.(conn.plasticities)
  return nothing
end

abstract type AbstractNonlinearity end
# this one is the default
struct NLRelu <: AbstractNonlinearity end

# for now, inputs are stationary and constant
# this can be generalized later
# abstract type ExternalInput end

abstract type AbstractPopulationState end
# includes the special case of ExpKernel, which has no unit type
# must still have n and label.

struct PopulationState{UT <:UnitType} <: AbstractPopulationState
  label::Symbol
  n::Int64
  unittype::UT
  trains::Vector{Vector{Float64}}
  trains_history::Vector{Vector{Float64}}
end
nneurons(ps::PopulationState) = ps.n

# constructor

function rand_label()
  return Symbol(randstring(3))
end
function PopulationState(unit_type::UnitType,n::Int64;
    label::Union{String,Nothing}=nothing)
  label = something(label,rand_label())
  trains = [Float64[] for _ in 1:n] # empty vectors
  trainsh = [Float64[] for _ in 1:n]
  PopulationState(label,n,unit_type,trains,trainsh)
end

function reset!(ps::PopulationState)
  for train in ps.trains
    if !isempty(train)
      deleteat!(train,1:length(train))
    end
  end
  for trainh in ps.trains_history
    if !isempty(trainh)
      deleteat!(trainh,1:length(trainh))
    end
  end
  return nothing
end

# this will cover Hawkes populations but also input populations
abstract type AbstractPopulation end
nneurons(p::AbstractPopulation) = nneurons(p.state)

struct PopulationHawkes{N,PS<:PopulationState,
    TC<:NTuple{N,Connection},
    TP<:NTuple{N,PopulationState},
    NL<:AbstractNonlinearity} <: AbstractPopulation
  state::PS
  connections::TC
  pre_states::TP
  input::Vector{Float64} # input only used for external currents
  nonlinearity::NL # nonlinearity if present
  spike_proposals::Vector{Float64} # memory allocation
end

function PopulationHawkes(state::PopulationState,input::Vector{Float64},
    (conn_pre::Tuple{C,PS} where {C<:Connection,PS<:PopulationState})...;
      nonlinearity::AbstractNonlinearity=NLRelu())
  connections = Tuple(getindex.(conn_pre,1))
  pre_states = Tuple(getindex.(conn_pre,2))
  spike_proposals = fill(-Inf,nneurons(state))
  return PopulationHawkes(state,connections,pre_states,input,nonlinearity,spike_proposals)
end

# one population only!
function PopulationHawkes(state::PopulationState,conn::Connection,input::Vector{Float64};
    nonlinearity::AbstractNonlinearity=NLRelu())
  return PopulationHawkes(state,input,(conn,state);nonlinearity=nonlinearity)
end

# recurrent network is just a tuple of (input) populations
struct RecurrentNetwork{N,TP<:NTuple{N,AbstractPopulation}}
  populations::TP
end
# base constructor
function RecurrentNetwork((pops::P where P<:AbstractPopulation)...)
  RecurrentNetwork(pops)
end
# reset
function reset!(rn::RecurrentNetwork)
  for pop in rn.populations
    reset!(pop.state)
    # reset plasticity operators if present
    if isdefined(pop,:connections) # because input has no connection field
      for conn in pop.connections
        for plast in conn.plasticities
          reset!(plast)
        end
      end
    end
  end
  return nothing
end
@inline function npopulations(ntw::RecurrentNetwork)
  return length(ntw.populations)
end

# more constructors
# one population simplification (assuming it is an Hawkes type!)
function RecurrentNetwork(state::PopulationState,conn::Connection,
    input::Vector{Float64};
    nonlinearity::AbstractNonlinearity=NLRelu())
  return RecurrentNetwork(PopulationHawkes(state,conn,input;nonlinearity=nonlinearity))
end
# or even simpler constructor
function RecurrentNetwork(state::PopulationState,weights::Matrix{Float64},
    input::Vector{Float64};
    nonlinearity=NLRelu())
  @assert size(weights,2) == length(input)
  return RecurrentNetwork(
    PopulationHawkes(state,ConnectionWeights(weights),input;
    nonlinearity=nonlinearity))
end

###################
# traces as structs
# (after all, why not?  ... why shouldn't I use a struct?)

# trace for dynamics : all receive an update
# as soon as the spike is produced. The update is normed
# so that the integral of kernel is one

# trace for plasticity : ignored by the dynamcs,
# plasticity rules should take care of it.

abstract type TracePurpose end
struct ForDynamics <: TracePurpose end
struct ForPlasticity <: TracePurpose end

struct Trace{P<:TracePurpose}
  val::Vector{Float64}
  τ::Float64
  t_last::Ref{Float64}
  purpose::P
  function Trace(τ::Real,n::Integer,
      (purpose::P)=ForPlasticity()) where {P<:TracePurpose}
    val = fill(0.0,n)
    t_last = Ref(0.0)
    return new{P}(val,τ,t_last,purpose)
  end
end

function nneurons(tra::Trace)
  return length(tra.val)
end

function propagate!(tnow::Float64,tra::Trace)
  Δt = tnow - tra.t_last[]
  tra.val .*= exp(-Δt/tra.τ)
  tra.t_last[]=tnow
  return nothing
end
# updates single element of trace, unless it is nothing
function update_now!(tra::Trace,idx_update::Int64,up_val::Float64=1.0)
  if idx_update != 0
    tra.val[idx_update] += up_val
  end
  return nothing
end

# updates only if the trace is specifically for dynamics
function update_for_dynamics!(::Trace,::Int64)
  return nothing
end
function update_for_dynamics!(tra::Trace{ForDynamics},idx_update::Int64)
  tra.val[idx_update] += inv(tra.τ)
  return nothing
end
# same with propagate
function propagate_for_dynamics!(::Float64,::Trace)
  return nothing
end
function propagate_for_dynamics!(tnow::Float64,tra::Trace{ForDynamics})
  Δt = tnow - tra.t_last[]
  tra.val .*= exp(-Δt/tra.τ)
  tra.t_last[]=tnow
  return nothing
end

function reset!(tra::Trace)
  fill!(tra.val,0.0)
  tra.t_last[]=0.0
  return nothing
end

# proposal of future trace. Useful to compute quantities without advancing the trace
function trace_proposal!(proposal::Vector{R},tnow::R,tra::Trace) where R
  Δt = tnow - tra.t_last[]
  @. proposal = tra.val * exp(-Δt/tra.τ)
  return proposal
end

#### Small general utility functions
@inline hardbounds(x::R,low::R,high::R) where R = min(high,max(x,low))

#############
# everything else

# general Hawkes for any type of kernel
include("spike_generation_hawkes.jl")
# the interaction_kernel and interaction_kernel_upper come from here :
include("kernels.jl")
# special case optimized for exponential kernels
include("spike_generation_expkernel.jl")
# apply_nonlinearity defined here
include("nonlinearities.jl")
# input-neurons
include("inputs.jl")
# apply_plasticities defined here
include("plasticity_rules.jl")
include("plasticity_heterosynaptic.jl")
# analyze the output: mean, covariance, etc
include("spike_analysis.jl")

end # of module


text:
@testset "Calculation Tools" begin
    @testset "Test radial_mean" begin
        a = rr((10,10)) .< 5.0
        (rm, ax) = MicroscopyTools.radial_mean(a)
        @test rm[1] ≈ 1.0
        @test rm[5] ≈ 1.0
        @test rm[6] ≈ 0.0
        @test ax[1] == 0.5
        @test ax[end] == 7.5
    end

    @testset "Δ_phase" begin
        @test MicroscopyTools.Δ_phase([1,1im,-1,-1im,1], 1) ≈ -pi/2 .* ones(4)
    end

    @testset "moment projection" begin
        @test MicroscopyTools.moment_proj([1 3; 2 4], 2, pdims = (1,)) == [0.25 0.25]
        @test MicroscopyTools.moment_proj_normed([1 3; 2 4], 2, pdims = (1,)) == [0.5 0.5]
    end
end

max_stars_repo_name: JuliaMicroscopy/MicroscopyTools

max_issues_repo_name: JuliaMicroscopy/MicroscopyTools


text:
# node_basics.jl
# ==============
#
# Types and methods for phylogenetic trees.
#
# Part of the Bio.Phylo module.
#
# This file contains the methods for updating or querying characteristics of
# nodes/clades in a phylogeny. In the Phylo representation these are
# vertices of a LightGraphs DiGraph, which are numbered higher than the
# number of taxa + 1.
#
# This file is a part of BioJulia. License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

"""
    isroot(tree::Phylogeny, vertex::Int)

Test whether a vertex is the root.

# Examples:
```julia
isroot(my_tree, 9)
```
"""
function isroot(tree::Phylogeny, vertex::Int)
    return vertex == root(tree)
end

"""
    isroot{T <: SymString}(tree::Phylogeny, vertex::T)

Test whether a vertex is the root.

# Examples:
```julia
isroot(my_tree, :Root)
isroot(my_tree, :A)
```
"""
function isroot(tree::Phylogeny, name::SymString)
    isroot(tree, tree.vertexindex[name])
end

"""
    isclade(tree::Phylogeny, vertex::Int)

Test whether a node is internal clade, i.e. has a parent and one or more children.

# Examples:
```julia
isclade(my_tree, 6)
```
"""
function isclade(tree::Phylogeny, vertex::Int)
    return vertex >= root(tree)
end

"""
    isclade(tree::Phylogeny, vertex::SymString)

Test whether a node is internal, i.e. has a parent and one or more children.

# Examples:
```julia
isclade(my_tree, :Node7)
isclade(my_tree, "Homo Sapiens")
isclade(my_tree, Symbol("Homo Sapiens"))
```
"""
function isclade(tree::Phylogeny, name::SymString)
    return isclade(tree, tree.vertexindex[name])
end

"""
    isleaf(tree::Phylogeny, vertex::Int)

Test whether a vertex is a leaf.

A vertex is a leaf, when it has a parent node, but no children nodes.

# Examples:
```julia
isleaf(my_phylogeny, 1)
```
"""
function isleaf(tree::Phylogeny, vertex::Int)
    return vertex <= tree.ntaxa
end

"""
    isleaf(tree::Phylogeny, name::SymString)

Test whether a vertex is a leaf.

A vertex is a leaf, when it has a parent node, but no children nodes.

# Examples:
```julia
isleaf(my_phylogeny, :A)
isleaf(my_phylogeny, "A")
```
"""
function isleaf(tree::Phylogeny, name::SymString)
    return isleaf(tree, tree.vertexindex[name])
end

"""
    hasparent(tree::Phylogeny, vertex::Int)

Check if a vertex in a phylogeny is linked to a parent.
"""
function hasparent(tree::Phylogeny, vertex::Int)
    return indegree(tree.graph, vertex) >= 1
end

"""
    hasparent(tree::Phylogeny, name::SymString)

Check if the vertex in a phylogeny is linked to a parent.
"""
function hasparent(tree::Phylogeny, name::SymString)
    return hasparent(tree, tree.vertexindex[name])
end

"""
    hasparent(tree::Phylogeny, vertices::AbstractVector{Int})

Check some vertices in a phylogeny and test if each each is linked to a parent.
"""
function hasparent(tree::Phylogeny, vertices::AbstractVector{Int})
    return indegree(tree.graph, vertices) .>= 1
end

"""
    parent(x::PhyNode, vertex::Int)

Get the parent of a node. In a phylogeny it is not possible
for any one node/vertex to have more than one parent.

# Examples:
```julia
parent(my_tree, 2)
```
"""
function parent(tree::Phylogeny, vertex::Int)
    @assert hasparent(tree, vertex) error("Vertex does not have a parent.")
    return in_neighbors(tree.graph, vertex)[1]
end

"""
    hasparent(tree::Phylogeny, child::Int, parent::Int)

Test if a vertex of a phylogeny is the parent of another vertex.
"""
function hasparent(tree::Phylogeny, child::Int, parent::Int)
    return has_edge(tree.graph, parent, child)
end

"""
    haschildren(tree::Phylogeny, vertex::Int)

Test whether a vertex in a phylogeny has children.
"""
function haschildren(tree::Phylogeny, vertex::Int)
    return nchildren(tree, vertex) > 0
end

"""
    haschildren(tree::Phylogeny, vertices::AbstractVector{Int})

Test vertices in a phylogeny to see if it has children.
"""
function haschildren(tree::Phylogeny, vertices::AbstractVector{Int})
    return nchildren(tree, vertices) .> 0
end

"""
    nchildren(tree::Phylogeny, vertex::Int)

Count the number of clades that are direct children of a clade.
"""
function nchildren(tree::Phylogeny, vertex::Int)
    return outdegree(tree.graph, vertex)
end

"""
    nchildren(tree::Phylogeny, vertices::AbstractVector{Int})

Count the number of clades that are direct children for a selection of clades.
"""
function nchildren(tree::Phylogeny, vertices::AbstractVector{Int})
    return outdegree(tree.graph, vertices)
end

"""
    nchildren(tree::Phylogeny, name::T)

Count the number of clades that are direct children of a clade.
"""
function nchildren(tree::Phylogeny, name::SymString)
    return nchildren(tree, tree.vertexindex[name])
end

"""
    haschild(tree::Phylogeny, child::Int, parent::Int)

Test if a vertex of a phylogeny is the child of another vertex.
"""
function haschild(tree::Phylogeny, parent::Int, child::Int)
    return has_edge(tree.graph, parent, child)
end

"""
    children(tree::Phylogeny, vertex::Int)

Get the children vertices of a vertex in the phylogeny.
"""
function children(tree::Phylogeny, vertex::Int)
    return copy(out_neighbors(tree.graph, vertex))
end

"""
    ispreterminal(tree::Phylogeny)

Test whether a node is preterminal i.e. It's children are all leaves.
"""
function ispreterminal(tree::Phylogeny, vertex::Int)
    if !haschildren(tree, vertex)
        return false
    end
    # We use out_neighbours below as it returns a reference, rather than a copy
    # as the children method above does, so we save unnessecery copy ops.
    return all([isleaf(tree, i) for i in out_neighbors(tree.graph, vertex)])
end

"""
    issemipreterminal(tree::Phylogeny, vertex::Int)

Test whether a node is semi-preterminal i.e. Some of it's children are leaves,
but not all are.
"""
function issemipreterminal(tree::Phylogeny, vertex::Int)
    if !haschildren(tree, vertex)
        return false
    end
    # We use out_neighbours below as it returns a reference, rather than a copy
    # as the children method above does, so we save unnessecery copy ops.
    areleaves = [isleaf(tree, i) for i in out_neighbors(tree.graph, vertex)]
    return any(areleaves) && !all(areleaves)
end

"""
    isredundant(tree::Phylogeny, vertex::Int)

Test if a vertex in a phylogeny is redundant or not.
A vertex is considered redundant is it has a parent and only one child.
"""
function isredundant(tree::Phylogeny, vertex::Int)
    return hasparent(tree, vertex) && nchildren(tree, vertex) == 1
end

"""
    isredundant(tree::Phylogeny, vertices::AbstractArray{Int})

Test if vertices in a phylogeny to see if each is redundant or not.
A vertex is considered redundant is it has a parent and only one child.
"""
function isredundant(tree::Phylogeny, vertices::AbstractArray{Int})
    return hasparent(tree, vertices) & nchildren(tree, vertices) .== 1
end


text:
# Methods to resolve ambiguities with `AbstractTriangular`
@inline *(tvec::TransposedVector, trimat::LinAlg.AbstractTriangular) = transpose(trimat * transpose(tvec))
@inline *(::LinAlg.AbstractTriangular, ::TransposedVector) = error("Cannot right-multiply matrix by transposed vector")

@inline A_mul_Bt(tvec::TransposedVector, mat::LinAlg.AbstractTriangular) = transpose(mat * transpose(tvec))
@inline A_mul_Bt(mat::LinAlg.AbstractTriangular, tvec::TransposedVector) = mat * transpose(tvec)

@inline At_mul_Bt(mat::LinAlg.AbstractTriangular, tvec::TransposedVector) = mat.' * transpose(tvec)
@inline At_mul_Bt(tvec::TransposedVector, mat::LinAlg.AbstractTriangular) = error("Cannot left-multiply matrix by vector")

@inline At_mul_B(mat::LinAlg.AbstractTriangular, tvec::TransposedVector) = error("Cannot right-multiply matrix by transposed vector")
@inline At_mul_B(tvec::TransposedVector, mat::LinAlg.AbstractTriangular) = error("Cannot left-multiply matrix by vector")

@inline A_mul_Bc(tvec::TransposedVector, mat::LinAlg.AbstractTriangular) = ctranspose(mat * ctranspose(tvec))
@inline A_mul_Bc(mat::LinAlg.AbstractTriangular, tvec::TransposedVector) = mat * ctranspose(tvec)

@inline Ac_mul_Bc(mat::LinAlg.AbstractTriangular, tvec::TransposedVector) = mat' * ctranspose(tvec)
@inline Ac_mul_Bc(tvec::TransposedVector, mat::LinAlg.AbstractTriangular) = error("Cannot left-multiply matrix by vector")

@inline Ac_mul_B(mat::LinAlg.AbstractTriangular, tvec::TransposedVector) = error("Cannot right-multiply matrix by transposed vector")
@inline Ac_mul_B(tvec::TransposedVector, mat::LinAlg.AbstractTriangular) = error("Cannot left-multiply matrix by vector")

@inline /(tvec::TransposedVector, tri::Union{UpperTriangular,LowerTriangular}) = transpose(transpose(tri) \ transpose(tvec))
@inline /(tvec::TransposedVector, tri::Union{LinAlg.UnitUpperTriangular,LinAlg.UnitLowerTriangular}) = transpose(transpose(tri) \ transpose(tvec))

@inline A_rdiv_Bt(tvec::TransposedVector, tri::Union{UpperTriangular,LowerTriangular}) = transpose(tri \ transpose(tvec))
@inline A_rdiv_Bt(tvec::TransposedVector, tri::Union{LinAlg.UnitUpperTriangular,LinAlg.UnitLowerTriangular}) = transpose(tri \ transpose(tvec))

@inline A_rdiv_Bc(tvec::TransposedVector, tri::Union{UpperTriangular,LowerTriangular}) = ctranspose(tri \ ctranspose(tvec))
@inline A_rdiv_Bc(tvec::TransposedVector, tri::Union{LinAlg.UnitUpperTriangular,LinAlg.UnitLowerTriangular}) = ctranspose(tri \ ctranspose(tvec))

@inline \(::Union{UpperTriangular,LowerTriangular}, ::TransposedVector) = error("Cannot left-divide matrix by transposed vector")
@inline \(::Union{LinAlg.UnitUpperTriangular,LinAlg.UnitLowerTriangular}, ::TransposedVector) = error("Cannot left-divide matrix by transposed vector")

@inline At_ldiv_B(::Union{UpperTriangular,LowerTriangular}, ::TransposedVector) = error("Cannot left-divide matrix by transposed vector")
@inline At_ldiv_B(::Union{LinAlg.UnitUpperTriangular,LinAlg.UnitLowerTriangular}, ::TransposedVector) = error("Cannot left-divide matrix by transposed vector")

@inline Ac_ldiv_B(::Union{UpperTriangular,LowerTriangular}, ::TransposedVector) = error("Cannot left-divide matrix by transposed vector")
@inline Ac_ldiv_B(::Union{LinAlg.UnitUpperTriangular,LinAlg.UnitLowerTriangular}, ::TransposedVector) = error("Cannot left-divide matrix by transposed vector")


text:
import Base.FloatRange

type CoordInterpGrid{T<:AbstractFloat, N, BC<:BoundaryCondition, IT<:InterpType, R} <: AbstractInterpGrid{T,N}
    coord::R
    grid::InterpGrid{T,N,BC,IT}
    function CoordInterpGrid(coord::NTuple{N,Range}, grid::InterpGrid{T,N,BC,IT})
        map(length,coord) == size(grid) || throw(DimensionMismatch("Coordinate lengths do not match grid size."))
        new(coord,grid)
    end
end

function CoordInterpGrid{T,N,BC,IT}(coord::NTuple{N,Range}, grid::InterpGrid{T,N,BC,IT})
    CoordInterpGrid{T,N,BC,IT,typeof(coord)}(coord,grid)
end

function CoordInterpGrid{T<:AbstractFloat,R<:Range}(coord::R, grid::InterpGrid{T,1})
    CoordInterpGrid((coord,),grid)
end

function CoordInterpGrid{N,T<:AbstractFloat}(coord::NTuple{N,Range}, A::Array{T,N}, args...)
    CoordInterpGrid(coord,InterpGrid(A,args...))
end
function CoordInterpGrid{R<:Range,T<:AbstractFloat}(coord::R,A::Array{T,1},args...)
    CoordInterpGrid((coord,),InterpGrid(A,args...))
end

size(C::CoordInterpGrid) = size(C.grid)

coordlookup(r::FloatRange,x::Real) = (r.divisor*x-r.start)/r.step + 1.0
coordlookup(r::StepRange,x::Real) = (x-r.start)/r.step + 1.0
coordlookup(r::UnitRange,x::Real) = x-r.start + 1
coordlookup{N}(r::NTuple{N},x::NTuple{N}) = map(coordlookup,r,x)

function getindex{T}(C::CoordInterpGrid{T,1}, x::Real)
    getindex(C.grid,coordlookup(C.coord[1],x))
end
function getindex{T}(C::CoordInterpGrid{T,2}, x::Real, y::Real)
    getindex(C.grid,coordlookup(C.coord[1],x),coordlookup(C.coord[2],y))
end
function getindex{T}(C::CoordInterpGrid{T,3}, x::Real, y::Real, z::Real)
    getindex(C.grid,coordlookup(C.coord[1],x),coordlookup(C.coord[2],y),coordlookup(C.coord[3],z))
end
function getindex{T}(C::CoordInterpGrid{T,4}, x::Real, y::Real, z::Real, q::Real)
    getindex(C.grid,coordlookup(C.coord[1],x),coordlookup(C.coord[2],y),coordlookup(C.coord[3],z),coordlookup(C.coord[4],q))
end
function getindex{T,N}(C::CoordInterpGrid{T,N}, x::Real...)
    getindex(C.grid,coordlookup(C.coord,x)...)
end

# gradients
# chain rule: dy/dx = (dy/du) / (dx/du)

coordiscale(r::FloatRange,x::Real) = x*r.divisor/r.step
coordiscale(r::StepRange,x::Real) = x/r.step
coordiscale(r::UnitRange,x::Real) = x

function valgrad{T,BC,IT,R}(C::CoordInterpGrid{T,1,BC,IT,R}, x::Real)
    val, g = valgrad(C.grid,coordlookup(C.coord[1],x))
    val, coordiscale(C.coord[1],g)
end

function valgrad{T,N,BC,IT,R}(g::Vector{T}, C::CoordInterpGrid{T,N,BC,IT,R}, x::Real...)
    val = valgrad(g,C.grid,coordlookup(C.coord,x)...)
    for i=1:length(g)
        g[i] = coordiscale(C.coord[i],g[i])
    end
    val
end


text:
using Test, CompGrids, ParallelStencil, MPI, ImplicitGlobalGrid
# Parallel tests

@init_backend(ParallelStencil, Threads, true, 1); # Backend with MPI and ImplicitGlobalGrid

# 1D
grid = RegularRectilinearCollocatedGrid(size=(16,), x=(-5,10))
if mpirank==1
    @test grid.N[1] == 16
    @test grid.Ng[1] == 18
    @test grid.Nl[1] == 6
    @test length(grid.Face[1]) == grid.Nl[1]
    @test length(grid.Face_global[1]) == grid.N[1]
    @test grid.Face[1][1] ≈ -2.0
end

# 1D with a wider stencil
grid = RegularRectilinearCollocatedGrid(size=(16,), x=(-5,10), stencilwidth=2)
if mpirank==1
    @test grid.N[1] == 16
    @test grid.Ng[1] == 20
    @test grid.Nl[1] == 8
    @test length(grid.Face[1]) == grid.Nl[1]
    @test length(grid.Face_global[1]) == grid.N[1]
    @test grid.Face[1][1] ≈ -3.0
end

# 2D
grid = RegularRectilinearCollocatedGrid(size=(16,8), extent=(10,12))
if mpirank==3
    @test grid.N  == (16, 8)
    @test grid.Ng == (18, 10)
    @test grid.Nl == (10, 6)
end

# 3D
grid = RegularRectilinearCollocatedGrid(size=(16,8, 8), extent=(10,12, 20.5))
if mpirank==3
    @test grid.N  == (16, 8, 8)
    @test grid.Ng == (18, 10, 10)
    @test grid.Nl == (10, 6, 10)
end

# 2D - non-default different processor layout
opts=Dict(:dimy=>1,:dimx=>4)
grid = RegularRectilinearCollocatedGrid(size=(16,8), extent=(10,12), opts=opts)
if mpirank==1
    @test grid.N  == (16, 8)
    @test grid.Ng == (18, 10)
    @test grid.Nl == (6, 10)
end

# 2D - with fields
ParallelStencil.@reset_parallel_stencil()
@init_backend(ParallelStencil, Threads, true, 2); # Backend with MPI and ImplicitGlobalGrid
grid = RegularRectilinearCollocatedGrid(size=(16,8), extent=(10,12), fields=(T=0, P=2))

# Loop over all local points (including ghost) and set values
ind = LinearIndices(grid.Nl)
#idx = LinearIndices(grid.ind_local)

idx(id::CartesianIndex) = LinearIndices(grid.Nl)[id]

for i in (grid.ghostcorners.lower:grid.ghostcorners.upper)
    grid.fields.T[i] = grid.Face[1][i[1]]   + grid.Face[2][i[2]]     # set it to the y-coordinate
    grid.fields.P[idx(i)] = grid.Face[1][i[1]] + grid.Face[2][i[2]]
end

# Loop over non-ghost local points
for i in (grid.corners.lower:grid.corners.upper)
    grid.fields.P[i] = grid.Face_global[1][i[1]]
    grid.fields.T[grid.idx(i)] = grid.Face_global[1][i[1]]
end
@test grid.fields.P[2,2] == grid.Face_global[1][2]

diff = sum(abs.(grid.fields.P - grid.fields.T))
@test diff==0.0

if mpirank==2
    #@show grid.fields.P grid.fields.T diff
end

# 2D - with fields and specifying local grid size
ParallelStencil.@reset_parallel_stencil()
grid = RegularRectilinearCollocatedGrid(local_size=(16,8), extent=(10,12), fields=(T=0, P=1))
if mpirank==1
    @test grid.N  == (28, 12)
    @test grid.Ng == (30, 14)
    @test grid.Nl == (16, 8)
end