text:
### Convenience access to Linear Operators
# Mostly copied over from LowRankApprox.jl
# This extends the linear operators from LowRankApprox to include custom getindex routines
# We call these structs LinearMaps

const AbstractLinOp = AbstractLinearOperator
const AbstractMatOrLinOp{T} = Union{AbstractMatrix{T}, AbstractLinOp{T}}

  mutable struct LinearMap{T, Mul<:Function, Mulc<:Function, GetIdx<:Function} <: AbstractLinOp{T}
    m::Int
    n::Int
    mul!::Mul
    mulc!::Mulc
    getidx::GetIdx
    LinearMap{T}(m, n, mul!, mulc!, getidx) where T = new{T, typeof(mul!), typeof(mulc!), typeof(getidx)}(m, n, mul!, mulc!, getidx)
  end
const LinMap = LinearMap

mutable struct HermitianLinearMap{T, Mul<:Function, GetIdx<:Function} <: AbstractLinOp{T}
  n::Int
  mul!::Mul
  getidx::GetIdx
  HermitianLinearMap{T}(m, n, mul!, getidx) where T = new{T, typeof(mul!), typeof(getidx)}(m, n, mul!, getidx)
end
const HermLinMap = HermitianLinearMap

function LinMap(A)
  ishermitian(A) && return HermLinMap(A)
  T = eltype(A)
  m, n = size(A)
  ml!  = (y, _, x) ->  mul!(y, A, x)
  mulc! = (y, _, x) -> mul!(y, A', x)
  getidx = (i, j) -> A[i,j]
  LinMap{T}(m, n, ml!, mulc!, getidx)
end

function HermLinMap(A)
  T = eltype(A)
  m, n = size(A)
  m == n || throw(DimensionMismatch)
  ml! = (y, _, x) ->  mul!(y, A, x)
  getidx = (i, j) -> A[i,j]
  HermLinMap{T}(n, ml!, getidx)
end

convert(::Type{LinMap}, A::HermLinMap) = convert(LinMap{eltype(A)}, A)
convert(::Type{LinMap{T}}, A::HermLinMap{T}) where T = LinMap{T}(A.n, A.n, A.mul!, A.mul!, A.getidx)

adjoint(A::LinMap{T}) where {T} = LinMap{T}(A.n, A.m, A.mulc!, A.mul!, (i,j) -> conj(A.getidx(j,i)))
adjoint(A::HermLinMap) = A

#Matrix(A::AbstractLinOp{T}) where {T} = A*Matrix{T}(I, size(A)...)

getindex(A::LinMap, ::Colon, ::Colon) = A.getidx(1:A.m, 1:A.n)
getindex(A::LinMap, rows, ::Colon) = A.getidx(rows, 1:A.n)
getindex(A::LinMap, ::Colon, cols) = A.getidx(1:A.m, cols)
getindex(A::LinMap, rows, cols) = A.getidx(rows, cols)

getindex(A::HermLinMap, ::Colon, ::Colon) = A.getidx(1:A.m, 1:A.n)
getindex(A::HermLinMap, rows, ::Colon) = A.getidx(rows, 1:A.n)
getindex(A::HermLinMap, ::Colon, cols) = A.getidx(1:A.m, cols)
getindex(A::HermLinMap, rows, cols) = A.getidx(rows, cols)

ishermitian(::LinMap) = false
ishermitian(::HermLinMap) = true
issymmetric(::LinMap) = false
issymmetric(A::HermLinMap) = isreal(A)

size(A::LinMap) = (A.m, A.n)
size(A::LinMap, dim::Integer) = dim == 1 ? A.m : (dim == 2 ? A.n : 1)
size(A::HermLinMap) = (A.n, A.n)
size(A::HermLinMap, dim::Integer) = (dim == 1 || dim == 2) ? A.n : 1

function transpose(A::LinMap{T}) where T
  n, m = size(A)
  mul!  = (y, L, x) -> (A.mulc!(y, L, conj(x)); conj!(y))
  mulc! = (y, L, x) -> ( A.mul!(y, L, conj(x)); conj!(y))
  getidx = (i,j) -> A.getidx(j,i)
  LinMap{T}(m, n, mul!, mulc!, getidx)
end
function transpose(A::HermLinMap{T}) where T
  n = size(A, 1)
  mul! = (y, L, x) -> (A.mul!(y, L, conj(x)); conj!(y))
  getidx = (i,j) -> A.getidx(j,i)
  HermLinMap{T}(n, mul!, getidx)
end

# # copied over from LowRankApprox.jl
# mul!(C, A::LinMap, B::AbstractVecOrMat) = A.mul!(C, A, B)
# mul!(C, A::Adjoint{<:Any,<:LinMap}, B::AbstractVecOrMat) = parent(A).mulc!(C, parent(A), B)
# mul!(C, A::Adjoint{<:Any,<:HermLinMap}, B::AbstractVecOrMat) = mul!(C, parent(A), B)
# mul!(C, A::AbstractMatrix, B::Adjoint{<:Any,<:AbstractLinMap}) = adjoint!(C, parent(B)*A')
# mul!(C, A::AbstractMatrix, B::AbstractLinMap) = adjoint!(C, B'*A')

# *(A::LinMap{T}, x::AbstractVector) where T = (y = Array{T}(undef, size(A,1)); mul!(y, A, x))
# *(A::LinMap{T}, B::AbstractMatrix) where T = (C = Array{T}(undef, size(A,1), size(B,2)); mul!(C, A, B))
# *(A::AbstractMatrix, B::LinMap{T}) where T = (C = Array{T}(undef, size(A,1), size(B,2)); mul!(C, A, B))

# scalar multiplication/division
for (f, g) in ((:(A::LinMap), :(c::Number)), (:(c::Number), :(A::LinMap)))
  @eval begin
    function *($f, $g)
      T = eltype(A)
      m, n = size(A)
      sc_mul!  = (y, _, x) -> ( mul!(y, A, x); lmul!(c, y))
      sc_mulc! = (y, _, x) -> (mul!(y, A', x); lmul!(c, y))
      sc_getidx = (i,j) -> c*A.getidx(i,j)
      LinMap{T}(m, n, sc_mul!, sc_mulc!, sc_getidx)
    end
  end
end

for (f, g) in ((:(A::HermLinMap), :(c::Number)), (:(c::Number), :(A::HermLinMap)))
  @eval begin
    function *($f, $g)
      T = eltype(A)
      n = size(A, 1)
      sc_mul! = (y, _, x) -> (mul!(y, A, x); lmul!(c, y))
      sc_getidx = (i,j) -> c*A.getidx(i,j)
      HermLinMap{T}(n, sc_mul!, sc_getidx)
    end
  end
end

# operator addition/subtraction
for (f, a) in ((:+, 1), (:-, -1))
  @eval begin
    function $f(A::HermLinMap{T}, B::HermLinMap{T}) where T
      size(A) == size(B) || throw(DimensionMismatch)
      n = size(A, 1)
      alpha = T($a)
      mul! = gen_linop_axpy(A, B, alpha)
      getidx = A.getidx(i,j) + B.getidx(i,j)
      HermLinMap{T}(n, mul!, getidx)
    end
  end
end


text:
using Test

include("definitions/arraypartition.jl")
include("definitions/compose.jl")

include("utilities/test_ad.jl")
include("utilities/test_iteration_tools.jl")
include("utilities/test_conjugate.jl")
include("utilities/test_fb_tools.jl")

include("accel/test_lbfgs.jl")
include("accel/test_anderson.jl")
include("accel/test_nesterov.jl")
include("accel/test_broyden.jl")

include("problems/test_equivalence.jl")
include("problems/test_elasticnet.jl")
include("problems/test_lasso_small.jl")
include("problems/test_lasso_small_strongly_convex.jl")
include("problems/test_lasso_small_v_split.jl")
include("problems/test_lasso_small_h_split.jl")
include("problems/test_linear_programs.jl")
include("problems/test_sparse_logistic_small.jl")
include("problems/test_nonconvex_qp.jl")
include("problems/test_verbose.jl")

ext: jl


text:
@doc raw"""
    CholeskySpace{N} <: AbstractManifold{ℝ}

The manifold of lower triangular matrices with positive diagonal and
a metric based on the cholesky decomposition. The formulae for this manifold
are for example summarized in Table 1 of [^Lin2019].

# Constructor

    CholeskySpace(n)

Generate the manifold of $n× n$ lower triangular matrices with positive diagonal.

[^Lin2019]:
    > Lin, Zenhua: "Riemannian Geometry of Symmetric Positive Definite Matrices via
    > Cholesky Decomposition", arXiv: [1908.09326](https://arxiv.org/abs/1908.09326).
"""
struct CholeskySpace{N} <: AbstractManifold{ℝ} end

CholeskySpace(n::Int) = CholeskySpace{n}()

@doc raw"""
    check_point(M::CholeskySpace, p; kwargs...)

Check whether the matrix `p` lies on the [`CholeskySpace`](@ref) `M`, i.e.
it's size fits the manifold, it is a lower triangular matrix and has positive
entries on the diagonal.
The tolerance for the tests can be set using the `kwargs...`.
"""
function check_point(M::CholeskySpace, p; kwargs...)
    if size(p) != representation_size(M)
        return DomainError(
            size(p),
            "The point $(p) does not lie on $(M), since its size is not $(representation_size(M)).",
        )
    end
    if !isapprox(norm(strictlyUpperTriangular(p)), 0.0; kwargs...)
        return DomainError(
            norm(UpperTriangular(p) - Diagonal(p)),
            "The point $(p) does not lie on $(M), since it strictly upper triangular nonzero entries",
        )
    end
    if any(diag(p) .<= 0)
        return DomainError(
            min(diag(p)...),
            "The point $(p) does not lie on $(M), since it hast nonpositive entries on the diagonal",
        )
    end
    return nothing
end

"""
    check_vector(M::CholeskySpace, p, X; kwargs... )

Check whether `v` is a tangent vector to `p` on the [`CholeskySpace`](@ref) `M`, i.e.
after [`check_point`](@ref)`(M,p)`, `X` has to have the same dimension as `p`
and a symmetric matrix.
The tolerance for the tests can be set using the `kwargs...`.
"""
function check_vector(M::CholeskySpace, p, X; kwargs...)
    if size(X) != representation_size(M)
        return DomainError(
            size(X),
            "The vector $(X) is not a tangent to a point on $(M) since its size does not match $(representation_size(M)).",
        )
    end
    if !isapprox(norm(strictlyUpperTriangular(X)), 0.0; kwargs...)
        return DomainError(
            norm(UpperTriangular(X) - Diagonal(X)),
            "The matrix $(X) is not a tangent vector at $(p) (represented as an element of the Lie algebra) since it is not lower triangular.",
        )
    end
    return nothing
end

@doc raw"""
    distance(M::CholeskySpace, p, q)

Compute the Riemannian distance on the [`CholeskySpace`](@ref) `M` between two
matrices `p`, `q` that are lower triangular with positive diagonal. The formula
reads

````math
d_{\mathcal M}(p,q) = \sqrt{\sum_{i>j} (p_{ij}-q_{ij})^2 +
\sum_{j=1}^m (\log p_{jj} - \log q_{jj})^2
}
````
"""
function distance(::CholeskySpace, p, q)
    return sqrt(
        sum((strictlyLowerTriangular(p) - strictlyLowerTriangular(q)) .^ 2) +
        sum((log.(diag(p)) - log.(diag(q))) .^ 2),
    )
end

@doc raw"""
    exp(M::CholeskySpace, p, X)

Compute the exponential map on the [`CholeskySpace`](@ref) `M` emanating from the lower
triangular matrix with positive diagonal `p` towards the lower triangular matrix `X`
The formula reads

````math
\exp_p X = ⌊ p ⌋ + ⌊ X ⌋ + \operatorname{diag}(p)
\operatorname{diag}(p)\exp\bigl( \operatorname{diag}(X)\operatorname{diag}(p)^{-1}\bigr),
````

where $⌊\cdot⌋$ denotes the strictly lower triangular matrix,
and $\operatorname{diag}$ extracts the diagonal matrix.
"""
exp(::CholeskySpace, ::Any...)

function exp!(::CholeskySpace, q, p, X)
    q .= (
        strictlyLowerTriangular(p) +
        strictlyLowerTriangular(X) +
        Diagonal(diag(p)) * Diagonal(exp.(diag(X) ./ diag(p)))
    )
    return q
end

@doc raw"""
    inner(M::CholeskySpace, p, X, Y)

Compute the inner product on the [`CholeskySpace`](@ref) `M` at the
lower triangular matric with positive diagonal `p` and the two tangent vectors
`X`,`Y`, i.e they are both lower triangular matrices with arbitrary diagonal.
The formula reads

````math
g_p(X,Y) = \sum_{i>j} X_{ij}Y_{ij} + \sum_{j=1}^m X_{ii}Y_{ii}p_{ii}^{-2}
````
"""
function inner(::CholeskySpace, p, X, Y)
    return (
        sum(strictlyLowerTriangular(X) .* strictlyLowerTriangular(Y)) +
        sum(diag(X) .* diag(Y) ./ (diag(p) .^ 2))
    )
end

@doc raw"""
    log(M::CholeskySpace, X, p, q)

Compute the logarithmic map on the [`CholeskySpace`](@ref) `M` for the geodesic emanating
from the lower triangular matrix with positive diagonal `p` towards `q`.
The formula reads

````math
\log_p q = ⌊ p ⌋ - ⌊ q ⌋ + \operatorname{diag}(p)\log\bigl(\operatorname{diag}(q)\operatorname{diag}(p)^{-1}\bigr),
````

where $⌊\cdot⌋$ denotes the strictly lower triangular matrix,
and $\operatorname{diag}$ extracts the diagonal matrix.
"""
log(::Cholesky, ::Any...)

function log!(::CholeskySpace, X, p, q)
    return copyto!(
        X,
        strictlyLowerTriangular(q) - strictlyLowerTriangular(p) +
        Diagonal(diag(p) .* log.(diag(q) ./ diag(p))),
    )
end

@doc raw"""
    manifold_dimension(M::CholeskySpace)

Return the manifold dimension for the [`CholeskySpace`](@ref) `M`, i.e.

````math
    \dim(\mathcal M) = \frac{N(N+1)}{2}.
````
"""
@generated manifold_dimension(::CholeskySpace{N}) where {N} = div(N * (N + 1), 2)

@doc raw"""
    reporesentation_size(M::CholeskySpace)

Return the representation size for the [`CholeskySpace`](@ref)`{N}` `M`, i.e. `(N,N)`.
"""
@generated representation_size(::CholeskySpace{N}) where {N} = (N, N)

Base.show(io::IO, ::CholeskySpace{N}) where {N} = print(io, "CholeskySpace($(N))")

# two small helpers for strictly lower and upper triangulars
strictlyLowerTriangular(p) = LowerTriangular(p) - Diagonal(diag(p))

strictlyUpperTriangular(p) = UpperTriangular(p) - Diagonal(diag(p))

@doc raw"""
    vector_transport_to(M::CholeskySpace, p, X, q, ::ParallelTransport)

Parallely transport the tangent vector `X` at `p` along the geodesic to `q`
on the [`CholeskySpace`](@ref) manifold `M`. The formula reads

````math
\mathcal P_{q←p}(X) = ⌊ X ⌋
+ \operatorname{diag}(q)\operatorname{diag}(p)^{-1}\operatorname{diag}(X),
````

where $⌊\cdot⌋$ denotes the strictly lower triangular matrix,
and $\operatorname{diag}$ extracts the diagonal matrix.
"""
vector_transport_to(::CholeskySpace, ::Any, ::Any, ::Any, ::ParallelTransport)

function vector_transport_to!(::CholeskySpace, Y, p, X, q, ::ParallelTransport)
    return copyto!(Y, strictlyLowerTriangular(p) + Diagonal(diag(q) .* diag(X) ./ diag(p)))
end

@doc raw"""
    zero_vector(M::CholeskySpace, p)

Return the zero tangent vector on the [`CholeskySpace`](@ref) `M` at `p`.
"""
zero_vector(::CholeskySpace, ::Any...)

zero_vector!(M::CholeskySpace, X, p) = fill!(X, 0)


text:
# test if fmd in getSample and factor eval is working right

using IncrementalInference
using Test

import IncrementalInference: getSample

##

struct MyFactor{T <: SamplableBelief} <: IIF.AbstractRelativeRoots
  Z::T
  # specialSampler approach will be deprecated
  # specialSampler::Function
end

function getSample(cf::CalcFactor{<:MyFactor}, N::Int=1)
  @warn "getSample(cf::CalcFactor{<:MyFactor},::Int) does not get hypo sub-selected FMD data"
  @show DFG.getLabel.(cf.metadata.fullvariables)
  # @assert DFG.getLabel.(fmd_[1].fullvariables) |> length < 3 "this factor is only between two variables"
  return (reshape(rand(cf.factor.Z, N),1,N),)
end

function (cf::CalcFactor{<:MyFactor})(z, X1, X2)
  @assert DFG.getLabel.(cf.metadata.fullvariables) |> length < 3 "this factor is only between two variables"

  # just a linear difference to complete the test
  return X2 - (X1 + z)
end

##

@testset "test FactorMetadata is properly populated" begin

##

fg = initfg()

addVariable!(fg, :x0, ContinuousScalar)
addVariable!(fg, :x1_a, ContinuousScalar)
addVariable!(fg, :x1_b, ContinuousScalar)

addFactor!(fg, [:x0], Prior(Normal()))

# create the object and add it to the graph
mf = MyFactor( Normal(10,1) )

##

# this sampling might error
addFactor!(fg, [:x0;:x1_a;:x1_b], mf, multihypo=[1;1/2;1/2])

##

meas = sampleFactor(fg, :x0x1_ax1_bf1, 10)

##

solveTree!(fg);

##

end

#


text:
"""
$(SIGNATURES)

Solve the time-dependent Redfield equation defined by `A` for a total evolution time `tf`.

...
# Arguments
- `A::Annealing`: the `Annealing`/`Evolution` object.
- `tf::Real`: the total evolution time.
- `unitary`: precomputed unitary operator of the corresponding closed-system evolution.
- `vectorize::Bool = false`: whether to vectorize the density matrix.
- `int_atol = 1e-8`: the absolute error tolerance for integration.
- `int_rtol = 1e-6`: the relative error tolerance for integration.
- `Ta = tf`: the timescale of the backward integration.
- `kwargs`: other keyword arguments supported by `DifferentialEquations`.
...
"""
function solve_redfield(
    A::Annealing,
    tf::Real,
    unitary;
    vectorize::Bool=false,
    int_atol=1e-8,
    int_rtol=1e-6,
    Ta=tf,
    kwargs...,
)
    u0 = build_u0(A.u0, :m, vectorize=vectorize)
    L = OpenQuantumBase.redfield_from_interactions(A.interactions, unitary, Ta, int_atol, int_rtol)
    R = DiffEqLiouvillian(A.H, [], L, size(A.H, 1))

    update_func = function (A, u, p, t)
        update_vectorized_cache!(A, p.L, p, t)
    end

    if vectorize == false
        ff = ODEFunction{true}(R)
    else
        cache = vectorize_cache(get_cache(A.H))
        diff_op = DiffEqArrayOperator(cache, update_func=update_func)
        jac_cache = similar(cache)
        jac_op = DiffEqArrayOperator(jac_cache, update_func=update_func)
        ff = ODEFunction(diff_op, jac_prototype=jac_op)
    end

    p = ODEParams(R, float(tf), A.annealing_parameter)
    prob = ODEProblem(ff, u0, (0.0, float(tf)), p)
    solve(prob; alg_hints=[:nonstiff], kwargs...)
end

"""
$(SIGNATURES)

Solve the time-dependent coarse-grained master equation (CGME) defined by `A` for a total evolution time `tf`.

...
# Arguments
- `A::Annealing`: the `Annealing`/`Evolution` object.
- `tf::Real`: the total evolution time.
- `unitary`: precomputed unitary operator of the corresponding closed-system evolution.
- `vectorize::Bool = false`: whether to vectorize the density matrix.
- `Ta = nothing`: coarse-graining time. If set to nothing, the solver will automatically choose the value.
- `int_atol = 1e-8`: the absolute error tolerance for integration.
- `int_rtol = 1e-6`: the relative error tolerance for integration.
- `kwargs`: other keyword arguments supported by `DifferentialEquations`.
...
"""
function solve_cgme(
    A::Annealing,
    tf::Real,
    unitary;
    vectorize::Bool=false,
    Ta=nothing,
    int_atol=1e-8,
    int_rtol=1e-6,
    kwargs...,
)
    u0 = build_u0(A.u0, :m, vectorize=vectorize)
    L = OpenQuantumBase.cg_from_interactions(A.interactions, unitary, tf, Ta,
            int_atol, int_rtol)
    R = DiffEqLiouvillian(A.H, [], L, size(A.H, 1))
    update_func = function (A, u, p, t)
        update_vectorized_cache!(A, p.L, p, t)
    end

    if vectorize == false
        ff = ODEFunction{true}(R)
    else
        cache = vectorize_cache(get_cache(A.H))
        diff_op = DiffEqArrayOperator(cache, update_func=update_func)
        jac_cache = similar(cache)
        jac_op = DiffEqArrayOperator(jac_cache, update_func=update_func)
        ff = ODEFunction(diff_op, jac_prototype=jac_op)
    end

    p = ODEParams(R, float(tf), A.annealing_parameter)
    prob = ODEProblem(ff, u0, (0.0, float(tf)), p)
    solve(prob; alg_hints=[:nonstiff], kwargs...)
end

"""
$(SIGNATURES)

Solve the time-dependent universal Lindblad equation (ULE) defined by `A` for a total evolution time `tf`.

...
# Arguments
- `A::Annealing`: the `Annealing`/`Evolution` object.
- `tf::Real`: the total evolution time.
- `unitary`: precomputed unitary operator of the corresponding closed-system evolution.
- `vectorize::Bool = false`: whether to vectorize the density matrix.
- `int_atol = 1e-8`: the absolute error tolerance for integration.
- `int_rtol = 1e-6`: the relative error tolerance for integration.
- `Ta = tf`: the timescale of the integration region.
- `kwargs`: other keyword arguments supported by `DifferentialEquations`.
...
"""
function solve_ule(
    A::Annealing,
    tf::Real,
    unitary;
    vectorize::Bool=false,
    int_atol=1e-8,
    int_rtol=1e-6,
    Ta=tf,
    kwargs...,
)
    u0 = build_u0(A.u0, :m, vectorize=vectorize)
    L = OpenQuantumBase.ule_from_interactions(A.interactions, unitary, Ta, int_atol, int_rtol)
    R = DiffEqLiouvillian(A.H, [], L, size(A.H, 1))

    update_func = function (A, u, p, t)
        update_vectorized_cache!(A, p.L, p, t)
        error("Vectorization of ULE is not currently supported.")
    end

    if vectorize == false
        ff = ODEFunction{true}(R)
    else
        cache = vectorize_cache(get_cache(A.H))
        diff_op = DiffEqArrayOperator(cache, update_func=update_func)
        jac_cache = similar(cache)
        jac_op = DiffEqArrayOperator(jac_cache, update_func=update_func)
        ff = ODEFunction(diff_op, jac_prototype=jac_op)
    end

    p = ODEParams(R, float(tf), A.annealing_parameter)
    prob = ODEProblem(ff, u0, (0.0, float(tf)), p)
    solve(prob; alg_hints=[:nonstiff], kwargs...)
end

function build_ensemble_redfield(
    A::Annealing,
    tf,
    unitary,
    output_func,
    reduction;
    vectorize::Bool=false,
    int_atol=1e-8,
    int_rtol=1e-6,
    Ta=tf,
    initializer=DEFAULT_INITIALIZER,
    save_positions=(false, false),
    kwargs...,
)
    u0 = build_u0(A.u0, :m, vectorize=vectorize)
    stocs = OpenQuantumBase.fluctuator_from_interactions(A.interactions)
    reds = OpenQuantumBase.redfield_from_interactions(A.interactions, unitary, Ta, int_atol, int_rtol)
    if isempty(stocs)
        error("No stochastic bath detected. Use the normal Redfeidl solver instead.")
    end
    # set the initializer
    initializer =
        initializer == DEFAULT_INITIALIZER ? (x, y) -> rand([-1, 1], x, y) :
        initializer

    if length(stocs) == 1
        cb = FluctuatorCallback(stocs[1], initializer, save_positions)
    else
        cb = CallbackSet([FluctuatorCallback(f, initializer, save_positions) for f in stocs]...)
    end
    R = DiffEqLiouvillian(A.H, [], [reds; stocs], size(A.H, 1))
    p = ODEParams(R, float(tf), A.annealing_parameter)

    update_func = function (A, u, p, t)
        update_vectorized_cache!(A, p.L, p, t)
    end

    if vectorize == false
        ff = ODEFunction{true}(R)
    else
        cache = vectorize_cache(get_cache(A.H))
        diff_op = DiffEqArrayOperator(cache, update_func=update_func)
        jac_cache = similar(cache)
        jac_op = DiffEqArrayOperator(jac_cache, update_func=update_func)
        ff = ODEFunction(diff_op, jac_prototype=jac_op)
    end

    prob = ODEProblem{true}(ff, u0, (0.0, float(tf)), p, callback=cb)

    ensemble_prob =
        EnsembleProblem(prob; output_func=output_func, reduction=reduction)

    ensemble_prob
end


text:
exists(path) = stat(path).inode != 0

"""this package assumes [x,y] indexing, with x varying fastest and
origin bottom left.  this is not how julia displays matrices (not at
all), so things can be confusing if you use "print" on raw structures.
this function converts for display."""
fix{T}(x::Array{T, 2}) = x'[size(x)[1]:-1:1,:]


text:
### These are generated by a) uncommenting import_sympy() in __init__ b) uncomment println parts in import_from
expj(ex::SymbolicObject, args...; kwargs...)=getproperty(mpmath, :expj)(ex, Sym.(args)...; kwargs...); export expj
fac(ex::SymbolicObject, args...; kwargs...)=getproperty(mpmath, :fac)(ex, Sym.(args)...; kwargs...); export fac
nint(ex::SymbolicObject, args...; kwargs...)=getproperty(mpmath, :nint)(ex, Sym.(args)...; kwargs...); export nint
Base.ceil(ex::SymbolicObject, args...; kwargs...)=getproperty(mpmath, :ceil)(ex, Sym.(args)...; kwargs...)
fib(ex::SymbolicObject, args...; kwargs...)=getproperty(mpmath, :fib)(ex, Sym.(args)...; kwargs...); export fib
monitor(ex::SymbolicObject, args...; kwargs...)=getproperty(mpmath, :monitor)(ex, Sym.(args)...; kwargs...); export monitor
Base.cospi(ex::SymbolicObject, args...; kwargs...)=getproperty(mpmath, :cospi)(ex, Sym.(args)...; kwargs...)
bernfrac(ex::SymbolicObject, args...; kwargs...)=getproperty(mpmath, :bernfrac)(ex, Sym.(args)...; kwargs...); export bernfrac
doctests(ex::SymbolicObject, args...; kwargs...)=getproperty(mpmath, :doctests)(ex, Sym.(args)...; kwargs...); export doctests
ei(ex::SymbolicObject, args...; kwargs...)=getproperty(mpmath, :ei)(ex, Sym.(args)...; kwargs...); export ei
Base.sinpi(ex::SymbolicObject, args...; kwargs...)=getproperty(mpmath, :sinpi)(ex, Sym.(args)...; kwargs...)
timing(ex::SymbolicObject, args...; kwargs...)=getproperty(mpmath, :timing)(ex, Sym.(args)...; kwargs...); export timing
rgamma(ex::SymbolicObject, args...; kwargs...)=getproperty(mpmath, :rgamma)(ex, Sym.(args)...; kwargs...); export rgamma
expjpi(ex::SymbolicObject, args...; kwargs...)=getproperty(mpmath, :expjpi)(ex, Sym.(args)...; kwargs...); export expjpi
ellipk(ex::SymbolicObject, args...; kwargs...)=getproperty(mpmath, :ellipk)(ex, Sym.(args)...; kwargs...); export ellipk
e1(ex::SymbolicObject, args...; kwargs...)=getproperty(mpmath, :e1)(ex, Sym.(args)...; kwargs...); export e1
SpecialFunctions.digamma(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :digamma)(ex, Sym.(args)...; kwargs...)
#Base.rem(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :rem)(ex, Sym.(args)...; kwargs...)
Base.asin(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :asin)(ex, Sym.(args)...; kwargs...)
Base.acosh(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :acosh)(ex, Sym.(args)...; kwargs...)
Base.im(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :im)(ex, Sym.(args)...; kwargs...)
Base.collect(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :collect)(ex, Sym.(args)...; kwargs...)
Base.Function(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :Function)(ex, Sym.(args)...; kwargs...)
SpecialFunctions.erfi(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :erfi)(ex, Sym.(args)...; kwargs...)
Base.floor(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :floor)(ex, Sym.(args)...; kwargs...)
Base.product(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :product)(ex, Sym.(args)...; kwargs...)
Base.gcd(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :gcd)(ex, Sym.(args)...; kwargs...)
Base.atan(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :atan)(ex, Sym.(args)...; kwargs...)
Base.sqrt(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :sqrt)(ex, Sym.(args)...; kwargs...)
Base.acsch(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :acsch)(ex, Sym.(args)...; kwargs...)
SpecialFunctions.besselj(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :besselj)(ex, Sym.(args)...; kwargs...)
Base.adjoint(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :adjoint)(ex, Sym.(args)...; kwargs...)
Base.asinh(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :asinh)(ex, Sym.(args)...; kwargs...)
SpecialFunctions.besselk(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :besselk)(ex, Sym.(args)...; kwargs...)
Base.binomial(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :binomial)(ex, Sym.(args)...; kwargs...)
Base.asec(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :asec)(ex, Sym.(args)...; kwargs...)
Base.exp(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :exp)(ex, Sym.(args)...; kwargs...)
Base.sech(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :sech)(ex, Sym.(args)...; kwargs...)
Base.cbrt(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :cbrt)(ex, Sym.(args)...; kwargs...)
Base.acsc(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :acsc)(ex, Sym.(args)...; kwargs...)
Base.factorial(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :factorial)(ex, Sym.(args)...; kwargs...)
Base.trunc(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :trunc)(ex, Sym.(args)...; kwargs...)
Base.acos(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :acos)(ex, Sym.(args)...; kwargs...)
SpecialFunctions.polygamma(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :polygamma)(ex, Sym.(args)...; kwargs...)
Base.tanh(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :tanh)(ex, Sym.(args)...; kwargs...)
SpecialFunctions.erfinv(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :erfinv)(ex, Sym.(args)...; kwargs...)
Base.sinh(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :sinh)(ex, Sym.(args)...; kwargs...)
SpecialFunctions.airybi(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :airybi)(ex, Sym.(args)...; kwargs...)
Base.asech(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :asech)(ex, Sym.(args)...; kwargs...)
SpecialFunctions.erfcinv(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :erfcinv)(ex, Sym.(args)...; kwargs...)
SpecialFunctions.besseli(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :besseli)(ex, Sym.(args)...; kwargs...)
Base.acot(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :acot)(ex, Sym.(args)...; kwargs...)
Base.coth(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :coth)(ex, Sym.(args)...; kwargs...)
SpecialFunctions.airyai(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :airyai)(ex, Sym.(args)...; kwargs...)
SpecialFunctions.erf(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :erf)(ex, Sym.(args)...; kwargs...)
Base.cosh(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :cosh)(ex, Sym.(args)...; kwargs...)
LinearAlgebra.diag(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :diag)(ex, Sym.(args)...; kwargs...)
Base.lcm(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :lcm)(ex, Sym.(args)...; kwargs...)
Base.zeros(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :zeros)(ex, Sym.(args)...; kwargs...)
Base.cot(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :cot)(ex, Sym.(args)...; kwargs...)
SpecialFunctions.zeta(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :zeta)(ex, Sym.(args)...; kwargs...)
Base.sign(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :sign)(ex, Sym.(args)...; kwargs...)
Base.permutedims(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :permutedims)(ex, Sym.(args)...; kwargs...)
Base.cos(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :cos)(ex, Sym.(args)...; kwargs...)
Base.transpose(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :transpose)(ex, Sym.(args)...; kwargs...)
Base.MathConstants.catalan(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :catalan)(ex, Sym.(args)...; kwargs...)
SpecialFunctions.erfc(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :erfc)(ex, Sym.(args)...; kwargs...)
SpecialFunctions.bessely(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :bessely)(ex, Sym.(args)...; kwargs...)
Base.diff(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :diff)(ex, Sym.(args)...; kwargs...)
Base.tan(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :tan)(ex, Sym.(args)...; kwargs...)
Base.decompose(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :decompose)(ex, Sym.(args)...; kwargs...)
SpecialFunctions.airyaiprime(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :airyaiprime)(ex, Sym.(args)...; kwargs...)
Base.csch(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :csch)(ex, Sym.(args)...; kwargs...)
Base.csc(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :csc)(ex, Sym.(args)...; kwargs...)
Base.ones(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :ones)(ex, Sym.(args)...; kwargs...)
SpecialFunctions.trigamma(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :trigamma)(ex, Sym.(args)...; kwargs...)
Base.prod(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :prod)(ex, Sym.(args)...; kwargs...)
SpecialFunctions.beta(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :beta)(ex, Sym.(args)...; kwargs...)
Base.sec(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :sec)(ex, Sym.(args)...; kwargs...)
Base.acoth(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :acoth)(ex, Sym.(args)...; kwargs...)
SpecialFunctions.airybiprime(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :airybiprime)(ex, Sym.(args)...; kwargs...)
Base.atanh(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :atanh)(ex, Sym.(args)...; kwargs...)
LinearAlgebra.det(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :det)(ex, Sym.(args)...; kwargs...)
SpecialFunctions.gamma(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :gamma)(ex, Sym.(args)...; kwargs...)
Base.reshape(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :reshape)(ex, Sym.(args)...; kwargs...)
Base.sin(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :sin)(ex, Sym.(args)...; kwargs...)
simplify(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :simplify)(ex, Sym.(args)...; kwargs...); export simplify
summation(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :summation)(ex, Sym.(args)...; kwargs...); export summation
solve(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :solve)(ex, Sym.(args)...; kwargs...); export solve
Max(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :Max)(ex, Sym.(args)...; kwargs...); export Max
unflatten(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :unflatten)(ex, Sym.(args)...; kwargs...); export unflatten
denom(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :denom)(ex, Sym.(args)...; kwargs...); export denom
nonlinsolve(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :nonlinsolve)(ex, Sym.(args)...; kwargs...); export nonlinsolve
cancel(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :cancel)(ex, Sym.(args)...; kwargs...); export cancel
solveset(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :solveset)(ex, Sym.(args)...; kwargs...); export solveset
DiracDelta(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :DiracDelta)(ex, Sym.(args)...; kwargs...); export DiracDelta
Or(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :Or)(ex, Sym.(args)...; kwargs...); export Or
conjugate(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :conjugate)(ex, Sym.(args)...; kwargs...); export conjugate
flatten(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :flatten)(ex, Sym.(args)...; kwargs...); export flatten
Not(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :Not)(ex, Sym.(args)...; kwargs...); export Not
integrate(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :integrate)(ex, Sym.(args)...; kwargs...); export integrate
roots(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :roots)(ex, Sym.(args)...; kwargs...); export roots
factor(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :factor)(ex, Sym.(args)...; kwargs...); export factor
pdsolve(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :pdsolve)(ex, Sym.(args)...; kwargs...); export pdsolve
Abs(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :Abs)(ex, Sym.(args)...; kwargs...); export Abs
Heaviside(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :Heaviside)(ex, Sym.(args)...; kwargs...); export Heaviside
nsimplify(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :nsimplify)(ex, Sym.(args)...; kwargs...); export nsimplify
isprime(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :isprime)(ex, Sym.(args)...; kwargs...); export isprime
apart(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :apart)(ex, Sym.(args)...; kwargs...); export apart
Min(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :Min)(ex, Sym.(args)...; kwargs...); export Min
intervals(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :intervals)(ex, Sym.(args)...; kwargs...); export intervals
intersection(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :intersection)(ex, Sym.(args)...; kwargs...); export intersection
line_integrate(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :line_integrate)(ex, Sym.(args)...; kwargs...); export line_integrate
real_roots(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :real_roots)(ex, Sym.(args)...; kwargs...); export real_roots
isolate(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :isolate)(ex, Sym.(args)...; kwargs...); export isolate
linsolve(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :linsolve)(ex, Sym.(args)...; kwargs...); export linsolve
Xor(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :Xor)(ex, Sym.(args)...; kwargs...); export Xor
real_root(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :real_root)(ex, Sym.(args)...; kwargs...); export real_root
nsolve(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :nsolve)(ex, Sym.(args)...; kwargs...); export nsolve
ln(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :ln)(ex, Sym.(args)...; kwargs...); export ln
rsolve(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :rsolve)(ex, Sym.(args)...; kwargs...); export rsolve
degree(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :degree)(ex, Sym.(args)...; kwargs...); export degree
prime(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :prime)(ex, Sym.(args)...; kwargs...); export prime
limit(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :limit)(ex, Sym.(args)...; kwargs...); export limit
And(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :And)(ex, Sym.(args)...; kwargs...); export And
root(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :root)(ex, Sym.(args)...; kwargs...); export root
rootof(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :rootof)(ex, Sym.(args)...; kwargs...); export rootof
ode_order(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :ode_order)(ex, Sym.(args)...; kwargs...); export ode_order
multiplicity(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :multiplicity)(ex, Sym.(args)...; kwargs...); export multiplicity
series(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :series)(ex, Sym.(args)...; kwargs...); export series
expand(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :expand)(ex, Sym.(args)...; kwargs...); export expand
hessian(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :hessian)(ex, Sym.(args)...; kwargs...); export hessian
srepr(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :srepr)(ex, Sym.(args)...; kwargs...); export srepr
nroots(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :nroots)(ex, Sym.(args)...; kwargs...); export nroots
interpolate(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :interpolate)(ex, Sym.(args)...; kwargs...); export interpolate
numer(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :numer)(ex, Sym.(args)...; kwargs...); export numer
cse(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :cse)(ex, Sym.(args)...; kwargs...); export cse
together(ex::SymbolicObject, args...; kwargs...)=getproperty(sympy, :together)(ex, Sym.(args)...; kwargs...); export together
Equality(ex::Number, args...; kwargs...)=getproperty(sympy, :Equality)(ex, Sym.(args)...; kwargs...); export Equality
Ne(ex::Number, args...; kwargs...)=getproperty(sympy, :Ne)(ex, Sym.(args)...; kwargs...); export Ne
LessThan(ex::Number, args...; kwargs...)=getproperty(sympy, :LessThan)(ex, Sym.(args)...; kwargs...); export LessThan
Gt(ex::Number, args...; kwargs...)=getproperty(sympy, :Gt)(ex, Sym.(args)...; kwargs...); export Gt
Eq(ex::Number, args...; kwargs...)=getproperty(sympy, :Eq)(ex, Sym.(args)...; kwargs...); export Eq
GreaterThan(ex::Number, args...; kwargs...)=getproperty(sympy, :GreaterThan)(ex, Sym.(args)...; kwargs...); export GreaterThan
Le(ex::Number, args...; kwargs...)=getproperty(sympy, :Le)(ex, Sym.(args)...; kwargs...); export Le
Lt(ex::Number, args...; kwargs...)=getproperty(sympy, :Lt)(ex, Sym.(args)...; kwargs...); export Lt
Unequality(ex::Number, args...; kwargs...)=getproperty(sympy, :Unequality)(ex, Sym.(args)...; kwargs...); export Unequality
Ge(ex::Number, args...; kwargs...)=getproperty(sympy, :Ge)(ex, Sym.(args)...; kwargs...); export Ge
StrictLessThan(ex::Number, args...; kwargs...)=getproperty(sympy, :StrictLessThan)(ex, Sym.(args)...; kwargs...); export StrictLessThan
StrictGreaterThan(ex::Number, args...; kwargs...)=getproperty(sympy, :StrictGreaterThan)(ex, Sym.(args)...; kwargs...); export StrictGreaterThan


text:
@testset "Test equations for Rayleigh resolution" begin
    f = 100e-3
    D = 50e-3

    @test rayleigh_criterion(calc_NA(f, D)) ≈ rayleigh_criterion(f, D)
    @test rayleigh_criterion(calc_NA(f, D), λ=500e-9) ≈ rayleigh_criterion(f, D)
    @test ~(rayleigh_criterion(calc_NA(f, D)) ≈ rayleigh_criterion(f, D, λ=200e-9))

end

@testset "Calc NA" begin
    @test calc_NA(100, 100) ≈ 0.5
    @test calc_NA(10e-9, 40e-9, 2) ≈ 4
end

max_stars_repo_name: JuliaPhysics/PhysicalOptics

max_issues_repo_name: JuliaPhysics/PhysicalOptics

max_forks_repo_name: JuliaPhysics/PhysicalOptics


text:
function stringtyped(x::Double64)
    str = string("Double64(", HI(x), ", ", LO(x), ")")
    return str
end
function stringtyped(x::Double32)
    str = string("Double32(", HI(x), ", ", LO(x), ")")
    return str
end
function stringtyped(x::Double16)
    str = string("Double16(", HI(x), ", ", LO(x), ")")
    return str
end

function stringtyped(x::Complex{DoubleFloat{Float64}})
    rea, ima = reim(x)
    str = string("ComplexDF64(", stringtyped(rea), ", ", stringtyped(imag(x)), ")")
    return str
end
function stringtyped(x::Complex{DoubleFloat{Float32}})
    rea, ima = reim(x)
    str = string("ComplexDF32(", stringtyped(rea), ", ", stringtyped(imag(x)), ")")
    return str
end
function stringtyped(x::Complex{DoubleFloat{Float16}})
    rea, ima = reim(x)
    str = string("ComplexDF16(", stringtyped(rea), ", ", stringtyped(imag(x)), ")")
    return str
end

"""
    BigFloatBits(::IEEEType)

significand bit precision for BigFloat when converting from DoubleFloat{IEEEType} into a string
""" BigFloatBits

"""
    BigFloatDigits(::IEEEType)

significant digits of BigFloat when converting to DoubleFloat{IEEEType} from a string
""" BigFloatDigits

# Float63 -> 106, Float32 -> 48, Float16 -> 22
minprecdblfloat(::Type{T}) where {T<:IEEEFloat} =
    (trailing_ones(Base.significand_mask(T))+1)*2

#=
  The effective precision of some DoubleFloat{T} may be larger than this
     where all bits betweeen the lsb of the HI oart and Lthe msb of the lo part
     are 0b0 and there are more than 22 zero bits interposed for eg Double64.
  When a DoubleFloat is the result of an arithmetic operation between DoubleFloats,
     the effective precision is almost never expected be larger than this.
=#
# Float63 -> 128, Float32 -> 58, Float16 -> 27
lrgprecdblfloat(::Type{T}) where {T<:IEEEFloat} =
    Int(ceil(Int, minprecdblfloat(T) * 1.2))

bigfloatbits(::Type{T}) where {T<:IEEEFloat} =
    nextpow(2, (trailing_ones(Base.significand_mask(T))+1)*2) - 1

BigFloatBits(::Type{Float64}) = 127 # bigfloatbits(Float64)
BigFloatBits(::Type{Float32}) =  63 # bigfloatbits(Float32)
BigFloatBits(::Type{Float16}) =  31 # bigfloatbits(Float16)

BigFloatDigits(::Type{Float64}) = 35 # bigfloatbits(Float64)
BigFloatDigits(::Type{Float32}) = 16 # bigfloatbits(Float32)
BigFloatDigits(::Type{Float16}) =  8 # bigfloatbits(Float16)

#= old
@inline BigFloatBits(::Type{Float64}) = 512
@inline BigFloatBits(::Type{Float32}) = 256
@inline BigFloatBits(::Type{Float16}) = 128
@inline BigFloatDigits(::Type{Float64}) = 35
@inline BigFloatDigits(::Type{Float32}) = 18
@inline BigFloatDigits(::Type{Float16}) =  9
=#

function string(x::DoubleFloat{T}) where {T<:IEEEFloat}
    !isfinite(HI(x)) && return string(HI(x))
    prec = precision(Base.BigFloat)
    setprecision(Base.BigFloat, BigFloatBits(T))
    bf = Base.BigFloat(x)
    bf = round(bf, digits=BigFloatBits(T), base=2)
    str = string(bf)
    if occursin('e', str)
       a, b = split(str, "e")
       n = min(length(a), BigFloatDigits(T)+1)
       str = string(a[1:n],"e",b)
    else
       n = min(length(str), BigFloatDigits(T)+1)
       str = str[1:n]
    end
    setprecision(Base.BigFloat, prec)
    return str
end

function string(x::Complex{DoubleFloat{T}}) where {T<:IEEEFloat}
    rea, ima = reim(x)
    s = signbit(ima) ? " - " : " + "
    imstr = isfinite(ima) ? "im" : "*im"
    str = string(string(rea), s, string(abs(ima)), imstr)
    return str
end


text:
### A Pluto.jl notebook ###
# v0.11.14

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ fafae38e-e852-11ea-1208-732b4744e4c2
md"_homework 0, version 2_"

# ╔═╡ 7308bc54-e6cd-11ea-0eab-83f7535edf25
# edit the code below to set your name and kerberos ID (i.e. email without @mit.edu)

student = (name = "Andreas Theodosiou", kerberos_id = "at694")

# press the ▶ button in the bottom right of this cell to run your edits
# or use Shift+Enter

# you might need to wait until all other cells in this notebook have completed running.
# scroll down the page to see what's up

# ╔═╡ cdff6730-e785-11ea-2546-4969521b33a7
md"""

Submission by: **_$(student.name)_** ($(student.kerberos_id)@mit.edu)
"""

# ╔═╡ a2181260-e6cd-11ea-2a69-8d9d31d1ef0e
md"""
# Homework 0: Getting up and running

First of all, **_welcome to the course!_** We are excited to teach you about real world applications of scientific computing, using the same tools that we work with ourselves.

Before we start next week, we'd like everyone to **submit this zeroth homework assignment**. It will not affect your grade, but it will help us get everything running smoothly when the course starts. If you're stuck or don't have much time, just fill in your name and ID and submit 🙂
"""

# ╔═╡ 094e39c8-e6ce-11ea-131b-07c4a1199edf

# ╔═╡ 31a8fbf8-e6ce-11ea-2c66-4b4d02b41995

# ╔═╡ 339c2d5c-e6ce-11ea-32f9-714b3628909c
md"## Exercise 1 - _Square root by Newton's method_

Computing the square of a number is easy -- you just multiply it with itself.

But how does one compute the square root of a number?

##### Algorithm:

Given: $x$

Output: $\sqrt{x}$

1. Take a guess `a`
1. Divide `x` by `a`
1. Set a = the average of `x/a` and `a`. (The square root must be between these two numbers. Why?)
1. Repeat until `x/a` is roughly equal to `a`. Return `a` as the square root.

In general, you will never get to the point where `x/a` is _exactly_ equal to `a`. So if our algorithm keeps going until `x/a == a`, then it will get stuck.

So instead, the algorithm takes a parameter `error_margin`, which is used to decide when `x/a` and `a` are close enough to halt.
"

# ╔═╡ 56866718-e6ce-11ea-0804-d108af4e5653
md"### Exercise 1.1

Step 3 in the algorithm sets the new guess to be the average of `x/a` and the old guess `a`.

This is because the square root must be between the numbers `x/a` and `a`. Why?
"

# ╔═╡ bccf0e88-e754-11ea-3ab8-0170c2d44628
ex_1_1 = md"""
The square root of a number $x$ is in the range $(1, x)$. Thus, for any value of $a$, where $a>0$, if $a>\sqrt{x}$ then the interval $(\frac{x}{a}, a)$ includes the square root and if $a<\sqrt{x}$ then the interval $(a, \frac{x}{a})$ includes the square root.
"""

# you might need to wait until all other cells in this notebook have completed running.
# scroll down the page to see what's up

# ╔═╡ e7abd366-e7a6-11ea-30d7-1b6194614d0a
if !(@isdefined ex_1_1)
	md"""Do not change the name of the variable - write you answer as `ex_1_1 = "..."`"""
end

# ╔═╡ d62f223c-e754-11ea-2470-e72a605a9d7e
md"### Exercise 1.2

Write a function newton_sqrt(x) which implements the above algorithm."

# ╔═╡ 4896bf0c-e754-11ea-19dc-1380bb356ab6
function newton_sqrt(x, error_margin=0.01, a=x / 2) # a=x/2 is the default value of `a`
	while (abs(x/a - a) >= error_margin)
		println("x is $(x)")
		println("a is $(a)")
		a = (a + x/a) / 2
	end
	return a # this is wrong, write your code here!
end

# ╔═╡ 7a01a508-e78a-11ea-11da-999d38785348
newton_sqrt(2)

# ╔═╡ 682db9f8-e7b1-11ea-3949-6b683ca8b47b
let
	result = newton_sqrt(2, 0.01)
	if !(result isa Number)
		md"""
!!! warning "Not a number"
    `newton_sqrt` did not return a number. Did you forget to write `return`?
		"""
	elseif abs(result - sqrt(2)) < 0.01
		md"""
!!! correct
    Well done!
		"""
	else
		md"""
!!! warning "Incorrect"
    Keep working on it!
		"""
	end
end

# ╔═╡ 088cc652-e7a8-11ea-0ca7-f744f6f3afdd
md"""
!!! hint
    `abs(r - s)` is the distance between `r` and `s`
"""

# ╔═╡ c18dce7a-e7a7-11ea-0a1a-f944d46754e5
md"""
!!! hint
    If you're stuck, feel free to cheat, this is homework 0 after all 🙃

    Julia has a function called `sqrt`
"""

# ╔═╡ 5e24d95c-e6ce-11ea-24be-bb19e1e14657
md"## Exercise 2 - _Sierpinksi's triangle_

Sierpinski's triangle is defined _recursively_:

- Sierpinski's triangle of complexity N is a figure in the form of a triangle which is made of 3 triangular figures which are themselves Sierpinski's triangles of complexity N-1.

- A Sierpinski's triangle of complexity 0 is a simple solid equilateral triangle
"

# ╔═╡ 6b8883f6-e7b3-11ea-155e-6f62117e123b
md"To draw Sierpinski's triangle, we are going to use an external package, [_Compose.jl_](https://giovineitalia.github.io/Compose.jl/latest/tutorial). Let's set up a package environment and add the package.

A package contains a coherent set of functionality that you can often use a black box according to its specification. There are [lots of Julia packages](https://juliahub.com/ui/Home).
"

# ╔═╡ 851c03a4-e7a4-11ea-1652-d59b7a6599f0
# setting up an empty package environment
begin
	import Pkg
	Pkg.activate(mktempdir())
	Pkg.Registry.update()
end

# ╔═╡ d6ee91ea-e750-11ea-1260-31ebf3ec6a9b
# add (ie install) a package to our environment
begin
	Pkg.add("Compose")
	# call `using` so that we can use it in our code
	using Compose
end

# ╔═╡ 5acd58e0-e856-11ea-2d3d-8329889fe16f
begin
	Pkg.add("PlutoUI")
	using PlutoUI
end

# ╔═╡ dbc4da6a-e7b4-11ea-3b70-6f2abfcab992
md"Just like the definition above, our `sierpinksi` function is _recursive_: it calls itself."

# ╔═╡ 02b9c9d6-e752-11ea-0f32-91b7b6481684
complexity = 3

# ╔═╡ 1eb79812-e7b5-11ea-1c10-63b24803dd8a
if complexity == 3
	md"""
Try changing the value of **`complexity` to `5`** in the cell above.

Hit `Shift+Enter` to affect the change.
	"""
else
	md"""
**Great!** As you can see, all the cells in this notebook are linked together by the variables they define and use. Just like a spreadsheet!
	"""
end

# ╔═╡ d7e8202c-e7b5-11ea-30d3-adcd6867d5f5
md"### Exercise 2.1

As you can see, the total area covered by triangles is lower when the complexity is higher."

# ╔═╡ f22222b4-e7b5-11ea-0ea0-8fa368d2a014
md"""
Can you write a function that computes the _area of `sierpinski(n)`_, as a fraction of the area of `sierpinski(0)`?

So:
```
area_sierpinski(0) = 1.0
area_sierpinski(1) = 0.??
...
```
"""

# ╔═╡ ca8d2f72-e7b6-11ea-1893-f1e6d0a20dc7
function area_sierpinski(n)
	return (3/4)^n
end

# ╔═╡ 71c78614-e7bc-11ea-0959-c7a91a10d481
if area_sierpinski(0) == 1.0 && area_sierpinski(1) == 3 / 4
	md"""
!!! correct
    Well done!
	"""
else
	md"""
!!! warning "Incorrect"
    Keep working on it!
	"""
end

# ╔═╡ c21096c0-e856-11ea-3dc5-a5b0cbf29335
md"**Let's try it out below:**"

# ╔═╡ 52533e00-e856-11ea-08a7-25e556fb1127
md"Complexity = $(@bind n Slider(0:6, show_value=true))"

# ╔═╡ c1ecad86-e7bc-11ea-1201-23ee380181a1
md"""
!!! hint
    Can you write `area_sierpinksi(n)` as a function of `area_sierpinski(n-1)`?
"""

# ╔═╡ c9bf4288-e6ce-11ea-0e13-a36b5e685998

# ╔═╡ a60a492a-e7bc-11ea-0f0b-75d81ce46a01
md"That's it for now, see you next week!"

# ╔═╡ b3c7a050-e855-11ea-3a22-3f514da746a4
if student.kerberos_id === "jazz"
	md"""
!!! danger "Oops!"
    **Before you submit**, remember to fill in your name and kerberos ID at the top of this notebook!
	"""
end

# ╔═╡ d3625d20-e6ce-11ea-394a-53208540d626

# ╔═╡ dfdeab34-e751-11ea-0f90-2fa9bbdccb1e
triangle() = compose(context(), polygon([(1, 1), (0, 1), (1 / 2, 0)]))

# ╔═╡ b923d394-e750-11ea-1971-595e09ab35b5
# It does not matter which order you define the building blocks (functions) of the
# program in. The best way to organize code is the one that promotes understanding.

function place_in_3_corners(t)
	# Uses the Compose library to place 3 copies of t
	# in the 3 corners of a triangle.
	# treat this function as a black box,
	# or learn how it works from the Compose documentation here https://giovineitalia.github.io/Compose.jl/latest/tutorial/#Compose-is-declarative-1
	compose(context(),
			(context(1 / 4,   0, 1 / 2, 1 / 2), t),
			(context(0, 1 / 2, 1 / 2, 1 / 2), t),
			(context(1 / 2, 1 / 2, 1 / 2, 1 / 2), t))
end

# ╔═╡ e2848b9a-e703-11ea-24f9-b9131434a84b
function sierpinski(n)
	if n == 0
		triangle()
	else
		t = sierpinski(n - 1) # recursively construct a smaller sierpinski's triangle
		place_in_3_corners(t) # place it in the 3 corners of a triangle
	end
end

# ╔═╡ 9664ac52-e750-11ea-171c-e7d57741a68c
sierpinski(complexity)

# ╔═╡ df0a4068-e7b2-11ea-2475-81b237d492b3
sierpinski.(0:6)

# ╔═╡ 147ed7b0-e856-11ea-0d0e-7ff0d527e352
md"""

Sierpinski's triangle of complexity $(n)

 $(sierpinski(n))

has area **$(area_sierpinski(n))**

"""

# ╔═╡ Cell order:
# ╠═fafae38e-e852-11ea-1208-732b4744e4c2
# ╟─cdff6730-e785-11ea-2546-4969521b33a7
# ╠═7308bc54-e6cd-11ea-0eab-83f7535edf25
# ╟─a2181260-e6cd-11ea-2a69-8d9d31d1ef0e
# ╟─094e39c8-e6ce-11ea-131b-07c4a1199edf
# ╟─31a8fbf8-e6ce-11ea-2c66-4b4d02b41995
# ╟─339c2d5c-e6ce-11ea-32f9-714b3628909c
# ╟─56866718-e6ce-11ea-0804-d108af4e5653
# ╠═bccf0e88-e754-11ea-3ab8-0170c2d44628
# ╟─e7abd366-e7a6-11ea-30d7-1b6194614d0a
# ╟─d62f223c-e754-11ea-2470-e72a605a9d7e
# ╠═4896bf0c-e754-11ea-19dc-1380bb356ab6
# ╠═7a01a508-e78a-11ea-11da-999d38785348
# ╟─682db9f8-e7b1-11ea-3949-6b683ca8b47b
# ╟─088cc652-e7a8-11ea-0ca7-f744f6f3afdd
# ╟─c18dce7a-e7a7-11ea-0a1a-f944d46754e5
# ╟─5e24d95c-e6ce-11ea-24be-bb19e1e14657
# ╟─6b8883f6-e7b3-11ea-155e-6f62117e123b
# ╠═851c03a4-e7a4-11ea-1652-d59b7a6599f0
# ╠═d6ee91ea-e750-11ea-1260-31ebf3ec6a9b
# ╠═5acd58e0-e856-11ea-2d3d-8329889fe16f
# ╟─dbc4da6a-e7b4-11ea-3b70-6f2abfcab992
# ╠═e2848b9a-e703-11ea-24f9-b9131434a84b
# ╠═9664ac52-e750-11ea-171c-e7d57741a68c
# ╠═02b9c9d6-e752-11ea-0f32-91b7b6481684
# ╟─1eb79812-e7b5-11ea-1c10-63b24803dd8a
# ╟─d7e8202c-e7b5-11ea-30d3-adcd6867d5f5
# ╠═df0a4068-e7b2-11ea-2475-81b237d492b3
# ╟─f22222b4-e7b5-11ea-0ea0-8fa368d2a014
# ╠═ca8d2f72-e7b6-11ea-1893-f1e6d0a20dc7
# ╟─71c78614-e7bc-11ea-0959-c7a91a10d481
# ╟─c21096c0-e856-11ea-3dc5-a5b0cbf29335
# ╟─52533e00-e856-11ea-08a7-25e556fb1127
# ╟─147ed7b0-e856-11ea-0d0e-7ff0d527e352
# ╟─c1ecad86-e7bc-11ea-1201-23ee380181a1
# ╟─c9bf4288-e6ce-11ea-0e13-a36b5e685998
# ╟─a60a492a-e7bc-11ea-0f0b-75d81ce46a01
# ╟─b3c7a050-e855-11ea-3a22-3f514da746a4
# ╟─d3625d20-e6ce-11ea-394a-53208540d626
# ╠═dfdeab34-e751-11ea-0f90-2fa9bbdccb1e
# ╠═b923d394-e750-11ea-1971-595e09ab35b5


text:
function producer(c::Channel)
    put!(c, "start")
    for i in 1:4
        put!(c, 2i)
    end
    put!(c, "stop")
end

function example1()
    ch = Channel(producer)
    println(take!(ch))

    println(take!(ch))
    println(take!(ch))
    println(take!(ch))
    println(take!(ch))

    println(take!(ch))
end

function example2()
    for x in Channel(producer)
        println(x)
    end
end

function main()
    example1()
    example2()
end

main()

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
# This is a Julia translation of the Matlab code PinfishSST.m
# Donglai Gong 2019-06-07
#
using Statistics
using NCDatasets, Missings, Dates
using FileIO
using Plotly

initflag = 1
loadflag = 1
analysisflag = 1
saveflag = 1

yrrangeall = [yyyy for yyyy=2003:2020];
yrrangehipinfish = [2012; 2013; 2015];
yrrangelopinfish = [2014; 2016; 2017; 2018; 2019];

yrrange1 = [2020];

ydayNlatest = 157
yearrange = yrrangeall; # change THIS to select different years!!!

if initflag == 1 #| !exist('latMAB')
    datapath = "https://podaac-opendap.jpl.nasa.gov:443/opendap/allData/ghrsst/data/GDS2/L4/GLOB/JPL/MUR/v4.1/2018/148/20180528090000-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.nc"
    ds = Dataset(datapath);
    lon = ds["lon"];
    lat = ds["lat"];
    lonarr = Float64.(collect(Missings.replace(lon[:],NaN)));
    latarr = Float64.(collect(Missings.replace(lat[:],NaN)));
    latind = findall( 44 .>= latarr .>= 33);
    lonind = findall(-80 .<= lonarr .<= -69);
    latMAB = latarr[latind];
    lonMAB = lonarr[lonind];
end #if

#startind = [lonind[1],latind[1],1];
#countind = [size(lonind,1),size(latind,1),1];
#lonind2 = [ind for ind=startind[1]:(startind[1]+countind[1]-1)];
#latind2 = [ind for ind=startind[2]:(startind[2]+countind[2]-1)];
#tind = [ind for ind=startind[3]:(startind[3]+countind[3]-1)];

# load all the SST data and compute composite for MAB
if loadflag == 1
    sstMAB = repeat([0],size(latMAB,1),size(lonMAB,1));
    sstMAB1yr = sstMAB;
    sstMAByr = repeat([NaN],size(latMAB,1),size(lonMAB,1),size(yearrange,1));

    t0 = Date(2002,1,1)+Day(152-1)
    tN = Date(2020,1,1)+Day(157-1) #dayofyear(Date(2019,5,5))

    yday0 = 1;
    global navg = 0;

    filenametail = "090000-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.nc";

    yyyy1 = yearrange[1];
    yyyy2 = yearrange[end];

    #df = DateFormat("yyyymmddTHMS")

    for ii = 1:size(yearrange,1)
        yyyy = yearrange[ii];
        println(yyyy)
        navg1 = 0;
        sstMAB1yr = repeat([0],size(latMAB,1),size(lonMAB,1));

        if Dates.isleapyear(Date(yyyy))
            yday0 = 1;
            ydayN = 366;
        else
            yday0 = 1;
            ydayN = 365;
        end #if

        yday1 = 120;
        yday2 = 160;

        if yyyy == 2002
            yday0 = 152;
            yday1 = yday0;
        elseif yyyy == 2020
            ydayN = ydayNlatest;
            yday2 = ydayN;
        else

        end

        yday1 = 90
        yday2 = ydayNlatest

        syday1000 = string(yday1+1000);
        syday2000 = string(yday2+1000);
        ydays = [yday for yday in yday1:yday2];
        for jj = 1:length(ydays)
            global navg
            global sstMAB
            yday = ydays[jj];
            println(yday)
            t = DateTime(yyyy,1,1)+Day(yday-1);
            st = string(t);
            filenamehead = string(st[1:4], st[6:7], st[9:10]);
            filename = string(filenamehead, filenametail);
            yday000 = string(yday+1000);
            datapath = string("https://podaac-opendap.jpl.nasa.gov:443/opendap/allData/ghrsst/data/GDS2/L4/GLOB/JPL/MUR/v4.1/", yyyy, "/", yday000[2:end], "/", filename);

            navg = navg + 1;
            navg1 = navg1 + 1;

            sstds = Dataset(datapath)["analysed_sst"];
            sstMABi = sstds[lonind[1]:lonind[end],latind[1]:latind[end],1];
            sstMABimx = collect(Missings.replace(sstMABi,NaN));

            sstMAB1yr = sstMAB1yr*(1.0-1.0/navg1) + sstMABimx*(1.0/navg1);
            sstMAB = sstMAB*(1.0-1.0/navg) + sstMABimx*(1.0/navg);
        end #for
        sstMAByr[:,:,ii] = sstMAB1yr;
    end #for
end #if

#latind2 = [ind for ind=startind[2]:(startind[2]+countind[2]-1)];
if analysisflag == 1
    LON = repeat(lonMAB',length(latMAB));
    LAT = repeat(reverse(latMAB,dims = 1),1,length(lonMAB));
    llind = findall((-73 .>= LON .>= -74.5) .& (39 .>= LAT .>= 37.5));
    nanmean(x) = mean(filter(!isnan,x)) # implementing nanmean similar to Matlab
    nanmean(x,y) = mapslices(nanmean,x,dims = y) # implementing nanmean by dim of array
    #sstMABts = nanmean(sstMAByr[llind],)
    sstMABts = [];
    tyr = yearrange;
end #if

if saveflag == 1
    workdir = "/Users/gong/Documents/Research/Projects/Pinfish/Pinfish/";
    filepath = string(workdir, "sstMAB_", yyyy1, "-", yyyy2, "_yd", syday1000[2:end], "-", syday2000[2:end], ".jld2")
    @save filepath latMAB lonMAB sstMAB sstMAByr navg tyr sstMABts
end #if

if plotflag == 1
    function sstmap()
        sst = heatmap(
            x=lonMAB,
            y=latMAB,
            z=sstMAB1yr
        )
        plot(sst)
    end
    sstmap()
end #if


text:
module Events
using ..StructTypes, ..Messages
import ..Optional, ..FriendId, ..GroupOrFriendId, ..GroupId, ..TimeStamp, ..DurationSeconds
import ..MessageId, ..EventId, ..DeviceKind, ..MessageChains.MessageChain, ..EventOrMessage
import ..Messages.SelfOrMember

abstract type AbstractEvent <: EventOrMessage end

abstract type BotLoginEvent <: AbstractEvent end

struct BotOnlineEvent <: BotLoginEvent
    qq::FriendId
end

abstract type BotLogoutEvent <: AbstractEvent end

struct BotOfflineEventActive <: BotLogoutEvent
    qq::FriendId
end

struct BotOfflineEventForce <: BotLogoutEvent
    qq::FriendId
end

struct BotOfflineEventDropped <: BotLogoutEvent
    qq::FriendId
end

struct BotReloginEvent <: BotLoginEvent
    qq::FriendId
end

abstract type FriendStatusEvent <: AbstractEvent end

struct FriendInputStatusChangedEvent <: FriendStatusEvent
    friend::Friend
    inputting::Bool
end

struct FriendNickChangedEvent <: FriendStatusEvent
    friend::Friend
    from::String
    to::String
end

abstract type RecallEvent <: AbstractEvent end

struct GroupRecallEvent <: RecallEvent
    authorId::FriendId
    messageId::MessageId
    time::TimeStamp
    group::Group
    operator::SelfOrMember
end

struct FriendRecallEvent <: RecallEvent
    authorId::FriendId
    messageId::MessageId
    time::TimeStamp
    operator::FriendId
end

module NudgeKinds
@enum NudgeKind begin
    GROUP
    FRIEND
    STRANGER
end
end
import .NudgeKinds.NudgeKind

function StructTypes.construct(::Type{NudgeKind}, s::String; kw...)
    s == "Group" && return NudgeKinds.GROUP
    s == "Friend" && return NudgeKinds.FRIEND
    s == "Stranger" && return NudgeKinds.STRANGER
    error("Unknown subject kind")
end

Base.string(x::NudgeKind) = x == NudgeKinds.GROUP ? "Group" : "Friend"

StructTypes.@Struct struct Subject
    id::GroupOrFriendId
    kind::NudgeKind
end

struct NudgeEvent <: AbstractEvent
    fromId::FriendId
    subject::Subject
    action::String
    suffix::String
    target::FriendId
end

abstract type GroupEvent <: AbstractEvent end

abstract type GroupStatusChangeEvent <: GroupEvent end

struct GroupNameChangeEvent <: GroupStatusChangeEvent
    origin::String  # original name
    current::String # current name
    group::Group
    operator::SelfOrMember
end

struct GroupEntranceAnnouncementChangeEvent <: GroupStatusChangeEvent
    origin::String
    current::String
    group::Group
    operator::SelfOrMember
end

struct GroupMuteAllEvent <: GroupStatusChangeEvent
    origin::Bool
    current::Bool
    group::Group
    operator::SelfOrMember
end

struct GroupAllowAnonymousChatEvent <: GroupStatusChangeEvent
    origin::Bool
    current::Bool
    group::Group
    operator::SelfOrMember
end

struct GroupAllowConfessTalkEvent <: GroupStatusChangeEvent
    origin::Bool
    current::Bool
    group::Group
    isByBot::Bool
end

struct GroupAllowMemberInviteEvent <: GroupStatusChangeEvent
    origin::Bool
    current::Bool
    group::Group
    operator::SelfOrMember
end

abstract type MemberEvent <: GroupEvent end

struct MemberJoinEvent <: MemberEvent
    member::Member
    invitor::SelfOrMember
end

struct BotJoinGroupEvent <: MemberEvent
    group::Group
    invitor::Union{Nothing,Member}
end

abstract type MemberLeaveEvent <: MemberEvent end

abstract type OtherMemberLeaveEvent <: MemberLeaveEvent end

struct MemberLeaveEventKick <: OtherMemberLeaveEvent
    member::Member
    operator::SelfOrMember
end

# If bot is leaving, see BotLeaveEvent
struct MemberLeaveEventQuit <: OtherMemberLeaveEvent
    member::Member
end

abstract type BotLeaveEvent <: MemberLeaveEvent end

struct BotLeaveEventActive <: BotLeaveEvent
    group::Group
end

struct BotLeaveEventKick <: BotLeaveEvent
    group::Group
    operator::Member
end

abstract type MemberStatusChangeEvent <: MemberEvent end

struct MemberCardChangeEvent <: MemberStatusChangeEvent
    origin::String
    current::String
    group::Group
    member::Member
end

struct MemberSpecialTitleChangeEvent <: MemberStatusChangeEvent
    origin::String
    current::String
    group::Group
    member::Member
end

struct MemberPermissionChangeEvent <: MemberStatusChangeEvent
    origin::Permission
    current::Permission
    group::Group
    member::Member
end

struct BotGroupPermissionChangeEvent <: MemberStatusChangeEvent
    origin::Permission
    current::Permission
    group::Group
end

struct MemberMuteEvent <: MemberStatusChangeEvent
    durationSeconds::DurationSeconds
    member::Member
    operator::SelfOrMember
end

struct MemberUnmuteEvent <: MemberStatusChangeEvent
    member::Member
    operator::SelfOrMember
end

struct BotMuteEvent <: MemberStatusChangeEvent
    durationSeconds::DurationSeconds
    operator::Member
end

struct BotUnmuteEvent <: MemberStatusChangeEvent
    oprerator::Member
end

module Actions
@enum Action achieve lose
end
import .Actions.Action

struct MemberHonorChangeEvent <: MemberStatusChangeEvent
    member::Member
    action::Action
    honor::String
end

abstract type RequestEvent <: AbstractEvent end

struct NewFriendRequestEvent <: RequestEvent
    eventId::EventId
    fromId::FriendId
    groupId::GroupId
    nick::String
    message::String
end

struct MemberJoinRequestEvent <: RequestEvent
    eventId::EventId
    fromId::FriendId # applyer
    groupId::GroupId
    groupName::String
    nick::String
    message::String
end

struct BotInvitedJoinGroupRequestEvent <: RequestEvent
    eventId::EventId
    fromId::FriendId # invitor
    groupId::GroupId
    groupName::String
    nick::String
    message::String
end

abstract type OtherClientEvent <: AbstractEvent end

struct OtherClientOnlineEvent <: OtherClientEvent
    client::Platform
    kind::Optional{DeviceKind}
end

struct OtherClientOfflineEvent <: OtherClientEvent
    client::Platform
end

struct CommandExecutedEvent <: AbstractEvent
    eventId::Optional{Int}   # Inconsistent documentation here
    name::String
    friend::Optional{Friend} # Nothing if sent from console
    member::Optional{Member} # Nothing if sent from console
    args::MessageChain
end

"""
    ProtocolAdapterConnected

This event is broadcasted by the `Broadcaster` when the protocol adapter is fully launched.
"""
struct ProtocolAdapterConnected <: AbstractEvent end

const event_types = (; BotOnlineEvent, BotReloginEvent, BotOfflineEventActive, BotOfflineEventDropped,
    BotOfflineEventForce, CommandExecutedEvent, FriendInputStatusChangedEvent, FriendNickChangedEvent,
    BotMuteEvent, BotUnmuteEvent, GroupAllowAnonymousChatEvent, GroupAllowConfessTalkEvent,
    GroupAllowMemberInviteEvent, GroupEntranceAnnouncementChangeEvent, GroupMuteAllEvent,
    GroupNameChangeEvent, BotJoinGroupEvent, MemberJoinEvent, BotLeaveEventActive, BotLeaveEventKick,
    MemberLeaveEventKick, MemberLeaveEventQuit, BotGroupPermissionChangeEvent, MemberCardChangeEvent,
    MemberHonorChangeEvent, MemberMuteEvent, MemberPermissionChangeEvent, MemberSpecialTitleChangeEvent,
    MemberUnmuteEvent, NudgeEvent, OtherClientOfflineEvent, OtherClientOnlineEvent, FriendRecallEvent,
    GroupRecallEvent, BotInvitedJoinGroupRequestEvent, MemberJoinRequestEvent, NewFriendRequestEvent,
    FriendMessage, GroupMessage, OtherClientMessage, StrangerMessage, TempMessage)

for n in keys(event_types)
    @eval export $n
end

export NudgeKinds, NudgeKind, Actions, Action

end


text:
__precompile__()

module Types

export IntP, FloatP
export DMDStruct

const IntP = Int64
const FloatP = Float64

mutable struct DMDStruct
    facecon::Array{IntP,3}; # faces-to-DG nodes connectivities (used to obtain DG nodes on faces)
    eblks::Array{IntP,2};  # blocks of elements for parallel computation
    fblks::Array{IntP,2};  # blocks of faces for parallel computation
    nbsd::Array{IntP,2};   # neighboring subdomains (neighbors)
    elemsend::Array{IntP,2};  # elements to be sent to neighbors
    elemrecv::Array{IntP,2};  # elements to be received from neighbors
    elemsendpts::Array{IntP,2}; # markers for elements to be sent to neighbors
    elemrecvpts::Array{IntP,2}; # markers for elements to be received from neighbors
    elempart::Array{IntP,2};    # element partitions
    elempartpts::Array{IntP,2}; # markers for element partitions
    facepart::Array{IntP,2};    # element partitions
    facepartpts::Array{IntP,2}; # markers for element partitions
    facepartbnd::Array{IntP,2};
    t2f::Array{IntP,2};
    f::Array{IntP,2};
    #elem2cpu::Array{IntP,2};    # element partitions
    DMDStruct() = new();
end

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
using ConstrainedDynamics

# Parameters
joint_axis = [.0;1.0;1.0]
v1 = -[0.5;0;0.]
v2 = -[0.;0.5;0.]
v3 = -[0.;0;0.5]

length1 = 1.
width, depth = 0.1, 0.1
box1 = Box(length1, width, depth, length1, color = RGBA(1., 1., 0.))
box2 = Box(width, length1, depth, length1, color = RGBA(1., 0., 0.))
box3 = Box(width, depth, length1, length1, color = RGBA(1., 0., 1.))

# Links
origin = Origin{Float64}()
link1 = Body(box1)
link2 = Body(box2)
link3 = Body(box3)

# Constraints
# joint1 = EqualityConstraint(Spherical(origin, link1; p2=v1))
# joint2 = EqualityConstraint(Spherical(origin, link2; p2=v2))
# joint3 = EqualityConstraint(Spherical(origin, link3; p2=v3))
joint1 = EqualityConstraint(Revolute(origin, link1, joint_axis; p2=v1))
joint2 = EqualityConstraint(Revolute(origin, link2, joint_axis; p2=v2))
joint3 = EqualityConstraint(Revolute(origin, link3, joint_axis; p2=v3))

links = [link1;link2;link3]
constraints = [joint1;joint2;joint3]
shapes = [box1;box2;box3]

mech = Mechanism(origin, links, constraints, g = 0., shapes = shapes)
setPosition!(origin,link1,Δq = UnitQuaternion(RotZ(0.)),p2=v1)
setPosition!(origin,link2,Δq = UnitQuaternion(RotX(0.)),p2=v2)
setPosition!(origin,link3,Δq = UnitQuaternion(RotX(0.)),p2=v3)
# setVelocity!(origin,link1,Δω = [0.;1;1],p2=v1)
# setVelocity!(origin,link2,Δω = [0.;1;1],p2=v2)
# setVelocity!(origin,link3,Δω = [0.;1;1],p2=v3)
# setVelocity!(link2,v = [-2.5;2.5;5],ω = [2;4;10.])
# setForce!(link1,τ=[0.1;0;0])


text:
# A slight re-write of simpson_vec from github.com/rodluger/limbdark.jl

struct IntegralArrays{T<:Real}
    g::Matrix{T}
    A::Matrix{T}
    S::Array{T,3}
    f_of_x::Vector{T}
    I_of_f::Vector{T}

    nf::Int64
    maxdepth::Int64
    tol::T

    function IntegralArrays(nf::Int64, maxdepth::Int64, tol::T) where T<:Real
        # Add one if derivatives are computed
        # length(dfdg) == length(dfdu) + 1
        if nf > 1; nf += 1; end
        g = zeros(T, nf, 5)
        A = zeros(T, nf, 3)
        S = zeros(T, nf, 3, maxdepth)
        f_of_x = zeros(T, nf)
        I_of_f = zeros(T, nf)
        return new{T}(g, A, S, f_of_x, I_of_f, nf, maxdepth, tol)
    end
end

function integrate_simpson!(a::T, b::T, f::Function, ia::IntegralArrays{T}) where T<:Real
    third = one(T) / 3
    i = zero(T)
    m = 0; n = 0

    # 'Unpack' IntegralArrays struct
    g = ia.g
    A = ia.A
    S = ia.S
    f_of_x = ia.f_of_x
    I_of_f = ia.I_of_f
    epsilon = ia.tol
    N = ia.maxdepth
    nf = ia.nf

    # Zero out integral array
    I_of_f .= zero(T)
    f_of_x .= zero(T)

    f(a,f_of_x)
    @inbounds for j = 1:nf
        g[j,1] = f_of_x[j]
    end
    f(0.5 * (a + b),f_of_x)
    @inbounds for j = 1:nf
        g[j,3] = f_of_x[j]
    end
    f(b,f_of_x)
    @inbounds for j = 1:nf
        g[j,5] = f_of_x[j]
    end
    bma = b - a
    @inbounds for j = 1:nf
        A[j,1] = 0.5 * bma * (g[j,1] + 4 * g[j,3] + g[j,5])
    end
    @label AA
    d = 2^n
    h = 0.25 * bma / d

    f(a + h * (4 * m + 1),f_of_x)
    @inbounds for j = 1:nf
        g[j,2] = f_of_x[j]
    end
    f(a + h * (4 * m + 3),f_of_x)
    @inbounds for j = 1:nf
        g[j,4] = f_of_x[j]
    end
    @inbounds for j = 1:nf
        A[j,2] = h * (g[j,1] + 4 * g[j,2] + g[j,3])
        A[j,3] = h * (g[j,3] + 4 * g[j,4] + g[j,5])
    end
    maxdiff = zero(T)
    @inbounds for j = 1:nf
        maxa231 = abs((A[j,2] + A[j,3]) - A[j,1])
        if maxa231 > maxdiff
            maxdiff = maxa231
        end
    end
    if maxdiff > 3 * epsilon
        m *= 2
        n += 1
        if n > N
            @goto CC
        end
        @inbounds for j = 1:nf
            A[j,1] = A[j,2]
            S[j,1,n] = A[j,3]
            S[j,2,n] = g[j,4]
            S[j,3,n] = g[j,5]
            g[j,5] = g[j,3]
            g[j,3] = g[j,2]
        end
        @goto AA
    else
        @inbounds for j = 1:nf
            I_of_f[j] += (A[j,2] + A[j,3]) * third
        end
        m += 1
        i = a + m * bma / d
        @label BB
        if m == 2 * div(m, 2)
            m = div(m, 2)
            n -= 1
            @goto BB
        end
        if (m != 1) || (n != 0)
            @inbounds for j = 1:nf
                A[j,1] = S[j,1,n]
                g[j,1] = g[j,5]
                g[j,3] = S[j,2,n]
                g[j,5] = S[j,3,n]
            end
            @goto AA
        end
    end
    @label CC
    return
end


text:
##
# Base hotfixes
##
module BaseHacks
    @info "activating BaseHacks..."

    import Base.IteratorsMD: inc
    import Base: tail

    function unsafe_inc(i::Int64)
        Base.llvmcall("""
        %i = add nsw i64 %0, 1
        ret i64 %i
        """, Int64, Tuple{Int64}, i)
    end

    @inline inc(state::Tuple{Int}, start::Tuple{Int}, stop::Tuple{Int}) = (unsafe_inc(state[1]),)
    @inline function inc(state, start, stop)
         if state[1] < stop[1]
           return (unsafe_inc(state[1]),tail(state)...)
         end
         newtail = inc(tail(state), tail(start), tail(stop))
         (start[1], newtail...)
    end

    @inline function Base.getindex(r::AbstractUnitRange, s::AbstractUnitRange{<:Integer})
        @boundscheck checkbounds(r, s)
        f = first(r)
        st = oftype(f, f + first(s)-1)
        Base._range(st, nothing, nothing, length(s)) #gotta use this version of range instead of const-propping through the kwarg version.
    end

    #FIXME This is a separate eachindex for scalar cases so that swizzles can use
    #a different thunk for the scalar case and avoid recursion depth limiting.
    Base.eachindex(args::AbstractArray{<:Any, 0}...) = CartesianIndices(())

    #=
    #Avoid checking shapes in eachindex when inbounds scopes
    Base.@propagate_inbounds function Base.eachindex(A::AbstractArray, B::AbstractArray)
        eachindex(IndexStyle(A,B), A, B)
    end
    Base.@propagate_inbounds function Base.eachindex(A::AbstractArray, B::AbstractArray...)
        eachindex(IndexStyle(A,B...), A, B...)
    end
    Base.@propagate_inbounds function Base.eachindex(::IndexLinear, A::AbstractArray, B::AbstractArray...)
        indsA = eachindex(IndexLinear(), A)
        @boundscheck Base._all_match_first(X->eachindex(IndexLinear(), X), indsA, B...) || Base.throw_eachindex_mismatch(IndexLinear(), A, B...)
        indsA
    end
    Base.@propagate_inbounds function Base.eachindex(::IndexCartesian, A::AbstractArray, B::AbstractArray...)
        axsA = axes(A)
        @boundscheck Base._all_match_first(axes, axsA, B...) || Base.throw_eachindex_mismatch(IndexCartesian(), A, B...)
        CartesianIndices(axsA)
    end
    =#

    #=
    Base.@propagate_inbounds function Base.getindex(iter::CartesianIndices{N,<:NTuple{N,Base.OneTo}}, I::Vararg{Int, N}) where {N}
        @boundscheck checkbounds(iter, I...)
        CartesianIndex(I)
    end

    Base.@propagate_inbounds function Base.getindex(iter::CartesianIndices{N,R}, I::Vararg{Int, N}) where {N,R}
        @boundscheck checkbounds(iter, I...)
        CartesianIndex(I .- first.(Base.axes1.(iter.indices)) .+ first.(iter.indices))
    end
    =#

    @info "BaseHacks is online. all your Base are belong to us."
end


text:
abstract type Node{T<:Real} end

abstract type Layer{T} <: Node{T} end
abstract type Activation{T} <: Node{T} end

const NameOrIndex = Union{Symbol,Int}

Base.eltype{T}(::Node{T}) = T

"""
    name(n::Node)
Return the Symbol name of this Node.
"""
@required_func name(n::Node)

"""
    output(n::Node)
Get the tensor output of the Node
"""
@required_func output(n::Node)

"""
    forward!(n::Node)
Compute the node output
"""
@required_func forward!(n::Node)


text:
module NaturalSort
export natural

function natural(x::AbstractString, y::AbstractString)
    iterx = iterate(x)
    itery = iterate(y)
    while iterx !== nothing && itery !== nothing
        cx, statex = iterx
        cy, statey = itery
        if isnumeric(cx) && isnumeric(cy)
            # Skip leading zeros
            while cx == '0' && (iterx = iterate(x, statex)) !== nothing
                cx, statex = iterx
            end
            while cy == '0' && (itery = iterate(y, statey)) !== nothing
                cy, statey = itery
            end
            if iterx === nothing || itery === nothing
                break
            end

            # Begin comparing numbers
            diff = false
            lt = false
            while true
                isnumx = isnumeric(cx)
                isnumy = isnumeric(cy)
                if isnumx && isnumy
                    if !diff && cx != cy
                        # Keep track of how numbers differ, in case the lengths match
                        diff = true
                        lt = cx < cy
                    end

                    iterx = iterate(x, statex)
                    itery = iterate(y, statey)
                    if iterx === nothing || itery === nothing
                        if iterx === nothing && itery !== nothing && isnumeric(itery[1])
                            # Number in y is longer than number in x
                            return true
                        elseif iterx !== nothing && itery === nothing && isnumeric(iterx[1])
                            # Number in x is longer than number in y
                            return false
                        end
                        # Both numbers ended and same length
                        return diff ? lt : iterx === nothing && itery !== nothing
                    end

                    cx, statex = iterx
                    cy, statey = itery
                elseif isnumx
                    # Number in x is longer than number in y
                    return false
                elseif isnumy
                    # Number in y is longer than number in x
                    return true
                elseif diff
                    # Numbers were same length but different
                    return lt
                else
                    # Numbers were the same
                    break
                end
            end
        end
        if cx != cy
            return cx < cy
        end
        iterx = iterate(x, statex)
        itery = iterate(y, statey)
    end

    return iterx == nothing && itery != nothing
end
end # module


text:
#########################
# Main types for meshes #
#########################
"""
A `Node` is a point in space.
"""
struct Node{dim,T}
    x::Vec{dim,T}
end
Node(x::NTuple{dim,T}) where {dim,T} = Node(Vec{dim,T}(x))
getcoordinates(n::Node) = n.x

"""
A `Cell` is a sub-domain defined by a collection of `Node`s as it's vertices.
"""
struct Cell{dim,N,M}
    nodes::NTuple{N,Int}
end
nfaces(c::Cell) = nfaces(typeof(c))
nfaces(::Type{Cell{dim,N,M}}) where {dim,N,M} = M

# Typealias for commonly used cells
const Line = Cell{1,2,2}
const QuadraticLine = Cell{1,3,2}

const Triangle = Cell{2,3,3}
const QuadraticTriangle = Cell{2,6,3}

const Quadrilateral = Cell{2,4,4}
const QuadraticQuadrilateral = Cell{2,9,4}

const Tetrahedron = Cell{3,4,4}
const QuadraticTetrahedron = Cell{3,10,4}

const Hexahedron = Cell{3,8,6}
const QuadraticHexahedron = Cell{3,20,6} # Function interpolation for this doesn't exist in JuAFEM yet

"""
A `CellIndex` wraps an Int and corresponds to a cell with that number in the mesh
"""
struct CellIndex
    idx::Int
end

"""
A `FaceIndex` wraps an (Int, Int) and defines a face by pointing to a (cell, face).
"""
struct FaceIndex
    idx::Tuple{Int,Int} # cell and side
end

"""
A `Grid` is a collection of `Cells` and `Node`s which covers the computational domain, together with Sets of cells, nodes and faces.
"""
mutable struct Grid{dim,N,T<:Real,M}
    cells::Vector{Cell{dim,N,M}}
    nodes::Vector{Node{dim,T}}
    # Sets
    cellsets::Dict{String,Set{Int}}
    nodesets::Dict{String,Set{Int}}
    facesets::Dict{String,Set{Tuple{Int,Int}}} # TODO: This could be Set{FaceIndex} which could result in nicer use later
    # Boundary matrix (faces per cell × cell)
    boundary_matrix::SparseMatrixCSC{Bool,Int}
end

function Grid(cells::Vector{Cell{dim,N,M}},
              nodes::Vector{Node{dim,T}};
              cellsets::Dict{String,Set{Int}}=Dict{String,Set{Int}}(),
              nodesets::Dict{String,Set{Int}}=Dict{String,Set{Int}}(),
              facesets::Dict{String,Set{Tuple{Int,Int}}}=Dict{String,Set{Tuple{Int,Int}}}(),
              boundary_matrix::SparseMatrixCSC{Bool,Int}=spzeros(Bool, 0, 0)) where {dim,N,M,T}
    return Grid(cells, nodes, cellsets, nodesets, facesets, boundary_matrix)
end

##########################
# Grid utility functions #
##########################
@inline getcells(grid::Grid) = grid.cells
@inline getcells(grid::Grid, v::Union{Int, Vector{Int}}) = grid.cells[v]
@inline getcells(grid::Grid, set::String) = grid.cells[grid.cellsets[set]]
@inline getncells(grid::Grid) = length(grid.cells)
@inline getcelltype(grid::Grid) = eltype(grid.cells)

@inline getnodes(grid::Grid) = grid.nodes
@inline getnodes(grid::Grid, v::Union{Int, Vector{Int}}) = grid.nodes[v]
@inline getnodes(grid::Grid, set::String) = grid.nodes[grid.nodesets[set]]
@inline getnnodes(grid::Grid) = length(grid.nodes)

@inline getcellset(grid::Grid, set::String) = grid.cellsets[set]
@inline getcellsets(grid::Grid) = grid.cellsets

@inline getnodeset(grid::Grid, set::String) = grid.nodesets[set]
@inline getnodesets(grid::Grid) = grid.nodesets

@inline getfaceset(grid::Grid, set::String) = grid.facesets[set]
@inline getfacesets(grid::Grid) = grid.facesets

n_faces_per_cell(grid::Grid) = nfaces(eltype(grid.cells))

# Transformations

function transform!(g::Grid, f::Function)
    c = similar(g.nodes)
    for i in 1:length(c)
        c[i] = Node(f(g.nodes[i].x))
    end
    copy!(g.nodes, c)
    g
end

# Sets

_check_setname(dict, name) = haskey(dict, name) && throw(ArgumentError("there already exists a set with the name: $name"))
_warn_emptyset(set) = length(set) == 0 && warn("no entities added to set")

function addcellset!(grid::Grid, name::String, cellid::Union{Set{Int},Vector{Int}})
    _check_setname(grid.cellsets,  name)
    cells = Set(cellid)
    _warn_emptyset(cells)
    grid.cellsets[name] = cells
    grid
end

function addcellset!(grid::Grid, name::String, f::Function; all::Bool=true)
    _check_setname(grid.cellsets,  name)
    cells = Set{Int}()
    for (i, cell) in enumerate(getcells(grid))
        pass = all
        for node_idx in cell.nodes
            node = grid.nodes[node_idx]
            v = f(node.x)
            all ? (!v && (pass = false; break)) : (v && break)
        end
        pass && push!(cells, i)
    end
    _warn_emptyset(cells)
    grid.cellsets[name] = cells
    grid
end

function addfaceset!(grid::Grid, name::String, faceid::Set{Tuple{Int,Int}})
    _check_setname(grid.facesets, name)
    faceset = Set(faceid)
    _warn_emptyset(faceset)
    grid.facesets[name] = faceset
    grid
end

function addfaceset!(grid::Grid, name::String, f::Function; all::Bool=true)
    _check_setname(grid.facesets, name)
    faceset = Set{Tuple{Int,Int}}()
    for (cell_idx, cell) in enumerate(getcells(grid))
        for (face_idx, face) in enumerate(faces(cell))
            pass = all
            for node_idx in face
                v = f(grid.nodes[node_idx].x)
                all ? (!v && (pass = false; break)) : (v && break)
            end
            pass && push!(faceset, (cell_idx, face_idx))
        end
    end
    _warn_emptyset(faceset)
    grid.facesets[name] = faceset
    grid
end

function addnodeset!(grid::Grid, name::String, nodeid::Union{Vector{Int},Set{Int}})
    _check_setname(grid.nodesets, name)
    grid.nodesets[name] = Set(nodeid)
    _warn_emptyset(grid.nodesets[name])
    grid
end

function addnodeset!(grid::Grid, name::String, f::Function)
    _check_setname(grid.nodesets, name)
    nodes = Set{Int}()
    for (i, n) in enumerate(getnodes(grid))
        f(n.x) && push!(nodes, i)
    end
    grid.nodesets[name] = nodes
    _warn_emptyset(grid.nodesets[name])
    grid
end

"""
    getcoordinates!(x::Vector, grid::Grid, cell::Int)

Update the coordinate vector `x` for cell number `cell`.
"""
@inline function getcoordinates!(x::Vector{Vec{dim,T}}, grid::Grid{dim,N,T}, cell::Int) where {dim,T,N}
    @assert length(x) == N
    @inbounds for i in 1:N
        x[i] = grid.nodes[grid.cells[cell].nodes[i]].x
    end
end
@inline getcoordinates!(x::Vector{Vec{dim,T}}, grid::Grid{dim,N,T}, cell::CellIndex) where {dim, T, N} = getcoordinates!(x, grid, cell.idx)
@inline getcoordinates!(x::Vector{Vec{dim,T}}, grid::Grid{dim,N,T}, face::FaceIndex) where {dim, T, N} = getcoordinates!(x, grid, face.idx[1])

"""
    getcoordinates(grid::Grid, cell)

Return a vector with the coordinates of the vertices of cell number `cell`.
"""
@inline function getcoordinates(grid::Grid{dim,N,T}, cell::Int) where {dim,N,T}
    nodeidx = grid.cells[cell].nodes
    return [grid.nodes[i].x for i in nodeidx]::Vector{Vec{dim,T}}
end
@inline getcoordinates(grid::Grid, cell::CellIndex) = getcoordinates(grid, cell.idx)
@inline getcoordinates(grid::Grid, face::FaceIndex) = getcoordinates(grid, face.idx[1])

# Iterate over cell vector
Base.start(c::Vector{Cell{dim,N}}) where {dim,N} = 1
Base.next(c::Vector{Cell{dim,N}}, state) where {dim,N} = (CellIndex(state), state + 1)
Base.done(c::Vector{Cell{dim,N}}, state) where {dim,N} = state > length(c)

function Base.show(io::IO, grid::Grid)
    print(io, "$(typeof(grid)) with $(getncells(grid)) $(celltypes[eltype(grid.cells)]) cells and $(getnnodes(grid)) nodes")
end

const celltypes = Dict{DataType, String}(Cell{1,2,2}  => "Line",
                                         Cell{1,3,2}  => "QuadraticLine",
                                         Cell{2,3,3}  => "Triangle",
                                         Cell{2,6,3}  => "QuadraticTriangle",
                                         Cell{2,4,4}  => "Quadrilateral",
                                         Cell{2,9,4}  => "QuadraticQuadrilateral",
                                         Cell{3,4,4}  => "Tetrahedron",
                                         Cell{3,10,4} => "QuadraticTetrahedron",
                                         Cell{3,8,6}  => "Hexahedron",
                                         Cell{3,20,6} => "QuadraticHexahedron")

# Functions to uniquely identify vertices, edges and faces, used when distributing
# dofs over a mesh. For this we can ignore the nodes on edged, faces and inside cells,
# we only need to use the nodes that are vertices.
# 1D: vertices
faces(c::Union{Line,QuadraticLine}) = (c.nodes[1], c.nodes[2])
vertices(c::Union{Line,QuadraticLine}) = (c.nodes[1], c.nodes[2])
# 2D: vertices, faces
vertices(c::Union{Triangle,QuadraticTriangle}) = (c.nodes[1], c.nodes[2], c.nodes[3])
faces(c::Union{Triangle,QuadraticTriangle}) = ((c.nodes[1],c.nodes[2]), (c.nodes[2],c.nodes[3]), (c.nodes[3],c.nodes[1]))
vertices(c::Union{Quadrilateral,QuadraticQuadrilateral}) = (c.nodes[1], c.nodes[2], c.nodes[3], c.nodes[4])
faces(c::Union{Quadrilateral,QuadraticQuadrilateral}) = ((c.nodes[1],c.nodes[2]), (c.nodes[2],c.nodes[3]), (c.nodes[3],c.nodes[4]), (c.nodes[4],c.nodes[1]))
# 3D: vertices, edges, faces
vertices(c::Union{Tetrahedron,QuadraticTetrahedron}) = (c.nodes[1], c.nodes[2], c.nodes[3], c.nodes[4])
edges(c::Union{Tetrahedron,QuadraticTetrahedron}) = ((c.nodes[1],c.nodes[2]), (c.nodes[2],c.nodes[3]), (c.nodes[3],c.nodes[1]), (c.nodes[1],c.nodes[4]), (c.nodes[2],c.nodes[4]), (c.nodes[3],c.nodes[4]))
faces(c::Union{Tetrahedron,QuadraticTetrahedron}) = ((c.nodes[1],c.nodes[2],c.nodes[3]), (c.nodes[1],c.nodes[2],c.nodes[4]), (c.nodes[2],c.nodes[3],c.nodes[4]), (c.nodes[1],c.nodes[4],c.nodes[3]))
vertices(c::Union{Hexahedron,QuadraticHexahedron}) = (c.nodes[1], c.nodes[2], c.nodes[3], c.nodes[4], c.nodes[5], c.nodes[6], c.nodes[7], c.nodes[8])
edges(c::Union{Hexahedron,QuadraticHexahedron}) = ((c.nodes[1],c.nodes[2]), (c.nodes[2],c.nodes[3]), (c.nodes[3],c.nodes[4]), (c.nodes[4],c.nodes[1]), (c.nodes[1],c.nodes[5]), (c.nodes[2],c.nodes[6]), (c.nodes[3],c.nodes[7]), (c.nodes[4],c.nodes[8]), (c.nodes[5],c.nodes[6]), (c.nodes[6],c.nodes[7]), (c.nodes[7],c.nodes[8]), (c.nodes[8],c.nodes[5]))
faces(c::Union{Hexahedron,QuadraticHexahedron}) = ((c.nodes[1],c.nodes[4],c.nodes[3],c.nodes[2]), (c.nodes[1],c.nodes[2],c.nodes[6],c.nodes[5]), (c.nodes[2],c.nodes[3],c.nodes[7],c.nodes[6]), (c.nodes[3],c.nodes[4],c.nodes[8],c.nodes[7]), (c.nodes[1],c.nodes[5],c.nodes[8],c.nodes[4]), (c.nodes[5],c.nodes[6],c.nodes[7],c.nodes[8]))

# random stuff
default_interpolation(::Type{Line}) = Lagrange{1,RefCube,1}()
default_interpolation(::Type{QuadraticLine}) = Lagrange{1,RefCube,2}()
default_interpolation(::Type{Triangle}) = Lagrange{2,RefTetrahedron,1}()
default_interpolation(::Type{QuadraticTriangle}) = Lagrange{2,RefTetrahedron,2}()
default_interpolation(::Type{Quadrilateral}) = Lagrange{2,RefCube,1}()
default_interpolation(::Type{QuadraticQuadrilateral}) = Lagrange{2,RefCube,2}()
default_interpolation(::Type{Tetrahedron}) = Lagrange{3,RefTetrahedron,1}()
default_interpolation(::Type{QuadraticTetrahedron}) = Lagrange{2,RefTetrahedron,2}()
default_interpolation(::Type{Hexahedron}) = Lagrange{3,RefCube,1}()
default_interpolation(::Type{QuadraticHexahedron}) = Lagrange{3,RefCube,2}()


text:
"""
Provides logging functionality for Genie apps.
"""
module Logger

using Lumberjack, Millboard, Genie

# color mappings for logging levels -- to be used in STDOUT printing
const colors = Dict{String,Symbol}("info" => :gray, "warn" => :yellow, "debug" => :green, "err" => :red, "error" => :red, "critical" => :magenta)

"""
    log(message, level = "info"; showst::Bool = true) :: Void
    log(message::Any, level::Any = "info"; showst::Bool = false) :: Void
    log(message::String, level::Symbol) :: Void

Logs `message` to all configured logs (STDOUT, FILE, etc) by delegating to `Lumberjack`.
Supported values for `level` are "info", "warn", "debug", "err" / "error", "critical".
If `level` is `error` or `critical` it will also dump the stacktrace onto STDOUT.

# Examples
```julia
julia> Logger.log("hello")

2016-12-21T18:38:09.105 - info: hello

julia> Logger.log("hello", "warn")

2016-12-21T18:38:22.461 - warn: hello

julia> Logger.log("hello", "debug")

2016-12-21T18:38:32.292 - debug: hello

julia> Logger.log("hello", "err")

2016-12-21T18:38:38.403 - err: hello
```
"""
function log(message, level::Any = "info"; showst = false) :: Void
  message = string(message)
  level = string(level)
  level == "err" && (level = "error")

  if ! isdefined(:Genie) || ! isdefined(Genie, :config) || ! Genie.config.suppress_output
    println()
    Lumberjack.log(string(level), string(message))
    println()
  end

  if (level == "critical" || level == "error") && showst
    println()
    stacktrace()
  end

  nothing
end
function log(message::String, level::Symbol; showst::Bool = false) :: Void
  log(message, level == :err ? "error" : string(level), showst = showst)
end

"""
    self_log(message, level::Union{String,Symbol}) :: Void

Basic logging function that does not rely on external logging modules (such as `Lumberjack`).

# Examples
```julia
julia> Logger.self_log("hello", :err)

err 2016-12-21T18:49:00.286
hello

julia> Logger.self_log("hello", :info)

info 2016-12-21T18:49:05.068
hello

julia> Logger.self_log("hello", :debug)

debug 2016-12-21T18:49:11.123
hello
```
"""
function self_log(message, level::Union{String,Symbol}) :: Void
  println()
  print_with_color(colors[string(level)], (string(level), " ", string(Dates.now()), "\n")...)
  print_with_color(colors[string(level)], string(message))
  println()

  nothing
end

"""
    truncate_logged_output(output::AbstractString) :: String

Truncates (shortens) output based on `output_length` settings and appends "..." -- to be used for limiting the output length when logging.

# Examples
```julia
julia> Genie.config.output_length
100

julia> Genie.config.output_length = 10
10

julia> Logger.truncate_logged_output("abc " ^ 10)
"abc abc ab..."
```
"""
function truncate_logged_output(output::AbstractString) :: String
  if length(output) > Genie.config.output_length
    output = output[1:Genie.config.output_length] * "..."
  end

  output
end

"""
    setup_loggers() :: Bool

Sets up default app loggers (STDOUT and per env file loggers) defferring to the `Lumberjack` module.
Automatically invoked.
"""
function setup_loggers() :: Bool
  configure(; modes=["debug", "info", "notice", "warn", "err", "critical", "alert", "emerg"])
  add_truck(LumberjackTruck(STDOUT, nothing, Dict{Any,Any}(:is_colorized => true)), "console")
  IS_IN_APP && add_truck(LumberjackTruck("$(joinpath(Genie.LOG_PATH, Genie.config.app_env)).log", nothing, Dict{Any,Any}(:is_colorized => true)), "file-logger")

  true
end

"""
    empty_log_queue() :: Vector{Tuple{String,Symbol}}

The Genie log queue is used to push log messages in the early phases of framework bootstrap,
when the logger itself is not available. Once the logger is ready, the queue is emptied and the
messages are logged.
Automatically invoked.
"""
function empty_log_queue() :: Vector{Tuple{String,Symbol}}
  for log_message in Genie.GENIE_LOG_QUEUE
    log(log_message...)
  end

  empty!(Genie.GENIE_LOG_QUEUE)
end

"""
    macro location()

Provides a macro that injects the FILE and the LINE where the logger was invoked.
"""
macro location()
  :(Logger.log(" in $(@__FILE__):$(@__LINE__)", :err))
end

setup_loggers()
empty_log_queue()

end