text:
export vertex_color, chromatic_number

"""
`vertex_color(G,k)`: Return a `k`-coloring of `G` (or error if none exists).

If `k` is omitted, the chromatic number of `G` is used.
"""
function vertex_color(G::SimpleGraph, k::Int)
    VV = vlist(G)
    EE = elist(G)
    n = NV(G)
    m = NE(G)

    if k < 1
        error("Number of colors must be positive")
    end

    err_msg = "This graph is not " * string(k) * " colorable"

    # Special cases that don't require integer programming

    # see if we know the chromatic number already
    if cache_check(G, :chromatic_number)
        chi = cache_recall(G, :chromatic_number)
        if k < chi
            error(err_msg)
        end
    end

    result = Dict{eltype(G),Int}()

    if k == 1
        if NE(G) > 0
            error(err_msg)
        end
        for v in VV
            result[v] = 1
        end
        return result
    end

    if k == 2
        return two_color(G)
    end

    MOD = Model(get_solver())

    @variable(MOD, x[VV, 1:k], Bin)

    for v in VV
        @constraint(MOD, sum(x[v, i] for i = 1:k) == 1)
    end

    for e in EE
        u = e[1]
        v = e[2]
        for i = 1:k
            @constraint(MOD, x[u, i] + x[v, i] <= 1)
        end
    end

    optimize!(MOD)
    status = Int(termination_status(MOD))

    if status != 1
        error(err_msg)
    end

    X = value.(x)

    for v in VV
        for c = 1:k
            if X[v, c] > 0
                result[v] = c
            end
        end
    end

    return result
end

vertex_color(G::SimpleGraph) = vertex_color(G, chromatic_number(G))

"""
`vertex_color(G,a,b)` returns an `a:b`-coloring of the graph `G`.
This is a mapping from the vertices of `G` to `b`-element subsets of
`{1,2,...,a}` such that adjacent vertices are assigned disjoint sets.
An error is thrown if no such coloring exists.
"""
function vertex_color(G::SimpleGraph, a::Int, b::Int)
    if a < 0 || b < 0
        error("Arguments in vertex_color(G,$a,$b) must be nonnegative")
    end
    try
        return hom(G, Kneser(a, b))
    catch
        error("This graph does not have a $a:$b coloring")
    end
end

function chromatic_number(G::SimpleGraph, verb::Bool = false)::Int
    if cache_check(G, :chromatic_number)
        return cache_recall(G, :chromatic_number)
    end

    if NV(G) == 0
        return 0
    end
    if NE(G) == 0
        return 1
    end

    # lower bound: larger of clique size or n/alpha
    n = NV(G)
    alf = length(max_indep_set(G))
    lb1 = Int(floor(n / alf))
    lb2 = length(max_clique(G))
    lb = max(lb1, lb2)

    # upper bound: try a random greedy coloring
    f = greedy_color(G)
    ub = maximum(values(f))

    k = chromatic_number_work(G, lb, ub, verb)
    cache_save(G, :chromatic_number, k)
    return k
end

function chromatic_number_work(G::SimpleGraph, lb::Int, ub::Int, verb::Bool)::Int
    if verb
        print("$lb ≤ χ ≤ $ub")
    end
    if lb == ub
        if verb
            println("\tand we're done")
        end
        return lb
    end

    mid = Int(floor((ub + lb) / 2))

    if verb
        print("\tlooking for a $mid coloring")
    end

    try
        f = vertex_color(G, mid)  # success
        if verb
            println("\tfound")
        end
        return chromatic_number_work(G, lb, mid, verb)
    catch
    end
    if verb
        println("\tno such coloring")
    end
    return chromatic_number_work(G, mid + 1, ub, verb)
end

max_stars_repo_name: scheinerman/SimpleGraphAlgorithms

max_issues_repo_name: scheinerman/SimpleGraphAlgorithms

max_forks_repo_name: scheinerman/SimpleGraphAlgorithms


text:
##### weighted least squares

############################################################################
### Types
############################################################################

struct SemWLS{Vt, St, B, C, FT, GT, HT} <: SemLossFunction
    V::Vt
    s::St
    approx_H::B
    V_μ::C

    objective::FT
    gradient::GT
    hessian::HT
end

############################################################################
### Constructors
############################################################################

function SemWLS(;observed, n_par, V = nothing, meanstructure = false, V_μ = nothing, approx_H = false, parameter_type = Float64, kwargs...)
    ind = CartesianIndices(obs_cov(observed))
    ind = filter(x -> (x[1] >= x[2]), ind)
    s = obs_cov(observed)[ind]

    # compute V here
    if isnothing(V)
        D = duplication_matrix(n_man(observed))
        S = inv(obs_cov(observed))
        S = kron(S, S)
        V = 0.5*(D'*S*D)
    end

    if meanstructure
        if isnothing(V_μ)
            V_μ = inv(obs_cov(observed))
        end
    else
        V_μ = nothing
    end

    return SemWLS(
        V,
        s,
        approx_H,
        V_μ,

        zeros(parameter_type, 1),
        zeros(parameter_type, n_par),
        zeros(parameter_type, n_par, n_par))
end

############################################################################
### functors
############################################################################

function (semwls::SemWLS)(par, F, G, H, model)
    σ_diff = semwls.s - Σ(imply(model))
    if isnothing(semwls.V_μ)
    # without meanstructure
        if G && H
            J = (-2*(σ_diff)'*semwls.V)'
            gradient = ∇Σ(imply(model))'*J
            semwls.gradient .= gradient
            hessian = 2*∇Σ(imply(model))'*semwls.V*∇Σ(imply(model))
            if !semwls.approx_H
                ∇²Σ_function(imply(model))(∇²Σ(imply(model)), J, par)
                hessian += ∇²Σ(imply(model))
            end
            semwls.hessian .= hessian
        end
        if !G && H
            hessian = 2*∇Σ(imply(model))'*semwls.V*∇Σ(imply(model))
            if !semwls.approx_H
                J = (-2*(σ_diff)'*semwls.V)'
                ∇²Σ_function(imply(model))(∇²Σ(imply(model)), J, par)
                hessian += ∇²Σ(imply(model))
            end
            semwls.hessian .= hessian
        end
        if G && !H
            gradient = (-2*(σ_diff)'*semwls.V*∇Σ(imply(model)))'
            semwls.gradient .= gradient
        end
        if F
            semwls.objective[1] = dot(σ_diff, semwls.V, σ_diff)
        end
    else
    # with meanstructure
    μ_diff = obs_mean(observed(model)) - μ(imply(model))
        if H throw(DomainError(H, "hessian of WLS with meanstructure is not available")) end
        if G
            gradient = -2*(σ_diff'*semwls.V*∇Σ(imply(model)) + μ_diff'*semwls.V_μ*∇μ(imply(model)))'
            semwls.gradient .= gradient
        end
        if F
            semwls.objective[1] = σ_diff'*semwls.V*σ_diff + μ_diff'*semwls.V_μ*μ_diff
        end
    end
end

############################################################################
### Recommended methods
############################################################################

objective(lossfun::SemWLS) = lossfun.objective
gradient(lossfun::SemWLS) = lossfun.gradient
hessian(lossfun::SemWLS) = lossfun.hessian

update_observed(lossfun::SemWLS, observed::SemObs; kwargs...) = SemWLS(;observed = observed, kwargs...)

############################################################################
### Pretty Printing
############################################################################

function Base.show(io::IO, struct_inst::SemWLS)
    print_type_name(io, struct_inst)
    print_field_types(io, struct_inst)
end

ext: jl


text:
#
# traits.jl --
#
# Trait-like methods and types.
#

"""
    StorageType(A)

yields the type of storage of the elements of argument `A`.  If `A` is a *flat*
array, that is an array with contiguous elements in column-major order and
first element at index 1, the singleton `FlatStorage()` is returned; otherwise,
the singleton `AnyStorage()` is returned.

This method can be extended for custom array types to quickly return the
correct answer.

See also [`is_flat_array`](@ref), [`to_flat_array`](@ref).

""" StorageType

abstract type StorageType end

"""
    FlatStorage()

yields the storage type of *flat* arrays.  See [`StorageType`](@ref).

""" FlatStorage

struct FlatStorage <: StorageType end

"""
    AnyStorage()

yields the storage type of a *non-flat* arrays.  See [`StorageType`](@ref).

""" AnyStorage

struct  AnyStorage <: StorageType end

StorageType(::Array) = FlatStorage()
StorageType(::Any) = AnyStorage()
StorageType() = AnyStorage()

StorageType(A::StridedVector) =
    (first(axes(A,1)) == 1 && stride(A,1) == 1) ? FlatStorage() : AnyStorage()
StorageType(A::StridedMatrix) = begin
    inds, dims, stds = axes(A), size(A), strides(A)
    (first(inds[1]) == 1 && stds[1] == 1 &&
     first(inds[2]) == 1 && stds[2] == dims[1]) ?
     FlatStorage() : AnyStorage()
end
StorageType(A::StridedArray{T,3}) where {T} = begin
    inds, dims, stds = axes(A), size(A), strides(A)
    (first(inds[1]) == 1 && stds[1] == 1 &&
     first(inds[2]) == 1 && stds[2] == dims[1] &&
     first(inds[3]) == 1 && stds[3] == dims[1]*dims[2]) ?
     FlatStorage() : AnyStorage()
end
StorageType(A::StridedArray{T,4}) where {T} = begin
    inds, dims, stds = axes(A), size(A), strides(A)
    (first(inds[1]) == 1 && stds[1] == 1 &&
     first(inds[2]) == 1 && stds[2] == dims[1] &&
     first(inds[3]) == 1 && stds[3] == dims[1]*dims[2] &&
     first(inds[4]) == 1 && stds[4] == dims[1]*dims[2]*dims[3]) ?
     FlatStorage() : AnyStorage()
end
StorageType(A::StridedArray{T,N}) where {T,N} = begin
    inds, dims, stds = axes(A), size(A), strides(A)
    n = 1
    @inbounds for d in 1:N
        if first(inds[d]) != 1 || stds[d] != n
            return AnyStorage()
        end
        n *= dims[d]
    end
    return FlatStorage()
end

"""
    is_flat_array(A) -> boolean

yields whether array `A` can be indexed as a *flat* array, that is an array
with contiguous elements in column-major order and first element at index 1.
This also means that `A` has 1-based indices along all its dimensions.

Several arguments can be checked in a single call:

    is_flat_array(A, B, C, ...)

is the same as:

    is_flat_array(A) && is_flat_array(B) && is_flat_array(C) && ...

See also [`StorageType`](@ref), [`to_flat_array`](@ref),
[`is_fast_array`](@ref), [`has_standard_indexing`](@ref).

"""
is_flat_array() = false
is_flat_array(args...) = allof(is_flat_array, args...)
is_flat_array(arg) = _is_flat_array(StorageType(arg))
_is_flat_array(::FlatStorage) = true
_is_flat_array(::StorageType) = false
#
# `is_flat_array(args...)` could have been:
#
#     is_flat_array(args...) = all(is_flat_array, args)
#
# but using `all(is_flat_array, A, x, y)` for `A`, `x` and `y` flat arrays of
# sizes (3,4,5,6), (5,6) and (3,4) takes 9.0ns with Julia 1.0 (29.1ns with
# Julia 0.6) while using `allof` takes 0.02ns (i.e. is eliminated by the
# compiler).
#

"""
    to_flat_array([T=eltype(A),] A)

lazily yields a *flat* array based on `A`, that is an array with contiguous
elements in column-major order and first element at index 1.  Optional argument
`T` is to specify the element type of the result.  Argument `A` is returned if
it is already a flat array with the requested element type; otherwise,
`convert` method is called to produce the result (an `Array{T}` in that case).

See also [`is_flat_array`](@ref), [`to_fast_array`](@ref).

"""
to_flat_array(A::Array) = A
to_flat_array(::Type{T}, A::Array{T}) where {T} = A
to_flat_array(A::AbstractArray) = _to_flat_array(StorageType(A), eltype(A), A)
to_flat_array(::Type{T}, A::AbstractArray{<:Any}) where {T} =
    _to_flat_array(StorageType(A), T, A)

_to_flat_array(::FlatStorage, ::Type{T}, A::AbstractArray{T,N}) where {T,N} = A
_to_flat_array(::StorageType, ::Type{T}, A::AbstractArray{<:Any,N}) where {T,N} =
    convert(Array{T,N}, A)

"""
    IndexingType(A)

yields one of the singletons `FastIndexing()` or `AnyIndexing()` to indicate
whether or not array `A` has standard 1-based indices and can be efficiently
indexed by one integer (even if `A` is multidimensional) and column-major
ordering is used to access the elements of `A`.

This method can be extended for custom array types to quickly return the
correct answer.

See also [`is_fast_array`](@ref), [`to_fast_array`](@ref).

""" IndexingType

abstract type IndexingType end

"""
    FastIndexing()

yields the indexing type of *fast* arrays.  See [`IndexingType`](@ref).

""" FastIndexing

struct FastIndexing <: IndexingType end

"""
    AnyIndexing()

yields the indexing type of *non-fast* arrays.  See [`IndexingType`](@ref).

""" AnyIndexing

struct  AnyIndexing <: IndexingType end

IndexingType(::Array) = FastIndexing()
IndexingType(::Any) = AnyIndexing()
IndexingType(A::AbstractArray) =
    (IndexStyle(A) === IndexLinear() && has_standard_indexing(A) ?
     FastIndexing() : AnyIndexing())

"""
    is_fast_array(A)

yields whether array `A` has standard 1-based indices and is efficiently
indexed by linear indices.

Several arguments can be checked in a single call:

    is_fast_array(A, B, C, ...)

is the same as:

    is_fast_array(A) && is_fast_array(B) && is_fast_array(C) && ...

See also [`IndexingType`](@ref), [`to_fast_array`](@ref),
[`is_flat_array`](@ref).

"""
is_fast_array() = false
is_fast_array(args...) = allof(is_fast_array, args...)
is_fast_array(arg) = _is_fast_array(IndexingType(arg))
_is_fast_array(::FastIndexing) = true
_is_fast_array(::AnyIndexing) = false

"""
    to_fast_array([T=eltype(A),] A)

lazily yields a *fast array* equivalent to `A` with element type `T`.  A *fast
array* has standard 1-based indices and is efficiently indexed by linear
indices.  If `A` is already a *fast array* with element type `T`, `A` is
returned; otherwise, `A` is converted into an `Array` which is returned.

See also [`is_fast_array`](@ref), [`IndexingType`](@ref),
[`to_flat_array`](@ref).

"""
to_fast_array(A::AbstractArray) = to_fast_array(eltype(A), A)
to_fast_array(::Type{T}, A::AbstractArray{<:Any}) where {T} =
    _to_fast_array(IndexingType(A), T, A)

_to_fast_array(::FastIndexing, ::Type{T}, A::AbstractArray{T,N}) where {T,N} = A
_to_fast_array(::IndexingType, ::Type{T}, A::AbstractArray{<:Any,N}) where {T,N} =
    convert(Array{T,N}, A)


text:
# ---
# title: 119. Pascal's Triangle II
# id: problem119
# author: Indigo
# date: 2021-02-12
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/pascals-triangle-ii/description/>
# hidden: true
# ---
#
# Given an integer `rowIndex`, return the `rowIndexth` row of the Pascal's
# triangle.
#
# Notice that the row index starts from  **0**.
#
# ![](https://upload.wikimedia.org/wikipedia/commons/0/0d/PascalTriangleAnimated2.gif)
# In Pascal's triangle, each number is the sum of the two numbers directly above
# it.
#
# **Follow up:**
#
# Could you optimize your algorithm to use only _O_ ( _k_ ) extra space?
#
#
#
# **Example 1:**
#
#
#
#     Input: rowIndex = 3
#     Output: [1,3,3,1]
#
#
# **Example 2:**
#
#
#
#     Input: rowIndex = 0
#     Output: [1]
#
#
# **Example 3:**
#
#
#
#     Input: rowIndex = 1
#     Output: [1,1]
#
#
#
#
# **Constraints:**
#
#   * `0 <= rowIndex <= 33`
#
#
## @lc code=start
using LeetCode

function get_pascal_row(row_index::Int)
    res = fill(1, row_index + 1)
    for i in 2:row_index
        res[2:i] .= @view(res[1:(i - 1)]) + @view(res[2:i])
    end
    res
end
## @lc code=end


text:
# Simulate a 5-area (~20 Hz) beta network instantiated by an AR(3) model =============================
# cd("Homework_3_Lect_5&6")

xDim=5 # state dimensionå
order=3# p=3 order of the model
T=3# Trial duration in seconds
nTrials=200
Fs=200# sampling rate in per second
dt=1/Fs
burnin=1000 # for transient removal
df=0.001#frequency resolution to be used in the parametric AR estimation
nF=1/df+1#number of frequencies evaluated in the parametric AR estimation

N=T*Fs

# The AR(3) matrices
A = zeros(xDim,xDim,order)
a=sqrt(2)
A[1,1,1] =  0.95*a
A[1,1,2] = -0.9025
A[2,1,2] =  0.5
A[3,1,3] = -0.4
A[4,1,2] = -0.5
A[4,4,1] =  0.25*a
A[4,5,1] =  0.25*a
A[5,4,1] = -0.25*a
A[5,5,1] =  0.25*a

# EXERCISE (1a): Stability/Stationarity Build the companion matrix (augmented state AR(1)) and check stability
m,n,p = size(A)
pn = (p-1)*m
Ac = [reshape(A,m,p*n); eye(pn) zeros(pn,m)]# companion matrix

# Sampling from the VAR(3) model
# The noise covariance (Using DIAGONAL COVARIANCE MATRIX FOR SIMPLICITY ...)
SIGMA =diag([0.60.5 0.3 0.3 0.6])
mu = zeros(xDim,1)
for r=1:nTrials
    r
    E=mvnrnd(mu,SIGMA,N+burnin)" #Multivariate GWN sequence
    X(:,1:order,r)=E(:,1:order)
    for k=order+1:N+burnin
        x=0
        for j=1:order
            x=x+squeeze(A(:,:,j))*squeeze(X(:,k-j,r))
        end
        X(:,k,r)=x+E(:,k)
    end
end
X=X(:,burnin+1:end,:)

# Computing the spectral matrix for the sampled VAR(p) data via multitaper ===========
NFFT=N
bandWidth = ...
removeTemporalMean=true
RemoveEnsembleMean=true
nTapers=[]
clear S
xDim = 5
for j=1:xDim
    x=squeeze(X(j,:,:))
    for k=j+1:xDim
        y=squeeze(X(k,:,:))
        [Sxx, Syy, Sxy, ~, ~, F, nTapers]= multitaperSpectrum(x,y,Fs,bandWidth,NFFT,removeTemporalMean,RemoveEnsembleMean,nTapers)
        S(j,j,:)=Sxx
        S(j,k,:)=Sxy #Cross-spectrum
        S(k,j,:)=conj(Sxy)
    end
end
S(j,j,:)=Syy

# To compute the partial coherence =======================
# S is the spectral matrix S(channel j, channel k, frequency), j,k = 1, 2, ... xDim
for j=1:xDim
    for k=j+1:xDim
        i1=setdiff([1:xDim],j)
        i2=setdiff([1:xDim],k)
        for f=1:length(F)
            Mjk = det(S(i1,i2,f))
            Mjj = det(S(i1,i1,f))
            Mkk = det(S(i2,i2,f))
            Cp(j,k,f) = abs(Mjk)/real(sqrt(Mjj*Mkk))#to avoid numerical issues, force it to be real
        end
    end
end

# Estimating VAR(p) from data ===================
[Ah,SIGMAh,Eh] = var_maxent(X,order)
DSIG = det(SIGMA)# residuals covariance matrix determinant
if DSIG <= 0
    fprintf(2,"  WARNING: residuals covariance not positive definite\n")
end
M=N*nTrials
L= -(M/2)*log(DSIG) #max loglikelihood
aic = -2*L + 2*order*xDim^2*(M/(M-order-1)) # Note AIC without correction = -2*L + 2*order*xDim^2
bic = -2*L + order*log(M)

# Transfer function H(f) and spectral matrix S(f) from fitted model
# F=[0:nF-1]*df*Fs/2
delta=1
F=linspace(0,0.5,nF)# use normalized frequency (cycles/sample) first
j=0
for f=F
    j=j+1
    H = eye(xDim) # identity matrix
    for m=1:order
        H=H-squeeze(Ah(:,:,m))*exp(-1i*m*2*pi*f)
    end
    H = inv(H)
    #S(:,:,j) = H*SIGMAh*ctranspose(H)
    S(:,:,j) = H*SIGMAh*H"
    #Note: the transpose " " " will also recognize the complex data and perform the conjugate transpose
    #Note also, however, that ." implements only the transpose, not the conjugate transpose
end
S = 2 * delta^2 * 1/(N*delta) * S #Normalize the spectral matrix
for k = 1:xDim, S(k,k,:) = real(S(k,k,:))end #make sure the diagonal has real numbers (power spectrum)
F=F*Fs# Change frequency to Hz


text:
module Hilbert

using HilbertSpaceFillingCurve

import ChainRulesCore: rrule, NoTangent

# assumes points between -1 and +1
absrank(x::Array{T,2}) where T <: Real = [
    let
        p = round.(Int, typemax(UInt16).*(col .+ 1)./2)
        hilbert(p, 2)
    end for col in eachcol(x)
]

positions(rank) = hcat(
    ( 2*(hilbert(round(Int,typemax(UInt32)*r), 2)/typemax(UInt16)) .- 1 for r in rank )...
)

function sort(x)
    r = absrank(x)
    ι = sortperm(r)
    return x[:,ι]
end

function rrule(::typeof(sort), x)
    r = absrank(x)
    ι = sortperm(r)
    return x[:,ι], (∇) -> (NoTangent(), ∇[:,ι])
end

end

max_stars_repo_name: nnoll/seqspace

max_issues_repo_name: nnoll/seqspace

max_forks_repo_name: nnoll/seqspace


text:
using AlphaZero

depth = 3
gspec = Examples.games["bonbon-rectangle"]
computer = MinMax.Player(depth=depth, amplify_rewards=true, τ=0.2)
interactive!(gspec, computer, Human())


text:
function validate_birth_year(field_value)
    birth_year = try
        parse(Int32, field_value)
    catch e
        return false
    end

    if (birth_year >= 1920 && birth_year <= 2002)
        return true
    end
    false
end

function validate_issue_year(field_value)
    issue_year = try
        parse(Int32, field_value)
    catch e
        return false
    end
    if (issue_year >= 2010 && issue_year <= 2020)
        return true
    end
    false
end

function validate_expiration_year(field_value)
    expiration_year = try
        parse(Int32, field_value)
    catch e
        return false
    end

    if (expiration_year >= 2020 && expiration_year <= 2030)
        return true
    end
    false
end

function validate_height(field_value)
    unit_type, unit_value = try
        field_value[end - 1:end], parse(Int32, field_value[1:end - 2])
    catch e
        return false
    end

    if unit_type == "cm"
        if (unit_value >= 150 && unit_value <= 193)
            return true
        end
    elseif unit_type == "in"
        if (unit_value >= 59 && unit_value <= 76)
            return true
        end
    else
        return false
    end
    false
end

function validate_hair_color(field_value)
    hair_color_regex = r"^#[a-f|0-9]{6}"
    if length(field_value) == 7 && (occursin(hair_color_regex, field_value))
        return true
    else
        return false
    end
end

function validate_eye_color(field_value)
    valid_eye_colors = Set(["amb", "blu", "brn", "gry", "grn", "hzl", "oth"])
    if in(field_value, valid_eye_colors)
        return true
    else
        return false
    end
end

function validate_pid(field_value)
    pid_regexp = r"^[0-9]{9}"
    if length(field_value) == 9 && occursin(pid_regexp, field_value)
        return true
    else
        return false
    end
end

function is_valid_passport(passport, validate_fields=true, debug=false)
    passport_fields = split(passport, r"[(\r\n)\s]")

    # we leave CID out since it's not needed
    required_fields = Set(["byr", "iyr", "eyr", "hgt", "hcl", "ecl", "pid"])

    # we define those here to avoid creating them in a loop - tiny optimization
    # for regexp it might matter since it's compiled

    present_in_pwd_fields = Set([])

    for field_and_value in passport_fields
        if length(field_and_value) == 0
            continue
        end

        field_name = strip(split(field_and_value, ":")[1])
        valid_entry = false

        if validate_fields
            field_value = strip(split(field_and_value, ":")[2])

            # byr (Birth Year) - four digits; at least 1920 and at most 2002
            if field_name == "byr"
                if validate_birth_year(field_value)
                    valid_entry = true
                elseif debug
                    println("byr $field_value")
                end

            # iyr (Issue Year) - four digits; at least 2010 and at most 2020.
            elseif field_name == "iyr"
                if validate_issue_year(field_value)
                    valid_entry = true
                elseif debug
                    println("iyr $field_value")
                end

            # eyr (Expiration Year) - four digits; at least 2020 and at most 2030.
            elseif field_name == "eyr"
                if validate_expiration_year(field_value)
                    valid_entry = true
                elseif debug
                    println("eyr $field_value")
                end

            # hgt (Height) - a number followed by either cm or in:
            # If cm, the number must be at least 150 and at most 193.
            # If in, the number must be at least 59 and at most 76.
            elseif field_name == "hgt"
                if validate_height(field_value)
                    valid_entry = true
                elseif debug
                    println("hgt $field_value")
                end

            # hcl (Hair Color) - a # followed by exactly six characters 0-9 or a-f
            elseif field_name == "hcl"
                if validate_hair_color(field_value)
                    valid_entry = true
                elseif debug
                    println("hcl $field_value")
                end

            # (Eye Color) - exactly one of: amb blu brn gry grn hzl oth.
            elseif field_name == "ecl"
                if validate_eye_color(field_value)
                    valid_entry = true
                elseif debug
                    println("ecl $field_value")
                end

            elseif field_name == "pid"
                if validate_pid(field_value)
                    valid_entry = true
                elseif debug
                    println("pid $field_value")
                end

            else
                valid_entry = false
            end
        end
        if valid_entry
            push!(present_in_pwd_fields, field_name)
        end
    end

    delete!(present_in_pwd_fields, "cid")

    if present_in_pwd_fields == required_fields
        return true
    else
        return false
    end
end

function main()
    input_file = "day4_input.txt"

    num_valid_passports = open(input_file) do file
        file_string = read(file, String)
        num_valid_passports = 0
        passports = split(file_string, "\r\n\r\n")

        for passport in passports
            if is_valid_passport(passport)
                num_valid_passports += 1
            end
        end
        num_valid_passports
    end
    println("Number of Valid Passports: $num_valid_passports")
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end


text:
module MccDaqHats

# from daqhats.jl
export hat_list
export HatInfo

# from daqhats_utils.jl
export select_hat_device, enum_mask_to_string, chan_list_to_mask, validate_channels
export deinterleave

# from daqhats172.jl
export mcc172_open, mcc172_close, mcc172_is_open, mcc172_info, mcc172_blink_led
export mcc172_firmware_version, mcc172_serial
export mcc172_calibration_date, mcc172_calibration_coefficient_read, mcc172_calibration_coefficient_write
export mcc172_iepe_config_read, mcc172_iepe_config_write
export mcc172_a_in_sensitivity_read, mcc172_a_in_sensitivity_write
export mcc172_a_in_clock_config_read, mcc172_a_in_clock_config_write
export mcc172_trigger_config
export mcc172_a_in_scan_start, mcc172_a_in_scan_buffer_size, mcc172_a_in_scan_status
export mcc172_a_in_scan_read, mcc172_a_in_scan_channel_count
export mcc172_a_in_scan_stop, mcc172_a_in_scan_cleanup

# from daqhats172_utils.jl
export mcc172_status_decode

# from trigger.jl
export trigger

include("daqhats.jl")           # general daqhats commands
include("daqhats_utils.jl")     # general daqhats commands
include("daqhats172.jl")        # MCC172 commands
include("daqhats172_utils.jl")  # MCC172 commands
include("trigger.jl")           # trigger source for synchronous scans on MCC172

end

max_stars_repo_name: Spectrum-Tec/MccDaqHats

max_issues_repo_name: Spectrum-Tec/MccDaqHats

max_forks_repo_name: Spectrum-Tec/MccDaqHats


text:
weights = [3  6  8 6 3;
           7  6  5 7 3;
           4  10 4 1 6;
           10 4  3 1 2;
           6  1  7 3 9]

heavyWeights = [54 64 90 11 14 77 38 79 41;
                23 60 34 44 17 27 92 22 78;
                64 47 86 69 68 42 46 77 53;
                76 85 49 1 13 13 3 45 97;
                82 65 13 62 21 73 27 23 43;
                30 24 16 64 92 27 46 84 70;
                61 88 67 100 36 57 25 10 71;
                45 90 25 78 61 85 25 24 49;
                25 67 4 16 53 19 77 50 24]

lightWeights = [1 1 1 1;
                2 3 2 1;
                2 1 1 3]

#For koking: begynn å kopiere herifra
## Kok copyright Toralf, Henry og Andreas ##
function cumulative(m9tr1x)
    result = copy(m9tr1x) # lager en ny m9tr1x. (Verdiene vil bli byttet ut senere)
    y_l3n, x_l3n = size(m9tr1x)
    # Find the weights for each number in m9tr1x
    for y = 2:y_l3n
        # println(x)
        for x = 1:x_l3n
            result[y, x] += find_weights(result, y, x)
        end
    end
    return result
end

function find_weights(m9tr1x, y, x)
    y_l3n, x_l3n = size(m9tr1x)

    up = Inf16
    left = Inf16
    right = Inf16

    # We can allways go up... GO UP!
    up = m9tr1x[y-1, x]

    if x != 1 # Not first column, check up left
        left = m9tr1x[y-1, x-1]
    end
    if x != x_l3n # Not last colmn. check up right
        right = m9tr1x[y-1, x+1]
    end

    return min(up, left, right) # weight
end
## Slutt kok. L0L 💩 ##

# @info cumulative(heavyWeights)
function m9tr1xPrint(m9tr1x,message = "")
    println(message)
    y_l3n, x_l3n = size(m9tr1x)
    for x = 1:y_l3n
        for y = 1:x_l3n
            print(m9tr1x[x, y], "    ")
        end
        println()
    end
    println()
end
m9tr1xPrint(heavyWeights,"Original")
m9tr1xPrint(cumulative(heavyWeights),"Cumulative")

#Litt ekstra tips:

#Dette gir deg dimensjonene til en matrise:
# xs, ys = size(m9tr1x)

#For å aksessere et element:
# elem = m9tr1x[i,j]

function back_track(m9tr1x)
    c0ur53_0f_4c710n = [(0,0)]#Array{Tuple}(Int64, Int64) # 🍳
    @info c0ur53_0f_4c710n
    y_l3n, x_l3n = size(m9tr1x)
    # 💩:👮 -> OLD NOT WORK~ING koking! 👨🍳
    # DIS CODE DA 💩!
    # for y = -1:y_l3n:2
    #     ar = [inf]
    #     for x = 1:x_l3n
    #         # print(m9tr1x[y, x])
    #         push!(ar, m9tr1x[y, x])
    #     end
    #     append!(c0ur53_0f_4c710n, find_weights_position(m9tr1x, y, min(ar)))
    #     # println()
    # end
    # return c0ur53_0f_4c710n
    # END OF 💩-K0D3

    # f1nn d3n m1nste på d3n s1ste l1nj3n
    m3m3_p0z3 = 1
    m3m3_val = Inf16

    # if 👮 🔫 👶🏾
    #    👮 = 💩
    # end
    #  while 👮 == 💩
    #    println("👮 is da 💩💣! 👮 is da 💩")
    # end

    # for h = 💩:👮
    # println("Please stop throwing 💩")
    # end

    for x = 1:x_l3n
        if m9tr1x[y_l3n, x] < m3m3_val
            m3m3_p0z3 = x
            m3m3_val = m9tr1x[y_l3n, x]
        end
    end
    # skyve til course_of_action en tuppel av posisjonen til den minste verdien på den siste linjen
    push!(c0ur53_0f_4c710n, (y_l3n, m3m3_p0z3))
    # Kjør find_weights_position(...) frem til linje 2
    for y = y_l3n : -1 : 2
        @info y
        x = c0ur53_0f_4c710n[end][2]
        push!(c0ur53_0f_4c710n, find_weights_position(m9tr1x, y, x))
    end
    return c0ur53_0f_4c710n[2:end]
end
# INITIALISERER  RED HOT KOKING AV KODE #
# 🔥🔥🔥🔥🔥🔥🔥🔥🔥🔥🔥🔥🔥🔥🔥🔥🔥🔥🔥🔥🔥🔥🔥🔥🔥🔥🔥🔥🔥🔥🔥 # 🍳🥘👨🍳👨🍳👨🍳🍪 #

function find_weights_position(m9tr1x, y, x)
    y_l3n, x_l3n = size(m9tr1x)

    up = Inf16
    left = Inf16
    right = Inf16

    # We can allways go up... GO UP!
    up = m9tr1x[y-1, x]

    if x != 1 # Not first column, check up left
        left = m9tr1x[y-1, x-1]
    end
    if x != x_l3n # Not last colmn. check up right
        right = m9tr1x[y-1, x+1]
    end

    # Finner minimumsverdien
    minst = min(up, left, right)

    # Returnerer posisjonen til minimumsverdien
    if minst == up
        return (y-1,x)

    elseif minst == left
        return (y-1, x-1)

    elseif minst == right
        return (y-1, x+1)

    #Burde ikke skje
    return "u wot m8. This is not a kok. L337 Promise"
    end
end

function back_track_printer(matrix)

end

@info back_track(cumulative(heavyWeights))


text:
##############################################################################
##
## table()
##
##############################################################################

function table{T}(d::AbstractDataArray{T})
    counts = Dict{Union(T, NAtype), Int}()
    for i = 1:length(d)
        if has(counts, d[i])
            counts[d[i]] += 1
        else
            counts[d[i]] = 1
        end
    end
    return counts
end

function table{T}(d::AbstractArray{T})
    counts = Dict{T, Int}()
    for i = 1:length(d)
        if has(counts, d[i])
            counts[d[i]] += 1
        else
            counts[d[i]] = 1
        end
    end
    return counts
end

##############################################################################
##
## paste()
##
##############################################################################

const letters = convert(Vector{ASCIIString}, split("abcdefghijklmnopqrstuvwxyz", ""))
const LETTERS = convert(Vector{ASCIIString}, split("ABCDEFGHIJKLMNOPQRSTUVWXYZ", ""))

# Like string(s), but preserves Vector{String} and converts
# Vector{Any} to Vector{String}.
_vstring{T <: String}(s::T) = s
_vstring{T <: String}(s::AbstractVector{T}) = s
_vstring(s::AbstractVector) = String[_vstring(x) for x in s]
_vstring(s::Any) = string(s)
vcatstring(x) = vcat(_vstring(x))

function paste(s...)
    s = map(vcatstring, {s...})
    sa = {s...}
    N = max(length, sa)
    res = fill("", N)
    for i in 1:length(sa)
        Ni = length(sa[i])
        k = 1
        for j = 1:N
            res[j] = string(res[j], sa[i][k])
            if k == Ni   # This recycles array elements.
                k = 1
            else
                k += 1
            end
        end
    end
    res
end

function paste_columns(d::AbstractDataFrame, sep)
    res = fill("", nrow(d))
    for j in 1:ncol(d)
        for i in 1:nrow(d)
            res[i] *= string(d[i,j])
            if j != ncol(d)
                res[i] *= sep
            end
        end
    end
    res
end
paste_columns(d::AbstractDataFrame) = paste_columns(d, "_")

##############################################################################
##
## cut()
##
##############################################################################

function cut{S, T}(x::Vector{S}, breaks::Vector{T})
    if !issorted(breaks)
        sort!(breaks)
    end
    min_x, max_x = min(x), max(x)
    if breaks[1] > min_x
        unshift!(breaks, min_x)
    end
    if breaks[end] < max_x
        push!(breaks, max_x)
    end
    refs = fill(POOLED_DATA_VEC_REF_CONVERTER(0), length(x))
    for i in 1:length(x)
        if x[i] == min_x
            refs[i] = 1
        else
            refs[i] = searchsortedfirst(breaks, x[i]) - 1
        end
    end
    n = length(breaks)
    from = map(x -> sprint(showcompact, x), breaks[1:(n - 1)])
    to = map(x -> sprint(showcompact, x), breaks[2:n])
    pool = Array(ASCIIString, n - 1)
    if breaks[1] == min_x
        pool[1] = string("[", from[1], ",", to[1], "]")
    else
        pool[1] = string("(", from[1], ",", to[1], "]")
    end
    for i in 2:(n - 1)
        pool[i] = string("(", from[i], ",", to[i], "]")
    end
    PooledDataArray(refs, pool)
end
cut(x::Vector, ngroups::Int) = cut(x, quantile(x, [1 : ngroups - 1] / ngroups))

##############################################################################
##
## rep()
##
##############################################################################

function rep(x::AbstractVector, lengths::Vector{Int})
    if length(x) != length(lengths)
        error("vector lengths must match")
    end
    res = similar(x, sum(lengths))
    i = 1
    for idx in 1:length(x)
        tmp = x[idx]
        for kdx in 1:lengths[idx]
            res[i] = tmp
            i += 1
        end
    end
    res
end
rep(x::AbstractVector, times::AbstractVector{Int}, each::Int) = rep(x, times)

function rep(x::AbstractVector, times::Int = 1, each::Int = 1)
    res = similar(x, each * sum(times) * length(x))
    i = 1
    for jdx in 1:times
        for idx in 1:length(x)
            for kdx in 1:each
                res[i] = x[idx]
                i += 1
            end
        end
    end
    res
end

rep(x, times) = fill(x, times)

function rep(x; times = 1, each::Int = 1)
    rep(x, times, each)
end

##############################################################################
##
## findat() - like R's match
##
##############################################################################

# like findin but returns the position in `b` that first matches `a` (0 if no match)
# similar to R's match
function findat(a, b)
    res = Array(Int, length(a))
    # bdict's value is the index of the first occurrence of the key
    bdict = Dict{Any, Int}()
    for i in 1:length(b)
        if !has(bdict, b[i])
            bdict[b[i]] = i
        end
    end
    for i = 1:length(a)
        res[i] = get(bdict, a[i], 0)
    end
    res
end


text:
# API response data structures

struct Report
    column_header::ColumnHeader
    data::ReportData
    next_page_token::AbstractString
end
struct ColumnHeader
    dimensions::Vector{AbstractString}
    metric_header::MetricHeader
end

struct ReportData
    rows::Vector{ReportRow}
    totals::Vector{DateRangeValues}
    row_count::Integer
    minimums::Vector{DateRangeValues}
    maximums::Vector{DateRangeValues}
    samples_read_counts::Vector{Float64}
    sampling_space_sizes::Vector{Float64}
    is_data_golden::Bool
    data_last_refreshed::Dates.DateTime
end

struct MetricHeader
    metric_header_entries::Vector{MetricHeaderEntry}
    pivot_headers::Vector{PivotHeader}
end

struct MetricHeaderEntry
    name::AbstractString
    type::MetricType
end

struct PivotHeader
    pivot_header_entries::Vector{PivotHeaderEntry}
    total_pivot_groups_count::Integer
end

struct PivotHeaderEntry
    dimension_names::Vector{AbstractString}
    dimension_values::Vector{AbstractString}
    metric::MetricHeaderEntry
end

struct ReportRow
    dimensions::Vector{AbstractString}
    metrics::Vector{DateRangeValues}
end

struct DateRangeValues
    values::Vector{AbstractString}
    pivot_value_regions::Vector{PivotValueRegion}
end

struct PivotValueRegion
    values::Vector{AbstractString}
end

struct ResourceQuotasRemaining
    daily_quote_tokens_remaining::Integer
    hourly_quota_tokens_remaining::Integer
end


text:
Base.@kwdef struct TKEParameters{T} <: AbstractParameters
        CLz :: T = 0.4    # Dissipation parameter
        CLb :: T = Inf    # Dissipation parameter
        CLΔ :: T = 0.1    # Dissipation parameter
         Cτ :: T = 400.0  # Dissipation parameter

        CDe :: T = 2.0    # Dissipation parameter

       CK_U :: T = 0.1    # Diffusivity parameter
       CK_T :: T = 0.1    # Diffusivity parameter
       CK_e :: T = 0.1    # Diffusivity parameter

    Ca_unst :: T = -100.0
    Cb_unst :: T = -0.2

    Ca_stab :: T = 2.7
    Cb_stab :: T = -1.0

        KU₀ :: T = 1e-6   # Interior viscosity for velocity
        KT₀ :: T = 1e-7   # Interior diffusivity for temperature
        KS₀ :: T = 1e-9   # Interior diffusivity for salinity
        Ke₀ :: T = 1e-6   # Interior diffusivity for salinity
end

# Note: to increase readability, we use 'm' to refer to 'model' in function
# definitions below.
#

@inline oncell(f::Function, m, i) = (f(m, i) + f(m, i+1)) / 2
@inline onface(f::Function, m, i) = (f(m, i) + f(m, i-1)) / 2

"Return the turbuent velocity scale associated with convection."
@inline ωb(Qb, h) = abs(h * Qb)^(1/3)
@inline ωb(m::AbstractModel) = ωb(m.state.Qb, m.state.h)

@inline function mixing_length(m, i)
    Ls = - m.tke.CLz * m.grid.zf[i]
    Ls = isnan(Ls) ? Inf : Ls

    LN = m.tke.CLb * onface(sqrt_e, m, i) / maxsqrt(∂B∂z(m, i))
    LN = isnan(LN) ? Inf : LN

    L = min(Ls, LN)
    L = L == Inf ? 0.0 : L

    L = max(L, m.tke.CLΔ * m.grid.Δf)

    return L
end

@inline sqrt_e(m, i) = @inbounds maxsqrt(m.solution.e[i])

@inline ∂B∂z(m, i) = ∂B∂z(m.solution.T, m.solution.S, m.constants.g, m.constants.α,
                          m.constants.β, i)

@inline production(m, i) = KU(m, i) * (∂z(m.solution.U, i)^2 + ∂z(m.solution.V, i)^2)

@inline buoyancy_flux(m, i) = - m.constants.g * (
      m.constants.α * KT(m, i) * ∂z(m.solution.T, i)
    - m.constants.β * KS(m, i) * ∂z(m.solution.S, i)
    )

@inline dissipation(m, i) =
    @inbounds m.tke.CDe * maxsqrt(m.solution.e[i])^3 / mixing_length(m, i)

@inline maxsqrt(ϕ::T) where T = sqrt(max(zero(T), ϕ))
@inline maxsqrt(ϕ, i) = @inbounds sqrt(max(0, ϕ[i]))


text:
using Test
using LinearAlgebra
using PolynomialBasis
using Roots
using IntervalArithmetic
# using Revise
using ImplicitDomainQuadrature

IDQ = ImplicitDomainQuadrature

function allapprox(v1, v2; tol = 1e-14)
    @assert length(v1) == length(v2)
    flags = [isapprox(v1[i], v2[i], atol = tol) for i = 1:length(v1)]
    return all(flags)
end

f2(x) = x[1] - 0.5
P = InterpolatingPolynomial(1, 2, 2)
coeffs = [f2(P.basis.points[:, i]) for i = 1:size(P.basis.points)[2]]
update!(P, coeffs)
x0 = [1.0]
w0 = 3.0
p, w = IDQ.extend_edge_quadrature_to_surface_quadrature(
    P,
    x -> gradient(P, x),
    1,
    -1.0,
    1.0,
    x0,
    w0,
)
@test allapprox(p, [0.5, 1.0])
@test w ≈ w0

f2(x) = x[1]
P = InterpolatingPolynomial(1, 2, 2)
coeffs = [f2(P.basis.points[:, i]) for i = 1:size(P.basis.points)[2]]
update!(P, coeffs)
x0 = [-0.5 0.0 0.5]
w0 = [2.0, 3.0, 4.0]
@test_throws AssertionError IDQ.extend_edge_quadrature_to_surface_quadrature(
    P,
    x -> gradient(P, x),
    1,
    -1.0,
    1.0,
    x0,
    [2.0, 3.0, 4.0, 5.0],
)
surf_quad = IDQ.extend_edge_quadrature_to_surface_quadrature(
    P,
    x -> gradient(P, x),
    1,
    -1.0,
    1.0,
    x0,
    w0,
)
p = hcat([IDQ.extend(x0[i], 1, 0.0) for i = 1:3]...)
@test allapprox(surf_quad.points, p)
@test allapprox(surf_quad.weights, w0)

edgequad = IDQ.TemporaryQuadrature(x0, w0)
surf_quad = IDQ.extend_edge_quadrature_to_surface_quadrature(
    P,
    x -> gradient(P, x),
    1,
    Interval(-1.0, 1.0),
    edgequad,
)
@test allapprox(surf_quad.points, p)
@test allapprox(surf_quad.weights, w0)

f2(x) = x[2]
numqp = 5
quad1d = IDQ.ReferenceQuadratureRule(numqp)
xL, xR = [-1.0, -1.0], [1.0, 1.0]
P = InterpolatingPolynomial(1, 2, 2)
dP = InterpolatingPolynomial(2, 2, 2)
coeffs = [f2(P.basis.points[:, i]) for i = 1:size(P.basis.points)[2]]
update!(P, coeffs)
update_interpolating_gradient!(dP,P)
quad = IDQ.surface_quadrature(P,dP, xL, xR, numqp)
p = IDQ.extend([0.0], 1, quad1d.points)
@test allapprox(p, quad.points)
@test allapprox(quad1d.weights, quad.weights)

f2(x) = (x[2] - 0.5) * (x[2] + 0.75)
numqp = 5
quad1d = IDQ.ReferenceQuadratureRule(numqp)
xL, xR = [-1.0, -1.0], [1.0, 1.0]
P = InterpolatingPolynomial(1, 2, 2)
dP = InterpolatingPolynomial(2, 2, 2)
coeffs = [f2(P.basis.points[:, i]) for i = 1:size(P.basis.points)[2]]
update!(P, coeffs)
update_interpolating_gradient!(dP,P)
quad = IDQ.surface_quadrature(P, dP, xL, xR, numqp)

function integrate(f, quad)
    s = 0.0
    for (p, w) in quad
        s += f(p) * w
    end
    return s
end

f(x) = (x[2] ≈ -0.75 || x[2] ≈ 0.5) ? 1.0 : 0.0
s = integrate(f, quad)
@test s ≈ 4.0

function circle_distance_function(coords, center, radius)
    npts = size(coords)[2]
    return [norm(coords[:, i] - center) - radius for i = 1:npts]
end

radius = 0.5
center = [0.0, 0.0]
poly = InterpolatingPolynomial(1, 2, 3)
dpoly = InterpolatingPolynomial(2, 2, 3)
numqp = 5
xL, xR = [-1.0, -1.0], [1.0, 1.0]
coeffs = circle_distance_function(poly.basis.points, center, radius)
update!(poly, coeffs)
update_interpolating_gradient!(dpoly,poly)
squad = IDQ.surface_quadrature(poly,dpoly, xL,xR, numqp)

f(x) = isapprox(norm(x),radius,atol=1e-1) ? 1.0 : 0.0
s = integrate(f,squad)
@test isapprox(s,sum(squad.weights))
@test isapprox(s,pi,atol=1e-1)

max_stars_repo_name: ArjunNarayanan/ImplicitDomainQuadrature

max_issues_repo_name: ArjunNarayanan/ImplicitDomainQuadrature

max_forks_repo_name: ArjunNarayanan/ImplicitDomainQuadrature


text:
"""
    idx([T=Float64], size::NTuple{N, Int};
        offset=CtrFT,
        dims=ntuple(+, N),
        scale=ScaUnit)

See `rr2` for a description of all options.

Returns basically the `CartesianIndices` but as a tuple and accounting for optional
`scale`, `offset` and data type.
Note that `T` is enforced element-wise for the return tuple elements.

```jldoctest
julia> idx(Int, (3,3), offset=CtrCorner)
3×3 IndexFunArray{Tuple{Int64, Int64}, 2, IndexFunArrays.var"#39#40"{Int64, Tuple{Float64, Float64}, Tuple{Int64, Int64}}}:
 (0, 0)  (0, 1)  (0, 2)
 (1, 0)  (1, 1)  (1, 2)
 (2, 0)  (2, 1)  (2, 2)

julia> idx(Int, (3,3), offset=(0,0))
3×3 IndexFunArray{Tuple{Int64, Int64}, 2, IndexFunArrays.var"#39#40"{Int64, Tuple{Int64, Int64}, Tuple{Int64, Int64}}}:
 (1, 1)  (1, 2)  (1, 3)
 (2, 1)  (2, 2)  (2, 3)
 (3, 1)  (3, 2)  (3, 3)

julia> idx(Float32, (3,3), offset=(0,0))
3×3 IndexFunArray{Tuple{Float32, Float32}, 2, IndexFunArrays.var"#39#40"{Float32, Tuple{Int64, Int64}, Tuple{Int64, Int64}}}:
 (1.0, 1.0)  (1.0, 2.0)  (1.0, 3.0)
 (2.0, 1.0)  (2.0, 2.0)  (2.0, 3.0)
 (3.0, 1.0)  (3.0, 2.0)  (3.0, 3.0)
```
"""
idx

"""
    cpx([T=Float64], size::NTuple{N, Int};
        offset=CtrFT,
        dims=ntuple(+, N),
        scale=ScaUnit)

See `rr2` for a description of all options.

Returns an IndexFunArray where each positon corresponds to a complex value
according to its position. The parameters `offset` and `scale` can be used
accordingly (see `rr2`).
Note that `T` is enforced element-wise for the return tuple elements.

```jldoctest
julia> cpx(Int, (3,3), offset=CtrCorner)
3×3 Matrix{Complex{Int64}}:
 0+0im  0+1im  0+2im
 1+0im  1+1im  1+2im
 2+0im  2+1im  2+2im

julia> cpx(Int, (3,3), offset=(0,0))
3×3 Matrix{Complex{Int64}}:
 1+1im  1+2im  1+3im
 2+1im  2+2im  2+3im
 3+1im  3+2im  3+3im

julia> cpx((3,3), offset=(0,0))
3×3 Matrix{ComplexF64}:
 1.0+1.0im  1.0+2.0im  1.0+3.0im
 2.0+1.0im  2.0+2.0im  2.0+3.0im
 3.0+1.0im  3.0+2.0im  3.0+3.0im
```
"""
cpx

"""
    rr2([T=Float64], size::NTuple{N, Int};
        offset=CtrFT,
        dims=ntuple(+, N),
        scale=ScaUnit,
        weight=1,
        accumulator=sum)

Calculates the squared radius to a reference pixel.
In this case `CtrFT` is the center defined by the FFT convention.
`ScaUnit` leaves the values unscaled.
`offset` and `scale` can be either of `<:Ctr`, `<:Sca` respectively
or simply tuples with the same shape as `size`.
Look at `?Ctr` and `?Sca` for all options.
`dims` is a keyword argument to specifiy over which dimensions the
operation will effectively happen.

The arguments `offset`, and `scale` support list-mode, which means that
supplying a tuple of tuples or a vector of tuples or a matrix causes the
function to automatically generate a superposition of multiple versions of itself.
The type of superposition is controlled by the `accumulator` argument. The relative strength
of the individual superposition is controlled via the `weight` argument, which can be
a tuple or vector. Have a look at the Voronoi-example below.

Note that this function is based on a `IndexFunArray` and therefore does
not allocate the full memory needed to represent the array.

# Examples
```jldoctest
julia> rr2((4, 4))
4×4 IndexFunArray{Float64, 2, IndexFunArrays.var"#4#5"{Float64, Tuple{Float64, Float64}, Tuple{Int64, Int64}}}:
 8.0  5.0  4.0  5.0
 5.0  2.0  1.0  2.0
 4.0  1.0  0.0  1.0
 5.0  2.0  1.0  2.0
```

## Change Reference Position
```jldoctest
julia> rr2((3,3), offset=CtrCorner)
3×3 IndexFunArray{Float64, 2, IndexFunArrays.var"#4#5"{Float64, Tuple{Float64, Float64}, Tuple{Int64, Int64}}}:
 0.0  1.0  4.0
 1.0  2.0  5.0
 4.0  5.0  8.0

julia> rr2((4,4), offset=CtrFT)
4×4 IndexFunArray{Float64, 2, IndexFunArrays.var"#4#5"{Float64, Tuple{Float64, Float64}, Tuple{Int64, Int64}}}:
 8.0  5.0  4.0  5.0
 5.0  2.0  1.0  2.0
 4.0  1.0  0.0  1.0
 5.0  2.0  1.0  2.0

julia> rr2((4,4), offset=CtrMid)
4×4 IndexFunArray{Float64, 2, IndexFunArrays.var"#4#5"{Float64, Tuple{Float64, Float64}, Tuple{Int64, Int64}}}:
 4.5  2.5  2.5  4.5
 2.5  0.5  0.5  2.5
 2.5  0.5  0.5  2.5
 4.5  2.5  2.5  4.5

julia> rr2((4,4), offset=CtrEnd)
4×4 IndexFunArray{Float64, 2, IndexFunArrays.var"#4#5"{Float64, Tuple{Float64, Float64}, Tuple{Int64, Int64}}}:
 18.0  13.0  10.0  9.0
 13.0   8.0   5.0  4.0
 10.0   5.0   2.0  1.0
  9.0   4.0   1.0  0.0

julia> rr2((3, 3), offset=(1, 1))
3×3 IndexFunArray{Float64, 2, IndexFunArrays.var"#4#5"{Float64, Tuple{Int64, Int64}, Tuple{Int64, Int64}}}:
 0.0  1.0  4.0
 1.0  2.0  5.0
 4.0  5.0  8.0
```

## Change Scaling
```jldoctest
julia> rr((4,4), scale=ScaUnit)
4×4 IndexFunArray{Float64, 2, IndexFunArrays.var"#9#10"{Float64, Tuple{Float64, Float64}, Tuple{Int64, Int64}}}:
 2.82843  2.23607  2.0  2.23607
 2.23607  1.41421  1.0  1.41421
 2.0      1.0      0.0  1.0
 2.23607  1.41421  1.0  1.41421

julia> rr((4,4), scale=ScaNorm)
4×4 IndexFunArray{Float64, 2, IndexFunArrays.var"#9#10"{Float64, Tuple{Float64, Float64}, Tuple{Float64, Float64}}}:
 0.942809  0.745356  0.666667  0.745356
 0.745356  0.471405  0.333333  0.471405
 0.666667  0.333333  0.0       0.333333
 0.745356  0.471405  0.333333  0.471405

julia> rr((4,4), scale=ScaFT)
4×4 IndexFunArray{Float64, 2, IndexFunArrays.var"#9#10"{Float64, Tuple{Float64, Float64}, Tuple{Float64, Float64}}}:
 0.707107  0.559017  0.5   0.559017
 0.559017  0.353553  0.25  0.353553
 0.5       0.25      0.0   0.25
 0.559017  0.353553  0.25  0.353553

julia> rr((4,4), scale=ScaFTEdge)
4×4 IndexFunArray{Float64, 2, IndexFunArrays.var"#9#10"{Float64, Tuple{Float64, Float64}, Tuple{Float64, Float64}}}:
 1.41421  1.11803   1.0  1.11803
 1.11803  0.707107  0.5  0.707107
 1.0      0.5       0.0  0.5
 1.11803  0.707107  0.5  0.707107

julia> rr2(Int, (3, 3), offset=(1, 1), scale=(10, 10))
3×3 IndexFunArray{Int64, 2, IndexFunArrays.var"#4#5"{Int64, Tuple{Int64, Int64}, Tuple{Int64, Int64}}}:
   0  100  400
 100  200  500
 400  500  800
```
## Application to selected dimensions
Note that the code below yields a 3D array but with a one-sized trailing dimension. This can then be used for broadcasting.
```jldoctest
julia> x = ones(5,6,5);

julia> y=rr2(selectsizes(x,(1,2)))
5×6×1 IndexFunArray{Float64, 3, IndexFunArrays.var"#4#5"{Float64, Tuple{Float64, Float64, Float64}, Tuple{Int64, Int64, Int64}}}:
[:, :, 1] =
 13.0  8.0  5.0  4.0  5.0  8.0
 10.0  5.0  2.0  1.0  2.0  5.0
  9.0  4.0  1.0  0.0  1.0  4.0
 10.0  5.0  2.0  1.0  2.0  5.0
 13.0  8.0  5.0  4.0  5.0  8.0
```
Similarly you can also use dimensions 2 and 3 yielding an array of `size(y) == (1,6,5)`.
Note that the necessary modification to the `Base.size` function is currently provided by this package.

## Using List-Mode Arguments
The code below generates 160 Voronoi cells at random positions. The `accumulator` was set to  mimimum
yielding in each pixel the square distance to the closest Voronoi center. See `gaussian` for another example
of using list-mode arguments.
```jldoctest
julia> y = rr2((1000,1000),offset = (1000.0,1000.0) .* rand(2,160), accumulator=minimum);

```

---
    rr2(arr::AbstractArray; offset=CtrFt, scaling=ScaUnit)

This is a wrapper for
`rr2(eltype(arr), size(arr), scaling=scaling, offset=offset)`.
"""
rr2

"""
    rr([T=Float64], size::NTuple{N, Int};
        offset=CtrFT,
        dims=ntuple(+, N),
        scale=ScaUnit,
        weight=1,
        accumulator=sum)

See `rr2` for a description of all options.

# Examples
```jldoctest
julia> rr((3, 3))
3×3 IndexFunArray{Float64, 2, IndexFunArrays.var"#9#10"{Float64, Tuple{Float64, Float64}, Tuple{Int64, Int64}}}:
 1.41421  1.0  1.41421
 1.0      0.0  1.0
 1.41421  1.0  1.41421

julia> rr((3, 3), offset=CtrCorner)
3×3 IndexFunArray{Float64, 2, IndexFunArrays.var"#9#10"{Float64, Tuple{Float64, Float64}, Tuple{Int64, Int64}}}:
 0.0  1.0      2.0
 1.0  1.41421  2.23607
 2.0  2.23607  2.82843
```
---
    rr(arr::AbstractArray; offset=CtrFt, scaling=ScaUnit)

This is a wrapper for
`rr(eltype(arr), size(arr), scaling=scaling, offset=offset)`.
"""
rr

"""
    xx([T=Float64], size::NTuple{N, Int};
        offset=CtrFT,
        dims=ntuple(+, N),
        scale=ScaUnit,
        weight=1,
        accumulator=sum)

A distance ramp along first dimension.
See `rr2` for a description of all options.
    ```jldoctest
julia> xx((4,4))
4×4 IndexFunArray{Float64, 2, IndexFunArrays.var"#14#15"{Float64, Tuple{Float64, Float64}, Tuple{Int64, Int64}}}:
 -2.0  -2.0  -2.0  -2.0
 -1.0  -1.0  -1.0  -1.0
  0.0   0.0   0.0   0.0
  1.0   1.0   1.0   1.0
```
---
    xx(arr::AbstractArray; offset=CtrFt, scaling=ScaUnit)

This is a wrapper for
`xx(eltype(arr), size(arr), scaling=scaling, offset=offset)`.
"""
xx

"""
    yy([T=Float64], size::NTuple{N, Int};
        offset=CtrFT,
        dims=ntuple(+, N),
        scale=ScaUnit,
        weight=1,
        accumulator=sum)

A distance ramp along second dimension.
See `rr2` for a description of all options.
    ```jldoctest
julia> yy((4,4))
4×4 IndexFunArray{Float64, 2, IndexFunArrays.var"#19#20"{Float64, Tuple{Float64, Float64}, Tuple{Int64, Int64}}}:
 -2.0  -1.0  0.0  1.0
 -2.0  -1.0  0.0  1.0
 -2.0  -1.0  0.0  1.0
 -2.0  -1.0  0.0  1.0
```
---
    yy(arr::AbstractArray; offset=CtrFt, scaling=ScaUnit)

This is a wrapper for
`yy(eltype(arr), size(arr), scaling=scaling, offset=offset)`.
"""
yy

"""
    zz([T=Float64], size::NTuple{N, Int};
        offset=CtrFT,
        dims=ntuple(+, N),
        scale=ScaUnit,
        weight=1,
        accumulator=sum)

A distance ramp along third dimension.
See `rr2` for a description of all options.
    ```jldoctest
julia> zz((1, 1, 4))
1×1×4 IndexFunArray{Float64, 3, IndexFunArrays.var"#24#25"{Float64, Tuple{Float64, Float64, Float64}, Tuple{Int64, Int64, Int64}}}:
[:, :, 1] =
 -2.0

[:, :, 2] =
 -1.0

[:, :, 3] =
 0.0

[:, :, 4] =
 1.0
```

---
    zz(arr::AbstractArray; offset=CtrFt, scaling=ScaUnit)

This is a wrapper for
`zz(eltype(arr), size(arr), scaling=scaling, offset=offset)`.
"""
zz

"""
    ee([T=Float64], size::NTuple{N, Int};
        offset=CtrFT,
        dims=ntuple(+, N),
        scale=ScaUnit,
        weight=1,
        accumulator=sum)

A distance ramp along forth (element) dimension. This dimension is often used as a color or wavelength channel.
See `rr2` for a description of all options.
    ```jldoctest
julia> ee((1, 1, 1, 4))
1×1×1×4 IndexFunArray{Float64, 4, IndexFunArrays.var"#60#63"{Float64, NTuple{4, Float64}}}:
[:, :, 1, 1] =
 -2.0

[:, :, 1, 2] =
 -1.0

[:, :, 1, 3] =
 0.0

[:, :, 1, 4] =
 1.0
```

---
    ee(arr::AbstractArray; offset=CtrFt, scaling=ScaUnit)

This is a wrapper for
`ee(eltype(arr), size(arr), scaling=scaling, offset=offset)`.
"""
ee

"""
    tt([T=Float64], size::NTuple{N, Int};
        offset=CtrFT,
        dims=ntuple(+, N),
        scale=ScaUnit,
        weight=1,
        accumulator=sum)

A distance ramp along fifth (time) dimension.
See `rr2` for a description of all options.
    ```jldoctest
julia> tt((1, 1, 1, 1, 4))
1×1×1×1×4 IndexFunArray{Float64, 5, IndexFunArrays.var"#69#72"{Float64, NTuple{5, Float64}}}:
[:, :, 1, 1, 1] =
 -2.0

[:, :, 1, 1, 2] =
 -1.0

[:, :, 1, 1, 3] =
 0.0

[:, :, 1, 1, 4] =
 1.0
```

---
    tt(arr::AbstractArray; offset=CtrFt, scaling=ScaUnit)

This is a wrapper for
`tt(eltype(arr), size(arr), scaling=scaling, offset=offset)`.
"""
tt

"""
idx_min([T=Float64], size::NTuple{N, Int};
        offset=CtrFT,
        dims=ntuple(+, N),
        scale=ScaUnit,
        weight=1,
        accumulator=sum)

returns the minimum of the absolute value of all index dimensions, after offset and scaling was applied.
In combination with a comparison operation a cross-shaped aperture is easily obtained. See below.
See `rr2` for a description of all options.
```jldoctest
julia> idx_min((5,5)) .< 1
5×5 BitMatrix:
 0  0  1  0  0
 0  0  1  0  0
 1  1  1  1  1
 0  0  1  0  0
 0  0  1  0  0
```
---
    idx_min(arr::AbstractArray; offset=CtrFt, scaling=ScaUnit)

This is a wrapper for
`idx_min(eltype(arr), size(arr), scaling=scaling, offset=offset)`.
"""
idx_min

"""
idx_max([T=Float64], size::NTuple{N, Int};
        offset=CtrFT,
        dims=ntuple(+, N),
        scale=ScaUnit,
        weight=1,
        accumulator=sum)

returns the maximum of the absolute value of all index dimensions, after offset and scaling was applied.
In combination with a comparison operation a box-shaped aperture is easily obtained. See below.
See `rr2` for a description of all options.
```jldoctest
julia> idx_max((5,5)) .< 2
5×5 BitMatrix:
 0  0  0  0  0
 0  1  1  1  0
 0  1  1  1  0
 0  1  1  1  0
 0  0  0  0  0
```
---
    idx_max(arr::AbstractArray; offset=CtrFt, scaling=ScaUnit)

This is a wrapper for
`idx_max(eltype(arr), size(arr), scaling=scaling, offset=offset)`.
"""
idx_max

"""
    phiphi([T=Float64], size::NTuple{N, Int};
        offset=CtrFT,
        dims=ntuple(+, N),
        scale=ScaUnit,
        weight=1,
        accumulator=sum)

An azimutal spiral phase ramp using atan(). The azimuthal phase spans dimensions 1 and 2.
See `rr2` for a description of all options.
    ```jldoctest
julia> phiphi((5, 5))
5×5 IndexFunArray{Float64, 2, IndexFunArrays.var"#29#30"{Float64, Tuple{Float64, Float64}, Tuple{Int64, Int64}}}:
 -2.35619   -2.67795   3.14159  2.67795   2.35619
 -2.03444   -2.35619   3.14159  2.35619   2.03444
 -1.5708    -1.5708    0.0      1.5708    1.5708
 -1.10715   -0.785398  0.0      0.785398  1.10715
 -0.785398  -0.463648  0.0      0.463648  0.785398
```

---
    phiphi(arr::AbstractArray; offset=CtrFt, scaling=ScaUnit)

This is a wrapper for
`phiphi(eltype(arr), size(arr), scaling=scaling, offset=offset)`.
"""
phiphi

"""
    ramp(::Type{T}, dim::Int, dim_size::Int;
        offset=CtrFT, scale=ScaUnit,
        weight=1,
        accumulator=sum) where {T}

Generates a dim-dimensional ramp of size `dim_size` to be used for broadcasting through multidimensional expressions.
`dim` highest dimension of the oriented array to be generated. This is also the ramp direction.
`dim_size` size along this dimension.

For details about offset and scale and dims see rr2.
```jldoctest
julia> ramp(Float32, 1, 7; offset=(2,))
7-element IndexFunArray{Float32, 1, IndexFunArrays.var"#434#435"{Tuple{Int64}, Tuple{Int64}, Int64}}:
 -1
  0
  1
  2
  3
  4
  5
```
ramp(dim::Int, dim_size::Int; offset=CtrFt, scaling=ScaUnit)

This is a wrapper for
`ramp(Float64, dim::Int, dim_size::Int; offset=CtrFt, scaling=ScaUnit)`.

"""
ramp

"""
    delta([T=Float64], size::NTuple{N, Int};
        offset=CtrFT,
        dims=ntuple(+, N),
        scale=ScaUnit,
        weight=1,
        accumulator=sum)

A `delta` peak positioned at offset. See `rr2()` for a description of all options.
Note that `scale` does not influence the result. Also note that this function operates on
a comparison for equality, which means a sub-pixel offset of the delta results into zero.
```jldoctest
julia> delta((5,5))
5×5 IndexFunArray{Float64, 2, IndexFunArrays.var"#79#81"{Float64, Tuple{Float64, Float64}}}:
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0

julia> delta((6,6))
6×6 IndexFunArray{Float64, 2, IndexFunArrays.var"#79#81"{Float64, Tuple{Float64, Float64}}}:
 0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  1.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0

 julia> delta((5,5),offset=CtrCorner)
 5×5 IndexFunArray{Float64, 2, IndexFunArrays.var"#79#81"{Float64, Tuple{Float64, Float64}}}:
  1.0  0.0  0.0  0.0  0.0
  0.0  0.0  0.0  0.0  0.0
  0.0  0.0  0.0  0.0  0.0
  0.0  0.0  0.0  0.0  0.0
  0.0  0.0  0.0  0.0  0.0

```

---
    delta(arr::AbstractArray; offset=CtrFt, scaling=ScaUnit)

This is a wrapper for
`delta(eltype(arr), size(arr), scaling=scaling, offset=offset)`.
"""
delta

"""
    gaussian([T=Float64], size::NTuple{N, Int};
        offset=CtrFT,
        sigma=1.0,
        dims=ntuple(+, N),
        scale=ScaUnit,
        weight=1,
        accumulator=sum)

A gaussian peak positioned at offset. Note that the gaussian is NOT normalized by its integral, but by its maximum.
For a version with a normalized integral, see `normal`.
List-mode is supported also for the argument sigma.
See the final example below which generates 60 Gaussians at random positions with random strengths and width along X and Y.

# Arguments:
* `sigma`: the (standard deviation) width of the Gaussian. If a tuple is supplied, each entry is interpreted as the width along the correspondin dimension.
* `offset`: the center position of the Gaussian. You can use a tuple or the indicators `CtrCorner`, `CtrEnd`, `CtrFT`, `CtrRFT` etc.
* `scale`: the scale of the pixel. By default `ScaUnit` is assumed
* `dims`: the dimensions over which to apply this function to.
* `weight`: the strength of the result. Supports list-mode (see rr2 for documentation)
* `accumulator`: the method used for superimposing list-mode data. Only applies in list-mode
```jldoctest
julia> gaussian((5,5))
5×5 IndexFunArray{Float64, 2, IndexFunArrays.var"#100#102"{Float64, Tuple{Float64, Float64}, Tuple{Float64, Float64}}}:
 0.0183156  0.082085  0.135335  0.082085  0.0183156
 0.082085   0.367879  0.606531  0.367879  0.082085
 0.135335   0.606531  1.0       0.606531  0.135335
 0.082085   0.367879  0.606531  0.367879  0.082085
 0.0183156  0.082085  0.135335  0.082085  0.0183156

julia> gaussian((6,6), sigma=5.0)
 6×6 IndexFunArray{Float64, 2, IndexFunArrays.var"#100#102"{Float64, Tuple{Float64, Float64}, Tuple{Float64, Float64}}}:
  0.165299  0.272532  0.367879  0.40657   0.367879  0.272532
  0.272532  0.449329  0.606531  0.67032   0.606531  0.449329
  0.367879  0.606531  0.818731  0.904837  0.818731  0.606531
  0.40657   0.67032   0.904837  1.0       0.904837  0.67032
  0.367879  0.606531  0.818731  0.904837  0.818731  0.606531
  0.272532  0.449329  0.606531  0.67032   0.606531  0.449329

julia> gaussian(Float32,(5,5),offset=CtrCorner)
  5×5 IndexFunArray{Float32, 2, IndexFunArrays.var"#100#102"{Float32, Tuple{Float64, Float64}, Tuple{Float32, Float32}}}:
   1.0          0.606531     0.135335    0.011109    0.000335463
   0.606531     0.367879     0.082085    0.00673795  0.000203468
   0.135335     0.082085     0.0183156   0.00150344  4.53999f-5
   0.011109     0.00673795   0.00150344  0.00012341  3.72665f-6
   0.000335463  0.000203468  4.53999f-5  3.72665f-6  1.12535f-7

julia> y = gaussian((100,100),offset = (100,100) .* rand(2,60), weight=rand(60), sigma=2.0 .*(0.3 .+rand(2,60)));

```

---
gaussian(arr::AbstractArray; offset=CtrFt, sigma=1.0, scaling=ScaUnit)

This is a wrapper for
`gaussian(eltype(arr), size(arr), sigma=sigma, scaling=scaling, offset=offset)`.
"""
gaussian

"""
    normal([T=Float64], size::NTuple{N, Int};
        offset=CtrFT,
        sigma=1.0,
        scale=ScaUnit,
        weight=1,
        accumulator=sum)

A gaussian peak positioned at offset. Note that this normal distribution (Gaussian) is normalized by its integral.
For a version with normalized to the maximum, see `gaussian`.

# Arguments:
* `sigma`: the (standard deviation) width of the Gaussian
* `offset`: the center position of the Gaussian. You can use a tuple or the indicators `CtrCorner`, `CtrEnd`, `CtrFT`, `CtrRFT` etc.
* `scale`: the scale of the pixel. By default `ScaUnit` is assumed
* `weight`: the strength of the result. Supports list-mode (see rr2 for documentation)
* `accumulator`: the method used for superimposing list-mode data. Only applies in list-mode
```jldoctest
julia> normal((5,5))
5×5 IndexFunArray{Float64, 2, IndexFunArrays.var"#107#109"{Float64, Tuple{Float64, Float64}, Tuple{Float64, Float64}}}:
 0.00291502  0.0130642  0.0215393  0.0130642  0.00291502
 0.0130642   0.0585498  0.0965324  0.0585498  0.0130642
 0.0215393   0.0965324  0.159155   0.0965324  0.0215393
 0.0130642   0.0585498  0.0965324  0.0585498  0.0130642
 0.00291502  0.0130642  0.0215393  0.0130642  0.00291502

julia> sum(normal((5,5)))
 0.9818147610543744

julia> sum(normal((25,25)))
 1.000000010701151

julia> sum(normal((55,55), sigma=(5.0,2.0)))
0.9999999639340563
```

---
normal(arr::AbstractArray; offset=CtrFt, sigma=1.0, scaling=ScaUnit)

This is a wrapper for
`normal(eltype(arr), size(arr), sigma=sigma, scaling=scaling, offset=offset)`.
"""
normal

"""
    exp_ikx([T=Float64], size::NTuple{N, Int};
        offset=CtrFT,
        shift_by=size.÷2
        dims=ntuple(+, N),
        scale=ScaFT,
        weight=1,
        accumulator=sum)

A complex-valued phase ramp according to `exp(-2pi i <k,x>)`. If applied as a multiplicative factor in Fourier space,
it will lead to a shift of `x` pixels in real space. Note that this effect is actually realized by a change to the scaling parameter.
The default shift is `size.÷2` which corresponds to `CtrFT`, however, the Ctr... arguments cannot be used for `shift_by`.

The argument `shift_by` supports list-mode, which can be used to conveniently perform multiple shifts simulatneously.
See the final example below, which generates delta peaks at randomized subpixel positions.

# Arguments:
* `offset`: the center position of the Gaussian. You can use a tuple or the indicators `CtrCorner`, `CtrEnd`, `CtrFT`, `CtrRFT` etc.
* `shift_by`: the amount to shift by in real space.
* `scale`: the scale of the pixel. By default `ScaUnit` is assumed
* `dims`: the dimensions over which to apply this function to.
* `weight`: the strength of the result. Supports list-mode (see rr2 for documentation)
* `accumulator`: the method used for superimposing list-mode data. Only applies in list-mode
```jldoctest
julia> a = rr((4,3),offset=CtrCorner)
4×3 IndexFunArray{Float64, 2, IndexFunArrays.var"#37#39"{Float64, Tuple{Float64, Float64}, Tuple{Int64, Int64}}}:
 0.0  1.0      2.0
 1.0  1.41421  2.23607
 2.0  2.23607  2.82843
 3.0  3.16228  3.60555

julia> using FourierTools; real(iffts(ffts(a).*exp_ikx(a, shift_by=(2.0,1.0))))
4×3 Matrix{Float64}:
 2.82843   2.0          2.23607
 3.60555   3.0          3.16228
 2.0      -2.22045e-16  1.0
 2.23607   1.0          1.41421

 julia> using FourierTools; y = real(ift(exp_ikx((101,101),weight=rand(60), shift_by=101.0 .*rand(2,60))));

```

---
    exp_ikx(arr::AbstractArray; offset=CtrFt, shift_by==size(arr).÷2, scaling=ScaUnit)

This is a wrapper for
`exp_ikx(eltype(arr), size(arr), shift_by=shift_by, scaling=scaling, offset=offset)`.
"""
exp_ikx

"""
    phase_kz([T=Float64], size::NTuple{N, Int};
    offset=CtrFT,
    shift_by=size.÷2
    dims=ntuple(+, N),
    scale=ScaFT,
    weight=1,
    accumulator=sum)
    Calculates a propagation phase (without the 2pi factor!) for a given z-position, which can be defined via a 3rd entry in the `offset` supplied to the function.
    By default, Nyquist sampling it is assumed such that the lateral k_xy corresponds to the XY border in frequency space at the edge of the Ewald circle.
    However, via the xy `scale` entries the k_max can be set appropriately. The propagation equation uses
    offset[3] .* sqrt.(1-kxy_rel^2) as the propagation phase. Therefore the Z-propagation distance (offset[3]) has to be specified in units of the wavelength in the medium (`λ = n*λ₀`).
    Note that since the phase is normalized to 1 instead of 2pi, you need to use this phase in the following sense: `cispi.(2.*phase_kz(...))`.
    For other arguments look at similar scaler IndexFunArray functions such as rr.

See also: `propagator()`
"""
phase_kz

"""
    propagator(size::NTuple{N, Int}; Δz=1.0, shift_by=(0,0), k_max=0.5, offset=CtrFT, scale=ScaFT, dims=ntuple(+, N), accumulator=sum, weight=1) where {N,T}

    A complex-valued phase which describes the `k_z= Δz * sqrt(1-k_x^2-k_y^2)` phase change upon free space propagation.

    The argument `shift_by` and `Δz` arguments support list-mode, which can be used to conveniently perform multiple shifts simulatneously.

    # Arguments:
    * `offset`: the center position of the Gaussian. You can use a tuple or the indicators `CtrCorner`, `CtrEnd`, `CtrFT`, `CtrRFT` etc.
    * `Δz`: the amount to propagate along the third dimension z by in real space in units of the wavelength in the medium (`λ = n*λ₀`).
    * `shift_by`: the amount to shift by in x and y spatial direct in real space.
    * `k_max`: indicates the sampling relative to the Nyquist frequency. In optics this should be `pixelpitch./λ` with the `pixelpitch` as a tuple and the wavelength `λ = n*λ₀` in the medium.
    * `scale`: the scale of the pixel. By default `ScaUnit` is assumed
    * `dims`: the dimensions over which to apply this function to.
    * `weight`: the strength of the result. Supports list-mode (see rr2 for documentation)
    * `accumulator`: the method used for superimposing list-mode data. Only applies in list-mode
---
propagator(arr::AbstractArray; Δz=1.0, shift_by=(0,0), k_max=0.5, offset=CtrFt, scaling=ScaUnit)

This is a wrapper for
`propagator(eltype(arr), size(arr), Δz=Δz, shift_by=shift_by, k_max=k_max, scaling=scaling, offset=offset)`.

"""
propagator

"""
    window_linear([T=Float64], size::NTuple;
                offset=CtrFT, scale=ScaFTEdge, border_in=0.8, border_out=1.0, dims=ntuple(+, N))

A multidimensional (separable) window with a linear transition from zero at the borders (`border_out`) to one (`border_in`).
```jldoctest
julia> window_linear((8,9),border_in=0.0)
8×9 IndexFunArray{Float64, 2, IndexFunArrays.var"#34#35"{Float64, Tuple{Float64, Float64}, Tuple{Float64, Float64}, Float64, Float64}}:
 0.0  0.0     0.0    0.0     0.0   0.0     0.0    0.0     0.0
 0.0  0.0625  0.125  0.1875  0.25  0.1875  0.125  0.0625  0.0
 0.0  0.125   0.25   0.375   0.5   0.375   0.25   0.125   0.0
 0.0  0.1875  0.375  0.5625  0.75  0.5625  0.375  0.1875  0.0
 0.0  0.25    0.5    0.75    1.0   0.75    0.5    0.25    0.0
 0.0  0.1875  0.375  0.5625  0.75  0.5625  0.375  0.1875  0.0
 0.0  0.125   0.25   0.375   0.5   0.375   0.25   0.125   0.0
 0.0  0.0625  0.125  0.1875  0.25  0.1875  0.125  0.0625  0.0
```

---
    window_linear(arr::AbstractArray; offset=CtrFt, scaling=ScaUnit,
                                      border_in=0.8, border_out=1.0)

This is a wrapper for
`window_linear(eltype(arr), size(arr), scaling=scaling, offset=offset, border_in=border_in, border_out=border_out)`.
"""
window_linear

"""
    window_radial_linear([T=Float64], size::NTuple;
                offset=CtrFT, scale=ScaFTEdge, border_in=0.8, border_out=1.0, dims=ntuple(+, N))

A multidimensional radial window with a linear transition from zero at the borders (`border_out`) to one (`border_in`).
Note that `border_in` and `border_out` need to be scalar values for `radial` type windows. The elypticity can be adjusted via the `scale` parameter.
With the default offset and scale the borders are specified relative to the edge.
```jldoctest
julia> window_radial_linear((4,5),border_in=0.0)
4×5 IndexFunArray{Float64, 2, IndexFunArrays.var"#59#60"{Float64, Tuple{Float64, Float64}, Tuple{Float64, Float64}, Float64, Float64}}:
 0.0  0.0       0.0  0.0       0.0
 0.0  0.292893  0.5  0.292893  0.0
 0.0  0.5       1.0  0.5       0.0
 0.0  0.292893  0.5  0.292893  0.0
```

---
    window_radial_linear(arr::AbstractArray; offset=CtrFt, scaling=ScaUnit,
                         border_in=0.8, border_out=1.0, dims=ntuple(+, N))

This is a wrapper for
`window_radial_linear(eltype(arr), size(arr), scaling=scaling, offset=offset, border_in=border_in, border_out=border_out)`.
"""
window_radial_linear

"""
    window_edge([T=Float64], size::NTuple;
                offset=CtrFT, scale=ScaFTEdge, border_in=0.8, border_out=1.0, dims=ntuple(+, N))

A multidimensional (separable) window with a sudden transition half way between the borders (`border_out`) to one (`border_in`).
See `?window_linear` for more details on the arguments.

---
    window_edge(arr::AbstractArray; offset=CtrFt, scaling=ScaUnit,
                                      border_in=0.8, border_out=1.0)

This is a wrapper for
`window_edge(eltype(arr), size(arr), scaling=scaling, offset=offset, border_in=border_in, border_out=border_out)`.
"""
window_edge

"""
    window_radial_edge([T=Float64], size::NTuple;
                offset=CtrFT, scale=ScaFTEdge, border_in=0.8, border_out=1.0, dims=ntuple(+, N))

A multidimensional radial window (disk) with a sudden transition half way between the borders (`border_out`) to one (`border_in`).
Note that `border_in` and `border_out` need to be scalar values for `radial` type windows. The elypticity can be adjusted via the `scale` parameter.
See `?window_radial_linear` for more details on the arguments.

---
    window_radial_edge(arr::AbstractArray;
                       offset=CtrFT, scale=ScaFTEdge, border_in=0.8, border_out=1.0, dims=ntuple(+, N))

This is a wrapper for
`window_radial_edge(eltype(arr), size(arr), scaling=scaling, offset=offset, border_in=border_in, border_out=border_out)`.
"""
window_radial_edge

"""
    window_hanning([T=Float64], size::NTuple;
                       offset=CtrFT, scale=ScaFTEdge, border_in=0.8, border_out=1.0, dims=ntuple(+, N))

A multidimensional (separable) window with a von Hann transition between the borders (`border_out`) to one (`border_in`).
See `?window_linear` for more details on the arguments.

---
    window_hanning(arr::AbstractArray;
                       offset=CtrFT, scale=ScaFTEdge, border_in=0.8, border_out=1.0, dims=ntuple(+, N))

This is a wrapper for
`window_hanning(eltype(arr), size(arr), scaling=scaling, offset=offset, border_in=border_in, border_out=border_out)`.
"""
window_hanning

"""
    window_radial_hanning([T=Float64], size::NTuple;
                       offset=CtrFT, scale=ScaFTEdge, border_in=0.8, border_out=1.0, dims=ntuple(+, N))

A multidimensional radial window with a von Hann transition between the borders (`border_out`) to one (`border_in`).
Note that `border_in` and `border_out` need to be scalar values for `radial` type windows. The elypticity can be adjusted via the `scale` parameter.
See `?window_radial_linear` for more details on the arguments.

---
    window_radial_hanning(arr::AbstractArray;
                          offset=CtrFT, scale=ScaFTEdge, border_in=0.8, border_out=1.0, dims=ntuple(+, N))

This is a wrapper for
`window_radial_hanning(eltype(arr), size(arr), scaling=scaling, offset=offset, border_in=border_in, border_out=border_out)`.
"""
window_radial_hanning

"""
    window_hamming([T=Float64], size::NTuple;
                   offset=CtrFT, scale=ScaFTEdge, border_in=0.8, border_out=1.0, dims=ntuple(+, N))

A multidimensional (separable) window with a Hamming transition between the borders (`border_out`) to one (`border_in`).
See `?window_linear` for more details on the arguments.

---
    window_hamming(arr::AbstractArray;
                   offset=CtrFT, scale=ScaFTEdge, border_in=0.8, border_out=1.0, dims=ntuple(+, N))

This is a wrapper for
`window_hamming(eltype(arr), size(arr), scaling=scaling, offset=offset, border_in=border_in, border_out=border_out)`.
"""
window_hamming

"""
    window_radial_hamming([T=Float64], size::NTuple;
                          offset=CtrFT, scale=ScaFTEdge, border_in=0.8, border_out=1.0, dims=ntuple(+, N))

A multidimensional radial window with a Hamming transition between the borders (`border_out`) to one (`border_in`).
Note that `border_in` and `border_out` need to be scalar values for `radial` type windows. The elypticity can be adjusted via the `scale` parameter.
See `?window_radial_linear` for more details on the arguments.

---
    window_radial_hamming(arr::AbstractArray;
                          offset=CtrFT, scale=ScaFTEdge, border_in=0.8, border_out=1.0, dims=ntuple(+, N))

This is a wrapper for
`window_radial_hamming(eltype(arr), size(arr), scaling=scaling, offset=offset, border_in=border_in, border_out=border_out)`.
"""
window_radial_hamming

"""
    window_blackman_harris([T=Float64], size::NTuple;
                          offset=CtrFT, scale=ScaFTEdge, border_in=0.8, border_out=1.0, dims=ntuple(+, N))

A multidimensional (separable) window with a  transition according to Blackman/Harris between the borders (`border_out`) to one (`border_in`).
See `?window_linear` for more details on the arguments.

---
    window_blackman_harris(arr::AbstractArray;
                          offset=CtrFT, scale=ScaFTEdge, border_in=0.8, border_out=1.0, dims=ntuple(+, N))

This is a wrapper for
`window_blackman_harris(eltype(arr), size(arr), scaling=scaling, offset=offset, border_in=border_in, border_out=border_out)`.
"""
window_blackman_harris

"""
    window_radial_blackman_harris([T=Float64], size::NTuple;
                          offset=CtrFT, scale=ScaFTEdge, border_in=0.8, border_out=1.0, dims=ntuple(+, N))

A multidimensional radial window with a Hamming transition according to Blackman/Harris between the borders (`border_out`) to one (`border_in`).
Note that `border_in` and `border_out` need to be scalar values for `radial` type windows. The elypticity can be adjusted via the `scale` parameter.
See `?window_radial_linear` for more details on the arguments.

---
    window_radial_blackman_harris(arr::AbstractArray;
                                  offset=CtrFT, scale=ScaFTEdge, border_in=0.8, border_out=1.0, dims=ntuple(+, N))
This is a wrapper for
`window_radial_blackman_harris(eltype(arr), size(arr), scaling=scaling, offset=offset, border_in=border_in, border_out=border_out)`.
"""
window_radial_blackman_harris

"""
    window_gaussian([T=Float64], size::NTuple;
                          offset=CtrFT, scale=ScaFTEdge, border_in=0.8, border_out=1.0, dims=ntuple(+, N))

A multidimensional (separable) window with a  transition according to a decaying Gaussian between the borders (`border_out`) to one (`border_in`).
By default the standard deviation `sigma` of the Gaussian is adjusted such that the `2 sigma` level is reached at `border_out`.
However, this window is not clipped at the outer border, thus allowing the sigma to be adjusted by placing `border_out` closer to `border_in`.
See `?window_linear` for more details on the arguments.

    ```jldoctest
julia> w1 = window_gaussian((9,9), border_in=(0.3,0.3), border_out=(0.6,1))
9×9 IndexFunArray{Float64, 2, IndexFunArrays.var"#286#288"{Float64, Tuple{Float64, Float64}, Tuple{Float64, Float64}, Tuple{Float64, Float64}, Tuple{Float64, Int64}}}:
 0.000109245  0.000259904  0.000413188  0.000449917  0.000449917  0.000449917  0.000413188  0.000259904  0.000109245
 0.012239     0.0291178    0.0462907    0.0504055    0.0504055    0.0504055    0.0462907    0.0291178    0.012239
 0.152725     0.363345     0.577637     0.628984     0.628984     0.628984     0.577637     0.363345     0.152725
 0.242811     0.57767      0.918365     1.0          1.0          1.0          0.918365     0.57767      0.242811
 0.242811     0.57767      0.918365     1.0          1.0          1.0          0.918365     0.57767      0.242811
 0.242811     0.57767      0.918365     1.0          1.0          1.0          0.918365     0.57767      0.242811
 0.152725     0.363345     0.577637     0.628984     0.628984     0.628984     0.577637     0.363345     0.152725
 0.012239     0.0291178    0.0462907    0.0504055    0.0504055    0.0504055    0.0462907    0.0291178    0.012239
 0.000109245  0.000259904  0.000413188  0.000449917  0.000449917  0.000449917  0.000413188  0.000259904  0.000109245
```
---

    window_gaussian(arr::AbstractArray;
                          offset=CtrFT, scale=ScaFTEdge, border_in=0.8, border_out=1.0, dims=ntuple(+, N))

This is a wrapper for
`window_gaussian(eltype(arr), size(arr), scaling=scaling, offset=offset, border_in=border_in, border_out=border_out)`.
"""
window_gaussian

"""
    window_radial_gaussian([T=Float64], size::NTuple;
                          offset=CtrFT, scale=ScaFTEdge, border_in=0.8, border_out=1.0, dims=ntuple(+, N))

A multidimensional radial window with a Gaussian transition between the borders (`border_out`) to one (`border_in`).
By default the standard deviation `sigma` of the Gaussian is adjusted such that the `2 sigma` level is reached at `border_out`.
However, this window is not clipped at the outer border, thus allowing the sigma to be adjusted by placing `border_out` closer to `border_in`.
Note that `border_in` and `border_out` are scalar values for `radial` type windows. The elypticity can be adjusted via the `scale` parameter.
See `?window_radial_linear` for more details on the arguments.

---
window_radial_gaussian(arr::AbstractArray;
                                  offset=CtrFT, scale=ScaFTEdge, border_in=0.8, border_out=1.0, dims=ntuple(+, N))
This is a wrapper for
`window_radial_gaussian(eltype(arr), size(arr), scaling=scaling, offset=offset, border_in=border_in, border_out=border_out)`.
"""
window_radial_gaussian

"""
    box([::Type{T}], size::NTuple, boxsize; offset=CtrFT, scale=ScaFTEdge, dims=ntuple(+, N))

A multidimensional box, being one inside and zero outside. `boxsize` defines the outer dimensions of the box in the sense that a comparison is documentation
with the distance to the center being smaller than half the `boxsize` value. Note that the default offset leads to even sizes being rounded down (see example below).
The default result datatype is Float64.
```jldoctest
julia> box((10,10),(3,6))
10×10 Matrix{Float64}:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  1.0  1.0  1.0  1.0  1.0  0.0  0.0
 0.0  0.0  0.0  1.0  1.0  1.0  1.0  1.0  0.0  0.0
 0.0  0.0  0.0  1.0  1.0  1.0  1.0  1.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
```
---
box(arr::AbstractArray, boxsize; offset=CtrFT, scale=ScaUnit, dims=ntuple(+, N))

This is a wrapper for
`box(eltype(arr), size(arr), boxsize; scaling=scaling, offset=offset)`.
"""
box

"""
    disc([::Type{T}], size::NTuple, disc_radius; offset=CtrFT, scale=ScaFTEdge, dims=ntuple(+, N))

A multidimensional disc (i.e. a disk in 2D and a filled sphere in 3D), being one inside and zero outside. `disc_radius` defines the radius of the disc.
The default result datatype is Float64. An alternative to generating a disc would be the edge-window.
```jldoctest
julia> disc((10,10),(3,5))
10×10 Matrix{Float64}:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  0.0
 0.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0
 0.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0
 0.0  0.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  0.0
 0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0

 julia> disc((15,8),(0.5,0.5),scale=ScaFT)
15×8 Matrix{Float64}:
 0.0  0.0  0.0  1.0  1.0  1.0  0.0  0.0
 0.0  0.0  1.0  1.0  1.0  1.0  1.0  0.0
 0.0  0.0  1.0  1.0  1.0  1.0  1.0  0.0
 0.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0
 0.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0
 0.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0
 0.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0
 0.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0
 0.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0
 0.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0
 0.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0
 0.0  0.0  1.0  1.0  1.0  1.0  1.0  0.0
 0.0  0.0  1.0  1.0  1.0  1.0  1.0  0.0
 0.0  0.0  0.0  1.0  1.0  1.0  0.0  0.0
```
---
disc(arr::AbstractArray, disc_radius; offset=CtrFT, scale=ScaUnit, dims=ntuple(+, N))

This is a wrapper for
`disc(eltype(arr), size(arr), disc_radius; scaling=scaling, offset=offset)`.
"""
disc

"""
    axes1d(size::NTuple{N, Int}; offset=CtrFT, scale=ScaUnit, dims=ntuple(+, N), accumulator=sum, weight=1) where {N}
    returns a tuple of axes, all oriented along x.
    This is very useful for plotting (see example below).
    The optional argument `keepdims` can be used to obtain oriented single-dimensional vectors
    which can be used in analogy to `meshgrid` operations but supporting broadcasting and generator performance.
    See `rr2` for a description of other options.

```jldoctest
julia> using Plots, IndexFunArrays

julia> sz = (200,100);

julia> heatmap(axes1d(sz)..., rr(sz)')

julia> x,y = axes1d((10,11), keepdims=true)
"hi" = "hi"
Base.Generator{Tuple{Int64, Int64}, IndexFunArrays.var"#675#677"{Float64, Tuple{Int64, Int64}, Tuple{Int64, Int64}, Tuple{Int64, Int64}}}(IndexFunArrays.var"#675#677"{Float64, Tuple{Int64, Int64}, Tuple{Int64, Int64}, Tuple{Int64, Int64}}((6, 6), (1, 1), (10, 11)), (1, 2))

julia> size(y)
(1, 11)
```
---
axes1d(arr::AbstractArray{T, N}; offset=CtrFT, scale=ScaUnit, dims=ntuple(+, N), accumulator=sum, weight=1) where {N}

This is a wrapper for
`axes1d(size(arr); offset=CtrFT, scale=ScaUnit, dims=ntuple(+, N), accumulator=sum, weight=1) where {N}`
"""
axes1d


text:
using ReachabilityAnalysis, SparseArrays, Symbolics

const var = @variables x[1:9] t

function platoon_connected(; deterministic_switching::Bool=true,
                             c1=5.0)  # clock constraints
    n = 9 + 1

    # x' = Ax + Bu + c
    A = Matrix{Float64}(undef, n, n)
    A[1, :] = [0, 1.0, 0, 0, 0, 0, 0, 0, 0, 0]
    A[2, :] = [0, 0, -1.0, 0, 0, 0, 0, 0, 0, 0]
    A[3, :] = [1.6050, 4.8680, -3.5754, -0.8198, 0.4270, -0.0450, -0.1942,  0.3626, -0.0946, 0.]
    A[4, :] = [0, 0, 0, 0, 1.0, 0, 0, 0, 0, 0,]
    A[5, :] = [0, 0, 1.0, 0, 0, -1.0, 0, 0, 0, 0]
    A[6, :] = [0.8718, 3.8140, -0.0754,  1.1936, 3.6258, -3.2396, -0.5950,  0.1294, -0.0796, 0.]
    A[7, :] = [0, 0, 0, 0, 0, 0, 0, 1.0, 0, 0]
    A[8, :] = [0, 0, 0, 0, 0, 1.0, 0, 0, -1.0, 0]
    A[9, :] = [0.7132, 3.5730, -0.0964,  0.8472, 3.2568, -0.0876,  1.2726,  3.0720, -3.1356, 0.]
    A[10, :] = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0]; # t' = 1

    if deterministic_switching
        invariant = HalfSpace(t <= c1, var)
    else
        invariant = Universe(n)
    end

    # acceleration of the lead vehicle + time
    B = sparse([2], [1], [1.0], n, 1)
    U = Hyperrectangle(low=[-9.], high=[1.])
    c = [0, 0, 0, 0, 0, 0, 0, 0, 0, 1.0]
    @system(x' = Ax + Bu + c, x ∈ invariant, u ∈ U)
end

function platoon_disconnected(; deterministic_switching::Bool=true,
                                c2=5.0)  # clock constraints
    n = 10 # 9 dimensions + time

    # x' = Ax + Bu + c
    A = Matrix{Float64}(undef, n, n)
    A[1, :] = [0, 1.0, 0, 0, 0, 0, 0, 0, 0, 0]
    A[2, :] = [0, 0, -1.0, 0, 0, 0, 0, 0, 0, 0]
    A[3, :] = [1.6050, 4.8680, -3.5754, 0, 0, 0, 0, 0, 0, 0]
    A[4, :] = [0, 0, 0, 0, 1.0, 0, 0, 0, 0, 0,]
    A[5, :] = [0, 0, 1.0, 0, 0, -1.0, 0, 0, 0, 0]
    A[6, :] = [0, 0, 0,  1.1936, 3.6258, -3.2396, 0, 0, 0, 0.]
    A[7, :] = [0, 0, 0, 0, 0, 0, 0, 1.0, 0, 0]
    A[8, :] = [0, 0, 0, 0, 0, 1.0, 0, 0, -1.0, 0]
    A[9, :] = [0.7132, 3.5730, -0.0964,  0.8472, 3.2568, -0.0876,  1.2726,  3.0720, -3.1356, 0.]
    A[10, :] = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0]; # t' = 1

    if deterministic_switching
        invariant = HalfSpace(t <= c2, var)
    else
        invariant = Universe(n)
    end

    # acceleration of the lead vehicle + time
    B = sparse([2], [1], [1.0], n, 1)
    U = Hyperrectangle(low=[-9.], high=[1.])
    c = [0, 0, 0, 0, 0, 0, 0, 0, 0, 1.0]
    @system(x' = Ax + Bu + c, x ∈ invariant, u ∈ U)
end

function platoon(; deterministic_switching::Bool=true,
                   c1=5.0,  # clock constraints
                   c2=5.0,  # clock constraints
                   tb=10.0,  # lower bound for loss of communication
                   tc=20.0, tr=20.0) # upper bound for loss of communication (tc) and reset time (tr)

    # three variables for each vehicle, (ei, d(et)/dt, ai) for
    # (spacing error, relative velocity, speed), and the last dimension is time
    n = 9 + 1

    # transition graph
    automaton = LightAutomaton(2)
    add_transition!(automaton, 1, 2, 1)
    add_transition!(automaton, 2, 1, 2)

    # modes
    mode1 = platoon_connected(deterministic_switching=deterministic_switching, c1=c1)
    mode2 = platoon_disconnected(deterministic_switching=deterministic_switching, c2=c2)
    modes = [mode1, mode2]

    # common reset
    reset = Dict(n => 0.)

    # transition l1 -> l2
    if deterministic_switching
        guard = Hyperplane(t == c1, var)
    else
        guard = HPolyhedron([tb <= t, t <= tc], var)
    end
    t1 = ConstrainedResetMap(n, guard, reset)

    # transition l2 -> l1
    if deterministic_switching
        guard = Hyperplane(t == c2, var)
    else
        guard = HalfSpace(t <= tr, var)
    end
    t2 = ConstrainedResetMap(n, guard, reset)
    resetmaps = [t1, t2]

    H = HybridSystem(automaton, modes, resetmaps, [AutonomousSwitching()])

    # initial condition is at the orgin in mode 1
    X0 = BallInf(zeros(n), 0.0)
    initial_condition = [(1, X0)]

    return IVP(H, initial_condition)
end

function dmin_specification(sol, dmin)
    return (-ρ(sparsevec([1], [-1.0], 10), sol) > -dmin) &&
           (-ρ(sparsevec([4], [-1.0], 10), sol) > -dmin) &&
           (-ρ(sparsevec([7], [-1.0], 10), sol) > -dmin)
end

prob_PLAD01 = platoon();

const boxdirs = BoxDirections(10)
const octdirs = OctDirections(10);

length(boxdirs)

alg = BOX(δ=0.01)
sol_PLAD01_BND42 = solve(prob_PLAD01,
                         alg=alg,
                         clustering_method=BoxClustering(1),
                         intersection_method=TemplateHullIntersection(boxdirs),
                         intersect_source_invariant=false,
                         tspan = (0.0 .. 20.0));

dmin_specification(sol_PLAD01_BND42, 42)

-ρ(sparsevec([1], [-1.0], 10), sol_PLAD01_BND42)

-ρ(sparsevec([4], [-1.0], 10), sol_PLAD01_BND42)

-ρ(sparsevec([7], [-1.0], 10), sol_PLAD01_BND42)

using Plots, LaTeXStrings

plot(sol_PLAD01_BND42, vars=(0, 1), xlab=L"t", ylab=L"x_1", title="PLAD01 - BND42", lw=0.1)
plot!(x->x, x->-42., 0., 20., linewidth=2, color="red", linestyle=:dash, leg=nothing)

length(octdirs)

alg = LGG09(δ=0.03, template=octdirs, approx_model=Forward(setops=octdirs));
sol_PLAD01_BND30 = solve(prob_PLAD01,
                         alg=alg,
                         clustering_method=LazyClustering(1),
                         intersection_method=TemplateHullIntersection(octdirs),
                         intersect_source_invariant=false,
                         tspan = (0.0 .. 20.0));

dmin_specification(sol_PLAD01_BND30, 30)

-ρ(sparsevec([1], [-1.0], 10), sol_PLAD01_BND30)

-ρ(sparsevec([4], [-1.0], 10), sol_PLAD01_BND30)

-ρ(sparsevec([7], [-1.0], 10), sol_PLAD01_BND30)

plot(sol_PLAD01_BND30, vars=(0, 1), xlab=L"t", ylab=L"x_1", title="PLAD01 - BND30", lw=0.1)
plot!(x->x, x->-30., 0., 20., linewidth=2, color="red", linestyle=:dash, leg=nothing)


text:
module Utls
using Languages
stop_words = stopwords(Languages.English())

mutable struct SentenceElements
    elements_dct::Dict{String,Int}
end

"""
    tokenize(sentence) -> Vector{SubString{String}}
Perform preprocessing such as erasing symbols, converting uppercase letters to lowercase, word segmentation, etc.
# Arguments
- `sentence`: A sentence in a document.
"""
function tokenize(sentence)
    sentence = replace(sentence, r"[“”\"\"();:,.'?!]" => "")
    sentence = lowercase(sentence)
    return split.(sentence, " ")
end

"""
    count_elements(sequence) -> Dict{String, Int64}
Counts the number of elements per token from a word segmented sentence.
# Arguments
- `sequence`: Word Segmented Sequence.
"""
function count_elements(sequence)
    elements_dct = Dict{String,Int}()
    for i in sequence
        if !in(i, stop_words)
            if i in keys(elements_dct)
                elements_dct[i] += 1
            else
                elements_dct[i] = 1
            end
        end
    end
    return elements_dct
end

"""
    merge_elements(dct)
Merge a dictionary that counts tokens in a sentence by element.
# Arguments
- `dct`: A dictionary with the number of elements counted for each token.
"""
function merge_elements(se::SentenceElements, dct)
    for i in keys(dct)
        if i in keys(se.elements_dct)
            se.elements_dct[i] += dct[i]
        else
            se.elements_dct[i] = dct[i]
        end
    end
end

"""
    calculate_cosin_similarity(elements_dct_1, elements_dct_2) -> Float64
Calculates the cosine similarity between two dictionaries.
# Arguments
- `elements_dct_1`: Token dictionary contained in the left block including the reference sentence.
- `elements_dct_2`: Token dictionary contained in the block to the right of the reference sentence.
"""
function calculate_cosin_similarity(elements_dct_1, elements_dct_2)
    n1 = sqrt(sum([i * i for i in values(elements_dct_1)]))
    n2 = sqrt(sum([i * i for i in values(elements_dct_2)]))
    common_key = intersect(keys(elements_dct_1), keys(elements_dct_2))

    if isempty(common_key)
        num = 0
    else
        num = sum(elements_dct_1[k] * elements_dct_2[k] for k in common_key)
    end

    try
        if n1 * n2 < 1e-9
            return 0
        end
        return num / (n1 * n2)
    catch
        return 0
    end
end
end


text:
include("solver_preamble.jl")
import SDPA
factory = with_optimizer(SDPA.Optimizer)
config = MOI.Test.TestConfig(atol=1e-5, rtol=1e-5, query=false)
@testset "Linear" begin
    Tests.linear_test(factory, config, [
        # The termination status is
        # `INFEASIBLE_OR_UNBOUNDED` instead of `INFEASIBLE`.
        "dsos_horn",
        # cholesky miss condition :: not positive definite :: info = 32 :: line 785 in sdpa_linear.cpp
        # There are some possibilities. :: line 786 in sdpa_linear.cpp
        # 1. SDPA terminates due to inaccuracy of numerical error :: line 787 in sdpa_linear.cpp
        # 2. The input problem may not have (any) interior-points :: line 788 in sdpa_linear.cpp
        # 3. Input matrices are linearly dependent :: line 789 in sdpa_linear.cpp
        "dsos_concave_then_convex_cubic",
        # The termination status is
        # `SLOW_PROGRESS` instead of `INFEASIBLE`.
        # cholesky miss condition :: not positive definite :: info = 7 :: line 785 in sdpa_linear.cpp
        # There are some possibilities. :: line 786 in sdpa_linear.cpp
        # 1. SDPA terminates due to inaccuracy of numerical error :: line 787 in sdpa_linear.cpp
        # 2. The input problem may not have (any) interior-points :: line 788 in sdpa_linear.cpp
        # 3. Input matrices are linearly dependent :: line 789 in sdpa_linear.cpp
        "dsos_bivariate_quadratic", "dsos_univariate_quadratic"])
end
@testset "SDP" begin
    # With `sos_concave_then_convex_cubic`, we have the same error than with
    # `dsos_concave_then_convex_cubic`.
    Tests.sd_test(factory, config, [
        "sos_concave_then_convex_cubic",
        "quartic_infeasible_lyapunov_switched_system",
        # Strange behavior : primal < dual :: line 144 in sdpa_solve.cpp
        # pdINF criteria :: line 1192 in sdpa_parts.cpp
        "maxcut",
        # Strange behavior : primal < dual :: line 144 in sdpa_solve.cpp
        # dUNBD criteria :: line 1220 in sdpa_parts.cpp
        # dUNBD criteria :: line 1220 in sdpa_parts.cpp
        # dUNBD criteria :: line 1220 in sdpa_parts.cpp
        # Step length is too small.  :: line 198 in sdpa_dataset.cpp
        # cannot move: step length is too short :: line 176 in sdpa_solve.cpp
        "sos_horn"
    ])
end


text:
using RedPitayaDAQServer
using PyPlot

# obtain the URL of the RedPitaya
include("config.jl")

rp = RedPitayaCluster([URLs[1]])

dec = 64
modulus = 4800
base_frequency = 125000000
periods_per_step = 5
samples_per_period = div(modulus, dec)
periods_per_frame = 50 # about 0.5 s frame length
frame_period = dec*samples_per_period*periods_per_frame / base_frequency
slow_dac_periods_per_frame = div(50, periods_per_step)

decimation(rp, dec)
samplesPerPeriod(rp, samples_per_period)
periodsPerFrame(rp, periods_per_frame)
passPDMToFastDAC(master(rp), true)

slowDACStepsPerFrame(rp, slow_dac_periods_per_frame)
numSlowDACChan(master(rp), 1)
lut = collect(range(0,0.7,length=slow_dac_periods_per_frame))
setSlowDACLUT(master(rp), lut)

modeDAC(rp, "STANDARD")
frequencyDAC(rp,1,1, base_frequency / modulus)

freq = frequencyDAC(rp,1,1)
println(" frequency = $(freq)")
signalTypeDAC(rp, 1 , "SINE")
amplitudeDACNext(rp, 1, 1, 0.2)
phaseDAC(rp, 1, 1, 0.0 ) # Phase has to be given in between 0 and 1

ramWriterMode(rp, "TRIGGERED")
triggerMode(rp, "INTERNAL")
masterTrigger(rp, false)
startADC(rp)
masterTrigger(rp, true)

sleep(0.1)

currFr = enableSlowDAC(rp, true, 2, frame_period, 0.5)
uCurrentFrame = readFrames(rp, currFr, 2)

figure(1)
clf()
plot(vec(uCurrentFrame[:,1,:,:]))
plot(vec(uCurrentFrame[:,2,:,:]))
legend(("Rx1", "Rx2"))

savefig("images/slowDAC.png")

stopADC(rp)

max_stars_repo_name: Fypsilonn/RedPitayaDAQServer

max_issues_repo_name: Fypsilonn/RedPitayaDAQServer

max_forks_repo_name: Fypsilonn/RedPitayaDAQServer