text:
include("../cli/cli_include.jl")

using GaPLAC
using Logging
using Random

Random.seed!(1)

global_logger(ConsoleLogger(stderr, Logging.Info))

gp = GaPLAC.parse_gp_formula("y : Gaussian(.01) ~| SExp(t; l=1.5)*Constant(1)", ["t", "y"], [false, false])

t = collect(-5.:0.1:5.)
f, y = GaPLAC.samplegp(gp.gp, [], [], [], [], t, t)

using DataFrames
using CSV
CSV.write("out.tsv", DataFrame(t=t, y=y), delim='\t')

using Plots
plot(t, y, size = (300, 200))
savefig("./gp.pdf")

gp = GaPLAC.parse_gp_formula("y*Reads/100 : Binomial(Reads) ~| Cat(Person) * SExp(Time) + Noise", ["Time", "y", "Reads", "Person"], [false, false, false, true])

function foo(f)
    f() + 1
end

function watest()
    f = eval(:(()->1))
    return Base.invokelatest(foo, f)
end

watest()

gp = GaPLAC.parse_gp_formula("y : Gaussian(.01) ~| Cat(person) * SExp(time)", ["person", "time"], [false, false])
x, z, y = GaPLAC.gp_inputs(gp, DataFrame(person=1:4, time=5:8))

cmd_sample(Dict(
    "data" => nothing,
    "atdata" => nothing,
    "mcmc" => nothing,
    "output" => "stdout",
    "formula" => "y :~| 1(1) * Cat(person) * SExp(time; l=1)",
    "at" => "person=1:3;time/person=range(-5,5,length=6)",
    "plot" => "sampleplot.png",
    "plotx" => "time:person"))

parsedgp = GaPLAC.parse_gp_formula("y : Gaussian(.1) ~| 1", ["y", "person", "time"], [false, false, false])

fw, ∇2lπ_fw, lπ = GaPLAC.laplace_approx(parsedgp.gp, GaPLAC.LowerTriangular(Matrix{Float64}(GaPLAC.I,(2,2))),
    [0.5, -0.5], zeros(2,0), [])

parsedgp = GaPLAC.parse_gp_formula("y ~| SExp(x)", ["y", "x"], [false, false])
x, z, y = GaPLAC.gp_inputs(parsedgp, DataFrame(x=-1:0.5:1, y=rand.([Uniform() for x in -1:0.5:1])))

chains = GaPLAC.mcmcgp(parsedgp.gp, x, y, z, [1., 1., 1.], 20)

max_stars_repo_name: biobakery/GaPLAC

max_issues_repo_name: biobakery/GaPLAC

max_forks_repo_name: biobakery/GaPLAC


text:
using FactorizationMachines
using Base.Test

include("test_read_libsvm.jl")
include("test_common.jl")
include("test_reco.jl")

include("matrix_factorization/test_identity.jl")
include("matrix_factorization/test_rank1.jl")

ext: jl


text:
module ReverseSearch

export reverse_search_producer, ReverseSearchParams

VERSION < v"0.4-" && using Docile

type ReverseSearchParams
    adj
    max_degere
    sink
    f
end

@doc doc"""
Equivalent to calling

    reverse_search_producer(p.adj, p.max_degree, p.sink, p.f)

Intended for use cases where `ReverseSearchParams` cleanly encapsulate some
problem domain.
""" ->
function reverse_search_producer(p::ReverseSearchParams)
    reverse_search_producer(p.adj, p.max_degere, p.sink, p.f)
end

@doc doc"""
# Input

- `adj` should be a function mapping `(vertex, n)` to `Union(child, Union{})`, where `child` is the `n`-th child of `vertex` in the traversal tree, or `Union{}` if `n` is greater than the degree of `vertex` in the tree
- `max_degree` should be a nonnegative integer denoting the maximum degree of a vertex in the tree
- `sink` should be a vertex
- `f` should map a vertex to its parent in the traversal tree

# Output
A producer for the nodes of the graph.
""" ->
function reverse_search_producer(adj, max_degree, sink, f)
    current_vertex = sink
    neighbor_counter = 0
    while true
        while neighbor_counter < max_degree
            neighbor_counter += 1
            next_vertex = adj(current_vertex, neighbor_counter)
            if (next_vertex != Union{}) && (f(next_vertex) == current_vertex)
                # reverse traverse (away from sink)
                current_vertex = next_vertex
                neighbor_counter = 0
            end
        end
        if current_vertex != sink
            produce(current_vertex)
            # forward traverse (towards sink)
            prev_vertex = current_vertex
            current_vertex = f(current_vertex)
            neighbor_counter = 0
            while true
                neighbor_counter += 1
                if adj(current_vertex, neighbor_counter) == prev_vertex
                    break
                end
            end
        end
        if current_vertex == sink && neighbor_counter == max_degree
            break
        end
    end
    produce(sink)
end

end # module


text:
#=――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
    NOTE
    This file provides a Quantity type to describe a number of a countable
    noun. Such objects are printed well, with correct pluralization.
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――=#

module Quantities

using ..Pluralize
using ..Semantics
import Base: isempty, length

export ItemQuantity, isnothing

"""
    ItemQuantity(n::Integer, item::AbstractString)

Represents a quantity of `n` occurrences of `item`. Although this is not a
collection, for ease of use, it implements some of the standard collection
methods `length` (for number of items) and `isempty` (for whether there are no
items).

```jldoctest
julia> using EnglishText

julia> ItemQuantity(2, "apple")
2 apples

julia> ItemQuantity(1, "standard canine")
1 standard canine
```
"""
struct ItemQuantity <: SemanticText
    count::Int
    noun::String
end

"""
    isempty(quantity::ItemQuantity)

Return `true` if the given `ItemQuantity` represents no items.

```jldoctest
julia> using EnglishText

julia> isempty(ItemQuantity(0, "orange"))
true

julia> isempty(ItemQuantity(4, "person"))
false
```
"""
isempty(quantity::ItemQuantity) = quantity.count == 0

"""
    length(quantity::ItemQuantity)

Return the number of items represented by this quantity.

```jldoctest
julia> using EnglishText

julia> length(ItemQuantity(7, "desk"))
7
```
"""
length(quantity::ItemQuantity) = quantity.count

function Base.show(io::IO, quantity::ItemQuantity)
    print(io, quantity.count)
    print(io, " ")
    if quantity.count == 1
        print(io, quantity.noun)
    else
        print(io, pluralize(quantity.noun))
    end
end

end

import .Quantities: ItemQuantity
export ItemQuantity


text:
function ASAF_actor_loss(πG, 𝒟_demo)
    (π, 𝒫, 𝒟,; info=Dict()) -> begin
        πsa_G = logpdf(π, 𝒟[:s], 𝒟[:a])
        πsa_E = logpdf(π, 𝒟_demo[:s], 𝒟_demo[:a])
        πGsa_G = logpdf(πG, 𝒟[:s], 𝒟[:a])
        πGsa_E = logpdf(πG, 𝒟_demo[:s], 𝒟_demo[:a])
        e = mean(entropy(π, 𝒟[:s]))

        ignore() do
            info[:entropy] = e
        end

        L = Flux.mean(log.(1 .+ exp.(πGsa_E - πsa_E))) + Flux.mean(log.(exp.(πsa_G - πGsa_G)  .+ 1))  - 0.1f0*e
        # if !isnothing(𝒟_nda)
        #     πsa_NDA = logpdf(π, 𝒟_nda[:s], 𝒟_nda[:a])
        #     πGsa_NDA = logpdf(πG, 𝒟_nda[:s], 𝒟_nda[:a])
        #     L += Flux.mean(log.(1 .+ exp.(πsa_NDA - πGsa_NDA)))
        # end
        L
    end
end

function ASAF(;π,
               S,
               𝒟_demo,
               normalize_demo::Bool=true,
               ΔN=50,
               λ_orth=1f-4,
               a_opt::NamedTuple=(;),
               c_opt::NamedTuple=(;),
               log::NamedTuple=(;),
               kwargs...)

    normalize_demo && (𝒟_demo = normalize!(deepcopy(𝒟_demo), S, action_space(π)))
    𝒟_demo = 𝒟_demo |> device(π)
    OnPolicySolver(;agent=PolicyParams(π),
                    S=S,
                    ΔN=ΔN,
                    post_batch_callback=(D; 𝒮, kwargs...) -> 𝒮.a_opt.loss = ASAF_actor_loss(deepcopy(𝒮.agent.π), 𝒟_demo),
                    log=LoggerParams(;dir="log/ASAF", period=100, log...),
                    a_opt=TrainingParams(;name="actor_", loss=nothing, a_opt...),
                    kwargs...)
end


text:
#
# coords.jl -
#
# Implement coordinate transforms.
#

module CoordinateTransforms

export
    CoordinateTransform,
    compose,
    jacobian,
    origin

"""
```julia
CoordinateTransform(stp, x0, y0) -> C
```

yields a simple 2D coordinate transform `C` such that:

```julia
C(u, v) -> (x0 + stp*u, y0 + stp*v)
```

The parameters of the coordinate transform `C` can be retrieved as follows:

```julia
step(C)   -> stp
origin(C) -> (x0, y0)
```

Input coordinates can be a 2-component tuple or vector (note that the
result is simlar to the argument):

```julia
C((u, v)) -> (x, y)
C([u, v]) -> [x, y]
```

Coordinate transforms may be inverted, pre/post-scaled, pre/post-translated:

```julia
A = inv(C)        # A(C(u,v)) = (u,v) and C(A(x,y)) = (x,y)
A = ρ*C           # A(u, v) = ρ*C(u, v)
A = C*ρ           # A(u, v) = C(ρ*u, ρ*v)
A = C + (tu, tv)  # A(u, v) = C(u + tu, v + tv)
A = C + [tu, tv]  # A(u, v) = C(u + tu, v + tv)
A = (tx, ty) + C  # A(u, v) = C(u, v) +. (tx, ty)
A = [tx, ty] + C  # A(u, v) = C(u, v) +. [tx, ty]
```

where `ρ` is some scalar factor.  The `-` operator can be used like `+` to
apply the opposite pre/post translation.

Two coordinate transforms `A` and `B` can be composed as follows:

```julia
C = A ∘ B
```

to form a new coordinate transform, `C` such that `C(u,v) = A(B(u,v))`.
Operators `\\` and  `/` can be used for left/right *division*:

```julia
A\\B = inv(A) ∘ B
A/B = A ∘ inv(B)
```

Two coordinate transforms `A` and `B` can be approximately compared:

```julia
A ≈ B
```

See also: [`compose`](@ref), [`jacobian`](@ref), [`origin`](@ref),
          [`step`](@ref).

"""
struct CoordinateTransform
    stp::Float64  # step size
    x0::Float64   # abscissa at (u=0,v)
    y0::Float64   # ordinate at (u,v=0)
end

# Yields the identity without arguments.
CoordinateTransform() = CoordinateTransform(1, 0, 0)

"""
```julia
step(C)
```

yields the step size the coordinate transform `C` (see
[`CoordinateTransform`](@ref)).

"""
Base.step(C::CoordinateTransform) = C.stp

"""
```julia
jacobian(C) -> jac
```

yields the Jacobian of the coordinate transform `C`.

See also: [`CoordinateTransform`](@ref), [`jacobian`](@ref), [`origin`](@ref).

"""
jacobian(C::CoordinateTransform) = step(C)^2

"""
```julia
origin(C) -> (x0, y0)
```

yields the origin of coordinates for the coordinate transform `C`; that is
the coordinates given by `C(0,0)`.

See also: [`CoordinateTransform`](@ref), [`jacobian`](@ref), [`step`](@ref).

"""
origin(C::CoordinateTransform) = (C.x0, C.y0)

compose(A::CoordinateTransform, B::CoordinateTransform) =
    CoordinateTransform(A.stp*B.stp,
                        A.x0 + A.stp*B.x0,
                        A.y0 + A.stp*B.y0)

Base.:∘(A::CoordinateTransform, B::CoordinateTransform) = compose(A, B)

Base.isapprox(A::CoordinateTransform, B::CoordinateTransform; kwds...) =
    (isapprox(A.stp, B.stp; kwds...) &&
     isapprox(A.x0,  B.x0;  kwds...) &&
     isapprox(A.y0,  B.y0;  kwds...))

# Union of possible coordinate types.
const Coordinates = Union{NTuple{2,Real},AbstractVector{<:Real}}

@inline __checkcoordinates(c::AbstractVector{Real}) =
    (length(c) == 2 || throw(ArgumentError("expecting 2 coordinates"));
     nothing)

@inline __checkcoordinates(c::NTuple{2,Real}) = nothing

(C::CoordinateTransform)(u::Real, v::Real) = (u*C.stp + C.x0,
                                              v*C.stp + C.y0)

(C::CoordinateTransform)(uv::NTuple{2,Real}) = C(uv[1], uv[2])

(C::CoordinateTransform)(uv::AbstractVector{Real}) =
    (__checkcoordinates(uv); return Float64[uv[1]*C.stp + C.x0,
                                            uv[2]*C.stp + C.y0])

Base.inv(C::CoordinateTransform) =
    (stp = 1/C.stp; return CoordinateTransform(stp, -stp*C.x0, -stp*C.y0))

Base.:\(A::CoordinateTransform, B::CoordinateTransform) =
    compose(inv(A), B)

Base.:/(A::CoordinateTransform, B::CoordinateTransform) =
    compose(A, inv(B))

#for op in (:(*), :(⋅), :(∘))
#    @eval begin
#
#        Base.$op(C::CoordinateTransform, uv::Coordinates) = C(uv)
#
#        Base.$op(A::CoordinateTransform, B::CoordinateTransform) =
#            compose(A, B)
#
#    end
#end

Base.:*(C::CoordinateTransform, ρ::Real) =
    CoordinateTransform(C.stp*ρ, C.x0, C.y0)

Base.:*(ρ::Real, C::CoordinateTransform) =
    CoordinateTransform(ρ*C.stp, ρ*C.x0, ρ*C.y0)

Base.:+(C::CoordinateTransform, t::Coordinates) =
    (__checkcoordinates(t);
     return CoordinateTransform(C.stp,
                                C.x0 + C.stp*t[1],
                                C.y0 + C.stp*t[2]))

Base.:-(C::CoordinateTransform) =
    CoordinateTransform(-C.stp, -C.x0, -C.y0)

Base.:-(C::CoordinateTransform, t::Coordinates) =
    (__checkcoordinates(t);
     return CoordinateTransform(C.stp,
                                C.x0 - C.stp*t[1],
                                C.y0 - C.stp*t[2]))

Base.:+(t::Coordinates, C::CoordinateTransform) =
    (__checkcoordinates(t);
     return CoordinateTransform(C.stp,
                                C.x0 + t[1],
                                C.y0 + t[2]))

Base.:-(t::Coordinates, C::CoordinateTransform) =
    (__checkcoordinates(t);
     return CoordinateTransform(-C.stp,
                                t[1] - C.x0,
                                t[2] - C.y0))

end # module


text:
using DrWatson
@quickactivate
using ArgParse
using CSV, DataFrames
using HierarchicalAD
using Flux, CUDA

arg_table = ArgParseSettings()
@add_arg_table arg_table begin
    "latent_dim"
        arg_type = Int
        default = 2
        help = "dimensionality of latent spaces"
    "hdim"
        arg_type = Int
        default = 8
        help = "width of the hidden layers in the discriminator"
    "channels"
        arg_type = Int
        nargs = '*'
        default = [16, 32, 64]
        help = "channel sizes"
    "--savepath"
        arg_type = String
        default = "test"
        help = "subdir of data/models"
    "--kernelsizes"
        arg_type = Int
        nargs = '*'
        default = [5, 3, 3]
        help = "kernelsizes"
    "--stride"
        help = "stride length, *to be implemented*"
        arg_type = Int
        default = 1
    "--layer_depth"
        arg_type = Int
        default = 1
        help = "depth of individual Conv blocks between latent spaces, *to be implemented*"
    "--seed"
        default = nothing
        help = "data split seed"
    "--nepochs"
        default = 20
        arg_type = Int
        help = "number of epochs"
    "--gamma"
        default = 0.1f0
        arg_type = Float32
        help = "scaling constant of total correlation in loss"
    "--lambda"
        default = 0f0
        arg_type = Float32
        help = "L2 regularization constant"
    "--test"
        action = :store_true
        help = "test run with limited data"
    "--batchsize"
        default = 128
        arg_type = Int
        help = "batchsize"
    "--activation"
        default = "relu"
        arg_type = String
        help = "activation function"
    "--xdist"
        default = :gaussian
        arg_type = Symbol
        help = "decoder distribution, one of gaussian/bernoulli"
    "--last_conv"
        action = :store_true
        help = "should the last layer of decoder be a dense or a conv layer"
    "--lr"
        help = "learning rate"
        arg_type = Float32
        default = 0.001f0
    "--gpu_id"
        help = "GPU device switch"
        arg_type = Int
        default = 0
    "--epochsize"
        help = "number of samples used in each epoch"
        default = nothing
    "--digit"
        help = "which digits to include"
        arg_type = Int
        nargs = '*'
        default = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
    "--area"
        help = "define the included area which is a number in [0,350] with mean of ~100, e.g. `--area >= 150`"
        arg_type = String
        nargs = 2
        default = [">=", "0"]
   "--length"
        help = "define the included digit length which is a number in [0,100] with mean of ~40"
        arg_type = String
        nargs = 2
        default = [">=", "0"]
   "--thickness"
        help = "define the included thickness which is a number in [0,10] with mean of ~3"
        arg_type = String
        nargs = 2
        default = [">=", "0"]
   "--slant"
        help = "define the included slant which is a number in [-1,1] with mean of ~0.1"
        arg_type = String
        nargs = 2
        default = [">=", "-1"]
   "--width"
        help = "define the included width which is a number in [0,25] with mean of ~14"
        arg_type = String
        nargs = 2
        default = [">=", "0"]
   "--height"
        help = "define the included height which is a number in [0,21] with mean of ~19.5"
        arg_type = String
        nargs = 2
        default = [">=", "0"]
end
args = parse_args(arg_table)
@unpack latent_dim, hdim, channels, kernelsizes, stride, layer_depth, last_conv, seed, test,
    lambda, batchsize, nepochs, gpu_id, epochsize, savepath, lr, activation, gamma, xdist = args
latent_count = length(channels)
(latent_count <= length(kernelsizes)) ? nothing : error("number of kernels and channels does not match.")
out_var = last_conv ? :conv : :dense
if seed != nothing
    seed = eval(Meta.parse(seed))
end
CUDA.device!(gpu_id)

# get filters
experiment_argnames = ["latent_dim", "hdim", "channels", "kernelsizes", "stride", "layer_depth", "last_conv",
    "seed", "lambda", "batchsize", "nepochs", "gpu_id", "epochsize", "savepath", "lr",
    "activation", "gamma", "xdist", "test"]
filter_dict, non_default_filters = HierarchicalAD.get_filter_info(experiment_argnames, args, arg_table)

# get the data
dataset = "morpho_mnist"
ratios = (0.8,0.199,0.001)
if test
    data, labels = HierarchicalAD.load_mnist("train")
    tr_x = HierarchicalAD.sample_tensor(data, 1000)
    val_x = HierarchicalAD.sample_tensor(data, 1000)
    tst_x = HierarchicalAD.sample_tensor(data, 100)
    tr_y, val_y, tst_y = nothing, nothing, nothing
    a_x, a_y = HierarchicalAD.sample_tensor(HierarchicalAD.load_mnist("test"), 1000), nothing
else
    (tr_x, tr_y), (val_x, val_y), (tst_x, tst_y), (a_x, a_y) =
        HierarchicalAD.load_train_val_test_data(dataset, filter_dict; ratios=ratios, seed=seed,
            categorical_key="digit")
end
if epochsize == nothing
    epochsize = size(tr_x, 4)
end

# now train the model
ncs = channels
ks = [(k,k) for k in kernelsizes][1:latent_count]
model, training_history, reconstructions, latent_representations =
    HierarchicalAD.train_fvlae(latent_dim, hdim, batchsize, ks, ncs, stride, nepochs, tr_x,
        val_x; λ=lambda, γ=gamma, epochsize=epochsize, layer_depth=layer_depth, lr=lr,
        var=out_var, activation=activation, xdist=xdist)

Flux.Zygote.ignore() do
    # compute scores
    tr_scores, val_scores, tst_scores, a_scores =
        map(x->HierarchicalAD.reconstruction_probability(gpu(model), x, 5, 16), (tr_x, val_x, tst_x, a_x))

    # compute encodings
    tr_encodings, val_encodings, tst_encodings, a_encodings =
        map(x->HierarchicalAD.encode_all(model,x,batchsize),(tr_x, val_x, tst_x, a_x))

    # now save everything
    model_id = HierarchicalAD.timetag()
    experiment_args = (model_id=model_id, data=dataset, latent_count=latent_count,
        latent_dim=latent_dim, channels=ncs, kernelsizes=ks, stride=stride, layer_depth=layer_depth,
        last_conv=last_conv, lr=lr, activation=activation, seed=seed, lambda=lambda,
        batchsize=batchsize, nepochs=nepochs, gpu_id=gpu_id, epochsize=epochsize, gamma=gamma,
        hdim=hdim, xdist=xdist, test=test)
    save_args = (model_id=model_id, data=dataset, model="fvlae", latent_dim=latent_dim,
        channels=channels,gamma=gamma, lambda=lambda, activation=activation, xdist=xdist)
    svn = HierarchicalAD.safe_savename(save_args, "bson", digits=5)
    svn = joinpath(datadir("models/$savepath"), svn)
    tagsave(svn, Dict(
            :model => cpu(model),
            :experiment_args => experiment_args,
            :filter_dict => filter_dict,
            :training_history => training_history,
            :reconstructions => reconstructions,
            :latent_representations => latent_representations,
            :ratios => ratios,
            :tr_scores => tr_scores,
            :val_scores => val_scores,
            :tst_scores => tst_scores,
            :a_scores => a_scores,
            :tr_encodings => tr_encodings,
            :val_encodings => val_encodings,
            :tst_encodings => tst_encodings,
            :a_encodings => a_encodings,
            :tr_labels => tr_y,
            :val_labels => val_y,
            :tst_labels => tst_y,
            :a_labels => a_y,
            :savepath => svn,
            :non_default_filters => non_default_filters
        ))
    @info "Results saved to $svn"
end

max_stars_repo_name: vitskvara/HierarchicalAD

max_issues_repo_name: vitskvara/HierarchicalAD

max_forks_repo_name: vitskvara/HierarchicalAD


text:
export
MarginalUpdateRule,
MarginalEntry,
MarginalTable,
marginalTable

"""
A `MarginalUpdateRule` specifies how a (joint) marginal is calculated from
incoming messages (and a node function).
"""
abstract type MarginalUpdateRule end
abstract type Product <: MarginalUpdateRule end

"""
A `MarginalEntry` defines a marginal computation.
The `marginal_update_rule <: MarginalUpdateRule` defines the rule that is used
to calculate the (joint) marginal over `target`.
"""
mutable struct MarginalEntry
    target::Region
    interfaces::Vector{Interface}
    marginal_update_rule::DataType

    # Fields for algorithm assembly
    marginal_id::Symbol # Specify the marginal identifier
    inbounds::Vector{Any} # Specify the inbounds required for the marginal update

    MarginalEntry() = new()
    MarginalEntry(target::Region, interfaces::Vector{Interface}, marginal_update_rule::DataType) = new(target, interfaces, marginal_update_rule)
end

"""
A `MarginalTable` defines the update order for marginal computations.
"""
const MarginalTable = Vector{MarginalEntry}

"""
Generate a `MarginalTable` that computes the marginals for target variables
"""
function marginalTable(targets::Vector{Variable})
    marginal_table = MarginalEntry[]
    for target in targets
        edge = first(target.edges) # For the sake of consistency, we always take the first edge as reference point for marginal computations.
        if edge.a == nothing # First handle cases where there is a `dangling` edge
            push!(marginal_table, MarginalEntry(target, [edge.b], Nothing))
        elseif edge.b == nothing
            push!(marginal_table, MarginalEntry(target, [edge.a], Nothing))
        elseif isClamped(edge)
            continue # Do not compute marginals for clamped edges
        else
            push!(marginal_table, MarginalEntry(target, [edge.a, edge.b], Product))
        end
    end

    return marginal_table
end

marginalTable(target::Variable) = marginalTable([target])

"""
Generate a mapping from target to marginal entry.
"""
function targetToMarginalEntry(table::MarginalTable)
    mapping = Dict{Region, MarginalEntry}()
    for entry in table
        target = entry.target
        mapping[target] = entry
    end

    return mapping
end


text:
function set_passive_or_contact_points(point_types::Array{PointType, 3}, potential::Array{T, 3},
                        grid::CylindricalGrid{T}, obj, pot::T, use_nthreads::Int = 1) where {T}
    if !isnan(pot)
        @onthreads 1:use_nthreads for iz in workpart(axes(potential, 3), 1:use_nthreads, Base.Threads.threadid())
            @inbounds for iφ in axes(potential, 2)
                for ir in axes(potential, 1)
                    pt::CylindricalPoint{T} = CylindricalPoint{T}( grid.axes[1].ticks[ir], grid.axes[2].ticks[iφ], grid.axes[3].ticks[iz] )
                    if pt in obj
                        potential[ ir, iφ, iz ] = pot
                        point_types[ ir, iφ, iz ] = zero(PointType)
                    end
                end
            end
        end
    end
    nothing
end

function set_point_types_and_fixed_potentials!(point_types::Array{PointType, 3}, potential::Array{T, 3},
        grid::CylindricalGrid{T}, det::SolidStateDetector{T};
        weighting_potential_contact_id::Union{Missing, Int} = missing,
        use_nthreads::Int = Base.Threads.nthreads(),
        not_only_paint_contacts::Val{NotOnlyPaintContacts} = Val{true}(),
        paint_contacts::Val{PaintContacts} = Val{true}())::Nothing where {T <: SSDFloat, NotOnlyPaintContacts, PaintContacts}

    @onthreads 1:use_nthreads for iz in workpart(axes(potential, 3), 1:use_nthreads, Base.Threads.threadid())
        @inbounds for iφ in axes(potential, 2)
            for ir in axes(potential, 1)
                pt::CylindricalPoint{T} = CylindricalPoint{T}( grid.axes[1].ticks[ir], grid.axes[2].ticks[iφ], grid.axes[3].ticks[iz] )
                if in(pt, det.semiconductor)
                    point_types[ ir, iφ, iz ] += pn_junction_bit
                end
            end
        end
    end
    isEP = ismissing(weighting_potential_contact_id)
    if !ismissing(det.passives)
        for passive in det.passives
            pot::T = isEP ? passive.potential : (isnan(passive.potential) ? passive.potential : zero(T))
            set_passive_or_contact_points(point_types, potential, grid, passive.geometry, pot, use_nthreads)
        end
    end
    if NotOnlyPaintContacts
        for contact in det.contacts
            pot::T = isEP ? contact.potential : contact.id == weighting_potential_contact_id
            set_passive_or_contact_points(point_types, potential, grid, contact.geometry, pot, use_nthreads)
        end
    end
    if PaintContacts
        for contact in det.contacts
            pot::T = isEP ? contact.potential : contact.id == weighting_potential_contact_id
            fs = ConstructiveSolidGeometry.surfaces(contact.geometry)
            for face in fs
                paint!(point_types, potential, face, contact.geometry, pot, grid)
            end
        end
    end
    nothing
end

function fill_ρimp_ϵ_ρfix(ρ_eff_imp_tmp::Array{T}, ϵ::Array{T}, ρ_eff_fix_tmp::Array{T},
    ::Type{Cylindrical}, mpz::Vector{T}, mpφ::Vector{T}, mpr::Vector{T}, axr::Vector{T}, use_nthreads::Int, obj) where {T}
    @inbounds begin
        @onthreads 1:use_nthreads for iz in workpart(axes(ϵ, 3), 1:use_nthreads, Base.Threads.threadid())
            pos_z::T = mpz[iz]
            for iφ in axes(ϵ, 2)
                pos_φ::T = mpφ[iφ]
                for ir in axes(ϵ, 1)
                    pos_r::T = mpr[ir]
                    if (ir == 1 && axr[1] == 0) pos_r = axr[2] * 0.5 end
                    pt::CylindricalPoint{T} = CylindricalPoint{T}(pos_r, pos_φ, pos_z)
                    if pt in obj
                        ρ_eff_imp_tmp[ir, iφ, iz]::T, ϵ[ir, iφ, iz]::T, ρ_eff_fix_tmp[ir, iφ, iz]::T = get_ρimp_ϵ_ρfix(pt, obj)
                    end
                end
            end
        end
    end
    nothing
end

function PotentialCalculationSetup(det::SolidStateDetector{T}, grid::CylindricalGrid{T},
                medium::NamedTuple = material_properties[materials["vacuum"]],
                potential_array::Union{Missing, Array{T, 3}} = missing,
                imp_scale::Union{Missing, Array{T, 3}} = missing;
                weighting_potential_contact_id::Union{Missing, Int} = missing,
                point_types = missing,
                use_nthreads::Int = Base.Threads.nthreads(),
                sor_consts = (1.0, 1.0),
                not_only_paint_contacts::Bool = true, paint_contacts::Bool = true)::PotentialCalculationSetup{T} where {T}
    r0_handling::Bool = typeof(grid.axes[1]).parameters[2] == :r0
    only_2d::Bool = length(grid.axes[2]) == 1 ? true : false
    @assert grid.axes[1][1] == 0 "Something is wrong. R-axis has `:r0`-boundary handling but first tick is $(axr[1]) and not 0."

    is_weighting_potential::Bool = !ismissing(weighting_potential_contact_id)
    depletion_handling::Bool = is_weighting_potential && !ismissing(point_types)

    @inbounds begin
        begin # Geometrical weights of the Axes
            nr = size(grid)[1]

            # R-axis
            axr::Vector{T} = collect(grid.axes[1]) # real grid points/ticks -> the potential at these ticks are going to be calculated
            r_inv::Vector{T} = inv.(axr)
            if r0_handling r_inv[1] = inv(axr[2] * 0.5) end
            r_ext::Vector{T} = get_extended_ticks(grid.axes[1])
            Δr_ext::Vector{T} = diff(r_ext)
            Δr_ext_inv::Vector{T} = inv.(Δr_ext)
            mpr::Vector{T} = midpoints(r_ext)
            Δmpr::Vector{T} = diff(mpr)
            Δmpr_inv::Vector{T} = inv.(Δmpr)
            Δmpr_squared::Vector{T} = T(0.5) .* ((mpr[2:end].^2) .- (mpr[1:end-1].^2))
            if r0_handling
                Δmpr_squared[1] = T(0.5) * (mpr[2]^2)
                Δr_ext_inv[1] = 0 # -> left weight @r=0 becomes 0 through this
            end
            Δhmprr::Vector{T} = mpr[2:end] - axr # distances between midpoints and real grid points (half distances -> h), needed for weights wrr, wrl, ... & dV (volume element)
            Δhmprl::Vector{T} = axr - mpr[1:end - 1]
            wrr::Vector{T} = Δmpr_inv .* Δhmprr # weights for epislon_r adding
            wrl::Vector{T} = Δmpr_inv .* Δhmprl
            if r0_handling wrl[1]::T = 1 - wrr[1] end
            wr::Array{T, 2} = zeros(T, 6, length(wrr))
            wr[1, :] = wrr
            wr[2, :] = wrl
            wr[3, :] = r_inv .* Δmpr
            wr[4, :] = Δr_ext_inv[2:end] .* mpr[2:end]
            wr[5, :] = Δr_ext_inv[1:length(wrr)] .* mpr[1:length(wrr)]
            wr[6, :] = Δmpr_squared
            gw_r::GeometricalRadialAxisWeights{T} = GeometricalRadialAxisWeights{T}( wr ) # Weights needed for Field Simulation loop
            r_ext_mid::T = (r_ext[end - 1] - r_ext[2]) / 2
            grid_boundary_factor_r_left::T = abs((r_ext[2] - r_ext_mid) / (r_ext[1] - r_ext_mid))
            grid_boundary_factor_r_right::T = r0_handling ? abs(r_ext[end - 1] / r_ext[end]) : abs((r_ext[end - 1] - r_ext_mid) / (r_ext[end] - r_ext_mid))

            # φ-axis
            axφ::Vector{T} = collect(grid.axes[2]) # real grid points/ticks -> the potential at these ticks are going to be calculated
            φ_ext::Vector{T} = get_extended_ticks(grid.axes[2])

            Δφ_ext::Vector{T} = diff(φ_ext)
            Δφ_ext_inv::Vector{T} = inv.(Δφ_ext)
            mpφ::Vector{T} = midpoints(φ_ext)
            Δmpφ::Vector{T} = diff(mpφ)
            Δmpφ_inv::Vector{T} = inv.(Δmpφ)
            Δhmpφr::Vector{T} = mpφ[2:end] - axφ # distances between midpoints and real grid points (half distances -> h), needed for weights wrr, wrl, ... & dV (volume element)
            Δhmpφl::Vector{T} = axφ - mpφ[1:end - 1]
            wφr::Vector{T} = Δmpφ_inv .* Δhmpφr # weights for epislon_r adding
            wφl::Vector{T} = Δmpφ_inv .* Δhmpφl
            wφ::Array{T, 2} = zeros(T, 4, length(wφr) + 1)
            wφ[1, 1:length(wφr)] = wφr
            wφ[2, 1:length(wφr)] = wφl
            wφ[3, 1:length(wφr)] = Δmpφ
            wφ[4, :] = Δφ_ext_inv
            gw_φ::GeometricalAzimutalAxisWeights{T} = GeometricalAzimutalAxisWeights{T}( wφ ) # Weights needed for Field Simulation loop
            φ_ext_mid::T = (φ_ext[end - 1] - φ_ext[2]) / 2
            grid_boundary_factor_φ_right::T = abs((φ_ext[end - 1] - φ_ext_mid) / (φ_ext[end] - φ_ext_mid))
            grid_boundary_factor_φ_left::T = abs((φ_ext[2] - φ_ext_mid) / (φ_ext[1] - φ_ext_mid))

            # Z-axis
            axz::Vector{T} = collect(grid.axes[3]) # real grid points/ticks -> the potential at these ticks are going to be calculated
            z_ext::Vector{T} = get_extended_ticks(grid.axes[3])
            Δz_ext::Vector{T} = diff(z_ext)
            Δz_ext_inv::Vector{T} = inv.(Δz_ext)
            mpz::Vector{T} = midpoints(z_ext)
            Δmpz::Vector{T} = diff(mpz)
            Δmpz_inv::Vector{T} = inv.(Δmpz)
            Δhmpzr::Vector{T} = mpz[2:end] - axz # distances between midpoints and real grid points (half distances -> h), needed for weights wrr, wrl, ... & dV (volume element)
            Δhmpzl::Vector{T} = axz - mpz[1:end - 1]
            wzr::Vector{T} = Δmpz_inv .* Δhmpzr # weights for epislon_r adding
            wzl::Vector{T} = Δmpz_inv .* Δhmpzl
            wz::Array{T, 2} = zeros(T, 4, length(wzr) + 1)
            wz[1, 1:length(wzr)] = wzr
            wz[2, 1:length(wzr)] = wzl
            wz[3, 1:length(wzr)] = Δmpz
            wz[4, :] = Δz_ext_inv
            gw_z::GeometricalCartesianAxisWeights{T} = GeometricalCartesianAxisWeights{T}( wz ) # Weights needed for Field Simulation loop
            z_ext_mid::T = (z_ext[end - 1] - z_ext[2]) / 2
            grid_boundary_factor_z_right::T = abs((z_ext[end - 1] - z_ext_mid) / (z_ext[end] - z_ext_mid))
            grid_boundary_factor_z_left::T = abs((z_ext[2] - z_ext_mid) / (z_ext[1] - z_ext_mid))

            geom_weights::NTuple{3, AbstractGeometricalAxisWeights{T}} = (gw_r, gw_φ, gw_z) # Weights needed for Field Simulation loop
            grid_boundary_factors::NTuple{3, NTuple{2, T}} = ((grid_boundary_factor_r_left, grid_boundary_factor_r_right), (grid_boundary_factor_φ_left, grid_boundary_factor_φ_right), (grid_boundary_factor_z_left, grid_boundary_factor_z_right))
        end

        bias_voltages::Vector{T} = if length(det.contacts) > 0
            [i.potential for i in det.contacts]
        else
            T[0]
        end
        minimum_applied_potential::T = minimum(bias_voltages)
        maximum_applied_potential::T = maximum(bias_voltages)
        bias_voltage::T = maximum_applied_potential - minimum_applied_potential
        sor_slope = (sor_consts[2] .- sor_consts[1]) / (nr - 1 )
        sor_const::Vector{T} = T[ sor_consts[1] + (i - 1) * sor_slope for i in 1:nr]

        medium_ϵ_r::T = medium.ϵ_r
        ϵ = fill(medium_ϵ_r, length(mpr), length(mpφ), length(mpz))
        ρ_eff_imp_tmp = zeros(T, length(mpr), length(mpφ), length(mpz))
        ρ_eff_fix_tmp = zeros(T, length(mpr), length(mpφ), length(mpz))
        fill_ρimp_ϵ_ρfix(ρ_eff_imp_tmp, ϵ, ρ_eff_fix_tmp, Cylindrical, mpz, mpφ, mpr, axr, use_nthreads, det.semiconductor)
        if !ismissing(det.passives)
            for passive in det.passives
                fill_ρimp_ϵ_ρfix(ρ_eff_imp_tmp, ϵ, ρ_eff_fix_tmp, Cylindrical, mpz, mpφ, mpr, axr, use_nthreads, passive)
            end
        end
        if depletion_handling
            for iz in axes(ϵ, 3)
                pos_z = mpz[iz]
                for iφ in axes(ϵ, 2)
                    pos_φ = mpφ[iφ]
                    for ir in axes(ϵ, 1)
                        pos_r = mpr[ir]
                        if (ir == 1 && axr[1] == 0) pos_r = axr[2] * 0.5 end
                        pt::CylindricalPoint{T} = CylindricalPoint{T}(pos_r, pos_φ, pos_z)
                        ig = find_closest_gridpoint(pt, point_types.grid)
                        if is_undepleted_point_type(point_types.data[ig...])
                            ϵ[ir, iφ, iz] *= scaling_factor_for_permittivity_in_undepleted_region(det.semiconductor) * (1 - imp_scale[ig...])
                        elseif is_fixed_point_type(point_types.data[ig...])
                            ϵ[ir, iφ, iz] *= scaling_factor_for_permittivity_in_undepleted_region(det.semiconductor)
                        end
                    end
                end
            end
        end

        ϵ0_inv::T = inv(ϵ0)
        ρ_eff_imp_tmp *= ϵ0_inv
        ρ_eff_fix_tmp *= ϵ0_inv

        volume_weights::Array{T, 4} = RBExtBy2Array(T, grid)
        q_eff_imp::Array{T, 4} = RBExtBy2Array(T, grid)
        q_eff_fix::Array{T, 4} = RBExtBy2Array(T, grid)
        for iz in range(2, stop = length(z_ext) - 1)
            inz::Int = iz - 1
            irbz::Int = rbidx(inz)
            for iφ in range(2, stop = length(φ_ext) - 1)
                inφ::Int = iφ - 1
                for ir in range(2, stop = length(r_ext) - 1)
                    inr::Int = ir - 1;

                    rbi::Int = iseven(inr + inφ + inz) ? rb_even::Int : rb_odd::Int
                    # rbinds = irbz, iφ, ir, rbi

                    ρ_imp_cell::T = 0
                    ρ_eff_fix_cell::T = 0
                    if !is_weighting_potential
                        if inr > 1
                            ρ_imp_cell += ρ_eff_imp_tmp[ ir,  iφ,  iz] * wzr[inz] * wrr[inr] * wφr[inφ]
                            ρ_imp_cell += ρ_eff_imp_tmp[ ir,  iφ, inz] * wzl[inz] * wrr[inr] * wφr[inφ]
                            ρ_imp_cell += ρ_eff_imp_tmp[ ir, inφ,  iz] * wzr[inz] * wrr[inr] * wφl[inφ]
                            ρ_imp_cell += ρ_eff_imp_tmp[ ir, inφ, inz] * wzl[inz] * wrr[inr] * wφl[inφ]

                            ρ_imp_cell += ρ_eff_imp_tmp[inr,  iφ,  iz] * wzr[inz] * wrl[inr] * wφr[inφ]
                            ρ_imp_cell += ρ_eff_imp_tmp[inr,  iφ, inz] * wzl[inz] * wrl[inr] * wφr[inφ]
                            ρ_imp_cell += ρ_eff_imp_tmp[inr, inφ,  iz] * wzr[inz] * wrl[inr] * wφl[inφ]
                            ρ_imp_cell += ρ_eff_imp_tmp[inr, inφ, inz] * wzl[inz] * wrl[inr] * wφl[inφ]

                            ρ_eff_fix_cell += ρ_eff_fix_tmp[ ir,  iφ,  iz] * wzr[inz] * wrr[inr] * wφr[inφ]
                            ρ_eff_fix_cell += ρ_eff_fix_tmp[ ir,  iφ, inz] * wzl[inz] * wrr[inr] * wφr[inφ]
                            ρ_eff_fix_cell += ρ_eff_fix_tmp[ ir, inφ,  iz] * wzr[inz] * wrr[inr] * wφl[inφ]
                            ρ_eff_fix_cell += ρ_eff_fix_tmp[ ir, inφ, inz] * wzl[inz] * wrr[inr] * wφl[inφ]

                            ρ_eff_fix_cell += ρ_eff_fix_tmp[inr,  iφ,  iz] * wzr[inz] * wrl[inr] * wφr[inφ]
                            ρ_eff_fix_cell += ρ_eff_fix_tmp[inr,  iφ, inz] * wzl[inz] * wrl[inr] * wφr[inφ]
                            ρ_eff_fix_cell += ρ_eff_fix_tmp[inr, inφ,  iz] * wzr[inz] * wrl[inr] * wφl[inφ]
                            ρ_eff_fix_cell += ρ_eff_fix_tmp[inr, inφ, inz] * wzl[inz] * wrl[inr] * wφl[inφ]
                        else
                            ρ_imp_cell += ρ_eff_imp_tmp[ ir,  iφ,  iz] * wzr[inz] * 0.5 #wφr[inφ]
                            ρ_imp_cell += ρ_eff_imp_tmp[ ir,  iφ, inz] * wzl[inz] * 0.5 #wφr[inφ]
                            ρ_imp_cell += ρ_eff_imp_tmp[ ir, inφ,  iz] * wzr[inz] * 0.5 #wφl[inφ]
                            ρ_imp_cell += ρ_eff_imp_tmp[ ir, inφ, inz] * wzl[inz] * 0.5 #wφl[inφ]

                            ρ_eff_fix_cell += ρ_eff_fix_tmp[ ir,  iφ,  iz] * wzr[inz] * 0.5 #wφr[inφ]
                            ρ_eff_fix_cell += ρ_eff_fix_tmp[ ir,  iφ, inz] * wzl[inz] * 0.5 #wφr[inφ]
                            ρ_eff_fix_cell += ρ_eff_fix_tmp[ ir, inφ,  iz] * wzr[inz] * 0.5 #wφl[inφ]
                            ρ_eff_fix_cell += ρ_eff_fix_tmp[ ir, inφ, inz] * wzl[inz] * 0.5 #wφl[inφ]
                        end
                    end
                    if inr > 1
                        wrr_eps::T = ϵ[  ir,  iφ, inz + 1] * wφr[inφ] * wzr[inz]
                        wrr_eps   += ϵ[  ir, inφ, inz + 1] * wφl[inφ] * wzr[inz]
                        wrr_eps   += ϵ[  ir,  iφ, inz ]    * wφr[inφ] * wzl[inz]
                        wrr_eps   += ϵ[  ir, inφ, inz ]    * wφl[inφ] * wzl[inz]
                        # # left weight in r: wrr
                        wrl_eps::T = ϵ[ inr,  iφ, inz + 1] * wφr[inφ] * wzr[inz]
                        wrl_eps   += ϵ[ inr, inφ, inz + 1] * wφl[inφ] * wzr[inz]
                        wrl_eps   += ϵ[ inr,  iφ, inz ]    * wφr[inφ] * wzl[inz]
                        wrl_eps   += ϵ[ inr, inφ, inz ]    * wφl[inφ] * wzl[inz]
                        # right weight in φ: wφr
                        wφr_eps::T = ϵ[ inr,  iφ, inz + 1] * wrl[inr] * wzr[inz]
                        wφr_eps   += ϵ[  ir,  iφ, inz + 1] * wrr[inr] * wzr[inz]
                        wφr_eps   += ϵ[ inr,  iφ, inz ]    * wrl[inr] * wzl[inz]
                        wφr_eps   += ϵ[  ir,  iφ, inz ]    * wrr[inr] * wzl[inz]
                        # left weight in φ: wφl
                        wφl_eps::T = ϵ[ inr, inφ, inz + 1] * wrl[inr] * wzr[inz]
                        wφl_eps   += ϵ[  ir, inφ, inz + 1] * wrr[inr] * wzr[inz]
                        wφl_eps   += ϵ[ inr, inφ, inz ]    * wrl[inr] * wzl[inz]
                        wφl_eps   += ϵ[  ir, inφ, inz ]    * wrr[inr] * wzl[inz]
                        # right weight in z: wzr
                        wzr_eps::T = ϵ[  ir,  iφ, inz + 1] * wrr[inr] * wφr[inφ]
                        wzr_eps   += ϵ[  ir, inφ, inz + 1] * wrr[inr] * wφl[inφ]
                        wzr_eps   += ϵ[ inr,  iφ, inz + 1] * wrl[inr] * wφr[inφ]
                        wzr_eps   += ϵ[ inr, inφ, inz + 1] * wrl[inr] * wφl[inφ]
                        # left weight in z: wzr
                        wzl_eps::T = ϵ[  ir,  iφ,    inz ] * wrr[inr] * wφr[inφ]
                        wzl_eps   += ϵ[  ir, inφ,    inz ] * wrr[inr] * wφl[inφ]
                        wzl_eps   += ϵ[ inr,  iφ,    inz ] * wrl[inr] * wφr[inφ]
                        wzl_eps   += ϵ[ inr, inφ,    inz ] * wrl[inr] * wφl[inφ]

                        volume_weight::T = wrr_eps * Δr_ext_inv[ ir] * mpr[ ir] * Δmpφ[inφ] * Δmpz[inz]
                        volume_weight += wrl_eps * Δr_ext_inv[inr] * mpr[inr] * Δmpφ[inφ] * Δmpz[inz]
                        volume_weight += wφr_eps * r_inv[inr] * Δφ_ext_inv[ iφ] * Δmpr[inr] * Δmpz[inz]
                        volume_weight += wφl_eps * r_inv[inr] * Δφ_ext_inv[inφ] * Δmpr[inr] * Δmpz[inz]
                        volume_weight += wzr_eps * Δz_ext_inv[ iz] * Δmpφ[inφ] * Δmpr_squared[inr]
                        volume_weight += wzl_eps * Δz_ext_inv[inz] * Δmpφ[inφ] * Δmpr_squared[inr]

                        volume_weights[ irbz, iφ, ir, rbi ] = inv(volume_weight)

                        dV::T = Δmpz[inz] * Δmpφ[inφ] * Δmpr_squared[inr]
                        q_eff_imp[ irbz, iφ, ir, rbi ] = dV * ρ_imp_cell
                        q_eff_fix[ irbz, iφ, ir, rbi ] = dV * ρ_eff_fix_cell
                    else
                        wrr_eps = ϵ[  ir,  iφ, inz + 1] * 0.5 * wzr[inz]
                        wrr_eps   += ϵ[  ir, inφ, inz + 1] * 0.5 * wzr[inz]
                        wrr_eps   += ϵ[  ir,  iφ, inz ]    * 0.5 * wzl[inz]
                        wrr_eps   += ϵ[  ir, inφ, inz ]    * 0.5 * wzl[inz]
                        # # left weight in r: wrr
                        wrl_eps = ϵ[ inr,  iφ, inz + 1] * 0.5 * wzr[inz]
                        wrl_eps   += ϵ[ inr, inφ, inz + 1] * 0.5 * wzr[inz]
                        wrl_eps   += ϵ[ inr,  iφ, inz ]    * 0.5 * wzl[inz]
                        wrl_eps   += ϵ[ inr, inφ, inz ]    * 0.5 * wzl[inz]
                        # right weight in φ: wφr
                        wφr_eps = ϵ[ inr,  iφ, inz + 1] * wrl[inr] * wzr[inz]
                        wφr_eps   += ϵ[  ir,  iφ, inz + 1] * wrr[inr] * wzr[inz]
                        wφr_eps   += ϵ[ inr,  iφ, inz ]    * wrl[inr] * wzl[inz]
                        wφr_eps   += ϵ[  ir,  iφ, inz ]    * wrr[inr] * wzl[inz]
                        # left weight in φ: wφl
                        wφl_eps = ϵ[ inr, inφ, inz + 1] * wrl[inr] * wzr[inz]
                        wφl_eps   += ϵ[  ir, inφ, inz + 1] * wrr[inr] * wzr[inz]
                        wφl_eps   += ϵ[ inr, inφ, inz ]    * wrl[inr] * wzl[inz]
                        wφl_eps   += ϵ[  ir, inφ, inz ]    * wrr[inr] * wzl[inz]
                        # right weight in z: wzr
                        wzr_eps = ϵ[  ir,  iφ, inz + 1] * wrr[inr] * 0.5
                        wzr_eps   += ϵ[  ir, inφ, inz + 1] * wrr[inr] * 0.5
                        wzr_eps   += ϵ[ inr,  iφ, inz + 1] * wrl[inr] * 0.5
                        wzr_eps   += ϵ[ inr, inφ, inz + 1] * wrl[inr] * 0.5
                        # left weight in z: wzr
                        wzl_eps = ϵ[  ir,  iφ,    inz ] * wrr[inr] * 0.5
                        wzl_eps   += ϵ[  ir, inφ,    inz ] * wrr[inr] * 0.5
                        wzl_eps   += ϵ[ inr,  iφ,    inz ] * wrl[inr] * 0.5
                        wzl_eps   += ϵ[ inr, inφ,    inz ] * wrl[inr] * 0.5

                        volume_weight = wrr_eps * Δr_ext_inv[ ir] * mpr[ ir] * 2π * Δmpz[inz]
                        volume_weight += wrl_eps * Δr_ext_inv[inr] * mpr[inr] * Δmpφ[inφ] * Δmpz[inz]
                        volume_weight += wφr_eps * r_inv[inr] * 0.15915494f0 * Δmpr[inr] * Δmpz[inz]
                        volume_weight += wφl_eps * r_inv[inr] * 0.15915494f0 * Δmpr[inr] * Δmpz[inz]
                        volume_weight += wzr_eps * Δz_ext_inv[ iz] * 2π * Δmpr_squared[inr]
                        volume_weight += wzl_eps * Δz_ext_inv[inz] * 2π * Δmpr_squared[inr]

                        volume_weights[ irbz, iφ, ir, rbi ] = inv(volume_weight)

                        dV = Δmpz[inz] * 2π * Δmpr_squared[inr]
                        q_eff_imp[ irbz, iφ, ir, rbi ] = dV * ρ_imp_cell
                        q_eff_fix[ irbz, iφ, ir, rbi ] = dV * ρ_eff_fix_cell
                    end
                end
            end
        end
        potential = ismissing(potential_array) ? zeros(T, size(grid)...) : potential_array
        point_types = ones(PointType, size(grid)...)
        set_point_types_and_fixed_potentials!( point_types, potential, grid, det,
                weighting_potential_contact_id = weighting_potential_contact_id,
                use_nthreads = use_nthreads,
                not_only_paint_contacts = Val(not_only_paint_contacts),
                paint_contacts = Val(paint_contacts)  )
        rbpotential  = RBExtBy2Array( potential, grid )
        rbpoint_types = RBExtBy2Array( point_types, grid )
    end # @inbounds

    pcs = PotentialCalculationSetup(
        grid,
        rbpotential,
        rbpoint_types,
        volume_weights,
        q_eff_imp,
        ones(T, size(rbpoint_types)),
        q_eff_fix,
        ϵ,
        broadcast(gw -> gw.weights, geom_weights),
        sor_const,
        bias_voltage,
        maximum_applied_potential,
        minimum_applied_potential,
        grid_boundary_factors
     )

    apply_boundary_conditions!(pcs, Val{ true}(), only_2d ? Val{true}() : Val{false}()) # even points
    apply_boundary_conditions!(pcs, Val{false}(), only_2d ? Val{true}() : Val{false}()) # odd
    return pcs
end

Grid(pcs::PotentialCalculationSetup) = pcs.grid

function ElectricPotentialArray(pcs::PotentialCalculationSetup{T, Cylindrical, 3, Array{T, 3}})::Array{T, 3} where {T}
    pot::Array{T, 3} = Array{T, 3}(undef, size(pcs.grid))
    for iz in axes(pot, 3)
        irbz::Int = rbidx(iz)
        for iφ in axes(pot, 2)
            irbφ::Int = iφ + 1
            idxsum::Int = iz + iφ
            for ir in axes(pot, 1)
                irbr::Int = ir + 1
                rbi::Int = iseven(idxsum + ir) ? rb_even::Int : rb_odd::Int
                pot[ir, iφ, iz] = pcs.potential[ irbz, irbφ, irbr, rbi ]
            end
        end
    end
    for iz in axes(pot, 3)
        p_r0 = mean(pot[1,:,iz])
        pot[1,:,iz] .= p_r0
    end
    return pot
end

function ImpurityScaleArray(pcs::PotentialCalculationSetup{T, Cylindrical, 3, Array{T, 3}})::Array{T, 3} where {T}
    s::Array{T, 3} = Array{T, 3}(undef, size(pcs.grid))
    for iz in axes(s, 3)
        irbz::Int = rbidx(iz)
        for iφ in axes(s, 2)
            irbφ::Int = iφ + 1
            idxsum::Int = iz + iφ
            for ir in axes(s, 1)
                irbr::Int = ir + 1
                rbi::Int = iseven(idxsum + ir) ? rb_even::Int : rb_odd::Int
                s[ir, iφ, iz] = pcs.imp_scale[ irbz, irbφ, irbr, rbi ]
            end
        end
    end
    return s
end

function PointTypeArray(pcs::PotentialCalculationSetup{T, Cylindrical, 3, Array{T, 3}})::Array{PointType, 3} where {T}
    point_types::Array{PointType, 3} = zeros(PointType, size(pcs.grid))
    for iz in axes(point_types, 3)
        irbz::Int = rbidx(iz)
        for iφ in axes(point_types, 2)
            irbφ::Int = iφ + 1
            idxsum::Int = iz + iφ
            for ir in axes(point_types, 1)
                irbr::Int = ir + 1
                rbi::Int = iseven(idxsum + ir) ? rb_even::Int : rb_odd::Int
                point_types[ir, iφ, iz] = pcs.point_types[irbz, irbφ, irbr, rbi ]
            end
        end
    end
    return point_types
end

function EffectiveChargeDensityArray(pcs::PotentialCalculationSetup{T, Cylindrical, 3, Array{T, 3}})::Array{T} where {T}
    ρ::Array{T, 3} = zeros(T, size(pcs.grid))
    for iz in axes(ρ, 3)
        irbz::Int = rbidx(iz)
        Δmpz::T = pcs.geom_weights[3][3, iz]
        for iφ in axes(ρ, 2)
            irbφ::Int = iφ + 1
            idxsum::Int = iz + iφ
            Δmpφ::T = pcs.geom_weights[2][3, iφ]
            Δmpzφ::T = Δmpz * Δmpφ
            for ir in axes(ρ, 1)
                irbr::Int = ir + 1
                rbi::Int = iseven(idxsum + ir) ? rb_even::Int : rb_odd::Int
                dV::T = pcs.geom_weights[1][6, ir] * Δmpzφ  #Δmpz[inz] * Δmpφ[inφ] * Δmpr_squared[inr]
                if ir == 1
                    dV = dV * 2π / Δmpφ
                end
                ρ[ir, iφ, iz] = pcs.q_eff_imp[irbz, irbφ, irbr, rbi ] / dV
            end
        end
    end
    return ρ
end

function FixedEffectiveChargeDensityArray(pcs::PotentialCalculationSetup{T, Cylindrical, 3, Array{T, 3}})::Array{T} where {T}
    ρ::Array{T, 3} = zeros(T, size(pcs.grid))
    for iz in axes(ρ, 3)
        irbz::Int = rbidx(iz)
        Δmpz::T = pcs.geom_weights[3][3, iz]
        for iφ in axes(ρ, 2)
            irbφ::Int = iφ + 1
            idxsum::Int = iz + iφ
            Δmpφ::T = pcs.geom_weights[2][3, iφ]
            Δmpzφ::T = Δmpz * Δmpφ
            for ir in axes(ρ, 1)
                irbr::Int = ir + 1
                rbi::Int = iseven(idxsum + ir) ? rb_even::Int : rb_odd::Int
                dV::T = pcs.geom_weights[1][6, ir] * Δmpzφ  #Δmpz[inz] * Δmpφ[inφ] * Δmpr_squared[inr]
                if ir == 1
                    dV = dV * 2π / Δmpφ
                end
                ρ[ir, iφ, iz] = pcs.q_eff_fix[irbz, irbφ, irbr, rbi ] / dV
            end
        end
    end
    return ρ
end

function DielectricDistributionArray(pcs::PotentialCalculationSetup{T, S, 3, Array{T, 3}})::Array{T, 3} where {T, S}
    return pcs.ϵ_r
end


text:
module PatientNurseAllocation

using JuMP
using Gurobi

using LinearAlgebra
using Distributions
using MathOptInterface

export patient_nurse_allocation, patient_nurse_block_allocation

function patient_nurse_allocation(
		beds::Array{<:Real,1},
		initial_patients::Array{<:Real,1},
		discharged_patients::Array{<:Real,2},
		admitted_patients::Array{<:Real,2},
		initial_nurses::Array{<:Real,1},
		adj_matrix::BitArray{2};
		los=11,
		nurse_days_per_patient_day::Real=2.0,
		sendreceive_switch_time::Int=0,
		min_send_amt::Real=0,
		smoothness_penalty::Real=0,
		setup_cost::Real=0,
		sent_penalty::Real=0,
		balancing_thresh_patients::Real=1.0,
		balancing_penalty_patients::Real=0,
		nurse_target_load::Real=1.25,
		nurse_target_load_gap::Real=0.25,
		nurse_load_penalty::Real=0,
		disallow_nurse_shortage_sent::Bool=false,
		disallow_nurse_shortage_newpatients::Bool=false,
		severity_weighting::Bool=false,
		no_artificial_overflow::Bool=false,
		no_artificial_shortage::Bool=false,
		verbose::Bool=false,
)
	N, T = size(admitted_patients)
	@assert(size(initial_patients, 1) == N)
	@assert(size(beds, 1) == N)
	@assert(size(initial_nurses, 1) == N)
	@assert(size(adj_matrix) == (N,N))
	@assert(size(discharged_patients) == (N, T))

	L = nothing
	if isa(los, Int)
		L = vcat(ones(Int, los), zeros(Int, T-los))
	elseif isa(los, Array{<:Real,1})
		if length(los) >= T
			L = los
		else
			L = vcat(los, zeros(Float64, T-length(los)))
		end
	elseif isa(los, Distribution)
		L = 1.0 .- cdf.(los, 0:T)
	else
		error("Invalid length of stay distribution")
	end

	model = Model(Gurobi.Optimizer)
	if !verbose set_silent(model) end

	@variable(model, sentpatients[1:N,1:N,1:T])
	@variable(model, sentnurses[1:N,1:N,1:T])

	@variable(model, obj_dummy_patients[1:N,1:T] >= 0)
	@variable(model, obj_dummy_nurses[1:N,1:T] >= 0)

	# enforce minimum transfer amount if enabled
	if min_send_amt <= 0
		@constraint(model, sentpatients .>= 0)
		@constraint(model, sentnurses .>= 0)
	else
		@constraint(model, [i=1:N,j=1:N,t=1:T], sentpatients[i,j,t] in MOI.Semicontinuous(Float64(min_send_amt), Inf))
		@constraint(model, [i=1:N,j=1:N,t=1:T], sentnurses[i,j,t] in MOI.Semicontinuous(Float64(min_send_amt), Inf))
	end

	objective = @expression(model, 100.0 * sum(obj_dummy_nurses))

	active_patients_null = [(
			initial_patients[i]
			- sum(discharged_patients[i,1:t])
			+ sum(L[t-t₁+1] * admitted_patients[i,t₁] for t₁ in 1:t)
		) for i in 1:N, t in 1:T
	]

	if severity_weighting
		max_load_null = [maximum(active_patients_null[i,:] / beds[i]) for i in 1:N]
		severity_weight = [max_load_null[i] > 1 ? 1.0 : 10.0 for i in 1:N]

		add_to_expression!(objective, dot(sum(obj_dummy_patients, dims=2), severity_weight))
	else
		add_to_expression!(objective, sum(obj_dummy_patients))
	end

	# penalize total sent if enabled
	if sent_penalty > 0
		add_to_expression!(objective, sent_penalty*sum(sentpatients))
		add_to_expression!(objective, sent_penalty*sum(sentnurses))
	end

	# penalize non-smoothness in sent patients if enabled
	if smoothness_penalty > 0
		@variable(model, smoothness_dummy[i=1:N,j=1:N,t=1:T-1] >= 0)
		@constraint(model, [t=1:T-1],  (sentpatients[:,:,t] - sentpatients[:,:,t+1]) .<= smoothness_dummy[:,:,t])
		@constraint(model, [t=1:T-1], -(sentpatients[:,:,t] - sentpatients[:,:,t+1]) .<= smoothness_dummy[:,:,t])

		add_to_expression!(objective, smoothness_penalty * sum(smoothness_dummy))
		add_to_expression!(objective, smoothness_penalty * sum(sentpatients[:,:,1]))
	end

	# add setup costs if enabled
	if setup_cost > 0
		@variable(model, setup_dummy_patients[i=1:N,j=i+1:N], Bin)
		@constraint(model, [i=1:N,j=i+1:N], [1-setup_dummy_patients[i,j], sum(sentpatients[i,j,:])+sum(sentpatients[j,i,:])] in MOI.SOS1([1.0, 1.0]))
		add_to_expression!(objective, setup_cost*sum(setup_dummy_patients))

		@variable(model, setup_dummy_nurses[i=1:N,j=i+1:N], Bin)
		@constraint(model, [i=1:N,j=i+1:N], [1-setup_dummy_nurses[i,j], sum(sentnurses[i,j,:])+sum(sentnurses[j,i,:])] in MOI.SOS1([1.0, 1.0]))
		add_to_expression!(objective, setup_cost*sum(setup_dummy_nurses))
	end

	# only send patients between connected locations
	for i = 1:N
		@constraint(model, sentnurses[i,i,:] .== 0)
		for j = 1:N
			if ~adj_matrix[i,j]
				@constraint(model, sentpatients[i,j,:] .== 0)
				# @constraint(model, sentnurses[i,j,:] .== 0)
			end
		end
	end

	# enforce a minimum time between sending and receiving
	if sendreceive_switch_time > 0
		@constraint(model, [i=1:N,t=1:T-1], [sum(sentpatients[:,i,t]), sum(sentpatients[i,:,t:min(t+sendreceive_switch_time,T)])] in MOI.SOS1([1.0, 1.0]))
		@constraint(model, [i=1:N,t=1:T-1], [sum(sentpatients[:,i,t:min(t+sendreceive_switch_time,T)]), sum(sentpatients[i,:,t])] in MOI.SOS1([1.0, 1.0]))
		@constraint(model, [i=1:N,t=1:T-1], [sum(sentnurses[:,i,t]), sum(sentnurses[i,:,t:min(t+sendreceive_switch_time,T)])] in MOI.SOS1([1.0, 1.0]))
		@constraint(model, [i=1:N,t=1:T-1], [sum(sentnurses[:,i,t:min(t+sendreceive_switch_time,T)]), sum(sentnurses[i,:,t])] in MOI.SOS1([1.0, 1.0]))
	end

	# send new patients only
	@constraint(model, [t=1:T], sum(sentpatients[:,:,t], dims=2) .<= admitted_patients[:,t])

	# expression for the number of active patients
	@expression(model, active_patients[i=1:N,t=1:T],
		initial_patients[i]
		- sum(discharged_patients[i,1:t])
		+ sum(L[t-t₁+1] * (
			admitted_patients[i,t₁]
			- sum(sentpatients[i,:,t₁])
			+ sum(sentpatients[:,i,t₁])
		) for t₁ in 1:t)
		+ sum(sentpatients[i,:,t])
	)

	# ensure the number of active patients is non-negative
	@constraint(model, [i=1:N,t=1:T], active_patients[i,t] >= 0)

	if no_artificial_overflow
		for i in 1:N, t in 1:T
			if active_patients_null[i,t] < beds[i]
				@constraint(model, active_patients[i,t] <= beds[i])
			end
		end
	end

	# load balancing for patients
	if balancing_penalty_patients > 0
		@variable(model, balancing_dummy_patients[1:N,1:T] >= 0)
		@constraint(model, [i=1:N,t=1:T], balancing_dummy_patients[i,t] >= (active_patients[i,t] / beds[i]) - balancing_thresh_patients)
		add_to_expression!(objective, balancing_penalty_patients * sum(balancing_dummy_patients))
	end

	# objective - patients
	@expression(model, patient_overflow[i=1:N,t=1:T], active_patients[i,t] - beds[i])
	@constraint(model, [i=1:N,t=1:T], obj_dummy_patients[i,t] >= patient_overflow[i,t])

	# @constraint(model, sentpatients .== 0)

	# compute active nurses
	@expression(model, active_nurses[i=1:N,t=0:T],
		initial_nurses[i]
		- sum(sentnurses[i,:,1:t])
		+ sum(sentnurses[:,i,1:t])
	)

	# compute nurse demand
	@expression(model, nurse_demand[i=1:N,t=1:T], active_patients[i,t] * nurse_days_per_patient_day)

	# sent nurses ≦ active nurses
	@constraint(model, [i=1:N,t=1:T], sum(sentnurses[i,:,t]) <= active_nurses[i,t-1])

	# active nurses ≧ 1/2 initial nurses
	@constraint(model, [i=1:N,t=1:T], active_nurses[i,t] >= 0.5 * initial_nurses[i])

	# nurses objective
	@constraint(model, [i=1:N,t=1:T], obj_dummy_nurses[i,t] >= nurse_demand[i,t] - active_nurses[i,t])

	nurse_demand_null = active_patients_null .* nurse_days_per_patient_day
	if no_artificial_shortage
		for i in 1:N, t in 1:T
			if nurse_demand_null[i,t] > initial_nurses[i]
				@constraint(model, active_nurses[i,t] >= initial_nurses[i])
			end
			if nurse_demand_null[i,t] <= initial_nurses[i]
				# @constraint(model, active_nurses[i,t] >= nurse_demand_null[i])
				@constraint(model, active_nurses[i,t] >= nurse_demand[i])
			end
		end
	end

	if disallow_nurse_shortage_sent
		# m = 1e-5
		# @variable(model, has_nurse_shortage[i=1:N,t=1:T], Bin)
		# @constraint(model, [i=1:N,t=1:T],     m*(nurse_demand[i,t] - active_nurses[i,t]) <= has_nurse_shortage[i,t])
		# @constraint(model, [i=1:N,t=1:T], 1 + m*(nurse_demand[i,t] - active_nurses[i,t]) >= has_nurse_shortage[i,t])
		# @constraint(model, [i=1:N,t=1:T], has_nurse_shortage[i,t] => {active_nurses[i,t] >= initial_nurses[i]})
		nurse_demand_null = active_patients_null .* nurse_days_per_patient_day
		for i in 1:N, t in 1:T
			if nurse_demand_null[i,t] >= initial_nurses[i]
				@constraint(model, sum(sentnurses[:,i,1:t]) >= sum(sentnurses[i,:,1:t]))
			end
		end
	end

	if disallow_nurse_shortage_newpatients
		m = 1e-5
		ts(t) = max(1,t-round(Int,mean(los))+1)
		@variable(model, has_outside_patients[i=1:N,t=1:T], Bin)
		@constraint(model, [i=1:N,t=1:T],     m*(sum(sentpatients[:,i,ts(t):t]) - sum(sentpatients[i,:,ts(t):t])) <= has_outside_patients[i,t])
		@constraint(model, [i=1:N,t=1:T], 1 + m*(sum(sentpatients[:,i,ts(t):t]) - sum(sentpatients[i,:,ts(t):t])) >= has_outside_patients[i,t])
		@constraint(model, [i=1:N,t=1:T], has_outside_patients[i,t] => {active_nurses[i,t] >= nurse_demand[i,t]})
	end

	# nurse load
	if nurse_load_penalty > 0
		@variable(model, load_dummy_nurses_abs[i=1:N,t=1:T] >= 0)
		@constraint(model, [i=1:N,t=1:T],  (active_nurses[i,t] - nurse_target_load*nurse_demand[i,t]) <= load_dummy_nurses_abs[i,t])
		@constraint(model, [i=1:N,t=1:T], -(active_nurses[i,t] - nurse_target_load*nurse_demand[i,t]) <= load_dummy_nurses_abs[i,t])

		@variable(model, load_dummy_nurses[i=1:N,t=1:T] >= 0)
		@constraint(model, [i=1:N,t=1:T], load_dummy_nurses[i,t] >= load_dummy_nurses_abs[i,t] - nurse_target_load_gap*nurse_demand[i,t])
		add_to_expression!(objective, nurse_load_penalty * sum(load_dummy_nurses))
	end

	@objective(model, Min, objective)

	optimize!(model)
	return model
end;

function patient_nurse_block_allocation(
		beds::Dict{Symbol,Array{Float32,1}},
		patient_blocks::Array,
		nurses::Array{Float32,1},
		adj_matrix::BitArray{2};
		nurse_hrs_per_week::Real=36,
		sendreceive_switch_time::Int=0,
		min_send_amt::Real=0,
		smoothness_penalty::Real=0,
		setup_cost::Real=0,
		sent_penalty::Real=0,
		balancing_thresh::Real=1.0,
		balancing_penalty::Real=0,
		nurse_target_load::Real=1.25,
		nurse_dead_zone::Real=0.25,
		nurse_load_penalty::Real=10.0,
		verbose::Bool=false,
)
	N, T = size(patient_blocks[1].admitted)
	G = length(patient_blocks)
	B = length(beds)

	bed_types = collect(keys(beds))

	nurse_hrs_per_day = nurse_hrs_per_week / 7
	nurses_days_per_day = 24 / nurse_hrs_per_day
	nurse_days_per_patient_day(g) = Float32(nurses_days_per_day / patient_blocks[g].patients_per_nurse)

	model = Model(Gurobi.Optimizer)
	if !verbose set_silent(model) end

	@variable(model, sent[1:G,1:N,1:N,1:T])
	@variable(model, obj_dummy[bed_types,1:N,1:T])

	@variable(model, sentnurses[1:N,1:N,1:T] >= 0)
	@variable(model, obj_dummy_nurses[1:N,1:T] >= 0)

	# enforce minimum transfer amount if enabled
	if min_send_amt <= 0
		@constraint(model, sent .>= 0)
		@constraint(model, sentnurses .>= 0)
	else
		@constraint(model, [i=1:N,j=1:N,t=1:T], sent[i,j,t] in MOI.Semicontinuous(Float64(min_send_amt), Inf))
		@constraint(model, [i=1:N,j=1:N,t=1:T], sentnurses[i,j,t] in MOI.Semicontinuous(Float64(min_send_amt), Inf))
	end

	objective = @expression(model, sum(obj_dummy) + sum(obj_dummy_nurses))

	# penalize total sent if enabled
	if sent_penalty > 0
		add_to_expression!(objective, sent_penalty*sum(sent))
		add_to_expression!(objective, sent_penalty*sum(sentnurses))
	end

	# penalize non-smoothness in sent patients if enabled
	if smoothness_penalty > 0
		@variable(model, smoothness_dummy[i=1:N,j=1:N,t=1:T-1] >= 0)
		@constraint(model, [t=1:T-1],  sum(sent[:,:,:,t] - sent[:,:,:,t+1], dims=1) .<= smoothness_dummy[:,:,t])
		@constraint(model, [t=1:T-1], -sum(sent[:,:,:,t] - sent[:,:,:,t+1], dims=1) .<= smoothness_dummy[:,:,t])

		add_to_expression!(objective, smoothness_penalty * sum(smoothness_dummy))
		add_to_expression!(objective, smoothness_penalty * sum(sent[:,:,:,1]))
	end

	# add setup costs if enabled
	if setup_cost > 0
		@variable(model, setup_dummy[i=1:N,j=i+1:N], Bin)
		@constraint(model, [i=1:N,j=i+1:N], [1-setup_dummy[i,j], sum(sent[:,i,j,:])+sum(sent[:,j,i,:])] in MOI.SOS1([1.0, 1.0]))
		add_to_expression!(objective, setup_cost*sum(setup_dummy))
	end

	# only send patients between connected locations
	for i = 1:N
		for j = 1:N
			if ~adj_matrix[i,j]
				@constraint(model, sent[:,i,j,:] .== 0)
				@constraint(model, sentnurses[i,j,:] .== 0)
			end
		end
	end

	# enforce a minimum time between sending and receiving
	if sendreceive_switch_time > 0
		@constraint(model, [i=1:N,t=1:T-1],
			[sum(sent[:,:,i,t]), sum(sent[:,i,:,t:min(t+sendreceive_switch_time,T)])] in MOI.SOS1([1.0, 1.0])
		)
		@constraint(model, [i=1:N,t=1:T-1],
			[sum(sent[:,:,i,t:min(t+sendreceive_switch_time,T)]), sum(sent[:,i,:,t])] in MOI.SOS1([1.0, 1.0])
		)
	end

	ts(t,g) = max(1,t-patient_blocks[g].hospitalized_days+1)

	@expression(model, transfer[g1=1:G,g2=1:G,i=1:N,t=1:T],
		((t > patient_blocks[g1].hospitalized_days) && (patient_blocks[g1].to == patient_blocks[g2].id)) ?
			patient_blocks[g1].admitted[i,t-patient_blocks[g1].hospitalized_days] : 0
	)

	# expression for the number of active patients
	@expression(model, active_patients[g=1:G,i=1:N,t=0:T],
		patient_blocks[g].initial[i]
		- sum(patient_blocks[g].discharged[i,1:t])
		+ sum(patient_blocks[g].admitted[i,ts(t,g):t])
		- sum(sent[g,i,:,ts(t,g):t])
		+ sum(sent[g,:,i,ts(t,g):t])
		+ sum(transfer[:,g,i,ts(t,g):t])
	)

	# expression for the patient overflow
	bed_groups = Dict(b => [i for (i,g) in enumerate(patient_blocks) if g.bed_type == b] for b in bed_types)
	@expression(model, overflow[b in bed_types,i=1:N,t=1:T],
		sum(active_patients[g,i,t] + sum(sent[g,i,:,t]) for g in bed_groups[b]) - beds[b][i]
	)

	# send new patients only
	@constraint(model, [g=1:G,i=1:N,t=1:T], sum(sent[g,i,:,t]) <= patient_blocks[g].admitted[i,t])

	# sent nurses <= active nurses
	@constraint(model, [i=1:N,t=1:T],
		sum(sentnurses[i,:,t]) <=
			nurses[i]
			- sum(sentnurses[i,:,1:t-1])
			+ sum(sentnurses[:,i,1:t-1])
	)

	# ensure the number of active patients is non-negative
	@constraint(model, [g=1:G,i=1:N,t=1:T], active_patients[g,i,t] >= 0)

	# load balancing
	if balancing_penalty > 0
		@variable(model, balancing_dummy[bed_types,1:N,1:T] >= 0)
		@constraint(model, [b in bed_types,i=1:N,t=1:T],
			balancing_dummy[b,i,t] >= (sum(active_patients[g,i,t] for g in bed_groups[b]) / beds[b][i]) - balancing_thresh)
		add_to_expression!(objective, balancing_penalty * sum(balancing_dummy))
	end

	# objective - patients
	@constraint(model, [b in bed_types,i=1:N,t=1:T], obj_dummy[b,i,t] >= overflow[b,i,t])

	@expression(model, active_nurses[i=1:N,t=1:T],
		nurses[i]
		- sum(sentnurses[i,:,1:t])
		+ sum(sentnurses[:,i,1:t-1])
	)
	@expression(model, total_nurse_demand[i=1:N,t=1:T],
		sum(active_patients[g,i,t] * nurse_days_per_patient_day(g) for g in 1:G)
	)
	@expression(model, nurse_demand[g=1:G,i=1:N,t=1:T],
		active_patients[g,i,t] * nurse_days_per_patient_day(g)
	)

	# objective - nurses
	@constraint(model, [i=1:N,t=1:T], obj_dummy_nurses[i,t] >= total_nurse_demand[i,t] - active_nurses[i,t])

	if nurse_load_penalty > 0
		@variable(model, nurse_load_dummy_abs[i=1:N,t=1:T] >= 0)
		@constraint(model, [i=1:N,t=1:T],  (active_nurses[i,t] - nurse_target_load*total_nurse_demand[i,t]) <= nurse_load_dummy_abs[i,t])
		@constraint(model, [i=1:N,t=1:T], -(active_nurses[i,t] - nurse_target_load*total_nurse_demand[i,t]) <= nurse_load_dummy_abs[i,t])

		@variable(model, nurse_load_dummy[i=1:N,t=1:T] >= 0)
		@constraint(model, [i=1:N,t=1:T], nurse_load_dummy[i,t] >= nurse_load_dummy_abs[i,t] - nurse_dead_zone*total_nurse_demand[i,t])
		add_to_expression!(objective, nurse_load_penalty * sum(nurse_load_dummy))
	end

	@objective(model, Min, objective)
	optimize!(model)

	return model
end;

end;

max_stars_repo_name: flixpar/covid-resource-allocation

max_issues_repo_name: flixpar/covid-resource-allocation

max_forks_repo_name: flixpar/covid-resource-allocation


text:
using Test

using problems.str

@testset "is unique 1" begin
    @test is_unique1("") == true
    @test is_unique2("fjdaslk") == true
    @test is_unique1("abcdefg")  == true
    @test is_unique1("abcda") == false
    @test is_unique1("is unique") == false
end

@testset "is unique 2" begin
    @test is_unique2("") == true
    @test is_unique2("fjdaslk") == true
    @test is_unique2("abcdefg")  == true
    @test is_unique2("abcda") == false
    @test is_unique2("is unique") == false
end

@testset "is permutation 1" begin
    @test is_permutation1("abcdef", "fedcba") == true
    @test is_permutation1("abcdef", "fedcbaa") == false
    @test is_permutation1("fabfcdef", "abcdefff")  == true
    @test is_permutation1("perm1", "1perm") == true
    @test is_permutation1("perm1", "perm2") == false
end

@testset "is permutation 2" begin
    @test is_permutation2("abcdef", "fedcba") == true
    @test is_permutation2("abcdef", "fedcbaa") == false
    @test is_permutation2("fabfcdef", "abcdefff")  == true
    @test is_permutation2("perm1", "1perm") == true
    @test is_permutation2("perm1", "perm2") == false
end

@testset "urlify 1" begin
    @test urlify1("", 0) == ""
    @test urlify1("NoSpace", 7) == "NoSpace"
    @test urlify1("Mr John Smith    ", 13) == "Mr%20John%20Smith"
    @test urlify1("Single Space  ", 12)  == "Single%20Space"
    @test urlify1("tr i pple Space      ", 15) == "tr%20i%20pple%20Space"
end

@testset "is palindrome permutation 1" begin
    @test is_pal_perm1("") == false
    @test is_pal_perm1("atco cta") == true
    @test is_pal_perm1("A nut for a jar of tuna") == true
    @test is_pal_perm1("Amore, Roma.") == true
    @test is_pal_perm1("Ed, I saw Harpo Marx ram Oprah W. aside.") == true
end

@testset "one away 1" begin
    @test one_away1("", "") == false
    @test one_away1("pale", "ple") == true
    @test one_away1("pales", "pale") == true
    @test one_away1("pale", "balep") == true
    @test one_away1("pale", "bake") == false
    @test one_away1("PALE", "baLe") == true
    @test one_away1("PALE    ", "b    aLe") == true
    @test one_away1("PALE", "PALEB") == true
end

@testset "compress 1" begin
    @test compress1("") == ""
    @test compress1("aa") == "aa"
    @test compress1("aaa") == "a3"
    @test compress1("aabcccccaaa") == "a2b1c5a3"
    @test compress1("aaa111bbb44ccccc") == "a313b342c5"
end

@testset "is rotation 1" begin
    @test is_rotation1("waterbottle", "") == false
    @test is_rotation1("waterbottle", "erbottlewat") == true
    @test is_rotation1("waterbottle", "bottlewater") == true
    @test is_rotation1("waterbottle", "ttlewaterbo") == true
    @test is_rotation1("waterbottle", "ttlewwaterbo") == false
end

max_stars_repo_name: bmoretz/Daily-Coding-Problem

max_issues_repo_name: bmoretz/Daily-Coding-Problem

max_forks_repo_name: bmoretz/Daily-Coding-Problem


text:
@testset "scomponents" begin
    A = load_matrix_network("cores_example")
    sizes = [15;1;5]
    components = 3
    @testset "MatrixNetwork" begin
        cc = scomponents(MatrixNetwork(A))
        @test cc.sizes == sizes
        @test length(cc.map) == 21
        @test cc.number == components
    end
    @testset "SparseMatrixCSC" begin
        cc = scomponents(A)
        @test cc.sizes == sizes
        @test length(cc.map) == 21
        @test cc.number == components
    end
    @testset "CSR" begin
        cc = scomponents(sparse_to_csr(A)...)
        @test cc.sizes == sizes
        @test length(cc.map) == 21
        @test cc.number == components
    end
    @testset "triplet" begin
        cc = scomponents(findnz(A)[1], findnz(A)[2])
        @test cc.sizes == sizes
        @test length(cc.map) == 21
        @test cc.number == components
    end
    @testset "strong_components_map" begin
        sci = strong_components_map(MatrixNetwork(A))
        @test strong_components_map(A) == sci
        @test strong_components_map(sparse_to_csr(A)...) == sci
        @test strong_components_map(findnz(A)[1], findnz(A)[2]) == sci
    end
    @testset "empty" begin
        @test scomponents(empty_graph(0)).sizes == []
    end
end


text:
using Signals
using Signals: state

using Test
include("push_pull.jl")
include("benchmark.jl")
include("operators.jl")
include("binding.jl")
include("time.jl")


text:
module DGMethods

using MPI
using StaticArrays
using DocStringExtensions
using KernelAbstractions
using KernelAbstractions.Extras: @unroll
using ..MPIStateArrays
using ..Mesh.Grids
using ..Mesh.Topologies
using ..Mesh.Filters
using ..VariableTemplates
using ..Courant
using ..BalanceLaws:
    BalanceLaw,
    AbstractStateType,
    Prognostic,
    Auxiliary,
    Gradient,
    GradientFlux,
    GradientLaplacian,
    Hyperdiffusive,
    UpwardIntegrals,
    DownwardIntegrals,
    vars_state,
    number_states

import ..BalanceLaws:
    BalanceLaw,
    init_state_prognostic_arr!,
    init_state_auxiliary!,
    flux_first_order_arr!,
    flux_second_order_arr!,
    compute_gradient_flux_arr!,
    compute_gradient_argument_arr!,
    source_arr!,
    transform_post_gradient_laplacian_arr!,
    wavespeed,
    boundary_conditions,
    boundary_state!,
    update_auxiliary_state!,
    update_auxiliary_state_gradient!,
    nodal_update_auxiliary_state!,
    nodal_init_state_auxiliary!,
    integral_load_auxiliary_state_arr!,
    integral_set_auxiliary_state_arr!,
    indefinite_stack_integral!,
    reverse_indefinite_stack_integral!,
    reverse_integral_load_auxiliary_state_arr!,
    reverse_integral_set_auxiliary_state_arr!

export DGModel,
    DGFVModel,
    SpaceDiscretization,
    init_ode_state,
    restart_ode_state,
    restart_auxiliary_state,
    basic_grid_info,
    init_state_auxiliary!,
    auxiliary_field_gradient!,
    courant

include("custom_filter.jl")
include("FVReconstructions.jl")
include("NumericalFluxes.jl")
include("DGModel.jl")
include("DGModel_kernels.jl")
include("DGFVModel_kernels.jl")
include("create_states.jl")

"""
    calculate_dt(dg, model, Q, Courant_number, direction, t)

For a given model, compute a time step satisying the nondiffusive Courant number
`Courant_number`
"""
function calculate_dt(dg, model, Q, Courant_number, t, direction)
    Δt = one(eltype(Q))
    CFL = courant(nondiffusive_courant, dg, model, Q, Δt, t, direction)
    return Courant_number / CFL
end

end