text:
module POMCPDPW

using POMDPs
using BasicPOMCP
using ParticleFilters
using Parameters
using MCTS
using CPUTime
using D3Trees
using Colors
using Random
using Printf
using POMDPPolicies

using BasicPOMCP: convert_estimator

import Base: insert!
import POMDPs: action, solve, mean, rand, updater, currentobs, history
import POMDPModelTools: action_info

import MCTS: n_children, next_action, isroot, node_tag, tooltip_tag

export
    POMCPOWSolver,
    POMCPOWPlanner,
    POMCPOWTree,
    POWNodeBelief,
    POWTreeObsNode,
    CategoricalVector,
    FORollout,
    FOValue,

    MaxUCB,
    MaxQ,
    MaxTries,

    init_N,
    init_V,

    n_children,
    belief,
    sr_belief,
    isroot,

    POMCPOWVisualizer,
    blink

const init_V = init_Q

include("categorical_vector.jl")
include("beliefs.jl")

include("tree.jl")
include("criteria.jl")

"""
    POMCPOWSolver

Partially observable Monte Carlo planning with observation widening.

Fields:

- `eps::Float64`:
    Rollouts and tree expansion will stop when discount^depth is less than this.
    default: `0.01`
- `max_depth::Int`:
    Rollouts and tree expension will stop when this depth is reached.
    default: `10`
- `criterion::Any`:
    Criterion to decide which action to take at each node. e.g. `MaxUCB(c)`, `MaxQ`, or `MaxTries`
    default: `MaxUCB(1.0)`
- `final_criterion::Any`:
    Criterion for choosing the action to take after the tree is constructed.
    default: `MaxQ()`
- `tree_queries::Int`:
    Number of iterations during each action() call.
    default: `1000`
- `max_time::Float64`:
    Time limit for planning at each steps (seconds).
    default: `Inf`
- `rng::AbstractRNG`:
    Random number generator.
    default: Base.GLOBAL_RNG
- `node_sr_belief_updater::Updater`:
    Updater for state-reward distribution at the nodes.
    default: `POWNodeFilter()`
- `estimate_value::Any`: (rollout policy can be specified by setting this to RolloutEstimator(policy))
    Function, object, or number used to estimate the value at the leaf nodes.
    If this is a function `f`, `f(pomdp, s, h::BeliefNode, steps)` will be called to estimate the value.
    If this is an object `o`, `estimate_value(o, pomdp, s, h::BeliefNode, steps)` will be called.
    If this is a number, the value will be set to that number
    default: `RolloutEstimator(RandomSolver(rng))`
- `enable_action_pw::Bool`:
    Controls whether progressive widening is done on actions; if `false`, the entire action space is used.
    default: `true`
- `check_repeat_obs::Bool`:
    Check if an observation was sampled multiple times. This has some dictionary maintenance overhead, but prevents multiple nodes with the same observation from being created. If the observation space is discrete, this should probably be used, but can be turned off for speed.
    default: `true`
- `check_repeat_act::Bool`:
    Check if an action was sampled multiple times. This has some dictionary maintenance overhead, but prevents multiple nodes with the same action from being created. If the action space is discrete, this should probably be used, but can be turned off for speed.
    default: `true`
- `tree_in_info::Bool`:
    If `true`, return the tree in the info dict when action_info is called, this can use a lot of memory if histories are being saved.
    default: `false`
- `k_action::Float64`, `alpha_action::Float64`, `k_observation::Float64`, `alpha_observation::Float64`:
        These constants control the double progressive widening. A new observation
        or action will be added if the number of children is less than or equal to kN^alpha.
        defaults: k: `10`, alpha: `0.5`
- `init_V::Any`:
    Function, object, or number used to set the initial V(h,a) value at a new node.
    If this is a function `f`, `f(pomdp, h, a)` will be called to set the value.
    If this is an object `o`, `init_V(o, pomdp, h, a)` will be called.
    If this is a number, V will be set to that number
    default: `0.0`
- `init_N::Any`:
    Function, object, or number used to set the initial N(s,a) value at a new node.
    If this is a function `f`, `f(pomdp, h, a)` will be called to set the value.
    If this is an object `o`, `init_N(o, pomdp, h, a)` will be called.
    If this is a number, N will be set to that number
    default: `0`
- `next_action::Any`
    Function or object used to choose the next action to be considered for progressive widening.
    The next action is determined based on the POMDP, the belief, `b`, and the current `BeliefNode`, `h`.
    If this is a function `f`, `f(pomdp, b, h)` will be called to set the value.
    If this is an object `o`, `next_action(o, pomdp, b, h)` will be called.
    default: `RandomActionGenerator(rng)`
- `default_action::Any`:
    Function, action, or Policy used to determine the action if POMCP fails with exception `ex`.
    If this is a Function `f`, `f(pomdp, belief, ex)` will be called.
    If this is a Policy `p`, `action(p, belief)` will be called.
    If it is an object `a`, `default_action(a, pomdp, belief, ex)` will be called, and
    if this method is not implemented, `a` will be returned directly.

"""
@with_kw mutable struct POMCPOWSolver{RNG<:AbstractRNG} <: AbstractPOMCPSolver
    eps::Float64                = 0.01
    max_depth::Int              = typemax(Int)
    criterion                   = MaxUCB(1.0)
    final_criterion             = MaxQ()
    tree_queries::Int           = 1000
    max_time::Float64           = Inf
    rng::RNG                    = Random.GLOBAL_RNG
    node_sr_belief_updater      = POWNodeFilter()

    estimate_value::Any         = RolloutEstimator(RandomSolver(rng))

    enable_action_pw::Bool      = true
    check_repeat_obs::Bool      = true
    check_repeat_act::Bool      = true
    tree_in_info::Bool          = false

    alpha_observation::Float64  = 0.5
    k_observation::Float64      = 10.0
    alpha_action::Float64       = 0.5
    k_action::Float64           = 10.0
    init_V::Any                 = 0.0
    init_N::Any                 = 0
    next_action::Any            = RandomActionGenerator(rng)
    default_action::Any         = ExceptionRethrow()
end

# unweighted ParticleCollections don't get anything pushed to them
function push_weighted!(::ParticleCollection, sp) end

include("planner2.jl")
include("solver2.jl")

function solve(solver::POMCPOWSolver, problem::POMDP)
    return POMCPOWPlanner(solver, problem)
end

include("updater.jl")
include("visualization.jl")

end # module


text:
## Cauchy

# pseudocauchy does not normalize at ∞
pseudostieltjes{LS,RR<:Arc}(S::Space{LS,RR},f,z,s...)=stieltjes(setcanonicaldomain(S),f,mobius(S,z),s...)
pseudohilbert{LS,RR<:Arc}(S::Space{LS,RR},f,z)=hilbert(setdomain(S,Interval()),f,mobius(S,z))

stieltjes{LS,RR<:Arc}(S::Space{LS,RR},f,z,s...)=stieltjes(setcanonicaldomain(S),f,mobius(S,z),s...)-stieltjes(setcanonicaldomain(S),f,mobius(S,Inf))
hilbert{LS,RR<:Arc}(S::Space{LS,RR},f,z)=hilbert(setcanonicaldomain(S),f,mobius(S,z))+(1/π)*stieltjes(setcanonicaldomain(S),f,mobius(S,Inf))

function PseudoHilbert{LS,RR<:Arc}(sp::JacobiWeight{LS,RR},k::Integer)
    @assert k==1
    csp=setcanonicaldomain(sp)
    H=Hilbert(csp)
    PseudoHilbertWrapper(SpaceOperator(H,setdomain(domainspace(H),domain(sp)),setdomain(rangespace(H),domain(sp))),k)
end

function Hilbert{LS,RR<:Arc}(sp::JacobiWeight{LS,RR},k::Integer)
    @assert k==1
    csp=setcanonicaldomain(sp)
    H=Hilbert(csp)+(1/π)*Stieltjes(csp,mobius(sp,Inf))
    HilbertWrapper(SpaceOperator(H,setdomain(domainspace(H),domain(sp)),setdomain(rangespace(H),domain(sp))),k)
end

## stieltjesintegral

function stieltjesintegral{LS,RR<:Arc}(sp::Space{LS,RR},w,z)
    g=Fun(w,setcanonicaldomain(sp))*fromcanonicalD(sp)
    stieltjesintegral(g,mobius(sp,z))-
        stieltjesintegral(g,mobius(sp,Inf))+
        sum(Fun(w,sp))*log(z-fromcanonical(sp,Inf))
end

function linesumstieltjesintegral{LS,RR<:Arc}(sp::Space{LS,RR},w,z)
    g=Fun(w,setcanonicaldomain(sp))*abs(fromcanonicalD(sp))
    stieltjesintegral(g,mobius(sp,z))-
        stieltjesintegral(g,mobius(sp,Inf))+
        linesum(Fun(w,sp))*log(z-fromcanonical(sp,Inf))
end

function logkernel{LS,RR<:Arc}(sp::Space{LS,RR},w,z)
    g=Fun(w,setcanonicaldomain(sp))*abs(fromcanonicalD(sp))
    logkernel(g,mobius(sp,z))-
        logkernel(g,mobius(sp,Inf))+
        linesum(Fun(w,sp))*log(abs(z-fromcanonical(sp,Inf)))/π
end

function SingularIntegral{JW,RR<:Arc}(S::JacobiWeight{JW,RR},k::Integer)
    @assert k==0
    tol=1E-15
    # the mapped logkernel
    d=domain(S)
    csp=setcanonicaldomain(S)
    Σ=SingularIntegral(csp,0)
    M=Multiplication(abs(fromcanonicalD(d,Fun(identity,csp))),csp)

    z∞=mobius(d,Inf)
    cnst=Array(Float64,0)
    for k=1:10000
        push!(cnst,logkernel(Fun([zeros(k-1);1.],csp),z∞))
        if k≥3&&norm(cnst[end-2:end])<tol
            break
        end
    end
    L∞=FiniteOperator(cnst.',csp,ConstantSpace())

    x=Fun(identity,S)
    SpaceOperator((Σ-L∞)*M,S,setdomain(rangespace(Σ),d))+(log(abs(x-fromcanonical(d,Inf)))/π)*DefiniteLineIntegral(S)
end

ext: jl


text:
mutable struct DVertex{T} <: Vertex{T}
  value::T
  inputs::Vector{DVertex{T}}
  outputs::OASet{DVertex{T}}

  DVertex{T}(x) where T = new(x, [], OASet{DVertex{T}}())
end

DVertex(x) = DVertex{typeof(x)}(x)

value(v::DVertex) = v.value
inputs(v::DVertex) = v.inputs
outputs(v::DVertex) = v.outputs

function thread!(to::DVertex, from::DVertex)
  push!(inputs(to), from)
  push!(outputs(from), to)
  return to
end

function prethread!(to::DVertex, from::DVertex)
  pushfirst!(inputs(to), from)
  push!(outputs(from), to)
  return to
end

dvertex(a...) = DVertex{Any}(a...)

dvertex(x::Vertex) = convert(DVertex{Any}, x)

dl(v::Vertex) = convert(DVertex, v)

nin(v::Vertex) = length(inputs(v))

function nout(v::Vertex)
  n = 0
  for o in outputs(v), i in inputs(o)
    i ≡ v && (n += 1)
  end
  return n
end

isfinal(v::Vertex) = nout(v) == 0

function equal(a::Vertex, b::Vertex, seen = OSet())
  (a, b) ∈ seen && return true
  (value(a) == value(b) &&
    length(inputs(a)) == length(inputs(b)) &&
    length(outputs(a)) == length(outputs(b))) || return false
  push!(seen, (a, b))
  for (i, j) in zip(inputs(a), inputs(b))
    equal(i, j, seen) || return false
  end
  each_val =  x -> Set((value(y) for y = x))
  @assert @>> a outputs each_val allunique
  @assert @>> b outputs each_val allunique
  for o in outputs(a)
    p = filter(p -> value(p) == value(o), collect(outputs(b)))
    isempty(p) && return false
    equal(o, first(p), seen) || return false
  end
  return true
end

import Base: ==

x::Vertex == y::Vertex = equal(x, y)

function mapv(f, v::Vertex; cache = ODict())
  haskey(cache, v) && return cache[v]
  node = cache[v] = typeof(v)(value(v))
  for out in outputs(v)
    push!(node.outputs, mapv(f, out, cache = cache))
  end
  for in in inputs(v)
    push!(node.inputs, mapv(f, in, cache = cache))
  end
  return f(node)
end

Base.map(f, v::DVertex) = mapv(v -> (v.value = f(v.value); v), v)

copy(v::DVertex) = mapv(identity, v)


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

#
# package typedb.client.api.answer;
#
# import javax.annotation.CheckReturnValue;
#
# public interface Numeric {
#
#     @CheckReturnValue
#     boolean isLong();
#
#     @CheckReturnValue
#     boolean isDouble();
#
#     @CheckReturnValue
#     boolean isNaN();
#
#     @CheckReturnValue
#     long asLong();
#
#     @CheckReturnValue
#     Double asDouble();
#
#     @CheckReturnValue
#     Number asNumber();
# }


text:
using FilesystemDatastructures, Test, Random

@testset "SizeConstrainedFileCache" begin
    include("scfc_tests.jl")
end

@testset "NFileCache" begin
    include("nfc_tests.jl")
end

@testset "Common FileCache tests" begin
    include("common_fc_tests.jl")
end

@testset "DiskUtils" begin
    include("diskutils_tests.jl")
end


text:
module GeoDatasets

using Printf

using RemoteFiles
using CodecZlib
using Shapefile
using GeoInterface
using ZipFile

"""
    lon,lat,data = GeoDatasets.landseamask(;resolution='l',grid=5)

Load the land-sea-lake raster from basemap: 0 is ocean, 1 is land and 2 is lake.
`grid` is the resolution in arc minutes and should be either 1.25, 2.5, 5 or 10.
The parameter `resolution` should be either `'c'`,`'l'`,`'i'`,`'h'` or `'f'`
(standing for crude, low, intermediate, high and full resolution)

The data is downloaded from [basemap](https://github.com/matplotlib/basemap/tree/master/lib/mpl_toolkits/basemap_data).

The data is originally from [GMT](http://gmt.soest.hawaii.edu/) and distributed under the terms of the
[GPL licences](https://github.com/matplotlib/basemap/blob/a551c2361314670dc7d95143190d3d48ba750d20/LICENSE_data).

```julia
using GeoDatasets
lon,lat,data = GeoDatasets.landseamask(;resolution='c',grid=5)
```

"""
function landseamask(;resolution='l',grid=5)
    if !(grid ∈ Set([1.25,2.5,5.,10.]))
        error("grid should be either 1.25, 2.5, 5 or 10")
    end
    if !(resolution ∈ Set(['c','l','i','h','f',"c","l","i","h","f"]))
        error("resolution should be either 'c','l','i','h' or 'f'")
    end

    url = "https://raw.githubusercontent.com/matplotlib/basemap/master/lib/mpl_toolkits/basemap_data/lsmask_$(grid)min_$(resolution).bin"
    @debug "url for landseamask: $url"

    # download remote file if it is not yet available
    @RemoteFile(lsmask, url)
    download(lsmask)

    # resolution in arc seconds
    res = round(Int,grid * 60)

    lon = (-180*3600 + res/2:res:180*3600 - res/2) / 3600
    lat = (-90*3600 + res/2:res:90*3600 - res/2) / 3600

    data = zeros(UInt8,length(lon)*length(lat));
    read!(GzipDecompressorStream(open(path(lsmask))),data);

    return lon,lat,reshape(data,length(lon),length(lat))
end

include("gshhg.jl")

end


text:
module Day3

function read_input(filepath)
    # Parse the file into a 1000x12 matrix of Ints
    lines = readlines(filepath)
    n_rows = length(lines)
    n_cols = length(lines[1])
    A = zeros(Int64, (n_rows, n_cols))
    for row = 1:n_rows
        for col = 1:n_cols
            A[row,col] = parse(Int64, lines[row][col])
        end
    end

    return A
end

function calculate_gamma_and_epsilon_rates(diagnostics::Matrix{Int64})
    column_sums = sum(diagnostics, dims=1)
    n_rows = size(diagnostics)[1]
    gamma_rate_binary = string("0b", join(Int(i >= n_rows/2) for i in column_sums))
    epsilon_rate_binary = string("0b", join(Int(i < n_rows/2) for i in column_sums))

    gamma_rate_decimal = parse(Int, gamma_rate_binary)
    epsilon_rate_decimal = parse(Int, epsilon_rate_binary)

    return (gamma_rate_decimal, epsilon_rate_decimal)
end

function calculate_oxygen_and_co2_ratings(diagnostics::Matrix{Int64})
    # To calculate the oxygen rating:
    # Start with the 1st column, determining the most common number in the column (tie goes to 1)
    # Discard all rows whose bit in this column are not that number.
    # Do we have only 1 row left? If not, move to the next column and repeat until we have only 1 row remaining.
    # if we cycle through all columns, start again

    n_cols = size(diagnostics)[2]

    # Calculate oxygen rating
    oxygen_rating_binary = begin
        candidates = diagnostics
        n_iter = 0
        while size(candidates)[1] > 1
            current_column = (n_iter % n_cols) + 1 # so that we can cycle through the columns indefinitely 1,2,...,n_cols,1,2,...

            n_rows = size(candidates)[1] # current number of rows remaining
            n_ones_in_column = sum(candidates[:, current_column])
            most_common_value_in_column = (n_ones_in_column / n_rows) >= 0.5 ? 1 : 0

            # Filter candidates matrix further to only include rows where current column contains most common bit
            candidates = candidates[candidates[:, current_column] .== most_common_value_in_column, :]

            n_iter += 1
        end
        join(candidates[1,:]) # collapse the lone remaining row into a string
    end

    # Calculate co2 rating
    co2_rating_binary = begin
        candidates = diagnostics
        n_iter = 0
        while size(candidates)[1] > 1
            current_column = (n_iter % n_cols) + 1 # so that we can cycle through the columns indefinitely 1,2,...,n_cols,1,2,...

            n_rows = size(candidates)[1] # current number of rows remaining
            n_ones_in_column = sum(candidates[:, current_column])
            least_common_value_in_column = (n_ones_in_column / n_rows) < 0.5 ? 1 : 0

            # Filter candidates matrix further to only include rows where current column contains most common bit
            candidates = candidates[candidates[:, current_column] .== least_common_value_in_column, :]

            n_iter += 1
        end
        join(candidates[1,:]) # collapse the lone remaining row into a string
    end

    oxygen_rating_decimal = parse(Int, string("0b", oxygen_rating_binary))
    co2_rating_decimal = parse(Int, string("0b", co2_rating_binary))
    return (oxygen_rating_decimal, co2_rating_decimal)
end

function part1_answer()
    diagnostics = read_input(joinpath(@__DIR__, "input.txt"))
    (γ, ϵ) = calculate_gamma_and_epsilon_rates(diagnostics)
    return γ * ϵ
end

function part2_answer()
    diagnostics = read_input(joinpath(@__DIR__, "input.txt"))
    (o2, co2) = calculate_oxygen_and_co2_ratings(diagnostics)
    return o2 * co2
end

end


text:
export read_dxf_file

include("line_numbering_io.jl")

struct DXFReadError <: Exception
    line::Int
    e::Exception
end

function Base.showerror(io::IO, e::DXFReadError)
    print(io, "Error while reading DXF at line $(e.line): $(e.e)")
end

function read_dxf_file(in::LineNumberingIO)
    groups = Vector{DXFGroup}()
    while true
        try
            l = readline(in)
            if l == ""
                break
            end
            code = parse(DXFGroupCode, strip(l))
            group_type = lookup_group_code(code)
            push!(groups,
                  # try
                  # read_value(group_type, in)
                  Base.invokelatest(read_value, group_type, in)
                  #=
                  catch e
                  if (e isa Base.MethodError &&
                      e.f == read_value &&
                      e.args[1] == group_type)
                  Base.invokelatest(read_value, group_type, in)
                  else
                  rethrow(e)
                  end
                  end=# )
        catch e
            throw(DXFReadError(in.line, e))
        end
    end
    groups
end

function read_dxf_file(path::String)
    open(io -> read_dxf_file(LineNumberingIO(io)),
         path, "r")
end


text:
using Base: KeySet
#=
      keys(::Dict{K,V})::KeySet{Symbol,Dict{Symbol,Any}}   # KeySet for a  Dict{Symbol,Any}
    values(::Dict{K,V})::ValueIterator{Dict{Symbol,Any}}   # ValueIterator for a  Dict{Symbol,Any}
=#

(getkeys(dict::Dict{K,V})::Array{Symbol,1}) where {K,V} = dict.keys[getslotidxs(dict)]

(getvalues(dict::Dict{K,V})::Array{Any,1}) where {K,V}  = dict.vals[getslotidxs(dict)]

(getslotidxs(dict::Dict{K,V})::Array{Int, 1}) where {K,V} = filter(i->isone(dict.slots[i]), 1:length(dict.slots))


text:
include("test_basics.jl")
include("test_approx.jl")
include("test_lininterp.jl")
include("test_FSpaceXD.jl")


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct Junction <: Topology
        name::String
        elevation::Float64
        head::Union{Nothing, Float64}
        coordinates::Union{Nothing, NamedTuple{(:lat, :lon),Tuple{Float64,Float64}}}
        pattern_name::Union{Nothing, String}
        _forecasts::InfrastructureSystems.Forecasts
        internal::InfrastructureSystemsInternal
    end

A water-system Junction.

# Arguments
- `name::String`: the name of the junction
- `elevation::Float64`: elevation of junction
- `head::Union{Nothing, Float64}`: pressure head at junction
- `coordinates::Union{Nothing, NamedTuple{(:lat, :lon),Tuple{Float64,Float64}}}`: latitude and longitude coordinates of junction
- `pattern_name::Union{Nothing, String}`: name of head pattern array, for reservoirs only
- `_forecasts::InfrastructureSystems.Forecasts`
- `internal::InfrastructureSystemsInternal`: internal reference, do not modify
"""
mutable struct Junction <: Topology
    "the name of the junction"
    name::String
    "elevation of junction"
    elevation::Float64
    "pressure head at junction"
    head::Union{Nothing, Float64}
    "latitude and longitude coordinates of junction"
    coordinates::Union{Nothing, NamedTuple{(:lat, :lon),Tuple{Float64,Float64}}}
    "name of head pattern array, for reservoirs only"
    pattern_name::Union{Nothing, String}
    _forecasts::InfrastructureSystems.Forecasts
    "internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function Junction(name, elevation, head, coordinates, pattern_name, _forecasts=InfrastructureSystems.Forecasts(), )
    Junction(name, elevation, head, coordinates, pattern_name, _forecasts, InfrastructureSystemsInternal(), )
end

function Junction(; name, elevation, head, coordinates, pattern_name, _forecasts=InfrastructureSystems.Forecasts(), )
    Junction(name, elevation, head, coordinates, pattern_name, _forecasts, )
end

# Constructor for demo purposes; non-functional.
function Junction(::Nothing)
    Junction(;
        name="init",
        elevation=0.0,
        head=nothing,
        coordinates=nothing,
        pattern_name=nothing,
        _forecasts=InfrastructureSystems.Forecasts(),
    )
end

"""Get Junction name."""
get_name(value::Junction) = value.name
"""Get Junction elevation."""
get_elevation(value::Junction) = value.elevation
"""Get Junction head."""
get_head(value::Junction) = value.head
"""Get Junction coordinates."""
get_coordinates(value::Junction) = value.coordinates
"""Get Junction pattern_name."""
get_pattern_name(value::Junction) = value.pattern_name
"""Get Junction _forecasts."""
get__forecasts(value::Junction) = value._forecasts
"""Get Junction internal."""
get_internal(value::Junction) = value.internal


text:
# interpolate_monotone.jl
# Ronald L. Rivest
# January 26, 2019

# Use monotone cubic splines to interpolate a curve fitting given data.

"""
    interpolate_monotone_curve(x, y)

    Input: x -- list of n distinct x-coords,
                e.g. [x_1, ..., x_n] = [2020, 2050, 2100]
                assumed to be in increasing order
           y -- list of n target values,
                e.g. [y_1, ..., y_n] = [1.1, 1.7, 2.0]

    Output: curve -- function (piecewise cubic spline) mapping from x to y

    Discussion:
           The desired curve is monotonic cubic interpolation
               https://en.wikipedia.org/wiki/Monotone_cubic_interpolation
           The "talk" section re correctness of step 5 is also relevant & necessary.

    Refs:
    https://en.wikipedia.org/wiki/Monotone_cubic_interpolation
    http://bl.ocks.org/niclasmattsson/7bceb05fba6c71c78d507adae3d29417

    Usage example:
	f = interpolate_monotone_curve([0, 1, 2, 3, 4], [0, 1, 4, 9, 16]);
	for x in 0:0.5:4.0
	    val = f(x)
            println(x, " squared is about ", val)
        end
"""

# Define basic cubic Hermite splines

h00(t) = (1+2t)*(1-t)^2
h10(t) = t*(1-t)^2
h01(t) = t^2 * (3-2t)
h11(t) = t^2 * (t-1)

function interpolate_monotone_curve(x, y)

    # deal with length issues
    @assert length(x)==length(y)
    n = length(x)
    if n == 0
        return function (x); return 0; end
    elseif n == 1
        return function (x); return y[1]; end
    end

    # 0. sort x and y similarly so x is in nondecreasing order
    p = sortperm(x)
    x = x[p]
    y = y[p]

    # 1. get consecutive differences and slopes
    dx = diff(x)
    dy = diff(y)
    delta = dy ./ dx      # note: assumes xs are distinct so no zerodiv!

    # 2. Compute m

    m = fill(0.0, n)      # slopes
    for k in 2:n-1
        if delta[k-1]*delta[k] <= 0.0
            m[k] = 0
        else
            m[k] = (delta[k-1]+delta[k])/2.0
        end
    end
    m[1] = delta[1]
    m[n] = delta[n-1]

    # 3. Check for flat spline segment
    alpha = fill(0.0, n)
    beta = fill(0.0, n)
    for k in 1:n-1
        if delta[k] == 0
            m[k] = 0
            m[k+1] = 0
        else
            # 4. local max or mins in data
            alpha[k] = m[k]/delta[k]
            beta[k] = m[k+1]/delta[k]
            if alpha[k] < 0
                m[k] = 0
            end
            if beta[k] < 0
                m[k+1] = 0
            end
            # 5. Prevent overshoot and ensure monotonicity
            # This is buggy, but it is what the wikipedia article seems to say was OK
            # if alpha[k] > 3 || beta[k] > 3
            #     m[k] = 3 * delta[k]
            # end
            # Here is fix, from "talk" post by Berland
            if alpha[k]^2 + beta[k]^2 > 9
                tau = 3 * (alpha[k]^2 + beta[k]^2)^(-0.5)
                m[k] = tau * alpha[k] * delta[k]
                m[k+1] = tau * beta[k] * delta[k]
            end
        end
    end

    # return interpolant function -- what y is predicted for input u?
    return function (u)
        # Give correct result if u too small or too large
        if u>=x[end]
            return y[end]
        elseif u<=x[1]
            return y[1]
        end
        # Search for the interval [x_k, x_{k+1}] that u is in
        low = 1
        high = n
        while low <= high
            mid = Int64(floor((low+high)/2))
            if x[mid] < u;   low = mid + 1  end
            if x[mid] > u;   high = mid - 1 end
            if x[mid] == u;  return y[mid] end
        end
        k = max(1, high)

        # Interpolate
        t = (u-x[k]) / dx[k]
        return y[k]*h00(t) + dx[k]*m[k]*h10(t) +
               y[k+1]*h01(t) + dx[k]*m[k+1]*h11(t)
    end # of function f(u)
end

################################################################

function test_interpolate_monotone_curve()
    f = interpolate_monotone_curve([1,2,3,4], [1,4,9,9])
    for u in 0:0.5:5.0
        val = f(u)
        println("Square of ", u, " is about ", val)
    end
end

# test_interpolate_monotone_curve()

################################################################
## The following illustrates non-monotonicity bug if
## step 5 is not implemented correctly.
## The correct implementation is from "talk" section of Wikipedia.

function demo_bug()
    xs = [0, 0.30, 0.5]
    ys = [0, 0.05, 0.5]
    f = interpolate_monotone_curve(xs, ys)
    ### BUGGY -- THIS f SHOULD BE MONOTONE BUT ISN'T!
    for i in 0.0:0.01:0.5
        println(i, " ",f(i))
    end
end

# demo_bug()

max_stars_repo_name: ron-rivest/ClimateControlTheory

max_issues_repo_name: ron-rivest/ClimateControlTheory

max_forks_repo_name: ron-rivest/ClimateControlTheory


text:
using DataStructures

type ISPCFile
    global_names::Set{ASCIIString}
    func_ids::Vector{Int64}
    compile_opts::Cmd
    compiled::Bool
    error::Bool
end
ISPCFile(opts::Cmd) = ISPCFile(Set{ASCIIString}(), Int64[], opts, false, false)

type ISPCFunction
    idx::Integer

    arg_names::Vector{Any}
    var_types::Dict{Any,Any}
    modified::Set
    ast::Expr

    ispc_name::ASCIIString
    ispc_code::ASCIIString

    file::ISPCFile
end

const ispc_fptr = Array(Ptr{Void},0)
const ispc_funcs = Array(ISPCFunction, 0)
const ispc_files = Dict{AbstractString, ISPCFile}()

const ispc_fragments = Dict()
const ispc_fragment_opts = Dict()

function run_typeinf(linfo, argtypes)
    linfo.ast = get_ast(linfo)
    arg_names = linfo.ast.args[1]
    captured = linfo.ast.args[2][2]
    captured_names = [varinfo[1] for varinfo in captured]
    kernel_argnames = Any[arg_names..., captured_names...]

    # Captured vars become arguments.
    # Unset the flags that applied to closure vars and do not
    # apply any more now that they're function arguments:
    # - 0x1 ('captured') -- does not seem to matter
    # - 0x4 ('assigned by inner function') -- would break type inference
    arg_varinfos = [Any[name, typ, flags & ~0x5]
                    for (name, typ, flags) in captured]
    local_vars = Any[linfo.ast.args[2][1]..., arg_varinfos...]

    # Transform captured vars into parameters:
    new_linfo = deepcopy(linfo)
    new_linfo.ast = get_ast(new_linfo)
    new_linfo.ast.args[1] = kernel_argnames
    new_linfo.ast.args[2][1] = local_vars
    new_linfo.ast.args[2][2] = Any[]

    typeinf(new_linfo, argtypes)
end

@generated function kernel_call{id}(::Type{Val{id}}, args...)
    println("Generating kernel call $id for argument types $args")

    opts = eval(ispc_fragment_opts[id])

    # Collect the modified arguments *before* type inference,
    # because we'll need to remove that information in order
    # for type inference to run correctly:
    linfo = ispc_fragments[id]
    modified = Set()
    for (name, typ, flags) in get_ast(linfo).args[2][2]
        if (flags & 0x4) != 0
            push!(modified, name)
        end
    end

    # println("Running type inference...")
    typed_ast = run_typeinf(linfo, args)
    # print_lambda(strip_lineno(typed_ast))

    # Register the new kernel:
    argexprs = [:(args[$i]) for i in 1:length(args)]
    func_call, func_idx = register_ispc_fragment!(id, argexprs,
                                            modified, typed_ast, opts)

    # Generate the main body:
    body = quote
        $(Expr(:meta, :inline))
#        @inbounds begin
            # First make sure that the function has been compiled:
            if UInt(ISPC.ispc_fptr[$(func_idx)]) == 0
                ISPC.compile_all()
                if UInt(ISPC.ispc_fptr[$(func_idx)]) == 0
                    error("Could not compile ISPC kernel $id")
                end
            end
            # Lazy compilation allows us to compile several fragments
            # with one call to the compiler. This is much faster.

            # Now ccall the function:
            $func_call
#        end
    end
    # println("@generated function kernel_call{$id}(::Type{Val{$id}}, args...)")
    # println(strip_lineno(body))
    body
end

function gen_code(funcs...)
    io = IOBuffer()
    write(io, ispc_includes)
    write(io, "\n\n")
    for func in funcs
        write(io, func.ispc_code)
        write(io, "\n")
    end
    takebuf_string(io)
end

function compile!(file)
    if file.compiled
        return
    end

    # Generate the code for all the functions in the file:
    funcs = [ispc_funcs[idx] for idx in file.func_ids]
    ispc_code = gen_code(funcs...)

#    println("Compiling ISPC file, compile options: ", file.compile_opts)
#    println(ispc_code)

    file.compiled = true
    try
        lib = load_ispc(ispc_code, file.compile_opts)
        # Now set all the function pointers:
        for func in funcs
            fptr = Libdl.dlsym(lib, func.ispc_name)
            ispc_fptr[func.idx] = fptr
            println("Loaded function $(func.ispc_name) at $fptr")
        end
        file.error = false
    catch e
        Base.display_error(e,catch_backtrace())
        file.error = true
    end
end

function compile_all()
    has_error = false

    # Compile all the open files:
    for opts in keys(ispc_files)
        file = ispc_files[opts]
        compile!(file)
        delete!(ispc_files, opts)
        has_error |= file.error
    end

    if has_error
        error("Some files did not compile successfully")
    end
end

function register_ispc_fragment!(id, call_args, modified, ast::Expr, opts::Cmd)

    ast_args = ast.args[1]
    ast_body = ast.args[3]

    # Collate all type information into a single dict.
    var_types = Dict()

    # Local variables:
    local_vars = ast.args[2][1]
    for (name, typ, flags) in local_vars
        var_types[name] = typ
    end

    # SSA vars:
    ssa_vars = ast.args[2][3]
    for (i, typ) in enumerate(ssa_vars)
        var_types[GenSym(i-1)] = typ # remember GenSym() starts at 0
    end

    arg_types = [] # ccall argument types
    arg_values = [] # ccall argument values
    arg_names = [] # function argument names

    pass_by_ref = []

    for (func_arg, call_arg) in zip(ast_args, call_args)
        T = var_types[func_arg]
        E = eltype(T)
        if !(isbits(E))
            error("Unsupported element type for $func_arg: $E")
        end
        if issubtype(T, DenseArray)
            if func_arg in modified
                error("Can't reassign arrays inside kernels: $func_arg")
            end

            # We have an array type, pass it as a ref:
            push!(arg_types, Ref{E})
            push!(arg_values, call_arg)
            argument = Any[func_arg]

            # Also pass the array size(s):
            for d in 1:ndims(T)
                push!(arg_types, Int64)
                push!(arg_values, :(size($call_arg, $d)))
                size_var = Symbol("$(func_arg)__len__$(d)")
                var_types[size_var] = Int64
                push!(argument, size_var)
            end
            push!(arg_names, argument)

        elseif E == T
            # Atomic isbits type:
            if func_arg in modified
                # pass-by-reference
                ref = gensym("ref")
                push!(pass_by_ref, (ref, T, call_arg))
                push!(arg_types, Ref{T})
                push!(arg_values, ref)
            else
                # pass-by-value
                push!(arg_types, T)
                push!(arg_values, call_arg)
            end
            push!(arg_names, func_arg)

        else
            error("Unsupported type for $func_arg: $T")
        end
    end

    # Get an ISPCFile to attach this function to.
    file = get(ispc_files, string(opts), ISPCFile(opts))
    if file.compiled
        # That one was already compiled, start a new one:
        file = ISPCFile(opts)
    end
    ispc_files[string(opts)] = file # use string() because of issue #14030

    # Create the function object:
    func_idx = length(ispc_funcs)+1
    func = ISPCFunction(func_idx, arg_names, var_types, modified,
                        ast_body, "", "", file)

    # Generate the kernel's ISPC code:
    ispc_codegen!(func, "kernel_$id")
    @assert func.ispc_name != ""
    @assert func.ispc_code != ""

    # Wrap all pass-by-reference arguments in Refs before
    # the call:
    wrap_refs = [:($ref = Ref{$T}($val))
                   for (ref, T, val) in pass_by_ref]

    # ... and return the modified value afterwards:
    unwrap_refs = [:($ref.x)
                   for (ref, T, val) in pass_by_ref]
    if isempty(unwrap_refs)
        ret = :(return)
    else
        ret = :(return $((unwrap_refs...,)))
    end

    # Ready! Let's put it all together:
    func_call = quote
        $(wrap_refs...)
        ccall(ISPC.ispc_fptr[$func_idx],
            Void, ($(arg_types...),), $(arg_values...))
        $ret
    end

    push!(ispc_funcs, func)
    push!(ispc_fptr, Ptr{Void}(0))
    push!(file.func_ids, func_idx)

    func_call, func_idx
end


text:
ages = parse.(Int, split(read("input.txt", String), ","))
# 🐟 = [count(==(i), ages) for i in 0:8] # this is costly
🐟 = zeros(Int, 9)
for i in eachindex(🐟); 🐟[i] = count(==(i-1), ages) end
🐡 = copy(🐟)

function evolve!(🐟, 🐡)
    copyto!(🐡, 🐟) # copy state of population
    circshift!(🐟, 🐡, -1) # reduce age by one, add newborn fish to age = 8
    🐟[7] += 🐡[1] # copy parent fish to age = 6
end

# Part 1
for t in 1:80; evolve!(🐟, 🐡); end
@info "Number of lanternfish after  80 days = $(sum(🐟))"

# Part 2
for t in 81:256; evolve!(🐟, 🐡); end
@info "Number of lanternfish after 256 days = $(sum(🐟))"


text:
# Example that performes a Linear 2-D Wavelet Approximation ,

#import libraries
using Shearlab
using PyPlot
using Images

n = 1024;
name = "./../../data_samples/boat.bmp";
f = load_image(name, n);
f = rescale(sum(f,3));
f = f[:,:,1];

Jmax = round(Int64,log2(n))-1;
Jmin = 1;

# Initialize the approximation.
eta = 4;
fWLin = zeros(n,n);
fWLin[1:round(Int64, n/eta),1:round(Int64, n/eta)] = fW[1:round(Int64, n/eta),1:round(Int64, n/eta)];

# Perform and display the approximation
min = 1;
# forward transform
fW = perform_wavortho_transf(f,Jmin,+1,h);
# linear approximation
eta = 4;
fWLin = zeros(n,n);
fWLin[1:round(Int64, n/eta),1:n/round(Int64, n/eta)] = fW[1:n/round(Int64, n/eta),1:round(Int64, n/eta)];
# backward transform
fLin = perform_wavortho_transf(fWLin,Jmin,-1,h);
elin = snr(f,fLin);
# display
clf;
imageplot(f, "Original", 1,2,1);
u = @sprintf("Linear, SNR=%.3f", elin);
imageplot(FWT.clamp(fLin), u, 1,2,2);


text:
using FDM: jvp, j′vp

const _fdm = central_fdm(5, 1)

"""
    frule_test(f, (x, ẋ)...; rtol=1e-9, atol=1e-9, fdm=central_fdm(5, 1), kwargs...)

# Arguments
- `f`: Function for which the `frule` should be tested.
- `x`: input at which to evaluate `f` (should generally be set randomly).
- `ẋ`: differential w.r.t. `x` (should generally be set randomly).

All keyword arguments except for `fdm` are passed to `isapprox`.
"""
function frule_test(f, (x, ẋ); rtol=1e-9, atol=1e-9, fdm=_fdm, kwargs...)
    return frule_test(f, ((x, ẋ),); rtol=rtol, atol=atol, fdm=fdm, kwargs...)
end

function frule_test(f, xẋs::Tuple{Any, Any}...; rtol=1e-9, atol=1e-9, fdm=_fdm, kwargs...)
    xs, ẋs = collect(zip(xẋs...))
    Ω, dΩ_rule = ChainRules.frule(f, xs...)
    @test f(xs...) == Ω

    dΩ_ad, dΩ_fd = dΩ_rule(ẋs...), jvp(fdm, xs->f(xs...), (xs, ẋs))
    @test isapprox(dΩ_ad, dΩ_fd; rtol=rtol, atol=atol, kwargs...)
end

"""
    rrule_test(f, ȳ, (x, x̄)...; rtol=1e-9, atol=1e-9, fdm=central_fdm(5, 1), kwargs...)

# Arguments
- `f`: Function to which rule should be applied.
- `ȳ`: adjoint w.r.t. output of `f` (should generally be set randomly).
- `x`: input at which to evaluate `f` (should generally be set randomly).
- `x̄`: currently accumulated adjoint (should generally be set randomly).

All keyword arguments except for `fdm` are passed to `isapprox`.
"""
function rrule_test(f, ȳ, (x, x̄)::Tuple{Any, Any}; rtol=1e-9, atol=1e-9, fdm=_fdm, kwargs...)
    # Check correctness of evaluation.
    fx, dx = ChainRules.rrule(f, x)
    @test fx ≈ f(x)

    # Correctness testing via finite differencing.
    x̄_ad, x̄_fd = dx(ȳ), j′vp(fdm, f, ȳ, x)
    @test isapprox(x̄_ad, x̄_fd; rtol=rtol, atol=atol, kwargs...)

    # Assuming x̄_ad to be correct, check that other ChainRules mechanisms are correct.
    test_accumulation(x̄, dx, ȳ, x̄_ad)
    test_accumulation(Zero(), dx, ȳ, x̄_ad)
end

function rrule_test(f, ȳ, xx̄s::Tuple{Any, Any}...; rtol=1e-9, atol=1e-9, fdm=_fdm, kwargs...)
    # Check correctness of evaluation.
    xs, x̄s = collect(zip(xx̄s...))
    Ω, Δx_rules = ChainRules.rrule(f, xs...)
    @test f(xs...) == Ω

    # Correctness testing via finite differencing.
    Δxs_ad = map(Δx_rule->Δx_rule(ȳ), Δx_rules)
    Δxs_fd = j′vp(fdm, f, ȳ, xs...)
    for (Δx_ad, Δx_fd) in zip(Δxs_ad, Δxs_fd)
        @test isapprox(Δx_ad, Δx_fd; rtol=rtol, atol=atol, kwargs...)
    end

    # Assuming the above to be correct, check that other ChainRules mechanisms are correct.
    for (x̄, Δx_rule, Δx_ad) in zip(x̄s, Δx_rules, Δxs_ad)
        test_accumulation(x̄, Δx_rule, ȳ, Δx_ad)
        test_accumulation(Zero(), Δx_rule, ȳ, Δx_ad)
    end
end

function Base.isapprox(ad::Wirtinger, fd; kwargs...)
    error("Finite differencing with Wirtinger rules not implemented")
end
function Base.isapprox(d_ad::Casted, d_fd; kwargs...)
    return all(isapprox.(extern(d_ad), d_fd; kwargs...))
end
function Base.isapprox(d_ad::DNE, d_fd; kwargs...)
    error("Tried to differentiate w.r.t. a DNE")
end
function Base.isapprox(d_ad::Thunk, d_fd; kwargs...)
    return isapprox(extern(d_ad), d_fd; kwargs...)
end

function test_accumulation(x̄, dx, ȳ, partial)
    @test all(extern(ChainRules.add(x̄, partial)) .== extern(x̄) .+ extern(partial))
    test_accumulate(x̄, dx, ȳ, partial)
    test_accumulate!(x̄, dx, ȳ, partial)
    test_store!(x̄, dx, ȳ, partial)
    return nothing
end

function test_accumulate(x̄::Zero, dx, ȳ, partial)
    @test extern(accumulate(x̄, dx, ȳ)) == extern(partial)
    return nothing
end

function test_accumulate(x̄::Number, dx, ȳ, partial)
    @test extern(accumulate(x̄, dx, ȳ)) == extern(x̄) + extern(partial)
    return nothing
end

function test_accumulate(x̄::AbstractArray, dx, ȳ, partial)
    x̄_old = copy(x̄)
    @test all(extern(accumulate(x̄, dx, ȳ)) .== (extern(x̄) .+ extern(partial)))
    @test x̄ == x̄_old
    return nothing
end

test_accumulate!(x̄::Zero, dx, ȳ, partial) = nothing

function test_accumulate!(x̄::Number, dx, ȳ, partial)
    @test accumulate!(x̄, dx, ȳ) == accumulate(x̄, dx, ȳ)
    return nothing
end

function test_accumulate!(x̄::AbstractArray, dx, ȳ, partial)
    x̄_copy = copy(x̄)
    accumulate!(x̄_copy, dx, ȳ)
    @test extern(x̄_copy) == (extern(x̄) .+ extern(partial))
    return nothing
end

test_store!(x̄::Zero, dx, ȳ, partial) = nothing
test_store!(x̄::Number, dx, ȳ, partial) = nothing

function test_store!(x̄::AbstractArray, dx, ȳ, partial)
    x̄_copy = copy(x̄)
    store!(x̄_copy, dx, ȳ)
    @test all(x̄_copy .== extern(partial))
    return nothing
end


text:
# To figure out what should be precompiled, run this script, then move
# precompile_Plots.jl in precompiles_path (see below) to src/precompile.jl

# This script works by using SnoopCompile to log compilations that take place
# while running the examples on the GR backend. So SnoopCompile must be
# installed, and StatsPlots, RDatasets, and FileIO are also required for
# certain examples.

# If precompilation fails with an UndefVarError for a module, probably what is
# happening is that the module appears in the precompile statements, but is
# only visible to one of Plots' dependencies, and not Plots itself. Adding the
# module to the blacklist below will remove these precompile statements.

# Anonymous functions may appear in precompile statements as functions with
# hashes in their name. Those of the form "#something##kw" have to do with
# compiling functions with keyword arguments, and are named reproducibly, so
# can be kept. Others generally will not work. Currently, SnoopCompile includes
# some anonymous functions that not reproducible, but SnoopCompile PR #30
# (which looks about to be merged) will ensure that anonymous functions are
# actually defined before attempting to precompile them. Alternatively, we can
# keep only the keyword argument related anonymous functions by changing the
# regex that SnoopCompile uses to detect anonymous functions to
# r"#{1,2}[^\"#]+#{1,2}\d+" (see anonrex in SnoopCompile.jl). To exclude all
# precompile statements involving anonymous functions, "#" can also be added to
# the blacklist below.

using SnoopCompile

log_path = joinpath(tempdir(), "compiles.log")
precompiles_path = joinpath(tempdir(), "precompile")

# run examples with GR backend, logging what needs to be compiled
SnoopCompile.@snoopc log_path begin
    using Plots
    Plots.test_examples(:gr, disp=true)
end

# precompile calls containing the following strings are dropped
blacklist = [
    # functions defined in examples
    "PlotExampleModule",
    # the following are not visible to Plots, only its dependencies
    "CategoricalArrays",
    "FixedPointNumbers",
    "SparseArrays",
    r"#{1,2}[^\"#]+#{1,2}\d+",
]

data = SnoopCompile.read(log_path)
pc = SnoopCompile.parcel(reverse!(data[2]), blacklist=blacklist)
SnoopCompile.write(precompiles_path, pc)


text:
 using Agents: vertices
 using MetaGraphs: add_vertex!, has_edge
 using Agents, AgentsPlots, Plots, LightGraphs, MetaGraphs, GraphPlot, GraphRecipes, NetworkLayout
 using Tables, DataFrames
 using CSV, JSON, Serialization, DelimitedFiles
 using BritishNationalGrid
 using ZipFile, Shapefile
 using Random
 using Match
 using LinearAlgebra
 using StatsBase
 using Distributions
 using StatsPlots
 using SparseArrays
 using Laplacians
 using DataStructures
 using RollingFunctions
 using BenchmarkTools
 using Statistics
 using Logging,LoggingExtras, LoggingFacilities
 using Dates,TimeZones

#Core structures
 include("core/core_structs.jl")
#graph-related
 include("core/graph_functions.jl")
#events
 include("eve/artificial_events.jl")
#various util functions
 include("utils/util_functions.jl")
#logging
 include("utils/logging_functions.jl")
#Plotting functions
 include("utils/plotting_functions.jl")
# time-to-event functions
 include("utils/tte_functions.jl")
#running functions
 include("utils/running_functions.jl")

#maintenance model
 include("phy/maintenance_model.jl")

 include("ntw/of_switch.jl")
 include("ntw/of_control.jl")
 include("ctl/agent_control.jl")

 include("phy/physical_model.jl")

 include("ntw/network_model.jl")
 include("phy/geo_model.jl")

#Main Functions
 include("model/netManFunctions.jl")

#Agents.jl function implementation for this model
 include("model/netManModel.jl")

#Basic queries
#include("model/queries_basic.jl")
#Multiple queries
 include("ctl/queries_multiple.jl")

export load_run_configs, single_run_with_logging

base_cfgs = load_base_cfgs("data/configs/configs_exp5_single.csv")
bcfg = base_cfgs[1]
#obtain services that provide this coverage of the network
coverage = 0.95
ntw_services = get_end_points(bcfg.seed,get_graph(bcfg.seed,bcfg.size,GraphModel(bcfg.ntw_topo_n);k=bcfg.k,adj_m_csv=bcfg.ntw_csv_adj_matrix),coverage)

cfg = config(bcfg,ntw_services)
single_run_with_logging(cfg)


text:
module CP

using Test

using Setfield
using StructArrays: StructArray

using QuantumESPRESSOBase
using QuantumESPRESSOBase.Inputs.CP

@testset "Constructing `AtomicVelocity`" begin
    # Data from https://gitlab.com/QEF/q-e/blob/master/CPV/examples/autopilot-example/reference/water.autopilot.out
    x = AtomicVelocity("H", [0.140374E-04, -0.333683E-04, 0.231834E-04])
    @test_throws AssertionError @set x.atom = "hydrogen"
    @test_throws InexactError @set x.velocity = [1im, 2im, 3im]
    @test_throws AssertionError x.velocity = [1, 2, 3, 4]
    @test_throws AssertionError x.velocity = [1, 2]
    x.atom = 'H'  # Setting `atom` with a `Char` still works
    @test x.atom == "H"
    @test_throws ErrorException x.velociti = [1, 2, 3] # Given a wrong field name
    @test x == AtomicVelocity("H", [0.140374E-04, -0.333683E-04, 0.231834E-04])
    y = AtomicVelocity('H')  # Incomplete initialization
    @test_throws UndefRefError y == AtomicVelocity("H")
    as = AtomicSpecies("H", 1.00794000, "H.pbe-rrkjus_psl.1.0.0.UPF")
    v = [0.140374E-04, -0.333683E-04, 0.231834E-04]
    @test AtomicVelocity(as).atom == "H"
    @test AtomicVelocity(as, v) == AtomicVelocity("H", v)
    ap = AtomicPosition('H', [0.500000000, 0.288675130, 1.974192764])
    @test AtomicVelocity(ap).atom == "H"
    @test AtomicVelocity(ap, v) == AtomicVelocity("H", v)
end # testset

@testset "Test constructing `AtomicVelocitiesCard` from `StructArray`" begin
    a = ["H", "O"]
    v = [
        [0.140374E-04, -0.333683E-04, 0.231834E-04],
        [-0.111125E-05, -0.466724E-04, 0.972745E-05],
    ]
    card = AtomicVelocitiesCard(StructArray{AtomicVelocity}((a, v)))
    @test card.data == [
        AtomicVelocity("H", [0.140374E-04, -0.333683E-04, 0.231834E-04]),
        AtomicVelocity("O", [-0.111125E-05, -0.466724E-04, 0.972745E-05]),
    ]
end # testset

@testset "Test `push_atom!`" begin
    @testset "`push_atom!` for `AtomicVelocity`" begin
        v = [AtomicVelocity("H", [0.140374E-04, -0.333683E-04, 0.231834E-04])]
        push_atom!(v, "S", "N")
        @test [x.atom for x in v] == ["H", "S", "N"]
    end # testset
    @testset "" begin
        a = ["H", "O"]
        v = [
            [0.140374E-04, -0.333683E-04, 0.231834E-04],
            [-0.111125E-05, -0.466724E-04, 0.972745E-05],
        ]
        card = AtomicVelocitiesCard(StructArray{AtomicVelocity}((a, v)))
        push_atom!(card, "S", "N")
        @test [x.atom for x in card.data] == ["H", "O", "S", "N"]
    end # testset
end # testset

@testset "Test `append_atom!`" begin
    @testset "`append_atom!` to `AtomicVelocity`" begin
        v = [AtomicVelocity("H", [0.140374E-04, -0.333683E-04, 0.231834E-04])]
        append_atom!(v, ["S", "N"])
        @test [x.atom for x in v] == ["H", "S", "N"]
    end # testset
    @testset "`append_atom!` to `AtomicVelocitiesCard`" begin
        a = ["H", "O"]
        v = [
            [0.140374E-04, -0.333683E-04, 0.231834E-04],
            [-0.111125E-05, -0.466724E-04, 0.972745E-05],
        ]
        card = AtomicVelocitiesCard(StructArray{AtomicVelocity}((a, v)))
        append_atom!(card, ["S", "N"])
        @test [x.atom for x in card.data] == ["H", "O", "S", "N"]
    end # testset
end # testset

@testset "Constructing `RefCellParametersCard`" begin
    option = "bohr"
    data = [
        12 0 0
        0 5 0
        0 0 5
    ]
    init = RefCellParametersCard(option, data)
    @test_throws AssertionError @set init.option = "alat"
    @test_throws AssertionError @set init.option = "crystal" # Allowed options are angstrom, bohr
    @test_throws AssertionError @set init.data = [ # Matrix size should be (3, 3)
        1 2
        3 4
    ]
    @test_throws AssertionError @set init.data = [
        1 2 3 4
        5 6 7 8
        4 3 2 1
        8 7 6 5
    ]
    @test RefCellParametersCard(data).option == "bohr"
end # testset

end # module CP


text:
using Luxor, Colors

L = 500
R = 130
r = 25
t = 18
h = 0.33

path_svg = joinpath(@__DIR__, "src", "assets", "logo.svg")
path_ico = joinpath(@__DIR__, "src", "assets", "favicon.ico")

@svg begin
	Drawing(L, L, path_svg)
	setlinecap("round")
	setlinejoin("round")
	p₊ = Point(cos(π/6), -sin(π/6)) * 2R/√3
	p₋ = -p₊/2
	angles = (-2π/3, -π/3, 0, π/3, 4π/3)

	for i in 1:3
		origin(L/2, L/2)
		rotate(4π/3*(i-1))
		sethue(Colors.JULIA_LOGO_COLORS[(:red, :purple, :green)[i]])
		circle(p₋, r, :fill)
		circle(p₊, r, :fill)
		setline(t)
		for j in 1:4
			arc(p₋, R, angles[j]+h, angles[j+1]-h, :stroke)
		end
		setline(t/2)
		poly(ngon(p₋*(1+√3), t, 3, π), :stroke, close=true)
	end
end

run(`convert -density 256x256 -background transparent $(path_svg) -define icon:auto-resize -colors 256 $(path_ico)`)


text:
@testset "Quantity Operations" begin
    At = CrossSectionalArea(collect(ts),randn(2))
    cvt = AlongChannelVelocity(collect(ts),randn(2))
    dvt = At*cvt
    @test quantity(dvt) == quantity(cvt).*quantity(At)

    ts2 = now():now()+Day(2)
    A2 = CrossSectionalArea(collect(ts2),randn(3))

    @test_throws ErrorException A2*cvt
end


text:
"""
    f_opt_eval(f,plan,samples;trace::Symbol=:silent)

Calculate the objective function value(s) and provide intermediate results printing
showing improvements over previous best iteration.
"""
function f_opt_eval(f,plan,samples;trace::Symbol=:silent)

    (trace == :verbose) && println("Evaluating function ",size(plan,2)," times")

    new_samples = mapslices(f,plan,dims=1)

    print_f_opt(trace,new_samples,samples,plan)

    return new_samples
end

function print_f_opt(trace,new_samples,samples,plan)

    new_min = minimum(new_samples)
    old_min = minimum(samples)
    new_max = maximum(new_samples)
    old_max = maximum(samples)

    if trace == :compact
        min_sample = minimum((old_min,new_min))
        if isless(new_min,old_min)
            print(@sprintf("%20.7g %20.7g\n", min_sample, old_min-new_min))
        else
            print(@sprintf("%20.7g %20.7s\n", min_sample, "N/A"))
        end

    elseif trace == :verbose
        println("    ---------------------------------------------------------------")
        print("    ")
        for sample in new_samples
            printstyled(@sprintf("%-15.7g",sample))
        end
        println("\t function evaluation\n")

        print("    Minimum sample value ")
        printstyled(@sprintf("%.7g",min(new_min,old_min)); color=:light_green, bold=true)
        print("\t(maximum = ")
        printstyled(@sprintf("%.7g",max(new_max,old_max));)
        print(")")

        if isless(new_min,old_min)
            print("\t (Improvement from last iteration ", @sprintf("%.7g",old_min-new_min),")\n")
        else
            print("\t (Improvement from last iteration N/A)\n")
        end
        print("\n")
    end
end

"""
    f_opt_eval(f,plan;trace::Symbol=:silent)

Calculate the objective function value(s) and plot value if trace.
"""
function f_opt_eval(f,plan;trace::Symbol=:silent)

    (trace == :verbose) && println("Evaluating function ",size(plan,2)," times")
    (trace == :compact) && print(@sprintf("%s %21s %20s\n", "   Iteration", "Function value", "Improvement"))

    new_samples = mapslices(f,plan,dims=1)

    if trace == :verbose
        new_min = minimum(new_samples)
        new_max = maximum(new_samples)

        print("    Minimum sample value ")
        printstyled(@sprintf("%.7g",new_min); color=:light_green, bold=true)
        print("\t(maximum = ")
        printstyled(@sprintf("%.7g",new_max);)
        println(")\n")
    end

    return new_samples
end

"""
    print_iteration(trace,i,iterations)

REPL printing of the current iteration if trace.
"""
function print_iteration(trace,i,iterations)
    if trace == :compact
        print(@sprintf("%4.4g%s%g",i, " out of ", iterations))
    end

    if trace == :verbose
        print("\n\nFunction eval. ")
        printstyled("iteration ";bold=true,color=:light_red)
        printstyled(i,bold=true)
        print(" out of ", iterations)
        print("\n")
    end
end

"""
    best_fitness(result::SurrogateResult)

Extract the best fitness from the `SurrogateResult` type.
"""
function best_fitness(result::SurrogateResult)
    minimum([result.lhc_samples result.infill_samples])
end

"""
    best_candidate(result::SurrogateResult)

Extract the design parameters yielding the best fitness from the `SurrogateResult` type.
"""
function best_candidate(result::SurrogateResult)
    min_ind = findmin([result.lhc_samples result.infill_samples])[2][2]
    return [result.lhc_plan result.infill_plan][:,min_ind:min_ind]
end

"""
    f_calls(result::SurrogateResult)

Extract the number of function calls used from the `SurrogateResult` type.
"""
function f_calls(result::SurrogateResult)
    return length([result.lhc_samples result.infill_samples])
end


text:
using Images

export fit_spike_shift

## these functions rely on the Ephys IO format for metadata & loading in spike times.

function fit!(model::PopulationModel, metadata::Dict, sessions,
		rates=(0.05,10.))
    spiketimes = loadspikes(SpikeTimes, metadata, sessions)
    cells2use = EPhys.ratefilt(metadata["cellIDs"], spiketimes, rates...)[2:end]
    Z = binZ(spiketimes)[:,cells2use]'
    fit!(model, Z)
    return model, Z
end

function fit!(model_type::Type, all_metadata::Array{<:Dict,1}, all_sessions;
		N=1000)
    Zs = []
    models = []
    btstrps = []

    for (metadata, sessions) in zip(all_metadata, all_sessions)
        model,Z = fit!(model_type(), metadata, sessions)
        btstrp = bootstrap(model, Z, N)
        push!(Zs     , Z)
        push!(models , model)
        push!(btstrps, btstrp)
    end
    return Zs,models,btstrps
end

function track_model(model_t, metadata::Dict, spiketimes::SpikeTimes;
                    ratefilt::Tuple{Float64,Float64} = (0.05,10.),
                    binsize::Float64 = 0.02,
                    tracking_resolution::Float64 = 0.001,
                    tracking_kernelwidth::Float64 = sqrt(binsize/
                                                         tracking_resolution),
                    tracking_type::Symbol = :full,
                    apply_tetrode_mask = true
                )
    cells2use = EPhys.ratefilt(metadata["cellIDs"], spiketimes,ratefilt[1],
                               ratefilt[2])[2:end]
    if apply_tetrode_mask
        # we don't want to track cells from the same tetrode together as
        # otherwise we might be biased by spike sorting errors
        mask   = [x==y for x in metadata["tetlist"][cells2use],
                     y in metadata["tetlist"][cells2use]]
    else
        mask   = [x==y for x in 1:sum(cells2use), y in 1:sum(cells2use)]
    end

    Z = binZ(spiketimes, binsize)[:,cells2use]'
    Zconv = binZ(spiketimes, tracking_resolution)[:,cells2use]
    Zconv = imfilter(Zconv, KernelFactors.IIRGaussian([tracking_kernelwidth;0]))

    model = model_t()
    PopulationDynamics.fit!(model, Z)

    if tracking_type == :full
        R_split = [track(Zconv', weights(model), mask)]
    elseif tracking_type == :regional
        idmatch = cat([contains.(metadata["cellIDs"][cells2use]) for id in regions]...,dims=2)
        multipliers = [BitArray(idmatch[:,i] * idmatch[:,j]')
            for i in indices(idmatch,2) for j in indices(idmatch,2)]
        R_split = track_partial(Zconv, weights(model), multipliers)
    elseif tracking_type == :diagvsoff
        idmatch = cat([occursin.(id,metadata["cellIDs"][cells2use]) for id in regions]..., dims=2)
        multiplier_ = BitArray(idmatch * idmatch')
        multipliers = [multiplier_, .!(multiplier_)]
        R_split = track_partial(Zconv, PopulationDynamics.weights(model), multipliers)
    else
        error("""
        No tracking type $tracking_type:
        Should be one of:
            - full
            - regional
            - diagvsoff
        """)
    end
    #@warn "This now returns both the model, the tracked trace, and the binned spike times"
    return model, R_split, Z
end

function fit_spike_shift(model_type::Type, metadata, sessions, jitters, nreps=10;
    verbose=false)
    # detect "real ensembles"
    spiketimes = loadspikes(SpikeTimes, metadata, sessions)
    cells2use = EPhys.ratefilt(metadata["cellIDs"], spiketimes, 0.05, 10.)[2:end]
    Z = binZ(spiketimes)[:,cells2use]'
    model = model_type()
    fit!(model, Z)

    # for each mean jitter jit
    models_alljitters = []
    for jitter in jitters
        verbose && println("Refitting with jitters in the range of $(jitter[1])-$(-jitter[1])")
        verbose && tic()
        models_jitter = pmap(_ -> shift_refit(model, spiketimes, jitter, cells2use), 1:nreps)
        push!(models_alljitters, models_jitter)
        verbose && toc()
    end
    model, models_alljitters
end

function shift_refit(model, spiketimes, jitter, cells2use)
    model_shifted = typeof(model)(model.k, pca(), ica(), srand())
    # shift spikes
    spiketimes_shifted = EPhys.addjitter(spiketimes, jitter)
    Z_shifted = binZ(spiketimes_shifted)[:,cells2use]' # this is hacky, "spiketimes should already be filtered"
    # re-detect
    fit!(model_shifted, Z_shifted)
    return model_shifted
end


text:
module Messaging

using Serialization, Base64, ZMQ, FileWatching

export sendrecv, recvsend, MessageSocket, RequestSocket, ReplySocket

function send2zmq(sock::Socket, obj)
    buff = IOBuffer()
    serialize(buff, obj)
    data = copy(buff.data)
    send(sock, base64encode(data))
    nothing
end

function recvfromzmq(sock::Socket)
    data = recv(sock, String)
    buff = IOBuffer(base64decode(data))
    deserialize(buff)
end

function request_socket(ep)
    sock = Socket(REQ)
    connect(sock, ep)
    return sock
end

function reply_socket(ep)
    sock = Socket(REP)
    bind(sock, ep)
    return sock
end

abstract type MessageSocket end

"""
A socket for making requests to a ZMQ service.

Parameters:
- `port::Int=5555`: The TCP port to connect to
"""
mutable struct RequestSocket <: MessageSocket
    sock::Socket
    ep::String

    function RequestSocket(port::Int=5555)
        ep = "tcp://localhost:$port"
        new(request_socket(ep), ep)
    end
end

"""
A socket for serving replies to ZMQ clients.

Parameters:
- `port::Int=5555`: The TCP port to connect to
"""
mutable struct ReplySocket <: MessageSocket
    sock::Socket
    ep::String

    function ReplySocket(port::Int=5555)
        ep = "tcp://*:$port"
        new(reply_socket(ep), ep)
    end
end

function reconnect!(sock::RequestSocket)
    close(sock.sock)
    sock.sock = request_socket(sock.ep)
    nothing
end

function reconnect!(sock::ReplySocket)
    close(sock.sock)
    sock.sock = reply_socket(sock.ep)
    nothing
end

"""
    sendrecv(sock::RequestSocket, obj; timeout::Real=60.0, nretry::Number=Inf)
Send an `::Any` object and wait for reply.

Parameters:
- `sock::RequestSocket`
- `obj::Any`: The object to be sent

Keyword Arguments:
- `timeout::Real=60.0`: Time in seconds to wait for reply
- `nretry::Number=Inf`: Number of retries after timeout

Returns: Returns `nothing` if `nretry`s are up else returns received `::Any` object
"""
function sendrecv(sock::RequestSocket, obj; timeout::Real=60.0, nretry::Number=Inf)
    send2zmq(sock.sock, obj)

    while nretry > 0
        rawfd = RawFD(sock.sock.fd)
        istimedout = false
        try
            event = FileWatching.poll_fd(rawfd, timeout; readable=true)
            istimedout = event.timedout
        catch ex
            @debug "Socket was closed when waiting for reply"
            ex isa EOFError && return nothing
        end

        if istimedout
            @debug "Timeout waiting for reply. Reconnecting..."
            reconnect!(sock)
            send2zmq(sock.sock, obj)
        else
            return recvfromzmq(sock.sock)
        end
        nretry -= 1
    end
    nothing
end

"""
    recvsend(f::Function, sock::ReplySocket)
Receive an ::Any object from a socket, process it and send a reply.

Example:
```
recvsend(repsock) do obj
    obj.x += 5
    obj
end
```

Parameters:
- `f::Function`: A function that accepts the received object and returns the object to send
- `sock::ReplySocket`

Returns:
- `true`: If send was successful after receive
- `false`: If socket was closed during receive
"""
function recvsend(f::Function, sock::ReplySocket)
    # TODO: This should be done in a better way. Reasons why this was not done in a better way:
    # 1) FileWatching.poll_fd does not poll reliable when waiting for a request.
    # 2) zmq_poll hangs the julia process so no other tasks can run
    while isopen(sock.sock) && sock.sock.events & ZMQ.POLLIN != ZMQ.POLLIN
        sleep(5)
    end
    !isopen(sock.sock) && return false
    obj = recvfromzmq(sock.sock)
    sobj = f(obj)
    send2zmq(sock.sock, sobj)
    true
end

end


text:
@testset "Equilibria" begin
    @testset "Nothing" begin
        feq = zeros(5,5,9)
        Swalbe.equilibrium!(feq, zeros(5,5), zeros(5,5), zeros(5,5), zeros(5,5), 0.0)
        @test all(feq[:,:,:] .== 0.0)
        Swalbe.equilibrium!(feq, zeros(5,5), zeros(5,5), zeros(5,5), zeros(5,5))
        @test all(feq[:,:,:] .== 0.0)
    end
    @testset "Density only" begin
        feq = zeros(5,5,9)
        Swalbe.equilibrium!(feq, ones(5,5), zeros(5,5), zeros(5,5), zeros(5,5), 0.0)
        @test all(feq[:,:,1] .== 1.0)
        @test all(feq[:,:,2:9] .== 0.0)
        Swalbe.equilibrium!(feq, ones(5,5), zeros(5,5), zeros(5,5), zeros(5,5))
        @test all(feq[:,:,1] .== 1.0)
        @test all(feq[:,:,2:9] .== 0.0)
    end
    @testset "Density and gravity" begin
        feq = zeros(5,5,9)
        Swalbe.equilibrium!(feq, ones(5,5), zeros(5,5), zeros(5,5), zeros(5,5), 0.1)
        @test all(feq[:,:,1] .== 1.0 - 1/12)
        @test all(feq[:,:,2:5] .≈ 1/9 * 1.5 * 0.1)
        @test all(feq[:,:,6:9] .≈ 1/36 * 1.5 * 0.1)
    end
    @testset "Density and velocity" begin
        feq = zeros(5,5,9)
        # with gravity
        Swalbe.equilibrium!(feq, ones(5,5), fill(0.1,5,5), fill(-0.1,5,5), zeros(5,5), 0.0)
        @test all(feq[:,:,1] .== 1.0 - 2/3 * 0.02)
        @test all(feq[:,:,2] .≈  1/9 * (3 * 0.1 + 4.5 * 0.01 - 3/2 * 0.02))
        @test all(feq[:,:,3] .≈  1/9 * (3 * -0.1 + 4.5 * 0.01 - 3/2 * 0.02))
        @test all(feq[:,:,4] .≈  1/9 * (3 * -0.1 + 4.5 * 0.01 - 3/2 * 0.02))
        @test all(feq[:,:,5] .≈  1/9 * (3 * 0.1 + 4.5 * 0.01 - 3/2 * 0.02))
        @test all(feq[:,:,6] .≈  1/36 * (-3/2 * 0.02))
        @test all(feq[:,:,7] .≈  1/36 * (3 * -0.2  + 4.5 * 0.2^2 - 3/2 * 0.02))
        @test all(feq[:,:,8] .≈  1/36 * (-3/2 * 0.02))
        @test all(feq[:,:,9] .≈  1/36 * (3 * 0.2 + 4.5 * 0.2^2 - 3/2 * 0.02))
        # without gravity
        Swalbe.equilibrium!(feq, ones(5,5), fill(0.1,5,5), fill(-0.1,5,5), zeros(5,5))
        @test all(feq[:,:,1] .== 1.0 - 2/3 * 0.02)
        @test all(feq[:,:,2] .≈  1/9 * (3 * 0.1 + 4.5 * 0.01 - 3/2 * 0.02))
        @test all(feq[:,:,3] .≈  1/9 * (3 * -0.1 + 4.5 * 0.01 - 3/2 * 0.02))
        @test all(feq[:,:,4] .≈  1/9 * (3 * -0.1 + 4.5 * 0.01 - 3/2 * 0.02))
        @test all(feq[:,:,5] .≈  1/9 * (3 * 0.1 + 4.5 * 0.01 - 3/2 * 0.02))
        @test all(feq[:,:,6] .≈  1/36 * (-3/2 * 0.02))
        @test all(feq[:,:,7] .≈  1/36 * (3 * -0.2  + 4.5 * 0.2^2 - 3/2 * 0.02))
        @test all(feq[:,:,8] .≈  1/36 * (-3/2 * 0.02))
        @test all(feq[:,:,9] .≈  1/36 * (3 * 0.2 + 4.5 * 0.2^2 - 3/2 * 0.02))
        # using the state struct
        state = Swalbe.Sys(Swalbe.SysConst(Lx=5, Ly=5), "CPU")
        state.feq .= 0.0
        state.height .= 1.0
        state.velx .= 0.1
        state.vely .= -0.1
        Swalbe.equilibrium!(state)
        @test all(state.feq[:,:,1] .== 1.0 - 2/3 * 0.02)
        @test all(state.feq[:,:,2] .≈  1/9 * (3 * 0.1 + 4.5 * 0.01 - 3/2 * 0.02))
        @test all(state.feq[:,:,3] .≈  1/9 * (3 * -0.1 + 4.5 * 0.01 - 3/2 * 0.02))
        @test all(state.feq[:,:,4] .≈  1/9 * (3 * -0.1 + 4.5 * 0.01 - 3/2 * 0.02))
        @test all(state.feq[:,:,5] .≈  1/9 * (3 * 0.1 + 4.5 * 0.01 - 3/2 * 0.02))
        @test all(state.feq[:,:,6] .≈  1/36 * (-3/2 * 0.02))
        @test all(state.feq[:,:,7] .≈  1/36 * (3 * -0.2  + 4.5 * 0.2^2 - 3/2 * 0.02))
        @test all(state.feq[:,:,8] .≈  1/36 * (-3/2 * 0.02))
        @test all(state.feq[:,:,9] .≈  1/36 * (3 * 0.2 + 4.5 * 0.2^2 - 3/2 * 0.02))

    end
    @testset "Density and gravity and velocity" begin
        feq = zeros(5,5,9)
        Swalbe.equilibrium!(feq, ones(5,5), fill(0.1,5,5), fill(-0.1,5,5), zeros(5,5), 0.1)
        @test all(feq[:,:,1] .== 1.0 - 1/12 - 2/3 * 0.02)
        @test all(feq[:,:,2] .≈  1/9 * (1.5 * 0.1 + 3 * 0.1 + 4.5 * 0.01 - 3/2 * 0.02))
        @test all(feq[:,:,3] .≈  1/9 * (1.5 * 0.1 + 3 * -0.1 + 4.5 * 0.01 - 3/2 * 0.02))
        @test all(feq[:,:,4] .≈  1/9 * (1.5 * 0.1 + 3 * -0.1 + 4.5 * 0.01 - 3/2 * 0.02))
        @test all(feq[:,:,5] .≈  1/9 * (1.5 * 0.1 + 3 * 0.1 + 4.5 * 0.01 - 3/2 * 0.02))
        @test all(feq[:,:,6] .≈  1/36 * (1.5 * 0.1 + -3/2 * 0.02))
        @test all(feq[:,:,7] .≈  1/36 * (1.5 * 0.1 + 3 * -0.2  + 4.5 * 0.2^2 - 3/2 * 0.02))
        @test all(feq[:,:,8] .≈  1/36 * (1.5 * 0.1 + -3/2 * 0.02))
        @test all(feq[:,:,9] .≈  1/36 * (1.5 * 0.1 + 3 * 0.2 + 4.5 * 0.2^2 - 3/2 * 0.02))
    end
end

@testset "Equilibria 1D" begin
    @testset "Nothing" begin
        feq = zeros(30,3)
        Swalbe.equilibrium!(feq, zeros(30), zeros(30), 0.0)
        @test all(feq[:,:] .== 0.0)
        Swalbe.equilibrium!(feq, zeros(30), zeros(30))
        @test all(feq[:,:] .== 0.0)
    end
    @testset "Density only" begin
        feq = zeros(30,3)
        Swalbe.equilibrium!(feq, ones(30), zeros(30), 0.0)
        @test all(feq[:,1] .== 1.0)
        @test all(feq[:,2:3] .== 0.0)
        Swalbe.equilibrium!(feq, ones(30), zeros(30))
        @test all(feq[:,1] .== 1.0)
        @test all(feq[:,2:3] .== 0.0)
    end
    @testset "Density and gravity" begin
        feq = zeros(30,3)
        Swalbe.equilibrium!(feq, ones(30), zeros(30), 0.1)
        @test all(feq[:,1] .== 1.0 - 0.05)
        @test all(feq[:,2:3] .== 0.025)

    end
    @testset "Density and velocity" begin
        feq = zeros(30,3)
        Swalbe.equilibrium!(feq, ones(30), fill(0.1,30), 0.0)
        @test all(feq[:,1] .== 1.0 - 0.01)
        @test all(feq[:,2] .≈  0.5 * 0.1 + 0.5 * 0.01)
        @test all(feq[:,3] .≈  -0.5 * 0.1 + 0.5 * 0.01)

        Swalbe.equilibrium!(feq, ones(30), fill(0.1,30))
        @test all(feq[:,1] .== 1.0 - 0.01)
        @test all(feq[:,2] .≈  0.5 * 0.1 + 0.5 * 0.01)
        @test all(feq[:,3] .≈  -0.5 * 0.1 + 0.5 * 0.01)

        state = Swalbe.Sys(Swalbe.SysConst_1D(L=30))
        state.feq .= 0.0
        state.height .= 1.0
        state.vel .= 0.1
        Swalbe.equilibrium!(state)
        @test all(state.feq[:,1] .== 1.0 - 0.01)
        @test all(state.feq[:,2] .≈  0.5 * 0.1 + 0.5 * 0.01)
        @test all(state.feq[:,3] .≈  -0.5 * 0.1 + 0.5 * 0.01)

    end
    @testset "Density and gravity and velocity" begin
        feq = zeros(30,3)
        Swalbe.equilibrium!(feq, ones(30), fill(0.1,30), 0.1)
        @test all(feq[:,1] .== 1.0 - 0.5 * 0.1 - 0.01)
        @test all(feq[:,2] .≈  0.025 + 0.5 * 0.1 + 0.5 * 0.01)
        @test all(feq[:,3] .≈  0.025 - 0.5 * 0.1 + 0.5 * 0.01)
    end
end


text:
function PotentialSimulationSetupRB(ssd::SolidStateDetector{T, :cylindrical}, grid::Grid{T, 3, :cylindrical} = Grid(ssd),
                potential_array::Union{Missing, Array{T, 3}} = missing; sor_consts = [1.0, 1.0],
                weighting_potential_contact_id::Union{Missing, Int} = missing
                )::PotentialSimulationSetupRB{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)

    @inbounds begin
        begin # Geometrical weights of the Axes
            nr, nφ, nz = size(grid)

            # R-axis
            axr::Vector{T} = collect(grid.axes[1]) # real grid points/ticks -> the potential at these ticks are going to be calculated
            Δr::Vector{T} = diff(axr) # difference between real ticks
            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_inv::Vector{T} = inv.(mpr)
            Δ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
            cyclic::T = grid.axes[2].interval.right
            Δφ::Vector{T} = diff(axφ)  # difference between real ticks
            φ_ext::Vector{T} = get_extended_ticks(grid.axes[2])
            # φ_ext::Vector{T} = if cyclic > 0
            #     get_extended_ticks(grid.axes[2])
            # else
            #     T[-2π, 0, 2π]
            # end
            Δφ_ext::Vector{T} = diff(φ_ext)
            Δφ_ext_inv::Vector{T} = inv.(Δφ_ext)
            mpφ::Vector{T} = midpoints(φ_ext)
            mpφ_inv::Vector{T} = inv.(mpφ)
            Δ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::Vector{T} = diff(axz) # difference between real ticks
            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_inv::Vector{T} = inv.(mpz)
            Δ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

        # detector_material_ϵ_r::T = ssd.material_detector.ϵ_r
        # detector_material_ϵ_r::T = ssd.semiconductors[1].material.ϵ_r
        # environment_material_ϵ_r::T = ssd.medium.ϵ_r

        # bulk_is_ptype::Bool = ssd.bulk_type == :ptype ? true : false
        bulk_is_ptype::Bool = ssd.semiconductors[1].bulk_type == :ptype ? true : false
        bias_voltages::Vector{T} = [i.potential for i in ssd.contacts]
        minimum_applied_potential::T = minimum(bias_voltages)
        maximum_applied_potential::T = maximum(bias_voltages)
        bias_voltage::T = maximum_applied_potential - minimum_applied_potential
        depletion_handling_potential_limit::T = -bias_voltage
        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]

        ϵ = Array{T, 3}(undef,    length(mpr), length(mpφ), length(mpz))
        ρ_tmp = Array{T, 3}(undef, length(mpr), length(mpφ), length(mpz))
        for iz in 1:size(ϵ, 3)
            pos_z::T = mpz[iz]
            for iφ in 1:size(ϵ, 2)
                pos_φ::T = mpφ[iφ]

                ir = 1
                pos_r::T = mpr[ir]
                if axr[1] == 0
                    pos_r = axr[2] * 0.5
                end
                pt::CylindricalPoint{T} = CylindricalPoint{T}(pos_r, pos_φ, pos_z)
                ρ_tmp[ir, iφ, iz]::T, ϵ[ir, iφ, iz]::T = get_ρ_and_ϵ(pt, ssd)

                for ir in 2:size(ϵ, 1)
                    pos_r = mpr[ir]
                    pt = CylindricalPoint{T}(pos_r, pos_φ, pos_z)
                    ρ_tmp[ir, iφ, iz]::T, ϵ[ir, iφ, iz]::T = get_ρ_and_ϵ(pt, ssd)
                end
            end
        end
        ϵ0_inv::T = inv(ϵ0)
        ρ_tmp *= ϵ0_inv

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

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

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

                            ρ_cell += ρ_tmp[inr,  iφ,  iz] * wzr[inz] * wrl[inr] * wφr[inφ]
                            ρ_cell += ρ_tmp[inr,  iφ, inz] * wzl[inz] * wrl[inr] * wφr[inφ]
                            ρ_cell += ρ_tmp[inr, inφ,  iz] * wzr[inz] * wrl[inr] * wφl[inφ]
                            ρ_cell += ρ_tmp[inr, inφ, inz] * wzl[inz] * wrl[inr] * wφl[inφ]
                        else
                            ρ_cell += ρ_tmp[ ir,  iφ,  iz] * wzr[inz] * 0.5 #wφr[inφ]
                            ρ_cell += ρ_tmp[ ir,  iφ, inz] * wzl[inz] * 0.5 #wφr[inφ]
                            ρ_cell += ρ_tmp[ ir, inφ,  iz] * wzr[inz] * 0.5 #wφl[inφ]
                            ρ_cell += ρ_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]
                        ρ[ irbz, iφ, ir, rbi ] = dV * ρ_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]
                        ρ[ irbz, iφ, ir, rbi ] = dV * ρ_cell
                    end
                end
            end
        end

        potential::Array{T, 3} = ismissing(potential_array) ? zeros(T, size(grid)...) : potential_array
        pointtypes::Array{PointType, 3} = ones(PointType, size(grid)...)
        set_pointtypes_and_fixed_potentials!( pointtypes, potential, grid, ssd, weighting_potential_contact_id = weighting_potential_contact_id  )
        rbpotential::Array{T, 4}  = RBExtBy2Array( potential, grid )
        rbpointtypes::Array{T, 4} = RBExtBy2Array( pointtypes, grid )
        potential = clear(potential)
        pointtypes = clear(pointtypes)
    end # @inbounds

    fssrb::PotentialSimulationSetupRB{T, 3, 4, :cylindrical} = PotentialSimulationSetupRB{T, 3, 4, :cylindrical}(
        grid,
        rbpotential,
        rbpointtypes,
        volume_weights,
        ρ,
        ϵ,
        geom_weights,
        sor_const,
        bias_voltage,
        maximum_applied_potential,
        minimum_applied_potential,
        depletion_handling_potential_limit,
        bulk_is_ptype,
        grid_boundary_factors
     )

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

function Grid(fssrb::PotentialSimulationSetupRB{T, N1, N2, S})::Grid{T, N1, S} where {T, N1, N2, S}
    return fssrb.grid
end

function ElectricPotentialArray(fssrb::PotentialSimulationSetupRB{T, 3, 4, :cylindrical})::Array{T, 3} where {T}
    pot::Array{T, 3} = Array{T, 3}(undef, size(fssrb.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] = fssrb.potential[ irbz, irbφ, irbr, rbi ]
            end
        end
    end
    return pot
end

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

function PointTypes(fss::PotentialSimulationSetup{T, N, S})::PointTypes{T, N, S} where {T, N, S}
    return PointTypes{T, N, S}( fss.pointtypes, fss.grid )
end

function ChargeDensityArray(fssrb::PotentialSimulationSetupRB{T, 3, 4, :cylindrical})::Array{T} where {T}
    ρ::Array{T, 3} = zeros(T, size(fssrb.grid))
    for iz in axes(ρ, 3)
        irbz::Int = rbidx(iz)
        Δmpz::T = fssrb.geom_weights[3].weights[3, iz]
        for iφ in axes(ρ, 2)
            irbφ::Int = iφ + 1
            idxsum::Int = iz + iφ
            Δmpφ::T = fssrb.geom_weights[2].weights[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 = fssrb.geom_weights[1].weights[6, ir] * Δmpzφ  #Δmpz[inz] * Δmpφ[inφ] * Δmpr_squared[inr]
                if ir == 1
                    dV = dV * 2π / Δmpφ
                end
                ρ[ir, iφ, iz] = fssrb.ρ[irbz, irbφ, irbr, rbi ] / dV
            end
        end
    end
    return ρ
end

function DielektrikumDistributionArray(fssrb::PotentialSimulationSetupRB{T, 3, 4, S})::Array{T, 3} where {T, S}
    return fssrb.ϵ
end


text:
# tests to be run on a computer where LTspice is installed.
using LTspice
using Test

function localtests()

  test1 = LTspiceSimulation("test1.asc",tempdir=true)
  v = 20.0
  r = 2.0
  i = test1(v,r)[1]
  @test i == 10.0

  testinc1 = LTspiceSimulation("testInc1.asc",tempdir=true)
  @test testinc1["incA"] == 1.0
  @test testinc1["incB"] == 2.0
  @test testinc1["incC"] == 3.0
  @test testinc1["incD"] == 4.0
  @test testinc1["incE"] == 5.0
  @test testinc1["incF"] == 6.0

  test14 = LTspiceSimulation("test14.asc",tempdir=true)
  testarray = [1.0,2.0]
  for a in testarray
    for b in 10*testarray
      for d in 100*testarray
        for e in 1000*testarray
          for this_is_a_long_name in 10000*testarray
            test14["a"] = a
            test14["b"] = b
            test14["d"] = d
            test14["e"] = e
            test14["this_is_a_long_name"] = this_is_a_long_name
            @test test14["m1"] == a + b + d + e + this_is_a_long_name
          end
        end
      end
    end
  end

end
localtests()


text:
using Blink
using CSV
using DataFrames
using Interact
using Plots

# specify backend
plotlyjs()

function main()
    df1 = DataFrame(Any[10*randn(500), 10*rand(500), 10*rand(500)], [:AAA, :BBB, :CCC])
    @show names(df1)
    CSV.write("sample.csv", df1)
    loadbutton = filepicker()
    columnbuttons = Observable{Any}(dom"div"())
    data = Observable{Any}(DataFrame)
    plt = Observable{Any}(plot())
    map!(CSV.read, data, loadbutton)

    function makebuttons(df)
        buttons = button.(string.(names(df)))
        for (btn, name) in zip(buttons, names(df))
            map!(t -> histogram(df[name]), plt, btn)
        end
        dom"div"(hbox(buttons))
    end

    map!(makebuttons, columnbuttons, data)

    ui = dom"div"(loadbutton, columnbuttons, plt)

    w = Window()
    body!(w, ui)
end

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
module InverseProblems

using OrderedCollections

using ..Observations: obs_str, AbstractObservation, OneDimensionalTimeSeries, initialize_simulation!, FieldTimeSeriesCollector,
                      observation_times

using ..TurbulenceClosureParameters: free_parameters_str, update_closure_ensemble_member!

using OffsetArrays

using Oceananigans: short_show, run!, fields, FieldTimeSeries, CPU
using Oceananigans.OutputReaders: InMemory
using Oceananigans.Fields: interior, location
using Oceananigans.Grids: Flat, Bounded,
                          Face, Center,
                          RegularRectilinearGrid, offset_data,
                          topology, halo_size,
                          interior_parent_indices

import ..Observations: normalize!

#####
##### Output maps (maps from simulation output to observation space)
#####

abstract type AbstractOutputMap end

output_map_type(fp) = output_map_str(fp)

struct ConcatenatedOutputMap end

output_map_str(::ConcatenatedOutputMap) = "ConcatenatedOutputMap"

"""
    ConcatenatedVectorNormMap()

Forward map transformation of simulation output to a scalar by
taking a naive `norm` of the difference between concatenated vectors of the
observations and simulation output.
"""
struct ConcatenatedVectorNormMap end

output_map_str(::ConcatenatedVectorNormMap) = "ConcatenatedVectorNormMap"

#####
##### InverseProblems
#####

struct InverseProblem{F, O, S, T, P}
    observations :: O
    simulation :: S
    time_series_collector :: T
    free_parameters :: P
    output_map :: F
end

"""
    InverseProblem(observations, simulation, free_parameters; output_map=ConcatenatedOutputMap())

Return an `InverseProblem`.
"""
function InverseProblem(observations, simulation, free_parameters; output_map=ConcatenatedOutputMap(), time_series_collector=nothing)

    if isnothing(time_series_collector) # attempt to construct automagically
        simulation_fields = fields(simulation.model)
        collected_fields = NamedTuple(name => simulation_fields[name] for name in keys(simulation_fields))
        time_series_collector = FieldTimeSeriesCollector(collected_fields, observation_times(observations))
    end

    return InverseProblem(observations, simulation, time_series_collector, free_parameters, output_map)
end

function Base.show(io::IO, ip::InverseProblem)
    sim_str = "Simulation on $(short_show(ip.simulation.model.grid)) with Δt=$(ip.simulation.Δt)"

    out_map_type = output_map_type(ip.output_map)
    out_map_str = output_map_str(ip.output_map)

    print(io, "InverseProblem{$out_map_type}", '\n',
              "├── observations: $(obs_str(ip.observations))", '\n',
              "├── simulation: $sim_str", '\n',
              "├── free_parameters: $(free_parameters_str(ip.free_parameters))", '\n',
              "└── output map: $out_map_str")

    return nothing
end

tupify_parameters(ip, θ) = NamedTuple{ip.free_parameters.names}(Tuple(θ))

#####
##### Forward map evaluation given vector-of-vector (one parameter vector for each ensemble member)
#####

const OneDimensionalEnsembleGrid = RegularRectilinearGrid{<:Any, Flat, Flat, Bounded}

n_ensemble(grid::OneDimensionalEnsembleGrid) = grid.Nx
n_observations(grid::OneDimensionalEnsembleGrid) = grid.Ny
n_z(grid::OneDimensionalEnsembleGrid) = grid.Nz

n_ensemble(ip::InverseProblem) = n_ensemble(ip.simulation.model.grid)

""" Transform and return `ip.observations` appropriate for `ip.output_map`. """
observation_map(ip::InverseProblem) = transform_observations(ip.output_map, ip.observations)

function forward_map(ip::InverseProblem, θ::Vector{<:NamedTuple})
    observations = ip.observations
    simulation = ip.simulation
    closures = simulation.model.closure

    for p in 1:length(θ)
        update_closure_ensemble_member!(closures, p, θ[p])
    end

    initialize_simulation!(simulation, observations, ip.time_series_collector)
    run!(simulation)

    return transform_output(ip.output_map, observations, ip.time_series_collector)
end

forward_map(ip, θ::Vector{<:Vector}) = forward_map(ip, [tupify_parameters(ip, p) for p in θ])
forward_map(ip, θ::Matrix) = forward_map(ip, [tupify_parameters(ip, θ[:, i]) for i in 1:size(θ, 2)])

(ip::InverseProblem)(θ) = forward_map(ip, θ)

function transform_observations(::ConcatenatedOutputMap, observation::OneDimensionalTimeSeries)
    flattened_normalized_data = []

    for field_name in keys(observation.field_time_serieses)
        field_time_series = observation.field_time_serieses[field_name]

        # *** FIXME ***
        # Here we hack an implementation of `interior` because
        # `field_time_series.grid` may be wrong (sometimes grid.Nx is wrong)
        grid = field_time_series.grid
        Hx, Hy, Hz = halo_size(grid)
        Nx, Ny, Nz = size(grid)
        topo = topology(grid)
        loc = location(field_time_series)
        y_indices = interior_parent_indices(loc[2], topo[2], Ny, Hy)
        z_indices = interior_parent_indices(loc[3], topo[3], Nz, Hz)
        field_time_series_data = Array(view(parent(field_time_series), :, y_indices, z_indices, :))

        Nx, Ny, Nz, Nt = size(field_time_series_data)
        field_time_series_data = reshape(field_time_series_data, Nx, Ny * Nz * Nt)

        normalize!(field_time_series_data, observation.normalization[field_name])

        push!(flattened_normalized_data, field_time_series_data)
    end

    transformed = hcat(flattened_normalized_data...)

    return Matrix(transpose(transformed))
end

transform_observations(map, observations::Vector) =
    hcat(Tuple(transform_observations(map, observation) for observation in observations)...)

function transform_output(map::ConcatenatedOutputMap,
                          observations::Union{OneDimensionalTimeSeries, Vector{<:OneDimensionalTimeSeries}},
                          time_series_collector)

    # transposed_output isa Vector{OneDimensionalTimeSeries} where OneDimensionalTimeSeries is Nx by Nz by Nt
    transposed_output = transpose_model_output(time_series_collector, observations)[1]

    return transform_observations(map, transposed_output)
end

vectorize(observation) = [observation]
vectorize(observations::Vector) = observations

"""
    transpose_model_output(time_series_collector, observations)

Transpose a `NamedTuple` of 4D `FieldTimeSeries` model output collected by `time_series_collector`
into a Vector of `OneDimensionalTimeSeries` for each member of the observation batch.

Return a 1-vector in the case of singleton observations.
"""
function transpose_model_output(time_series_collector, observations)
    observations = vectorize(observations)
    times = time_series_collector.times

    transposed_output = []

    n_ensemble = time_series_collector.grid.Nx
    n_batch = time_series_collector.grid.Ny
    Nz = time_series_collector.grid.Nz
    Hz = time_series_collector.grid.Hz
    Nt = length(times)

    for j = 1:n_batch
        observation = observations[j]
        grid = observation.grid
        time_serieses = OrderedDict{Any, Any}()

        for name in keys(observation.field_time_serieses)
            loc = LX, LY, LZ = location(observation.field_time_serieses[name])
            topo = topology(grid)

            field_time_series = time_series_collector.field_time_serieses[name]

            raw_data = parent(field_time_series.data)
            data = OffsetArray(view(raw_data, :, j:j, :, :), 0, 0, -Hz, 0)

            # Note: FieldTimeSeries.grid.Nx is in general incorrect
            time_series = FieldTimeSeries{LX, LY, LZ, InMemory}(data, CPU(), grid, nothing, times)
            time_serieses[name] = time_series
        end

        # Convert to NamedTuple
        time_serieses = NamedTuple(name => time_series for (name, time_series) in time_serieses)

        batch_output = OneDimensionalTimeSeries(time_serieses,
                                                grid, # this grid has the wrong grid.Nx --- forgive us our sins
                                                times,
                                                nothing,
                                                nothing,
                                                observation.normalization)

        push!(transposed_output, batch_output)
    end

    return transposed_output
end

end # module