text:
"""
    Interval{T, L <: Bound, R <: Bound}

An `Interval` represents a non-iterable range or span of values (non-interable because,
unlike a `StepRange`, no step is defined).

An `Interval` can be closed (both `first` and `last` are included in the interval), open
(neither `first` nor `last` are included), or half-open. This openness is defined by the
bounds information which is stored as the type parameters `L` and `R`.

### Example

```julia
julia> interval = Interval{Closed,Open}(0, 100)
Interval{Int64,Closed,Open}}(0, 100)

julia> 0 in interval
true

julia> 50 in interval
true

julia> 100 in interval
false

julia> intersect(Interval{Open,Open}(0, 25), Interval{Closed,Closed}(20, 50)
Interval{Int64,Closed,Open}(20, 25)
```

### Infix Constructor: `..`

A closed `Interval` can be constructed with the `..` infix constructor:

```julia
julia> Dates.today() - Dates.Week(1) .. Dates.today()
Interval{Date,Closed,Closed}(2018-01-24, 2018-01-31)
```

### Note on Ordering

The `Interval` constructor will compare `first` and `last`; if it finds that
`first > last`, they will be reversed to ensure that `first < last`. This simplifies
calls to `in` and `intersect`:

```julia
julia> i = Interval{Open,Closed}(Date(2016, 8, 11), Date(2013, 2, 13))
Interval{Date,Closed,Open}(2013-02-13, 2016-08-11)
```

Note that the bounds are also reversed in this case.

See also: [`AnchoredInterval`](@ref)
"""
struct Interval{T, L <: Bound, R <: Bound} <: AbstractInterval{T,L,R}
    first::T
    last::T

    function Interval{T,L,R}(f::T, l::T) where {T, L <: Bounded, R <: Bounded}
        # Ensure that `first` preceeds `last`.
        if f ≤ l
            return new{T,L,R}(f, l)
        elseif l ≤ f
            # Note: Most calls to this inner constructor will be from other constructors
            # which may make it hard to identify the source of this deprecation. Use
            # `--depwarn=error` to see a full stack trace.
            Base.depwarn(
                "Constructing an `Interval{T,X,Y}(x, y)` " *
                "where `x > y` is deprecated, use `Interval{T,Y,X}(y, x)` instead.",
                :Interval,
            )
            return new{T,R,L}(l, f)
        else
            throw(ArgumentError("Unable to determine an ordering between: $f and $l"))
        end
    end

    function Interval{T,L,R}(f::Nothing, l::T) where {T, L <: Unbounded, R <: Bounded}
        # Note: Using `<` enforces that the type `T` defines `isless`
        if !(l ≤ l)
            throw(ArgumentError(
                "Unable to determine an ordering between $l and other values of type $T"
            ))
        end
        return new{T,L,R}(l, l)
    end

    function Interval{T,L,R}(f::T, l::Nothing) where {T, L <: Bounded, R <: Unbounded}
        # Note: Using `<` enforces that the type `T` defines `isless`
        if !(f ≤ f)
            throw(ArgumentError(
                "Unable to determine an ordering between $f and other values of type $T"
            ))
        end
        return new{T,L,R}(f, f)
    end

    function Interval{T,L,R}(f::Nothing, l::Nothing) where {T, L <: Unbounded, R <: Unbounded}
        return new{T,L,R}()
    end
end

function Interval{T,L,R}(f, l) where {T, L <: Bounded, R <: Bounded}
    return Interval{T,L,R}(convert(T, f), convert(T, l))
end
function Interval{T,L,R}(f, l::Nothing) where {T, L <: Bounded, R <: Unbounded}
    return Interval{T,L,R}(convert(T, f), l)
end
function Interval{T,L,R}(f::Nothing, l) where {T, L <: Unbounded, R <: Bounded}
    return Interval{T,L,R}(f, convert(T, l))
end

Interval{L,R}(f::T, l::T) where {T,L,R} = Interval{T,L,R}(f, l)
Interval{L,R}(f, l) where {L,R} = Interval{promote_type(typeof(f), typeof(l)), L, R}(f, l)
Interval{L,R}(f::Nothing, l::T) where {T,L,R} = Interval{T,L,R}(f, l)
Interval{L,R}(f::T, l::Nothing) where {T,L,R} = Interval{T,L,R}(f, l)
Interval{L,R}(f::Nothing, l::Nothing) where {L,R} = Interval{Nothing,L,R}(f, l)

Interval{T}(f, l) where T = Interval{T, Closed, Closed}(f, l)
Interval{T}(f::Nothing, l) where T = Interval{T, Unbounded, Closed}(f, l)
Interval{T}(f, l::Nothing) where T = Interval{T, Closed, Unbounded}(f, l)
Interval{T}(f::Nothing, l::Nothing) where T = Interval{T, Unbounded, Unbounded}(f, l)

Interval(f::T, l::T) where T = Interval{T}(f, l)
Interval(f, l) = Interval(promote(f, l)...)
Interval(f::Nothing, l::T) where T = Interval{T}(f, l)
Interval(f::T, l::Nothing) where T = Interval{T}(f, l)
Interval(f::Nothing, l::Nothing) = Interval{Nothing}(f, l)

(..)(first, last) = Interval(first, last)

# In Julia 0.7 constructors no longer automatically fall back to using `convert`
Interval(interval::AbstractInterval) = convert(Interval, interval)
Interval{T}(interval::AbstractInterval) where T = convert(Interval{T}, interval)

# Endpoint constructors
function Interval{T}(left::LeftEndpoint{T,L}, right::RightEndpoint{T,R}) where {T,L,R}
    Interval{T,L,R}(endpoint(left), endpoint(right))
end

function Interval{T}(left::LeftEndpoint, right::RightEndpoint) where T
    Interval{T, bound_type(left), bound_type(right)}(endpoint(left), endpoint(right))
end

function Interval(left::LeftEndpoint{S}, right::RightEndpoint{T}) where {S,T}
    Interval{promote_type(S, T)}(left, right)
end

# Empty Intervals
Interval{T}() where T = Interval{T, Open, Open}(zero(T), zero(T))
Interval{T}() where T <: TimeType = Interval{T, Open, Open}(T(0), T(0))

function Interval{T}() where T <: ZonedDateTime
    return Interval{T, Open, Open}(T(0, tz"UTC"), T(0, tz"UTC"))
end

Base.copy(x::T) where T <: Interval = T(x.first, x.last)

function Base.hash(interval::AbstractInterval, h::UInt)
    h = hash(LeftEndpoint(interval), h)
    h = hash(RightEndpoint(interval), h)
    return h
end

##### ACCESSORS #####

function Base.first(interval::Interval{T,L,R}) where {T,L,R}
    return L !== Unbounded ? interval.first : nothing
end

function Base.last(interval::Interval{T,L,R}) where {T,L,R}
    return R !== Unbounded ? interval.last : nothing
end

function span(interval::Interval)
    if isbounded(interval)
        interval.last - interval.first
    else
        throw(DomainError(
            "unbounded endpoint(s)",
            "Unable to determine the span of an non-bounded interval",
        ))
    end
end

isclosed(interval::AbstractInterval{T,L,R}) where {T,L,R} = L === Closed && R === Closed
Base.isopen(interval::AbstractInterval{T,L,R}) where {T,L,R} = L === Open && R === Open
isunbounded(interval::AbstractInterval{T,L,R}) where {T,L,R} = L === Unbounded && R === Unbounded
isbounded(interval::AbstractInterval{T,L,R}) where {T,L,R} = L !== Unbounded && R !== Unbounded

function Base.minimum(interval::AbstractInterval{T,L,R}; increment=nothing) where {T,L,R}
    return L === Unbounded ? typemin(T) : first(interval)
end

function Base.minimum(interval::AbstractInterval{T,Open,R}; increment=eps(T)) where {T,R}
    isempty(interval) && throw(BoundsError(interval, 0))
    min_val = first(interval) + increment
    # Since intervals can't have NaN, we can just use !isfinite to check if infinite
    !isfinite(min_val) && return typemin(T)
    min_val ∈ interval && return min_val
    throw(BoundsError(interval, min_val))
end

function Base.minimum(interval::AbstractInterval{T,Open,R}) where {T<:Integer,R}
    return minimum(interval, increment=one(T))
end

function Base.minimum(interval::AbstractInterval{T,Open,R}; increment=nothing) where {T<:AbstractFloat,R}
    isempty(interval) && throw(BoundsError(interval, 0))
    min_val = first(interval)
    # Since intervals can't have NaN, we can just use !isfinite to check if infinite
    next_val = if !isfinite(min_val) || increment === nothing
        nextfloat(min_val)
    else
        min_val + increment
    end
    next_val ∈ interval && return next_val
    throw(BoundsError(interval, next_val))
end

function Base.maximum(interval::AbstractInterval{T,L,R}; increment=nothing) where {T,L,R}
    return R === Unbounded ? typemax(T) : last(interval)
end

function Base.maximum(interval::AbstractInterval{T,L,Open}; increment=eps(T)) where {T,L}
    isempty(interval) && throw(BoundsError(interval, 0))
    max_val = last(interval) - increment
    # Since intervals can't have NaN, we can just use !isfinite to check if infinite
    !isfinite(max_val) && return typemax(T)
    max_val ∈ interval && return max_val
    throw(BoundsError(interval, max_val))
end

function Base.maximum(interval::AbstractInterval{T,L,Open}) where {T<:Integer,L}
    return maximum(interval, increment=one(T))
end

function Base.maximum(interval::AbstractInterval{T,L,Open}; increment=nothing) where {T<:AbstractFloat,L}
    isempty(interval) && throw(BoundsError(interval, 0))
    max_val = last(interval)
    next_val = if !isfinite(max_val) || increment === nothing
        prevfloat(max_val)
    else
        max_val - increment
    end
    next_val ∈ interval && return next_val
    throw(BoundsError(interval, next_val))
end

##### CONVERSION #####

# Allows an interval to be converted to a scalar when the set contained by the interval only
# contains a single element.
function Base.convert(::Type{T}, interval::Interval{T}) where T
    if first(interval) == last(interval) && isclosed(interval)
        return first(interval)
    else
        throw(DomainError(interval, "The interval is not closed with coinciding endpoints"))
    end
end

##### DISPLAY #####

function Base.show(io::IO, interval::Interval{T,L,R}) where {T,L,R}
    if get(io, :compact, false)
        print(io, interval)
    else
        print(io, "$(typeof(interval))(")
        L === Unbounded ? print(io, "nothing") : show(io, interval.first)
        print(io, ", ")
        R === Unbounded ? print(io, "nothing") : show(io, interval.last)
        print(io, ")")
    end
end

function Base.print(io::IO, interval::AbstractInterval{T,L,R}) where {T,L,R}
    # Print to io in order to keep properties like :limit and :compact
    if get(io, :compact, false)
        io = IOContext(io, :limit=>true)
    end

    print(
        io,
        L === Closed ? "[" : "(",
        L === Unbounded ? "" : first(interval),
        " .. ",
        R === Unbounded ? "" : last(interval),
        R === Closed ? "]" : ")",
    )
end

##### ARITHMETIC #####

Base.:+(a::T, b) where {T <: Interval} = T(first(a) + b, last(a) + b)

Base.:+(a, b::Interval) = b + a
Base.:-(a::Interval, b) = a + -b
Base.:-(a, b::Interval) = a + -b
Base.:-(a::Interval{T,L,R}) where {T,L,R} = Interval{T,R,L}(-last(a), -first(a))

##### EQUALITY #####

function Base.:(==)(a::AbstractInterval, b::AbstractInterval)
    return LeftEndpoint(a) == LeftEndpoint(b) && RightEndpoint(a) == RightEndpoint(b)
end

function Base.isequal(a::AbstractInterval, b::AbstractInterval)
    le = isequal(LeftEndpoint(a), LeftEndpoint(b))
    re = isequal(RightEndpoint(a), RightEndpoint(b))
    return le && re
end

# While it might be convincingly argued that this should define < instead of isless (see
# https://github.com/invenia/Intervals.jl/issues/14), this breaks sort.
Base.isless(a::AbstractInterval, b) = LeftEndpoint(a) < b
Base.isless(a, b::AbstractInterval) = a < LeftEndpoint(b)

less_than_disjoint(a::AbstractInterval, b) = RightEndpoint(a) < b
less_than_disjoint(a, b::AbstractInterval) = a < LeftEndpoint(b)

function Base.:isless(a::AbstractInterval, b::AbstractInterval)
    return LeftEndpoint(a) < LeftEndpoint(b)
end

function less_than_disjoint(a::AbstractInterval, b::AbstractInterval)
    return RightEndpoint(a) < LeftEndpoint(b)
end

greater_than_disjoint(a, b) = less_than_disjoint(b, a)

"""
    ≪(a::AbstractInterval, b::AbstractInterval) -> Bool
    less_than_disjoint(a::AbstractInterval, b::AbstractInterval) -> Bool

Less-than-and-disjoint comparison operator. Returns `true` if `a` is less than `b` and they
are disjoint (they do not overlap).

```
julia> 0..10 ≪ 10..20
false

julia> 0..10 ≪ 11..20
true
```
"""
≪(a, b) = less_than_disjoint(a, b)
# ≪̸(a, b) = !≪(a, b)

"""
    ≫(a::AbstractInterval, b::AbstractInterval) -> Bool
    greater_than_disjoint(a::AbstractInterval, b::AbstractInterval) -> Bool

Greater-than-and-disjoint comparison operator. Returns `true` if `a` is greater than `b` and
they are disjoint (they do not overlap).

```
julia> 10..20 ≫ 0..10
false

julia> 11..20 ≫ 0..10
true
```
"""
≫(a, b) = greater_than_disjoint(a, b)
# ≫̸(a, b) = !≫(a, b)

##### SET OPERATIONS #####

Base.isempty(i::AbstractInterval) = LeftEndpoint(i) > RightEndpoint(i)
Base.in(a, b::AbstractInterval) = !(a ≫ b || a ≪ b)

function Base.in(a::AbstractInterval, b::AbstractInterval)
    # Intervals should be compared with set operations
    throw(ArgumentError("Intervals can not be compared with `in`. Use `issubset` instead."))
end

function Base.issubset(a::AbstractInterval, b::AbstractInterval)
    return LeftEndpoint(a) ≥ LeftEndpoint(b) && RightEndpoint(a) ≤ RightEndpoint(b)
end

function isdisjoint(a::AbstractInterval, b::AbstractInterval)
    return RightEndpoint(a) < LeftEndpoint(b) || LeftEndpoint(a) > RightEndpoint(b)
end

Base.:⊈(a::AbstractInterval, b::AbstractInterval) = !issubset(a, b)
Base.:⊉(a::AbstractInterval, b::AbstractInterval) = !issubset(b, a)

function overlaps(a::AbstractInterval, b::AbstractInterval)
    left = max(LeftEndpoint(a), LeftEndpoint(b))
    right = min(RightEndpoint(a), RightEndpoint(b))

    return left <= right
end

function contiguous(a::AbstractInterval, b::AbstractInterval)
    left = max(LeftEndpoint(a), LeftEndpoint(b))
    right = min(RightEndpoint(a), RightEndpoint(b))

    return (
        !isunbounded(right) && !isunbounded(left) &&
        right.endpoint == left.endpoint && isclosed(left) != isclosed(right)
    )
end

function Base.intersect(a::AbstractInterval{T}, b::AbstractInterval{T}) where T
    !overlaps(a,b) && return Interval{T}()
    left = max(LeftEndpoint(a), LeftEndpoint(b))
    right = min(RightEndpoint(a), RightEndpoint(b))

    return Interval{T}(left, right)
end

function Base.intersect(a::AbstractInterval{S}, b::AbstractInterval{T}) where {S,T}
    !overlaps(a, b) && return Interval{promote_type(S, T)}()
    left = max(LeftEndpoint(a), LeftEndpoint(b))
    right = min(RightEndpoint(a), RightEndpoint(b))

    return Interval(left, right)
end

# There is power in a union.
"""
    union(intervals::AbstractVector{<:AbstractInterval})

Flattens a vector of overlapping intervals into a new, smaller vector containing only
non-overlapping intervals.
"""
function Base.union(intervals::AbstractVector{<:AbstractInterval})
    return union!(convert(Vector{AbstractInterval}, intervals))
end

"""
    union!(intervals::AbstractVector{<:Union{Interval, AbstractInterval}})

Flattens a vector of overlapping intervals in-place to be a smaller vector containing only
non-overlapping intervals.
"""
function Base.union!(intervals::Union{AbstractVector{<:Interval}, AbstractVector{AbstractInterval}})
    sort!(intervals)

    i = 2
    n = length(intervals)
    while i <= n
        prev = intervals[i - 1]
        curr = intervals[i]

        # If the current and previous intervals don't meet then move along
        if !overlaps(prev, curr) && !contiguous(prev, curr)
            i = i + 1

        # If the two intervals meet then we absorb the current interval into
        # the previous one.
        else
            intervals[i - 1] = merge(prev, curr)
            deleteat!(intervals, i)
            n -= 1
        end
    end

    return intervals
end

"""
    superset(intervals::AbstractArray{<:AbstractInterval}) -> Interval

Create the smallest single interval which encompasses all of the provided intervals.
"""
function superset(intervals::AbstractArray{<:AbstractInterval})
    left = minimum(LeftEndpoint.(intervals))
    right = maximum(RightEndpoint.(intervals))

    return Interval(left, right)
end

function Base.merge(a::AbstractInterval, b::AbstractInterval)
    if !overlaps(a, b) && !contiguous(a, b)
        throw(ArgumentError("$a and $b are neither overlapping or contiguous."))
    end

    left = min(LeftEndpoint(a), LeftEndpoint(b))
    right = max(RightEndpoint(a), RightEndpoint(b))
    return Interval(left, right)
end

##### ROUNDING #####
const RoundingFunctionTypes = Union{typeof(floor), typeof(ceil), typeof(round)}

for f in (:floor, :ceil, :round)
    @eval begin
        """
           $($f)(interval::Interval, args...; on::Symbol)

        Round the interval by applying `$($f)` to a single endpoint, then shifting the
        interval so that the span remains the same. The `on` keyword determines which
        endpoint the rounding will be applied to. Valid options are `:left` or `:right`.
        """
        function Base.$f(interval::Interval, args...; on::Symbol)
            return _round($f, interval, Val(on), args...)
        end
    end
end

function _round(f::RoundingFunctionTypes, interval::Interval, on::Val{:anchor}, args...)
    throw(ArgumentError(":anchor is only usable with an AnchoredInterval."))
end

function _round(
    f::RoundingFunctionTypes, interval::Interval{T,L,R}, on::Val{:left}, args...
) where {T, L <: Bounded, R <: Bounded}
    left_val = f(first(interval), args...)
    return Interval{T,L,R}(left_val, left_val + span(interval))
end

function _round(
    f::RoundingFunctionTypes, interval::Interval{T,L,R}, on::Val{:left}, args...
) where {T, L <: Bounded, R <: Unbounded}
    left_val = f(first(interval), args...)
    return Interval{T,L,R}(left_val, nothing)
end

function _round(
    f::RoundingFunctionTypes, interval::Interval{T,L,R}, on::Val{:left}, args...
) where {T, L <: Unbounded, R <: Bound}
    return interval
end

function _round(
    f::RoundingFunctionTypes, interval::Interval{T,L,R}, on::Val{:right}, args...
) where {T, L <: Bounded, R <: Bounded}
    right_val = f(last(interval), args...)
    return Interval{T,L,R}(right_val - span(interval), right_val)
end

function _round(
    f::RoundingFunctionTypes, interval::Interval{T,L,R}, on::Val{:right}, args...
) where {T, L <: Unbounded, R <: Bounded}
    right_val = f(last(interval), args...)
    return Interval{T,L,R}(nothing, right_val)
end

function _round(
    f::RoundingFunctionTypes, interval::Interval{T,L,R}, on::Val{:right}, args...
) where {T, L <: Bound, R <: Unbounded}
    return interval
end

##### TIME ZONES #####

function TimeZones.astimezone(i::Interval{ZonedDateTime, L, R}, tz::TimeZone) where {L,R}
    return Interval{ZonedDateTime, L, R}(astimezone(first(i), tz), astimezone(last(i), tz))
end

function TimeZones.timezone(i::Interval{ZonedDateTime})
    if timezone(first(i)) != timezone(last(i))
        throw(ArgumentError("Interval $i contains mixed timezones."))
    end
    return timezone(first(i))
end


text:
using Distributions
using Base.LinAlg.BLAS

function sample(::Type{Bernoulli}, means::Mat{Float64})
    s = zeros(means)
    r = rand(size(means))
    @simd for i=1:length(means)
        @inbounds s[i] = r[i] < means[i] ? 1.0 : 0.0
    end

    return s
end

function sample(::Type{Gaussian}, means::Mat{Float64})
    sigma2 = 1                   # using fixed standard diviation
    samples = zeros(size(means))
    for j=1:size(means, 2), i=1:size(means, 1)
        samples[i, j] = rand(Normal(means[i, j], sigma2))
    end
    return samples
end

function sample_hiddens{V,H}(rbm::RBM{V,H}, vis::Mat{Float64})
    means = ProbHidCondOnVis(rbm, vis)
    return sample(H, means), means
end

function sample_visibles{V,H}(rbm::RBM{V,H}, hid::Mat{Float64})
    means = ProbVisCondOnHid(rbm, hid)
    return sample(V, means), means
end

function gibbs(rbm::RBM, vis::Mat{Float64}; n_times=1)
    v_pos = vis
    h_samp, h_pos = sample_hiddens(rbm, v_pos)
    h_neg = Array(Float64,0,0)::Mat{Float64}
    v_neg = Array(Float64,0,0)::Mat{Float64}
    if n_times > 0
    # Save computation by setting `n_times=0` in the case
    # of persistent CD.
        v_neg = sample_visibles(rbm, h_samp)
        h_samp, h_neg = sample_hiddens(rbm, v_neg)
        for i=1:n_times-1
            v_neg = sample_visibles(rbm, h_samp)
            h_samp, h_neg = sample_hiddens(rbm, v_neg)
        end
    end
    return v_pos, h_pos, v_neg, h_neg
end

function MCMC(rbm::RBM, init::Mat{Float64}; iterations=1, StartMode="visible")
    if StartMode == "visible"
    # In this first mode we assume that we are starting from the visible samples. E.g. in
    # the case of binary RBM, we should be starting with binary samples.
        vis_samples = copy(init)                                    # Start from the visible samples
        vis_means   = copy(init)                                    # Giving a starting point for the means
        hid_samples, hid_means = sample_hiddens(rbm,vis_samples)     # Get the first hidden means [NMF-ish]
    end

    if StartMode == "hidden"
    # In this second mode we assume that we are starting from a set of hidden
    # samples. Because of this, we increment the iteration count by 1
        hid_samples = copy(init)
        hid_means = copy(init)
        iterations+=1
    end

    for i=1:iterations-1
        vis_samples, vis_means = sample_visibles(rbm,hid_samples)          # Sample the visible units from true distribution
        hid_samples, hid_means = sample_hiddens(rbm,vis_samples)           # Update the hidden unit means, a NMF-ish approach
    end

    return vis_samples, hid_samples, vis_means, hid_means
end

ext: jl


text:
module Stats

import Base

export logsumexp

"""
    logsumexp(x,y)
    logsumexp(x)

Return `log(exp(x) + exp(y))` but stably without overflow.

Single-argument version returns `log(exp(x[1]) + exp(x[2]) + ...)`.
"""
function logsumexp(x, y)
    if x == -Inf
        y
    elseif y == -Inf
        x
    elseif x > y
        x + log1p(exp(y-x))
    else
        y + log1p(exp(x-y))
    end
end

function logsumexp(x::AbstractArray{T}) where {T <: Number}
    sum = -Inf
    for i in eachindex(x)
        sum = logsumexp(x[i], sum)
    end
    sum
end

end


text:
# Collect Time To Finish Task (TTFT)

# Using `@eval` to avoid missing compile time, see the `@time` docstring for more info.
@timeit TOUT "Pluto.Cell" cell = @eval Pluto.Cell("1 + 1")

@timeit TOUT "Pluto.Notebook" nb = @eval Pluto.Notebook([cell])

module Foo end
@timeit TOUT "PlutoRunner.run_expression" @eval Pluto.PlutoRunner.run_expression(Foo, :(1 + 1), Pluto.uuid1(), nothing);

function wait_for_ready(notebook::Pluto.Notebook)
    while notebook.process_status != Pluto.ProcessStatus.ready
        sleep(0.1)
    end
end

🍭 = Pluto.ServerSession()
🍭.options.server.disable_writing_notebook_files = true
🍭.options.evaluation.workspace_use_distributed = false

path = joinpath(pkgdir(Pluto), "sample", "Basic.jl")

@timeit TOUT "SessionActions.open" nb = @eval Pluto.SessionActions.open(🍭, path; run_async=true)

wait_for_ready(nb)

@timeit TOUT "SessionActions.shutdown" @eval Pluto.SessionActions.shutdown(🍭, nb; async=true)

# According to SnoopCompile, this is a big part of the time for `Pluto.run()`.
# However, it's very tricky to measure this via the `Pluto.run` below.
@timeit TOUT "Configuration.from_flat_kwargs" @eval Pluto.Configuration.from_flat_kwargs()

# Compile HTTP get.
HTTP.get("http://github.com")

@timeit TOUT "Pluto.run" server_task = @eval let
    port = 13435
    options = Pluto.Configuration.from_flat_kwargs(; port, launch_browser=false, workspace_use_distributed=false, require_secret_for_access=false, require_secret_for_open_links=false)
    🍭 = Pluto.ServerSession(; options)
    server_task = @async Pluto.run(🍭)

    # Give the async task time to start.
    sleep(1)

    HTTP.get("http://localhost:$port/edit").status == 200
    server_task
end


text:
#=
This model is taken from [1]. The theoretical derivations of the analytic solution can be found in [2].

[1] Malakiyeh, Mohammad Mahdi, Saeed Shojaee, and Klaus-Jürgen Bathe. "The Bathe time integration method revisited for prescribing desired numerical dissipation." Computers & Structures 212 (2019): 289-298.

[2] Mechanical Vibrations, Gerardin et al, page 250-251.
=#

using ReachabilityAnalysis, LinearAlgebra, LazySets
using ReachabilityAnalysis: solve, discretize
using LazySets.Arrays
using SparseArrays

LazySets.set_ztol(Float64, 1e-15)
LazySets.set_atol(Float64, 1e-15)
LazySets.set_rtol(Float64, 1e-15)

function clamped_matrices(; N=1000, E=30e6, ρ=7.3e-4, A=1, L=200)
    ℓ = L / N
    K = (E * A / ℓ) * SymTridiagonal(fill(2.0, N), fill(-1.0, N))
    K[end, end] = E * A / ℓ

    M = (ρ * A * ℓ / 2) * Diagonal(vcat(fill(2.0, N-1), 1.0))
    M[end, end] = ρ * A * ℓ / 2

    return M, K
end

function clamped_free(; N)
    M, K = clamped_matrices(N=N)
    C = zeros(N, N) # no damping
    sys = SecondOrderLinearContinuousSystem(M, C, K)
end

function clamped_forced(; N, F=10e3, E=30e6, A=1)
    M, K = clamped_matrices(N=N)
    C = zeros(N, N) # no damping
    F = vcat(zeros(N-1), F) # the right-most node is excited
    sys = SecondOrderAffineContinuousSystem(M, C, K, F)
end

function clamped(; a=0.0, b=0.0, # ignored if damped = false
                   constant=true,
                   N,            # number of elements
                   homogeneize,  # flag to homogeneize the system
                   damped)       # flag to consider damping C = a*K + b*M

    sys = clamped_forced(N=N)
    M = sys.M; K = sys.K; F = sys.b

    invM = inv(M)
    ZN = zeros(N, N)
    IN = Matrix(1.0I, N, N)

    if damped
        C = a*K + b*M
        A = [ZN           IN       ;
            -invM*K       -invM*C]

    else
        A = [ZN           IN ;
            -invM*K       ZN]
    end

    f0 = vcat(zeros(N), invM * F)

    if homogeneize
        n = 2N
        Aext = zeros(n+1, n+1)
        Aext[1:n, 1:n] .= A
        Aext[1:n, n+1] .= f0
        Aext = sparse(Aext)
        S = @system(x' = Aext*x)

    else
        S = @system(x' = A*x + f0)
    end

    if !constant
        @assert homogeneize == false

        # model time-varying forcing
        X = Universe(statedim(S))
        ΔF = Interval(0.99, 1.01)
        S = @system(x' = S.A * x + S.c * u, u ∈ ΔF, x ∈ X)
    end

    return S
end

# "nominal" step size
const dtn = 9.88e-7

#=
# problem instances
s5 = clamped_forced(N=5);
s10 = clamped_forced(N=10);
s20 = clamped_forced(N=20);
s30 = clamped_forced(N=30);
s50 = clamped_forced(N=50);
s100 = clamped_forced(N=100);
s200 = clamped_forced(N=200);
s500 = clamped_forced(N=500);
s1000 = clamped_forced(N=1000);
=#

#=

# singleton initial condition
X0(N) = Singleton(zeros(2N))
prob(s) = InitialValueProblem(s, X0(statedim(s)))

# box initial condition
prob_box(s, e) = InitialValueProblem(s, X0(statedim(s)) + BallInf(2*statedim(s), e))

# solve with GLGM06
function solve_glg(s; δ=dtn, NSTEPS=10, model=Forward(inv=false))
    ivp = prob(s)
    alg = GLGM06(δ=δ, approx_model=model)
    sol = solve(ivp, NSTEPS=NSTEPS, alg=alg, homogeneize=true)
end

# solve with LGG09 for nodeidx
function solve_lgg(s; δ=dtn, NSTEPS=100, nodeidx, model=Forward(inv=false))
    ivp = prob(s)
    n = statedim(s)
    eplus = SingleEntryVector(nodeidx, 2n+1, 1.0)
    eminus = SingleEntryVector(nodeidx, 2n+1, -1.0)

    dirs = [eminus, eplus]

    alg = LGG09(δ=δ, template=dirs, approx_model=model)
    sol = solve(ivp, NSTEPS=NSTEPS, alg=alg, homogeneize=true)
end

# solve with BOX
# usage: sol = solve_box(s10, NSTEPS=100)
function solve_box(s; δ=dtn, NSTEPS=100, model=Forward(inv=false))

    # p = prob(s)
    # ph = homogeneize(normalize(p))
    # @time pd = discretize(ph, dt, model) # ~35 sec for 1000 nodes

    alg = BOX(δ=δ, approx_model=model)
    p = prob(s)
    sbox = solve(p, alg=alg, NSTEPS=NSTEPS, homogeneize=true)
    return sbox
end
=#

# Ec. (4.181) ref Geradin & Rixen (2015)
function _sol_analytic_clamped(x, t, Nmax)

    F = 10E3
    E = 30E6
    dens = 7.3E-4
    L = 200
    A = 1

    # reescalar
    nMasses = 1000
    x = x*L/nMasses

    m = A*dens
    c = sqrt(E*A/m)

    α = 8*F*L/π^2/E/A

    β = π/2/L
    u = sin.( β * x ) * ( 1.0 .- cos.( β * c * t ) )
    v = sin.( β * x ) * β * c * sin.( β * c * t )

    for s in 2:Nmax
        u .+= (-1)^(s-1)/(2*s-1)^2 * sin.( (2*s-1)* β * x ) * ( 1.0 .- cos.( (2*s-1)* β * c * t ) )
        v .+= (-1)^(s-1)/(2*s-1)^2 * sin.( (2*s-1)* β * x ) * ( (2*s-1)* β * c * sin.( (2*s-1)* β * c * t ) )
    end

    return u*α, v*α

end

max_stars_repo_name: JuliaReach/ClampedFreeBar

max_issues_repo_name: JuliaReach/ClampedFreeBar

max_forks_repo_name: JuliaReach/ClampedFreeBar


text:
import Pkg
Pkg.activate("./")

using PDMesh

c = Cuboid([0 5; -10 10; 0 3])

# function rand_(a, b)
#     return a + rand()*(b-a)
# end

# for i in 1:40
#     global obj
#     center = [rand_(-10, 10), rand_(-10, 10), rand_(0, 3)]
#     radius = 0.2 + 1.0*rand()
#     obj = changetype(obj, out -> begin x=out[:x]; sum((x .- vec(center)).^2, dims=1) .< radius^2 end, 4)
#     obj = changetype(obj, out -> begin x=out[:x]; sum((x .- vec(center)).^2, dims=1) .< (radius - 0.1)^2 end, 5)
# end

# obj = delete!(obj, out -> begin x=out[:x]; sum(x[1:2, :].^2, dims=1) .> 4.0^2 end)
# obj = velocity(obj, out -> begin type=out[5]; type.==1 end, [1.0, 0.0, 0.0])

obj = copy(c)

for i in 1:1000
    global obj
    obj = changetype(obj, out -> begin x=out[:x]; sum(x[1:2, :].^2, dims=1) .< 3.0^2 end, 2)
    obj = changetype(obj, out -> begin x=out[:x]; sum(x[1:2, :].^2, dims=1) .< 2.0^2 end, 3)
    obj = move(obj, by=[5.0, 0.0, 0.0])
    obj = rotate(obj, angle=2, point=[0.0, 0.0, 0.0], vector_=[1.0, 1.0, 0.0])
    obj = combine(obj, c)
end

x, v, y, vol, type = create(obj, resolution=1, rand_=0.0, type=1)

write_data("./examples/ovito.data", y, v, type, vol)
write_data("./examples/ovito1.data", 1.1y, v, type, vol)

max_stars_repo_name: sajidmannan/PDMesh

max_issues_repo_name: sajidmannan/PDMesh

max_forks_repo_name: sajidmannan/PDMesh


text:
using Base.Test

require("../src/StatsdClient.jl")
using StatsdClient

# server socket
function setup_statsd()
    svr = UdpSocket()
    bind(svr,ip"127.0.0.1",2376)
    # configure client with server address
    server = StatsdClient.Statsd("127.0.0.1",2376)
    (svr, server)
end

c = Condition()
function test_statsd(svr,val)
    @async begin
        arg = bytestring(recv(svr))
        println(arg,":",val)
        @test arg == val
        notify(c)
    end
end

function test_increment()
    svr, server = setup_statsd()
    test_statsd(svr,"test.increment:1|c")
    StatsdClient.increment(server,"test.increment")
    wait(c)
    close(svr)
    true
end

function test_decrement()
    svr, server = setup_statsd()
    test_statsd(svr,"test.decrement:-1|c")
    StatsdClient.decrement(server,"test.decrement")
    wait(c)
    close(svr)
    true
end

function test_count()
    svr, server = setup_statsd()
    test_statsd(svr,"test.count:412|c")
    StatsdClient.count(server,"test.count",412)
    wait(c)
    close(svr)
    true
end

function test_timing()
    svr, server = setup_statsd()
    test_statsd(svr,"test.timing:38|ms")
    StatsdClient.timing(server,"test.timing",38)
    wait(c)
    close(svr)
    true
end

function test_gauge()
    svr, server = setup_statsd()
    test_statsd(svr,"test.gauge:123|g")
    StatsdClient.gauge(server,"test.gauge",123)
    wait(c)
    close(svr)
    true
end

function test_set()
    svr, server = setup_statsd()
    test_statsd(svr,"test.set:7623|s")
    StatsdClient.set(server,"test.set","7623")
    wait(c)
    close(svr)
    true
end

function test_sample_rate()
    svr, server = setup_statsd()
    test_statsd(svr,"test.sample_rate:3|c|@1.0")
    StatsdClient.count(server,"test.sample_rate",3,1.0)
    wait(c)
    close(svr)
    true
end

@test test_increment()
@test test_decrement()
@test test_count()
@test test_timing()
@test test_gauge()
@test test_set()
@test test_sample_rate()


text:
using SlackSDK
using Documenter

DocMeta.setdocmeta!(SlackSDK, :DocTestSetup, :(using SlackSDK); recursive=true)

indexio = IOBuffer()
bodyio = IOBuffer()

for x in SlackSDK.methodNames
    parts =  split(x, "_")
    for i in 2:length(parts)
        z = parts[i]
        parts[i] = uppercase(z[1]) * z[2:end]
    end
    fname = join(parts, "")
    fn = getfield(SlackSDK, Symbol(fname))
    api_object_name = join(map(x->titlecase(x, strict=false),parts[1:end-1])) * "Api"
    for m in methods(fn)
        tv, decls, file, line = Base.arg_decl_parts(m)
        if decls[2][2] != "SlackSDK.Web.$api_object_name" || (length(decls)>2 && decls[3][1] == "response_stream")
            continue
        end
        fndocs = string(Base.Docs.doc(Base.Docs.Binding(SlackSDK.Web, Symbol(fname))))
        u = length(fndocs) ÷ 2    #We get duplicate docstrings
        fndocs = fndocs[1:u+1]
        fndocs = replace(fndocs, "*"=>"_")
        fndocs = replace(fndocs, "Param:"=>"\n\tParam:")
        fndocs = replace(fndocs, "Param:"=>"\n\tParam:", count=1)
        fndocs = replace(fndocs, "Return:"=>"\n\n\tReturn:")
        write(bodyio, """
        ### `$x`

        $fndocs

        [https://api.slack.com/methods/$(replace(x,"_"=>"."))](https://api.slack.com/methods/$(replace(x,"_"=>".")))

        """)
        write(indexio, "- [`$x`](#$x)\n")
    end
end

f = open(joinpath(@__DIR__, "src", "methods.md"), "w+")
write(f, """
# Slack API methods

This package exposes the following API methods. Each of these methods can be called as member functions of the `WebClient` stuct.

Required arguments should be passed as positional arguments to these functions, and optional arguments are passed as keyword arguments.
There are some documented [examples](./index.md#Examples) of how to call these methods.

There are some [discrepancies](https://github.com/slackapi/slack-api-specs/issues/70) between the online documentation and the api
specs about which arguments are required. Refer to this page to double check the parameter names and types.
Required parameters will be explicitly marked as such below.

Each methods documentation below also links to Slack's online documentation about that method.

""")
write(f, String(take!(indexio)))
write(f, "\n\n\n")
write(f, String(take!(bodyio)))
close(f)

makedocs(;
    modules=[SlackSDK],
    authors="Avik Sengupta",
    repo="https://github.com/aviks/SlackSDK.jl/blob/{commit}{path}#{line}",
    sitename="SlackSDK.jl",
    doctest=false,
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Methods" => "methods.md",
        "Types" => "types.md"
    ],
)


text:
import TextParse: tryparsenext

function _iterate(parser::AbstractToken, sequence, till, before_i, next_i, state,opts=TextParse.default_opts)
    if parser isa CombinedParser
        @warn "define _iterate(parser::$(typeof(parser)), sequence, till, start_i, next_i, state::$(typeof(state)))"
        return nothing
    end
    if state === nothing
        r,next_i_ = tryparsenext(parser, sequence, next_i, till,opts)
        if isnull(r)
            nothing
        else
            NCodeunitsState(next_i,next_i_,get(r))
        end
    else
        nothing
    end
end

"""
    TextParse.tryparsenext(x::CombinedParser,str,i,till,opts=TextParse.default_opts)

TextParse.jl integrates with CombinedParsers.jl both ways.

```jldoctest
julia> using TextParse

julia> p = ("Number:" * Repeat(' ') * TextParse.Numeric(Int))[3]
🗄 Sequence[3]
├─ Number\\:
├─ \\ *  |> Repeat
└─ <Int64>
::Int64

julia> parse(p, "Number:    42")
42

julia> TextParse.tryparsenext(p, "Number:    42")
(Nullable{Int64}(42), 14)
```

"""
function TextParse.tryparsenext(x::CombinedParser,str,i,till,opts=TextParse.default_opts)
    s = _iterate(x,str,till,i,nothing)
    if s === nothing
        Nullable{result_type(x)}(),i
    else
        Nullable(get(x,str,till,tuple_pos(s),i,tuple_state(s))),tuple_pos(s)
    end
end


text:
export simulate

"""
    sample_locus(locus::Dict, n::Int, ploidy::Signed)
Internal function used by `simulate` to take a `Dict` of alleles => frequencies of a locus and return
`n` number of genotypes (n_alleles = `ploidy`) by using weighted sampling of the
allele-frequency pairs.

**Example**
```julia
d = Dict(133 => 0.125,135 => 0.5625,143 => 0.25,137 => 0.0625)

julia> sample_locus(d, 3, 2)
5-element Array{Tuple{Int16,Int16},1}:
 (133, 135)
 (135, 135)
 (143, 137)

julia> sample_locus(d, 3, 3)
5-element Array{Tuple{Int16,Int16,Int16},1}:
 (135, 135, 133)
 (143, 135, 133)
 (137, 135, 135)
```
"""
function sample_locus(locus::Dict, n::Int, ploidy::Signed)
    isempty(locus) && return fill(missing, n)
    k,v = collect(keys(locus)), collect(values(locus))
    alleles = [sample(k, Weights(v), n) for i in 1:ploidy]
    Tuple.(sort.(eachrow(hcat(alleles...))))
end

"""
    simulate(data::PopData; n::Int = 100)
Simulate `n` number of individuals per population (default: `100`) using per-population
allele frequencies derived from a `PopData` object. Returns a new `PopData` object.

**Example**
```julia
cats = @nanycats;

julia> sims = simulate(cats , n = 100)
PopData{Diploid, 9 Microsatellite Loci}
  Samples: 1700
  Populations: 17

julia> sims.sampleinfo

  1700×5 DataFrame
  Row │ name      population  ploidy
      │ String    String      Int8
──────┼───────────────────────────────
    1 │ sim_1     1                2
    2 │ sim_2     1                2
    3 │ sim_3     1                2
    4 │ sim_4     1                2
    5 │ sim_5     1                2
  ⋮   │    ⋮          ⋮         ⋮
 1697 │ sim_1697  17               2
 1698 │ sim_1698  17               2
 1699 │ sim_1699  17               2
 1700 │ sim_1700  17               2
                                         1691 rows omitted

julia> sims.genodata
15300×4 DataFrame
   Row │ name      population  locus   genotype
       │ String    String      String  Tuple…?
───────┼──────────────────────────────────────────
     1 │ sim_1     1           fca8    (135, 143)
     2 │ sim_1     1           fca23   (136, 146)
     3 │ sim_1     1           fca43   (141, 145)
     4 │ sim_1     1           fca45   (120, 126)
     5 │ sim_1     1           fca77   (156, 156)
   ⋮   │    ⋮          ⋮         ⋮         ⋮
 15297 │ sim_1700  17          fca78   (150, 150)
 15298 │ sim_1700  17          fca90   (197, 197)
 15299 │ sim_1700  17          fca96   (113, 113)
 15300 │ sim_1700  17          fca37   (208, 208)
                                15291 rows omitted
```
"""
function simulate(data::PopData; n::Int = 100)
    data.metadata.ploidy != 2 && error("Simulations do not work on mixed-ploidy data (yet)")
    ploidy = data.metadata.ploidy
    pops = unique(data.sampleinfo.population)
    npops = data.metadata.populations
    nloci = data.metadata.loci

    # instantiate output df
    simnames = repeat(["sim_" * "$i" for i in 1:(n*npops)], inner = nloci)
    popnames = repeat(pops, inner = (nloci * n))
    locinames = repeat(unique(data.genodata.locus), outer = (n * npops))
    geno_out = DataFrame(:name => simnames, :population => popnames, :locus => locinames, :genotype => similar(data.genodata.genotype, length(locinames)))

    # generate allele freqs per population
    gdf = groupby(data.genodata, [:population, :locus])
    freqs = DataFrames.combine(
        gdf,
        :genotype => allele_freq => :frq
    )
    # create new genotypes
    transform!(freqs, :frq => (i -> sample_locus.(i,n,ploidy)) => :frq)

    # populate out df
    out_gdf = groupby(geno_out, :population)
    geno_gdf = groupby(freqs, :population)
    for pop in pops
        out_gdf[(population = pop,)][:,:genotype] = reduce(hcat, geno_gdf[(population = pop,)].frq) |> permutedims |> vec
    end
    transform!(
        geno_out,
        :name => (i -> PooledArray(i, compress = true)) => :name,
        :population => (i -> PooledArray(i, compress = true)) => :population,
        :locus => (i -> PooledArray(i, compress = true)) => :locus,
        :genotype
    )
    PopData(geno_out)
end


text:
"""`lqr(A, B, Q, R)`

Calculate the optimal gain matrix `K` for the state-feedback law `u = K*x` that
minimizes the cost function:

J = integral(x'Qx + u'Ru, 0, inf).

For the continuous time model `dx = Ax + Bu`.

`lqr(sys, Q, R)`

Solve the LQR problem for state-space system `sys`. Works for both discrete
and continuous time systems.

See also `LQG`

Usage example:
```julia
using LinearAlgebra # For identity matrix I
A = [0 1; 0 0]
B = [0;1]
C = [1 0]
sys = ss(A,B,C,0)
Q = I
R = I
L = lqr(sys,Q,R)

u(x,t) = -L*x # Form control law,
t=0:0.1:5
x0 = [1,0]
y, t, x, uout = lsim(sys,u,t,x0=x0)
plot(t,x, lab=["Position" "Velocity"], xlabel="Time [s]")
```
"""
function lqr(A, B, Q, R)
    S = care(A, B, Q, R)
    K = R\B'*S
    return K
end

"""`kalman(A, C, R1, R2)` kalman(sys, R1, R2)`

Calculate the optimal Kalman gain

See also `LQG`

"""
kalman(A, C, R1,R2) = Matrix(lqr(A',C',R1,R2)')

function lqr(sys::StateSpace, Q, R)
    if iscontinuous(sys)
        return lqr(sys.A, sys.B, Q, R)
    else
        return dlqr(sys.A, sys.B, Q, R)
    end
end

function kalman(sys::StateSpace, R1,R2)
    if iscontinuous(sys)
        return Matrix(lqr(sys.A', sys.C', R1,R2)')
    else
        return Matrix(dlqr(sys.A', sys.C', R1,R2)')
    end
end

"""`dlqr(A, B, Q, R)`, `dlqr(sys, Q, R)`

Calculate the optimal gain matrix `K` for the state-feedback law `u[k] = K*x[k]` that
minimizes the cost function:

J = sum(x'Qx + u'Ru, 0, inf).

For the discrte time model `x[k+1] = Ax[k] + Bu[k]`.

See also `lqg`

Usage example:
```julia
using LinearAlgebra # For identity matrix I
h = 0.1
A = [1 h; 0 1]
B = [0;1]
C = [1 0]
sys = ss(A,B,C,0, h)
Q = I
R = I
L = dlqr(A,B,Q,R) # lqr(sys,Q,R) can also be used

u(x,t) = -L*x # Form control law,
t=0:h:5
x0 = [1,0]
y, t, x, uout = lsim(sys,u,t,x0=x0)
plot(t,x, lab=["Position"  "Velocity"], xlabel="Time [s]")
```
"""
function dlqr(A, B, Q, R)
    S = dare(A, B, Q, R)
    K = (B'*S*B + R)\(B'S*A)
    return K
end

"""`dkalman(A, C, R1, R2)` kalman(sys, R1, R2)`

Calculate the optimal Kalman gain for discrete time systems

"""
dkalman(A, C, R1,R2) = Matrix(dlqr(A',C',R1,R2)')

"""`place(A, B, p)`, `place(sys::StateSpace, p)`

Calculate gain matrix `K` such that
the poles of `(A-BK)` in are in `p`.

Uses Ackermann's formula."""
function place(A, B, p)
    n = length(p)
    n != size(A,1) && error("Must define as many poles as states")
    n != size(B,1) && error("A and B must have same number of rows")
    if size(B,2) == 1
        acker(A,B,p)
    else
        error("place only implemented for SISO systems")
    end
end

function place(sys::StateSpace, p)
    return place(sys.A, sys.B, p)
end

#Implements Ackermann's formula for placing poles of (A-BK) in p
function acker(A,B,P)
    n = length(P)
    #Calculate characteristic polynomial
    poly = mapreduce(p -> Polynomial([1, -p]), *, P, init=Polynomial(one(eltype(P))))
    q = zero(Array{promote_type(eltype(A),Float64),2}(undef, n,n))
    for i = n:-1:0
        q += A^(n-i)*poly[i]
    end
    S = Array{promote_type(eltype(A),eltype(B),Float64),2}(undef, n,n)
    for i = 0:(n-1)
        S[:,i+1] = A^i*B
    end
    return [zeros(1,n-1) 1]*(S\q)
end


text:
"""
    MaxwellOutput
Return type for Maxwell computations.
Fields are `nothing` on return unless requested.
See individual FMM/direct computation function
documentation for specifics.
"""
mutable struct MaxwellOutput <: FMMVals
    e
    h
    grade
    gradh

    etarg
    htarg
    gradetarg
    gradhtarg

    ifehtarg
    ifeh
end
function MaxwellOutput(ifehtarg,ifeh=nothing)
    MaxwellOutput(nothing,nothing,nothing,nothing,
                  nothing,nothing,nothing,nothing,
                  ifehtarg, ifeh)
end


text:
function apply_laser_over_atoms(laser::Laser{PlaneWave3D}, atoms)
    @debug "start: apply_laser_over_atoms - Laser{PlaneWave3D}"

    s, Δ = laser.s, laser.Δ

    v_atoms = view(Float64.(atoms.r), :, :)
    direction = view(laser.pump.direction, :)
    E = zeros(Complex{eltype(atoms.r)}, atoms.N)
    E₀ = estimate_E₀(laser)

    Threads.@threads for n in 1:(atoms.N)
        oneAtom = v_atoms[:,n]
        @inbounds E[n] = _core_PlaneWave3D(oneAtom, E₀, direction)
    end

    @debug "end  : apply_laser_over_atoms - Laser{PlaneWave3D}"
    return E
end

function apply_laser_over_oneAtom(laser::Laser{PlaneWave3D}, oneAtom::AbstractArray)
    E₀ = estimate_E₀(laser)
    direction = laser.pump.direction
    E = _core_PlaneWave3D(oneAtom, E₀, direction)

    return E
end

"""
    E₀ * exp(+im * k₀ *dot(oneAtom, direction))
"""
function _core_PlaneWave3D(oneAtom, E₀, direction)
    E₀ * cis(+k₀ * (oneAtom[1]*direction[1] + oneAtom[2]*direction[2] + oneAtom[3]*direction[3]))
end


text:
using LinAlgTools
using Test

@testset "LinAlgTools.jl" begin
  @testset "In-place calls" begin
    A = [(i+j)//j for i in 1:3, j in 1:4]
    rref!(A)
    @test A == [1 0 -1//3 -1//2; 0 1 4//3 3//2; 0 0 0 0]

    newrows = row_swap!(A, 1, 3)
    @test A == [0 0 0 0; 0 1 4//3 3//2; 1 0 -1//3 -1//2]
    @test newrows == [0 0 0 0; 1 0 -1//3 -1//2]
    @test !Base.mightalias(newrows, A)

    mulrow = row_mul!(A, 3 => 10)
    @test A == [0 0 0 0; 0 1 4//3 3//2; 10 0 -10//3 -5]
    @test mulrow == [10, 0, -10//3, -5]
    @test !Base.mightalias(mulrow, A)

    addrow = row_add!(A, 3 => 10 => 2)
    @test A == [0 0 0 0; 100 1 -32 -97//2; 10 0 -10//3 -5]
    @test addrow == [100, 1, -32, -97//2]
    @test !Base.mightalias(addrow, A)
  end

  @testset "Non-mutating calls" begin
    A = [(i+j)//j for i in 1:3, j in 1:4]
    Aref = ref(A)
    @test Aref ==  [2 3//2 4//3 5//4; 0 -1//4 -1//3 -3//8; 0 0 0 0]
    @test !Base.mightalias(Aref, A)

    Arref = rref(A)
    @test Arref == [1 0 -1//3 -1//2; 0 1 4//3 3//2; 0 0 0 0]
    @test !Base.mightalias(Arref, A)

    Aswap = row_swap(A, 1, 3)
    @test Aswap == [4 5//2 2 7//4; 3 2 5//3 3//2; 2 3//2 4//3 5//4]
    @test !Base.mightalias(Aswap, A)

    Amul = row_mul(A, 3 => 10)
    @test Amul == [2 3//2 4//3 5//4; 3 2 5//3 3//2; 40 25 20 35//2]
    @test !Base.mightalias(Amul, A)

    Aadd = row_add(A, 3 => 10 => 2)
    @test Aadd == [2 3//2 4//3 5//4; 43 27 65//3 19; 4 5//2 2 7//4]
    @test !Base.mightalias(Aadd, A)
  end
end


text:
using DiffEqBase, Catalyst, Random, Test
using ModelingToolkit: get_states, get_ps

using StableRNGs
rng = StableRNG(12345)

### Tests various cutom made functions ###
new_hill(x, v, k, n) = v*x^n/(k^n+x^n)
new_poly(x,p1,p2) = .5*p1*x^2
new_exp(x,p) = exp(-p*x)

custom_function_network_1 = @reaction_network begin
    hill(X1,v1,K1,2), X1 + Y1--> Z1
    mm(X2,v2,K2), X2 + Y2 --> Z2
    p1*X3^2+p2, X3 + Y3 --> Z3
    exp(-p3*Y4), X4 + Y4 --> Z4
    hillr(X5,v3,K3,2), X5 + Y5 --> Z5
    mmr(X6,v4,K4), X6 + Y6 --> Z6
    hillar(X7,Y7,v5,K5,2), X7 + Y7 --> Z7
end v1 K1 v2 K2 p1 p2 p3 v3 K3 v4 K4 v5 K5

custom_function_network_2 = @reaction_network begin
    new_hill(X1,v1,K1,2), X1 + Y1 --> Z1
    v2*X2/(X2+K2), X2 + Y2 --> Z2
    2*new_poly(X3,p1,p2)+p2, X3 + Y3 --> Z3
    new_exp(Y4,p3), X4 + Y4--> Z4
    v3*(K3^2)/(K3^2+X5^2), X5 + Y5 --> Z5
    v4*K4/(X6+K4), X6 + Y6 --> Z6
    v5*(X7^2)/(K5^2+X7^2+Y7^2), X7 + Y7 --> Z7
end v1 K1 v2 K2 p1 p2 p3 v3 K3 v4 K4 v5 K5

function permute_ps(pvals, rn1, rn2)
    ps1    = parameters(rn1)
    ps2    = parameters(rn2)
    pvals2 = similar(pvals)
    for (i,p) in enumerate(ps2)
        pidx      = findfirst(isequal(p), ps1)
        pvals2[i] = pvals[pidx]
    end
    pvals2
end

f1 = ODEFunction(convert(ODESystem,custom_function_network_1),jac=true)
f2 = ODEFunction(convert(ODESystem,custom_function_network_2),jac=true)
g1 = SDEFunction(convert(SDESystem,custom_function_network_1))
g2 = SDEFunction(convert(SDESystem,custom_function_network_2))
for factor in [1e-2, 1e-1, 1e0, 1e1, 1e2]
    u0 = factor*rand(rng,length(get_states(custom_function_network_1)))
    p = factor*rand(rng,length(get_ps(custom_function_network_2)))

    # needed as this code assumes an ordering of the parameters and species...
    p2 = permute_ps(p, custom_function_network_1, custom_function_network_2)

    t = rand(rng)
    @test all(abs.(f1(u0,p,t) .- f2(u0,p2,t)) .< 10e-10)
    @test all(abs.(f1.jac(u0,p,t) .- f2.jac(u0,p2,t)) .< 10e-10)
    @test all(abs.(g1(u0,p,t) .- g2(u0,p2,t)) .< 10e-10)
end

### Tests that the various notations gives identical results ###

# Michaelis-Menten function.
mm_network = @reaction_network begin
    (1.,1.), 0 ↔ X
    mm(X,v,K), 0 --> X1
    mm(X,v,K), 0 --> X2
    mm(X,v,K), 0 --> X3
end v K
f_mm = ODEFunction(convert(ODESystem,mm_network),jac=true)

u0 = 10*rand(rng,length(get_states(mm_network)))
p = 10*rand(rng,length(get_ps(mm_network)))
t = 10*rand(rng)

f_mm_output = f_mm(u0,p,t)[2:end]
f_mm_jac_output = f_mm.jac(u0,p,t)[2:end,1]
@test (maximum(f_mm_output) - minimum(f_mm_output)) .< 100*eps()
@test (maximum(f_mm_jac_output) - minimum(f_mm_jac_output)) .< 100*eps()

# Repressing Michaelis-Menten function.
mmr_network = @reaction_network begin
    (1.,1.), 0 ↔ X
    mmr(X,v,K), 0 --> X1
    mmr(X,v,K), 0 --> X2
    mmr(X,v,K), 0 --> X3
end v K
f_mmr = ODEFunction(convert(ODESystem,mmr_network),jac=true)

u0 = 10*rand(rng,length(get_states(mmr_network)))
p = 10*rand(rng,length(get_ps(mmr_network)))
t = 10*rand(rng)

f_mmr_output = f_mmr(u0,p,t)[2:end]
f_mmr_jac_output = f_mmr.jac(u0,p,t)[2:end,1]
@test (maximum(f_mmr_output) - minimum(f_mmr_output)) .< 100*eps()
@test (maximum(f_mmr_jac_output) - minimum(f_mmr_jac_output)) .< 100*eps()

# Hill function.
hill_network = @reaction_network begin
    (1.,1.), 0 ↔ X
    hill(X,v,K,2), 0 --> X1
    hill(X,v,K,2), 0 --> X2
end v K
f_hill = ODEFunction(convert(ODESystem,hill_network),jac=true)

u0 = 10*rand(rng,length(get_states(hill_network)))
p = 10*rand(rng,length(get_ps(hill_network)))
t = 10*rand(rng)

f_hill_output = f_hill(u0,p,t)[2:end]
f_hill_jac_output = f_hill.jac(u0,p,t)[2:end,1]
@test (maximum(f_hill_output) - minimum(f_hill_output)) .< 100*eps()
@test (maximum(f_hill_jac_output) - minimum(f_hill_jac_output)) .< 100*eps()

# Repressing Hill function.
hillr_network = @reaction_network begin
    (1.,1.), 0 ↔ X
    hillr(X,v,K,2), 0 --> X1
    hillr(X,v,K,2), 0 --> X2
end v K
f_hillr = ODEFunction(convert(ODESystem,hillr_network),jac=true)

u0 = 10*rand(rng,length(get_states(hillr_network)))
p = 10*rand(rng,length(get_ps(hillr_network)))
t = 10*rand(rng)

f_hillr_output = f_hillr(u0,p,t)[2:end]
f_hillr_jac_output = f_hillr.jac(u0,p,t)[2:end,1]
@test (maximum(f_hillr_output) - minimum(f_hillr_output)) .< 100*eps()
@test (maximum(f_hillr_jac_output) - minimum(f_hillr_jac_output)) .< 100*eps()

# Activation/repressing Hill function.
hillar_network = @reaction_network begin
    (1.,1.), 0 ↔ (X,Y)
    hillar(X,Y,v,K,2), 0 --> X1
    hillar(X,Y,v,K,2), 0 --> X2
    hillar(X,Y,v,K,2), 0 --> X3
    hillar(X,Y,v,K,2), 0 --> X4
end v K
f_hillar = ODEFunction(convert(ODESystem,hillar_network),jac=true)

u0 = 10*rand(rng,length(get_states(hillar_network)))
p = 10*rand(rng,length(get_ps(hillar_network)))
t = 10*rand(rng)

f_hillar_output = f_hillar(u0,p,t)[3:end]
f_hillar_jac_output = f_hillar.jac(u0,p,t)[3:end,1]
@test (maximum(f_hillar_output) - minimum(f_hillar_output)) .< 100*eps()
@test (maximum(f_hillar_jac_output) - minimum(f_hillar_jac_output)) .< 100*eps()


text:
#
module Spectral

using Reexport

@reexport using SciMLBase

using LinearAlgebra
using LinearSolve

import UnPack: @unpack
import Setfield: @set!
import Base.ReshapedArray
import SciMLBase: AbstractDiffEqOperator
import Lazy: @forward
import SparseArrays: sparse
import NNlib: gather, scatter
import FastGaussQuadrature: gausslobatto, gausslegendre, gausschebyshev
import FFTW: plan_rfft, plan_irfft

# AbstractVector subtyping
import Base: summary, show, similar, zero
import Base: size, getindex, setindex!, IndexStyle
import Base.Broadcast: BroadcastStyle

# overload maths
import Base: +, -, *, /, \, adjoint, ∘, inv, one
import LinearAlgebra: mul!, ldiv!, lmul!, rmul!

###
# Abstract Supertypes
###

""" Scalar function field in D-Dimensional space """
abstract type AbstractField{T,D} <: AbstractVector{T} end
""" Operators acting on fields in D-Dimensional space """
abstract type AbstractOperator{T,D} <: AbstractDiffEqOperator{T} end
""" D-Dimensional physical domain """
abstract type AbstractDomain{T,D} end
""" Function space in D-Dimensional space """
abstract type AbstractSpace{T,D} end
""" Boundary condition on domain in D-Dimensional space """
abstract type AbstractBonudaryCondition{T,D} end

AbstractSupertypes{T,D} = Union{
                                AbstractField{T,D},
                                AbstractOperator{T,D},
                                AbstractSpace{T,D},
                                AbstractDomain{T,D},
                                AbstractBonudaryCondition{T,D}
                               }

Base.eltype(::AbstractSupertypes{T,D}) where{T,D} = T
dims(::AbstractSupertypes{T,D}) where{T,D} = D

###
# AbstractOperator subtypes
###

""" Tensor product operator in D-Dimensional space """
abstract type AbstractTensorProductOperator{T,D} <: AbstractOperator{T,D} end

""" Gather-Scatter operator in D-Dimensional space """
abstract type AbstractGatherScatterOperator{T,D} <: AbstractOperator{T,D} end

###
# AbstractSpace subtypes
###

""" D-Dimensional spectral space """
abstract type AbstractSpectralSpace{T,D} <: AbstractSpace{T,D} end

""" D-Dimensional tensor-product space """
abstract type AbstractTensorProductSpace{T,D} <: AbstractSpace{T,D} end

###
# other abstract types
###

""" Boundary Condition on D-Dimensional domain """
abstract type AbstractBoundaryCondition{T,D} end

#""" Deform D-Dimensional domain """
#abstract type AbstractDeformation{D} end

# misc
include("utils.jl")
include("NDgrid.jl")

# field
include("Field.jl")

# operators
include("OperatorBasics.jl")
include("Operators.jl")

# domain
include("Domain.jl")

# space
include("Space.jl")
include("DeformSpace.jl")

# polynomial spaces
include("LagrangeMatrices.jl")
include("LagrangePolynomialSpace.jl")
#include("GalerkinOperators.jl")

# fourier spaces
#include("FourierSpace.jl")

# boundary value problem
#include("BoundaryValueProblem.jl")

export
       # Domains
       IntervalDomain, BoxDomain,
       deform, end_points, isperiodic,
       unit_box, annulus_2D,

       # space functionality
       grid, gradOp, massOp, laplaceOp, advectionOp, interpOp,
#      BoundaryCondition,

       # spaces
       LagrangePolynomialSpace,
       GaussLobattoLegendre1D, GaussLegendre1D, GaussChebychev1D,
       GaussLobattoLegendre2D, GaussLegendre2D, GaussChebychev2D

       # boundary value problem
#      BoundaryValuePDEProblem

end # module


text:
# # Histogram of the most significant digit (MSD)

using LiterateTest                                                     #src
using Test                                                             #src

# ## An allocation-free function to compute MSD

function msd(x::Real)
    x = abs(x)
    d = x
    while true
        x < 1 && return floor(Int, d)
        d = x
        x ÷= 10
    end
end
nothing  # hide
#-

@test begin
    msd(34513)
end == 3
#-

@test begin
    msd(-51334)
end == 5
#-

@assert 2.76e19 > typemax(Int64)  #src
@test begin
    msd(2.76e19)
end == 2
#-

# ## MSD of random numbers

using CUDA
using FLoops
using FoldsCUDA
using Setfield

function histogram_msd(xs, ex = nothing)
    zs = ntuple(_ -> 0, 9)  # a tuple of 9 zeros
    @floop ex for x in xs
        d = msd(x)
        1 <= d <= 9 || continue  # skip it if `msd` returns 0
        h2 = @set zs[d] = 1      # set `d`-th position of the tuple to 1
        @reduce(h1 = zs .+ h2)   # point-wise addition merges the histogram
    end
    return h1
end
nothing  # hide

# Generate some random numbers

xs = let randn = has_cuda_gpu() ? CUDA.randn : randn
    exp.(10.0 .* (randn(10^8) .+ 6))
end
@assert all(isfinite, xs)
collect(view(xs, 1:10))  # preview

# Pass an array of (real) numbers to `histogram_msd` to compute the
# histogram of MSD:

hist1 = histogram_msd(xs)

# Frequency in percentage:

aspercentage(h) = pairs(round.((collect(h) ./ length(xs) .* 100); digits = 1))
aspercentage(hist1)
#-

# ## MSD of lazily computed numbers

# FoldsCUDA supports running reduction over non-`CuArray` containers
# such as `UnitRange` and iterator transformation wrapping it.
# However, you need to explicitly specify to use CUDA by, e.g.,
# passing `CUDAEx` to `@floop`:

hist2 = histogram_msd(x^2 for x in 1:10^8)

# Frequency in percentage:
aspercentage(hist2)
#-

hist3 = histogram_msd(exp(x) for x in range(1, 35, length=10^8))

# Frequency in percentage:
aspercentage(hist3)


text:
using LinearAlgebra
using Nystrom
using StaticArrays
using Random
using Plots
plotlyjs()
Random.seed!(1)

## general definitions
xs = SVector(3,3,3)
ops = (
    Laplace(;dim=3),
    Helmholtz(;dim=3,k=1.2),
    Elastostatic(;dim=3,μ=2,λ=3),
    Maxwell(;k=2)
)
pde = ops[1] # chose your pde here
T    = Nystrom.default_density_eltype(pde)
c    = rand(T)

## place to store dofs and error
ndofs = []
ee = []

npatches = [2,4,8]
qorder   =  6
for n in npatches
    ## generate mesh
    Geometry.clear_entities!()
    Ω   = ParametricSurfaces.Sphere() |> Geometry.Domain
    Γ   = boundary(Ω)
    M   = meshgen(Γ,(n,n))
    mesh = NystromMesh(view(M,Γ),order=qorder)
    push!(ndofs,length(Nystrom.dofs(mesh)))
    ## create exact solution
    u    = (qnode) -> SingleLayerKernel(pde)(xs,qnode)*c
    dudn = (qnode) -> transpose(DoubleLayerKernel(pde)(xs,qnode))*c
    γ₀u   = Density(u,mesh)
    γ₁u   = Density(dudn,mesh)
    γ₀u_norm = norm(norm.(γ₀u,Inf),Inf)

    ## compute error
    S,D  = Nystrom.single_doublelayer_dim(pde,mesh)
    e1    = Nystrom.error_interior_green_identity(S,D,γ₀u,γ₁u)/γ₀u_norm
    push!(ee,norm(e1,Inf))
    @info ndofs[end],ee[end]
end

xx = sqrt.(ndofs)
fig = plot(xx,ee,xscale=:log10,yscale=:log10,m=:o,label="error",lc=:black)
plot!(xlabel="√n",ylabel="error",title="$pde")
for p in 1:5
    c = ee[end]*xx[end]^p
    plot!(fig,xx,c ./ xx.^(p),xscale=:log10,yscale=:log10,label="h^$p",ls=:dash)
end
display(fig)

max_stars_repo_name: WaveProp/Nystrom

max_issues_repo_name: WaveProp/Nystrom

max_forks_repo_name: WaveProp/Nystrom


text:
using LinearMaps: LinearMap

"""
    run(parameters_file; verbose=false, refinement_level_increment=0, parameters...)

Run a Trixi simulation with the parameters in `parameters_file`.
Parameters can be overriden by specifying them as keyword arguments (see examples).

If `verbose` is `true`, additional output will be generated on the terminal
that may help with debugging. If a value for `refinement_level_increment` is given,
`initial_refinement_level` will be increased by this value before running the simulation (mostly
used by EOC analysis).

# Examples
```julia
julia> Trixi.run("examples/parameters.toml", verbose=true)
[...]
```

Without changing the parameters file we can start a simulation with `polydeg = 1` and
`t_end = 0.5` as follows:
```julia
julia> Trixi.run("examples/parameters.toml", polydeg=1, t_end=0.5)
[...]
```
"""
function run(parameters_file; verbose=false, refinement_level_increment=0, parameters...)
  # Reset timer
  reset_timer!(timer())

  # Read command line or keyword arguments and parse parameters file
  init_parameters(parameters_file; verbose=verbose,
      refinement_level_increment=refinement_level_increment, parameters...)

  # Separate initialization and execution into two functions such that Julia can specialize
  # the code in `run_simulation` for the actual type of `solver` and `mesh`
  if parameter("equations") == "euler_gravity"
    globals[:euler_gravity] = true
    mesh, solver, time_parameters, time_integration_function = init_simulation_euler_gravity()
    run_simulation_euler_gravity(mesh, solver, time_parameters, time_integration_function)
  else
    globals[:euler_gravity] = false
    mesh, solver, time_parameters, time_integration_function = init_simulation()
    run_simulation(mesh, solver, time_parameters, time_integration_function)
  end
end

function init_parameters(parameters_file=nothing; verbose=false, refinement_level_increment=0, parameters...)
  # Set global verbosity
  globals[:verbose] = verbose

  # Parse parameters file
  @timeit timer() "read parameter file" parse_parameters_file(parameters_file)

  # Override specified parameters
  for (parameter, value) in parameters
    setparameter(string(parameter), value)
  end

  # Start simulation with an increased initial refinement level if specified
  # for convergence analysis
  if refinement_level_increment != 0
    setparameter("initial_refinement_level",
      parameter("initial_refinement_level") + refinement_level_increment)
  end
end

function init_simulation()
  # Print starup message
  print_startup_message()

  # Get number of dimensions
  ndims_ = parameter("ndims")::Int

  # Check if this is a restart from a previous result or a new simulation
  restart = parameter("restart", false)
  if restart
    restart_filename = parameter("restart_filename")
  end

  # Initialize mesh
  if restart
    print("Loading mesh... ")
    @timeit timer() "mesh loading" mesh = load_mesh(restart_filename)
    println("done")
  else
    print("Creating mesh... ")
    @timeit timer() "mesh creation" mesh = generate_mesh()
    mesh.current_filename = save_mesh_file(mesh)
    mesh.unsaved_changes = false
    println("done")
  end

  # Initialize system of equations
  print("Initializing system of equations... ")
  equations_name = parameter("equations")
  equations = make_equations(equations_name, ndims_)
  println("done")

  # Initialize solver
  print("Initializing solver... ")
  solver_name = parameter("solver", valid=["dg"])
  solver = make_solver(solver_name, equations, mesh)
  println("done")

  # Sanity checks
  # If DG volume integral type is weak form, volume flux type must be flux_central,
  # as everything else does not make sense
  if solver.volume_integral_type === Val(:weak_form) && solver.volume_flux_function !== flux_central
    error("using the weak formulation with a volume flux other than 'flux_central' does not make sense")
  end

  if equations isa AbstractIdealGlmMhdEquations && solver.volume_integral_type === Val(:weak_form)
    error("The weak form is not implemented for $equations.")
  end

  # Initialize solution
  amr_interval = parameter("amr_interval", 0)
  adapt_initial_conditions = parameter("adapt_initial_conditions", true)
  adapt_initial_conditions_only_refine = parameter("adapt_initial_conditions_only_refine", true)
  if restart
    print("Loading restart file...")
    time, step = load_restart_file!(solver, restart_filename)
    println("done")
  else
    print("Applying initial conditions... ")
    t_start = parameter("t_start")
    time = t_start
    step = 0
    set_initial_conditions!(solver, time)
    println("done")

    # If AMR is enabled, adapt mesh and re-apply ICs
    if amr_interval > 0 && adapt_initial_conditions
      @timeit timer() "initial condition AMR" has_changed = adapt!(mesh, solver, time,
          only_refine=adapt_initial_conditions_only_refine)

      # Iterate until mesh does not change anymore
      while has_changed
        set_initial_conditions!(solver, time)
        @timeit timer() "initial condition AMR" has_changed = adapt!(mesh, solver, time,
            only_refine=adapt_initial_conditions_only_refine)
      end

      # Save mesh file
      mesh.current_filename = save_mesh_file(mesh)
      mesh.unsaved_changes = false
    end
  end
  t_end = parameter("t_end")

  # Init time integration
  time_integration_scheme = Symbol(parameter("time_integration_scheme", "timestep_carpenter_kennedy_erk54_2N!"))
  time_integration_function = eval(time_integration_scheme)

  # Print setup information
  solution_interval = parameter("solution_interval", 0)
  restart_interval = parameter("restart_interval", 0)
  polydeg = parameter("polydeg") # FIXME: This is currently the only DG-specific code in here
  n_steps_max = parameter("n_steps_max")
  cfl = parameter("cfl")
  sources = parameter("sources", "none")
  n_leaf_cells = count_leaf_cells(mesh.tree)
  min_level = minimum_level(mesh.tree)
  max_level = maximum_level(mesh.tree)
  domain_center = mesh.tree.center_level_0
  domain_length = mesh.tree.length_level_0
  min_dx = domain_length / 2^max_level
  max_dx = domain_length / 2^min_level
  s = ""
  s *= """| Simulation setup
          | ----------------
          | working directory:  $(pwd())
          | parameters file:    $(parameter("parameters_file"))
          | equations:          $(get_name(equations))
          | | #variables:       $(nvariables(equations))
          | | variable names:   $(join(varnames_cons(equations), ", "))
          | sources:            $sources
          | restart:            $(restart ? "yes" : "no")
          """
  if restart
    s *= "| | restart timestep: $step\n"
    s *= "| | restart time:     $time\n"
  else
    s *= "| initial conditions: $(get_name(solver.initial_conditions))\n"
    s *= "| t_start:            $t_start\n"
  end
  s *= """| t_end:              $t_end
          | AMR:                $(amr_interval > 0 ? "yes" : "no")
          """
  if amr_interval > 0
    s *= "| | AMR interval:     $amr_interval\n"
    s *= "| | adapt ICs:        $(adapt_initial_conditions ? "yes" : "no")\n"
  end
  s *= """| n_steps_max:        $n_steps_max
          | time integration:   $(get_name(time_integration_function))
          | restart interval:   $restart_interval
          | solution interval:  $solution_interval
          | #parallel threads:  $(Threads.nthreads())
          |
          | Solver
          | | solver:           $solver_name
          | | polydeg:          $polydeg
          | | CFL:              $cfl
          | | volume integral:  $(get_name(solver.volume_integral_type))
          | | volume flux:      $(get_name(solver.volume_flux_function))
          | | surface flux:     $(get_name(solver.surface_flux_function))
          | | #elements:        $(solver.n_elements)
          | | #interfaces:      $(solver.n_interfaces)
          | | #boundaries:      $(solver.n_boundaries)
          | | #l2mortars:       $(solver.n_l2mortars)
          | | #DOFs:            $(ndofs(solver))
          |
          | Mesh
          | | #cells:           $(length(mesh.tree))
          | | #leaf cells:      $n_leaf_cells
          | | minimum level:    $min_level
          | | maximum level:    $max_level
          | | domain center:    $(join(domain_center, ", "))
          | | domain length:    $domain_length
          | | minimum dx:       $min_dx
          | | maximum dx:       $max_dx
          """
  println()
  println(s)

  # Set up main loop
  save_final_solution = parameter("save_final_solution", true)
  save_final_restart = parameter("save_final_restart", true)
  analysis_interval = parameter("analysis_interval", 0)
  if analysis_interval > 0
    alive_interval = parameter("alive_interval", div(analysis_interval, 10))
  else
    alive_interval = 0
  end

  # Save initial conditions if desired
  if !restart && parameter("save_initial_solution", true)
    # we need to make sure, that derived quantities, such as e.g. blending
    # factor is already computed for the initial condition
    @notimeit timer() rhs!(solver, time)
    save_solution_file(solver, mesh, time, 0, step)
  end

  # Print initial solution analysis and initialize solution analysis
  if analysis_interval > 0
    analyze_solution(solver, mesh, time, 0, step, 0, 0)
  end

  time_parameters = (time=time, step=step, t_end=t_end, cfl=cfl,
                    n_steps_max=n_steps_max,
                    save_final_solution=save_final_solution,
                    save_final_restart=save_final_restart,
                    analysis_interval=analysis_interval,
                    alive_interval=alive_interval,
                    solution_interval=solution_interval,
                    amr_interval=amr_interval,
                    restart_interval=restart_interval)
  return mesh, solver, time_parameters, time_integration_function
end

function run_simulation(mesh, solver, time_parameters, time_integration_function)
  @unpack time, step, t_end, cfl, n_steps_max,
          save_final_solution, save_final_restart,
          analysis_interval, alive_interval,
          solution_interval, amr_interval,
          restart_interval = time_parameters

  loop_start_time = time_ns()
  analysis_start_time = time_ns()
  output_time = 0.0
  n_analysis_timesteps = 0

  # Declare variables to return error norms calculated in main loop
  # (needed such that the values are accessible outside of the main loop)
  local l2_error, linf_error

  # Start main loop (loop until final time step is reached)
  finalstep = false
  first_loop_iteration = true
  @timeit timer() "main loop" while !finalstep
    # Calculate time step size
    @timeit timer() "calc_dt" dt = calc_dt(solver, cfl)

    # Abort if time step size is NaN
    if isnan(dt)
      error("time step size `dt` is NaN")
    end

    # If the next iteration would push the simulation beyond the end time, set dt accordingly
    if time + dt > t_end
      dt = t_end - time
      finalstep = true
    end

    # Evolve solution by one time step
    time_integration_function(solver, time, dt)
    step += 1
    time += dt
    n_analysis_timesteps += 1

    # Check if we reached the maximum number of time steps
    if step == n_steps_max
      finalstep = true
    end

    # Check steady-state integration residual
    if solver.equations isa AbstractHyperbolicDiffusionEquations
      if maximum(abs, view(solver.elements.u_t, 1, .., :)) <= solver.equations.resid_tol
        println()
        println("-"^80)
        println("  Steady state tolerance of ",solver.equations.resid_tol," reached at time ",time)
        println("-"^80)
        println()
        finalstep = true
      end
    end

    # Analyze solution errors
    if analysis_interval > 0 && (step % analysis_interval == 0 || finalstep)
      # Calculate absolute and relative runtime
      runtime_absolute = (time_ns() - loop_start_time) / 10^9
      runtime_relative = ((time_ns() - analysis_start_time - output_time) / 10^9 /
                          (n_analysis_timesteps * ndofs(solver)))

      # Analyze solution
      l2_error, linf_error = @timeit timer() "analyze solution" analyze_solution(
          solver, mesh, time, dt, step, runtime_absolute, runtime_relative)

      # Reset time and counters
      analysis_start_time = time_ns()
      output_time = 0.0
      n_analysis_timesteps = 0
      if finalstep
        println("-"^80)
        println("Trixi simulation run finished.    Final time: $time    Time steps: $step")
        println("-"^80)
        println()
      end
    elseif alive_interval > 0 && step % alive_interval == 0
      runtime_absolute = (time_ns() - loop_start_time) / 10^9
      @printf("#t/s: %6d | dt: %.4e | Sim. time: %.4e | Run time: %.4e s\n",
              step, dt, time, runtime_absolute)
    end

    # Write solution file
    if solution_interval > 0 && (
        step % solution_interval == 0 || (finalstep && save_final_solution))
      output_start_time = time_ns()
      @timeit timer() "I/O" begin
        # Compute current AMR indicator values such that it can be written to
        # the solution file for the current number of elements
        if amr_interval > 0
          calc_amr_indicator(solver, mesh, time)
        end

        # If mesh has changed, write a new mesh file name
        if mesh.unsaved_changes
          mesh.current_filename = save_mesh_file(mesh, step)
          mesh.unsaved_changes = false
        end

        # Then write solution file
        save_solution_file(solver, mesh, time, dt, step)
      end
      output_time += time_ns() - output_start_time
    end

    # Write restart file
    if restart_interval > 0 && (
        step % restart_interval == 0 || (finalstep && save_final_restart))
      output_start_time = time_ns()
      @timeit timer() "I/O" begin
        # If mesh has changed, write a new mesh file
        if mesh.unsaved_changes
          mesh.current_filename = save_mesh_file(mesh, step)
          mesh.unsaved_changes = false
        end

        # Then write restart file
        save_restart_file(solver, mesh, time, dt, step)
      end
      output_time += time_ns() - output_start_time
    end

    # Perform adaptive mesh refinement
    if amr_interval > 0 && (step % amr_interval == 0) && !finalstep
      @timeit timer() "AMR" has_changed = adapt!(mesh, solver, time)

      # Store if mesh has changed to write changed mesh file before next restart/solution output
      if has_changed
        mesh.unsaved_changes = true
      end
    end

    # The following call ensures that when doing memory allocation
    # measurements, the memory allocations for JIT compilation are discarded
    # (since virtually all relevant methods have already been called by now)
    if first_loop_iteration
      clear_malloc_data()
      first_loop_iteration = false
    end
  end

  # Print timer information
  print_timer(timer(), title="Trixi.jl", allocations=true, linechars=:ascii, compact=false)
  println()

  # Return error norms for EOC calculation
  return l2_error, linf_error, varnames_cons(solver.equations)
end

"""
    convtest(parameters_file, iterations; parameters...)

Run multiple Trixi simulations with the parameters in `parameters_file` and compute
the experimental order of convergence (EOC) in the ``L^2`` and ``L^\\infty`` norm.
The number of runs is specified by `iterations` and in each run the initial
refinement level will be increased by 1. Parameters can be overriden by specifying them as
additional keyword arguments, which are passed to the respective call to `run`..
"""
function convtest(parameters_file, iterations; parameters...)
  @assert(iterations > 1, "Number of iterations must be bigger than 1 for a convergence analysis")

  # Types of errors to be calcuated
  errors = Dict(:L2 => Float64[], :Linf => Float64[])

  # Declare variable to access variable names after for loop
  local variablenames

  # Run trixi and extract errors
  for i = 1:iterations
    println(string("Running convtest iteration ", i, "/", iterations))
    l2_error, linf_error, variablenames = run(parameters_file; refinement_level_increment = i - 1,
                                              parameters...)

    # Collect errors as one vector to reshape later
    append!(errors[:L2], l2_error)
    append!(errors[:Linf], linf_error)
  end

  # Number of variables
  nvariables = length(variablenames)

  # Reshape errors to get a matrix where the i-th row represents the i-th iteration
  # and the j-th column represents the j-th variable
  errorsmatrix = Dict(kind => transpose(reshape(error, (nvariables, iterations))) for (kind, error) in errors)

  # Calculate EOCs where the columns represent the variables
  # As dx halves in every iteration the denominator needs to be log(1/2)
  eocs = Dict(kind => log.(error[2:end, :] ./ error[1:end-1, :]) ./ log(1 / 2) for (kind, error) in errorsmatrix)

  for (kind, error) in errorsmatrix
    println(kind)

    for v in variablenames
      @printf("%-20s", v)
    end
    println("")

    for k = 1:nvariables
      @printf("%-10s", "error")
      @printf("%-10s", "EOC")
    end
    println("")

    # Print errors for the first iteration
    for k = 1:nvariables
      @printf("%-10.2e", error[1, k])
      @printf("%-10s", "-")
    end
    println("")

    # For the following iterations print errors and EOCs
    for j = 2:iterations
      for k = 1:nvariables
        @printf("%-10.2e", error[j, k])
        @printf("%-10.2f", eocs[kind][j-1, k])
      end
      println("")
    end
    println("")

    # Print mean EOCs
    for k = 1:nvariables
      @printf("%-10s", "mean")
      @printf("%-10.2f", sum(eocs[kind][:, k]) ./ length(eocs[kind][:, k]))
    end
    println("")
    println("-"^80)
  end
end

"""
    compute_linear_structure(parameters_file, source_terms=nothing; verbose=false, parameters...)

Computes the exact Jacobian `A` of a linear DG operator wrapped as a `LinearMap` and the right hand side `b`.
Returns `A, b`.
"""
function compute_linear_structure(parameters_file, source_terms=nothing; verbose=false, parameters...)
  # Read command line or keyword arguments and parse parameters file
  init_parameters(parameters_file; verbose=verbose, parameters...)
  globals[:euler_gravity] = false
  mesh, solver, time_parameters = init_simulation()

  equations(solver) isa Union{AbstractLinearScalarAdvectionEquation, AbstractHyperbolicDiffusionEquations} ||
    throw(ArgumentError("Only linear problems are supported."))

  # get the right hand side from the source terms
  solver.elements.u .= 0
  rhs!(solver, 0)
  b = vec(-solver.elements.u_t) |> copy

  # set the source terms to zero to extract the linear operator
  if solver isa Dg2D
    solver = Dg2D(solver.equations, solver.surface_flux_function, solver.volume_flux_function, solver.initial_conditions,
                  source_terms, mesh, polydeg(solver))
  elseif solver isa Dg3D
    solver = Dg3D(solver.equations, solver.surface_flux_function, solver.volume_flux_function, solver.initial_conditions,
                  source_terms, mesh, polydeg(solver))
  else
    error("not implemented")
  end
  A = LinearMap(length(solver.elements.u), ismutating=true) do dest,src
    vec(solver.elements.u) .= src
    rhs!(solver, 0)
    dest .= vec(solver.elements.u_t)
  end

  A, b
end

"""
    compute_jacobian_dg(parameters_file=nothing; verbose=false, parameters...)

Uses DG right hand side operator and simple second order finite difference to compute the Jacobian `J` of the operator.
The linearisation state is the initial condition from the parameter file.
Returns `J`.
"""
function compute_jacobian_dg(parameters_file; verbose=false, parameters...)
  # Read command line or keyword arguments and parse parameters file
  init_parameters(parameters_file; verbose=verbose, parameters...)
  # function does not support multi-physics
  if parameter("equations") == "euler_gravity"
    throw(ArgumentError("Multi-physics such as Euler-gravity is not supported"))
  end
  globals[:euler_gravity] = false

  # linearisation state is initial condition
  mesh, dg, time_parameters = init_simulation()
  # store initial state
  u0 = dg.elements.u |> copy

  #compute residual of linearisation state
  rhs!(dg, 0)
  res0 = vec(dg.elements.u_t) |> copy

  # initialize Jacobian matrix
  J = zeros(length(dg.elements.u),length(dg.elements.u))

  # use second order finite difference to estimate Jacobian matrix
  for idx in eachindex(dg.elements.u)
    # determine size of fluctuation
    epsilon = sqrt(eps(u0[idx]))
    # plus fluctuation
    dg.elements.u[idx] = u0[idx] + epsilon
    rhs!(dg, 0)
    # stores the right hand side with plus epsilon fluctuation
    res_p = vec(dg.elements.u_t) |> copy
    # minus fluctuation
    dg.elements.u[idx] = u0[idx] - epsilon
    rhs!(dg, 0)
    # stores the right hand side with minus epsilon fluctuation
    res_m = vec(dg.elements.u_t) |> copy
    # restore linearisation state
    dg.elements.u[idx] = u0[idx]
    # central second order finite difference
    J[:,idx] = (res_p - res_m) / (2 * epsilon)
  end

  return J
end

# Include source file with init and run methods for coupled Euler-gravity simulations
include("run_euler_gravity.jl")


text:
using CurveProximityQueries
using IntervalArithmetic
using Plots

pgfplots()

include("SymbolicCurves.jl")
using .SymbolicCurves

function heartcurve(a::Real, θ, xcenter, ycenter)
    x = :(($a*0.02*16)*(sin(t))^3)
    y = :(($a*0.02)*(13*cos(t) - 5*cos(t*2) - 2*cos(t*3) - cos(t*4)))
    SymCurve(x, y, 0..2π, θ, xcenter, ycenter)
end

function circinvolute(a::Real, θ, xcenter, ycenter)
    x = :($a*(cos(t) + t*sin(t)))
    y = :($a*(sin(t) - t*cos(t)))
    SymCurve(x, y, 0..10, θ, xcenter, ycenter)
end

a = heartcurve(1.0, 0.0, 0.0, 0.0)
b = circinvolute(0.035, -π, -1.0, 0.0)
pts = closest_points(a,b)
plot(a,
     linewidth=2.0,
     linecolor=:black,
    )
plot!(b,
      linewidth=2.0,
      linecolor=:black,
     )
plot!(pts,
      linestyle=:dash,
      markersize=8.0,
      linewidth=2.0,
     )
#  savefig("out/fig1b.pdf")


text:
module TestUtils

using Bootstrap
using Test
using Random

@testset "Statistic functions" begin

    x = collect(100:-1:0) ## avoid sorting

    ## reference values, computed with boot:::norm_inter (R)
    qr = (
        (1.0,      Inf),
        (0.0,     -Inf),
        (0.5,     50.0),
        (0.4713,  47.07265),
        (0.1,      9.205641),
        (0.052,    4.319606)
    )

    @testset "Normal-interpolated quantiles" begin

        @testset "Reference values: alpha scalar" begin
            for (alpha, ref) in qr
                @test isapprox( Bootstrap.iquantile(x, alpha), ref, atol = 1e-5 )
            end
        end

        @testset "Reference values: alpha vector" begin
            a = Float64[x[1] for x in qr]
            r = Float64[x[2] for x in qr]
            y = Bootstrap.iquantile(x, a)
            @test isapprox(y, r, atol = 1e-4 )
        end

        @testset "Monotone quantiles" begin
            x = collect(0:100)
            ## check that quantiles are monotone
            y0 = -Inf
            for alpha in 0:0.001:1
                y1 = Bootstrap.iquantile(x, alpha)
                @test y1 >= y0
                y0 = y1
            end
        end

    end

end

@testset "Utility functions" begin

    @testset "Exact sampling" begin

        mss = collect(Bootstrap.exact(4))
        @test length(mss) == 35
        @test reduce(&, [issorted(c) for c in mss])

        for n in 1:6
            @test length(Bootstrap.exact(n)) == Bootstrap.nrun_exact(n)
        end

    end

    @testset "draw! for views" begin

        x = rand(100)
        xv = view(x, 1:2:100)
        xvc = collect(xv)

        Random.seed!(1)
        d1 = draw!(xv, copy(xv))
        @test d1 isa Vector{Float64}

        Random.seed!(1)
        d2 = draw!(xvc, copy(xvc))
        @test d2 == d1

    end

end

end


text:
# The mesh data structure is based on the paper:
# Efficient Representation of Computational Meshes. A. Logg [1]
# This describes the general data structure as well as the algorithms to compute
# mesh connectivities

using FixedSizeArrays
import GeometryTypes: Point, Simplex

abstract BoundaryCondition
"""
The ambient dimension of the problem. These types don't hold data and only serve
as a flag for function dispatch.
"""
abstract Dimension
immutable Dimension2 <: Dimension end
immutable Dimension3 <: Dimension end

type Dirichlet{T<:Function} <: BoundaryCondition
    val ::T
end

type Neumann{T<:Function} <: BoundaryCondition
    val ::T
end

type MeshBoundaryCondition{T<:BoundaryCondition}
    bcs     ::Array{T}
    markers ::Array{Int64}
end

"""
A MeshEntity is a unique tuple (d,i), where d is the topological dimension of
the entitiy and i is the index for the mesh entity.
"""
type MeshEntity{N,T} <: FixedVector{N,T}
    _::NTuple{N,T}
end
function call{T<:MeshEntity, F<:FixedVector}(::Type{T}, f::F...)
    MeshEntitiy{length(f),F}(f)
end

"""
The mesh connectivity describes how different mesh entities are connected to
each other. For example, it describes which vertices belong to a mesh face or
which edges belong to a mesh vertex. From [1], all the connectivity combinations
can be derived from knowing only the connectivity between the mesh cells and the
mesh vertices.

The notation for the connectivity of dimension d to dimension d' is d -> d'.

The indices array contains the indices for the dimension d' entities that are
connected to the dimension d entities. The kth element in the offsets array
gives the index into the indices array for the entities connected to the kth
entity of dimension d. See [1] for details.
"""
type MeshConnectivity
    indices ::Array{Int64}
    offsets ::Array{Int64}
end

type MeshTopology{D<:Dimension}
    dim          ::D # Ambient dimension of the model
    connectivity ::Dict{AbstractString, MeshConnectivity}
end

type MeshGeometry
    vertices     ::Array{Point}
    nodes        ::Array{Point}
    normals      ::Array{Point} # One normal for every mesh cell
end

type Mesh{D<:Dimension}
    topology ::MeshTopology{D}
    geometry ::MeshGeometry
end


text:
using stoosc
rc("font", family="arial")

MODEL = "model02b"  # model01, model02, model02b

if MODEL == "model01"
    I_strings = ["0.1", "0.25", "0.3"]
    σ_strings = ["0.003", "0.01", "0.03"]
else
    I_strings = ["3.0", "3.5", "5.0"]
    σ_strings = ["0.01", "0.08", "0.2"]
end

fig, axarr = subplots(3, 3, figsize=(6.5, 4.5))
for (i, I_string) in enumerate(I_strings)
    for (j, σ_string) in enumerate(σ_strings)
        println("I = $(I_string), σ = $(σ_string); ")
        phase_shift_arr = []
        t = nothing
        for r = 1:10  # 10
            filename = joinpath("prcs", MODEL, "sde_I=$(I_string)_σ=$(σ_string)_r=$(r).csv")
            df = load_data(filename)
            push!(phase_shift_arr, df[:, "phase_shift"])
            t = df[:, "time"]
        end

        phase_shift = [cmean([2π*x[i] for x in phase_shift_arr])/2π for i in 1:length(phase_shift_arr[1])]
        phase_shift_std = [cstd([2π*x[i] for x in phase_shift_arr])/2π for i in 1:length(phase_shift_arr[1])]
        axarr[i, j].errorbar(t, phase_shift,
            yerr=phase_shift_std, fmt=".", color="black")

        # Draw a grid
        grid_color = fill(0.8, 3)
        axarr[i, j].hlines(0.0, 0, 1, color=grid_color, zorder=-1000)
        axarr[i, j].vlines(0.5, -0.5, 0.5, color=grid_color, zorder=-1000)
        axarr[i, j].set_xlim(0, 1)
        axarr[i, j].set_xticks(0:1.0:1)
        axarr[i, j].set_ylim(-0.5, 0.5)
        axarr[i, j].set_yticks(-0.5:1.0:0.5)
        plot_events([0 0.5], ax=axarr[i, j], ylims=[-0.5, 0.5], zorder=-2000)
        axarr[i, j].set_title("I = $(I_string), σ = $(σ_string)", pad=0.1)
        axarr[i, j].set_xlabel("Pulse time", labelpad=0.0)
        axarr[i, j].set_ylabel("Phase shift", labelpad=0.0)
    end
end

fig.tight_layout(pad=0)
save_figure(fig, "prcs", MODEL, "PRC.svg")

max_stars_repo_name: vkumpost/stoosc

max_issues_repo_name: vkumpost/stoosc

max_forks_repo_name: vkumpost/stoosc


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Tests related to ECEF to ECEF transformations.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

# Get the current EOP Data.
#
# TODO: The EOP data obtained from IERS website does not match the values in the
# examples in [1]. However, this should be enough, because 1) the individual
# functions at the low level are tested using the same values of [1], and 2) the
# difference is smaller than 30 cm.

eop_iau1980  = read_iers_eop("./eop_IAU1980.txt",  :IAU1980)
eop_iau2000a = read_iers_eop("./eop_IAU2000A.txt", :IAU2000A)

# File: ./src/transformations/ecef_to_ecef.jl
# ===========================================

# Functions: rECEFtoECEF
# ----------------------

# ITRF <=> PEF
# ============

################################################################################
#                                  IAU-76/FK5
################################################################################

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_itrf = -1033.4793830    i + 7901.2952754    j + 6380.3565958    k [km]
#   v_itrf =    -3.225636520  i -    2.872451450  j +    5.531924446  k [km/s]
#
# one gets:
#
#   r_pef  = -1033.47503130   i + 7901.30558560   j + 6380.34453270   k [km]
#   v_pef  =    -3.2256327470 i -    2.8724425110 j +    5.5319312880 k [km/s]
#
################################################################################

@testset "Function rECEFtoECEF ITRF <=> PEF" begin
    JD_UTC = DatetoJD(2004, 4, 6, 7, 51, 28.386009)

    ## ITRF => PEF
    ## ===========

    r_itrf = [-1033.4793830; 7901.2952754; 6380.3565958]
    v_itrf = [-3.225636520; -2.872451450; +5.531924446]

    ## DCM
    ## ---

    D_PEF_ITRF = rECEFtoECEF(ITRF(), PEF(), JD_UTC, eop_iau1980)

    r_pef = D_PEF_ITRF*r_itrf
    v_pef = D_PEF_ITRF*v_itrf

    @test r_pef[1] ≈ -1033.47503130 atol=1e-4
    @test r_pef[2] ≈ +7901.30558560 atol=1e-4
    @test r_pef[3] ≈ +6380.34453270 atol=1e-4

    @test v_pef[1] ≈ -3.2256327470  atol=1e-7
    @test v_pef[2] ≈ -2.8724425110  atol=1e-7
    @test v_pef[3] ≈ +5.5319312880  atol=1e-7

    ## Quaternion
    ## ----------

    q_PEF_ITRF = rECEFtoECEF(Quaternion, ITRF(), PEF(), JD_UTC, eop_iau1980)

    r_pef = vect(conj(q_PEF_ITRF)*r_itrf*q_PEF_ITRF)
    v_pef = vect(conj(q_PEF_ITRF)*v_itrf*q_PEF_ITRF)

    @test r_pef[1] ≈ -1033.47503130 atol=1e-4
    @test r_pef[2] ≈ +7901.30558560 atol=1e-4
    @test r_pef[3] ≈ +6380.34453270 atol=1e-4

    @test v_pef[1] ≈ -3.2256327470  atol=1e-7
    @test v_pef[2] ≈ -2.8724425110  atol=1e-7
    @test v_pef[3] ≈ +5.5319312880  atol=1e-7

    ## PEF => ITRF
    ## ===========

    r_pef  = [-1033.47503130; 7901.30558560; 6380.34453270]
    v_pef  = [-3.2256327470; -2.8724425110; +5.5319312880]

    ## DCM
    ## ---

    D_ITRF_PEF = rECEFtoECEF(PEF(), ITRF(), JD_UTC, eop_iau1980)

    r_itrf = D_ITRF_PEF*r_pef
    v_itrf = D_ITRF_PEF*v_pef

    @test r_itrf[1] ≈ -1033.4793830 atol=1e-4
    @test r_itrf[2] ≈ +7901.2952754 atol=1e-4
    @test r_itrf[3] ≈ +6380.3565958 atol=1e-4

    @test v_itrf[1] ≈ -3.225636520  atol=1e-7
    @test v_itrf[2] ≈ -2.872451450  atol=1e-7
    @test v_itrf[3] ≈ +5.531924446  atol=1e-7

    ## Quaternion
    ## ----------

    q_ITRF_PEF = rECEFtoECEF(Quaternion, PEF(), ITRF(), JD_UTC, eop_iau1980)

    r_itrf = vect(conj(q_ITRF_PEF)*r_pef*q_ITRF_PEF)
    v_itrf = vect(conj(q_ITRF_PEF)*v_pef*q_ITRF_PEF)

    @test r_itrf[1] ≈ -1033.4793830 atol=1e-4
    @test r_itrf[2] ≈ +7901.2952754 atol=1e-4
    @test r_itrf[3] ≈ +6380.3565958 atol=1e-4

    @test v_itrf[1] ≈ -3.225636520  atol=1e-7
    @test v_itrf[2] ≈ -2.872451450  atol=1e-7
    @test v_itrf[3] ≈ +5.531924446  atol=1e-7
end

################################################################################
#                                IAU-2006/2010
################################################################################

## ITRF <=> TIRS
## =============

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-14: Performing an IAU-2000 reduction [1, p. 220]
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_itrf = -1033.4793830    i + 7901.2952754    j + 6380.3565958    k [km]
#   v_itrf =    -3.225636520  i -    2.872451450  j +    5.531924446  k [km/s]
#
# one gets the following:
#
#   r_tirs  = -1033.47503120   i + 7901.30558560   j + 6380.34453270   k [km]
#   v_tirs  =    -3.2256327470 i -    2.8724425110 j +    5.5319312880 k [km/s]
#
################################################################################

@testset "Function rECEFtoECEF ITRF <=> TIRS" begin
    JD_UTC = DatetoJD(2004, 4, 6, 7, 51, 28.386009)

    ## ITRF => TIRS
    ## ===========

    r_itrf = [-1033.4793830; 7901.2952754; 6380.3565958]
    v_itrf = [-3.225636520; -2.872451450; +5.531924446]

    ## DCM
    ## ---

    D_TIRS_ITRF = rECEFtoECEF(ITRF(), TIRS(), JD_UTC, eop_iau2000a)

    r_tirs = D_TIRS_ITRF*r_itrf
    v_tirs = D_TIRS_ITRF*v_itrf

    @test r_tirs[1] ≈ -1033.47503120 atol=1e-4
    @test r_tirs[2] ≈ +7901.30558560 atol=1e-4
    @test r_tirs[3] ≈ +6380.34453270 atol=1e-4

    @test v_tirs[1] ≈ -3.2256327470  atol=1e-7
    @test v_tirs[2] ≈ -2.8724425110  atol=1e-7
    @test v_tirs[3] ≈ +5.5319312880  atol=1e-7

    ## Quaternion
    ## ----------

    q_TIRS_ITRF = rECEFtoECEF(Quaternion, ITRF(), TIRS(), JD_UTC, eop_iau2000a)

    r_tirs = vect(conj(q_TIRS_ITRF)*r_itrf*q_TIRS_ITRF)
    v_tirs = vect(conj(q_TIRS_ITRF)*v_itrf*q_TIRS_ITRF)

    @test r_tirs[1] ≈ -1033.47503120 atol=1e-4
    @test r_tirs[2] ≈ +7901.30558560 atol=1e-4
    @test r_tirs[3] ≈ +6380.34453270 atol=1e-4

    @test v_tirs[1] ≈ -3.2256327470  atol=1e-7
    @test v_tirs[2] ≈ -2.8724425110  atol=1e-7
    @test v_tirs[3] ≈ +5.5319312880  atol=1e-7

    ## TIRS => ITRF
    ## ===========

    r_tirs  = [-1033.47503120; 7901.30558560; 6380.34453270]
    v_tirs  = [-3.2256327470; -2.8724425110; +5.5319312880]

    ## DCM
    ## ---

    D_ITRF_TIRS = rECEFtoECEF(TIRS(), ITRF(), JD_UTC, eop_iau2000a)

    r_itrf = D_ITRF_TIRS*r_tirs
    v_itrf = D_ITRF_TIRS*v_tirs

    @test r_itrf[1] ≈ -1033.4793830 atol=1e-4
    @test r_itrf[2] ≈ +7901.2952754 atol=1e-4
    @test r_itrf[3] ≈ +6380.3565958 atol=1e-4

    @test v_itrf[1] ≈ -3.225636520  atol=1e-7
    @test v_itrf[2] ≈ -2.872451450  atol=1e-7
    @test v_itrf[3] ≈ +5.531924446  atol=1e-7

    ## Quaternion
    ## ----------

    q_ITRF_TIRS = rECEFtoECEF(Quaternion, TIRS(), ITRF(), JD_UTC, eop_iau2000a)

    r_itrf = vect(conj(q_ITRF_TIRS)*r_tirs*q_ITRF_TIRS)
    v_itrf = vect(conj(q_ITRF_TIRS)*v_tirs*q_ITRF_TIRS)

    @test r_itrf[1] ≈ -1033.4793830 atol=1e-4
    @test r_itrf[2] ≈ +7901.2952754 atol=1e-4
    @test r_itrf[3] ≈ +6380.3565958 atol=1e-4

    @test v_itrf[1] ≈ -3.225636520  atol=1e-7
    @test v_itrf[2] ≈ -2.872451450  atol=1e-7
    @test v_itrf[3] ≈ +5.531924446  atol=1e-7
end


text:
using PackageCompiler, Libdl

PackageCompiler.create_sysimage(Symbol[:Diffusion];
                                project=joinpath(@__DIR__, ".."),
                                precompile_execution_file=[joinpath(@__DIR__, "generate_precompile.jl")],
                                precompile_statements_file=[joinpath(@__DIR__, "additional_precompile.jl")],
                                sysimage_path=joinpath(@__DIR__, "..", "libdiffusion.$(Libdl.dlext)"),
                                incremental=false,
                                filter_stdlibs=true)

max_stars_repo_name: omlins/libdiffusion

max_issues_repo_name: omlins/libdiffusion

max_forks_repo_name: omlins/libdiffusion