text:
module TESOptimalFiltering

"Return the autocorrelation of vector data with length `nlags`. Generated by averaging together
the autocorrelations of sequenes of length `seqlen` of points chosen from `data`.
Rejects sequences that contain a point more than `max_excursion` from the mean of `data`.

The preferred way to call is `autocorrelation(data::Vector, nlags::Integer, max_excursion=250)` where `seqlen` is calculated by huersitic."
function autocorrelation(data::Vector, nlags::Integer, seqlen::Integer, max_excursion)
    @assert seqlen >= 2*nlags
    @assert seqlen < length(data)
    nseq = div( length(data), seqlen)
    d = zeros(Float64, seqlen+nlags)
    acorr = zeros(Float64, nlags)
    segments_used = 0

    for i = 1:nseq
        d[1:seqlen] = float(data[(i-1)*seqlen+1:i*seqlen]) - mean(data[(i-1)*seqlen+1:i*seqlen])
        if maximum(d) > max_excursion || minimum(d) < -max_excursion
            continue
        end
        acorr += xcorr(d, d)[seqlen+nlags:seqlen+2*nlags-1]
        segments_used += 1
    end
    if segments_used < max(div(nseq,2), 20)
        error("Rejected too many sequences. Used $segments_used segments out of $nseq available.")
    end
    acorr ./ (segments_used*collect(seqlen:-1:seqlen-nlags+1))
end

"Choose a reasonble sequence length for use in `autocorrelation`."
function autocorr_seqlen_heuristic(nlags::Integer)
    # We'll use sequences at least 7x the # of lags, then
    # round up to one of (1,3,5) times a power of 2
    total =8*nlags
    total = max(total,8)
    lowpow2 = 2^floor(Int, log2(total))
    pow2ratio = total/lowpow2
    if pow2ratio > 1.5
        total = lowpow2*2
    elseif pow2ratio > 1.25
        total = div(lowpow2*3, 2)
    else
        total = div(lowpow2*5, 4)
    end
    @assert total >= 8*nlags
    total-nlags
end

autocorrelation(data::Vector, nlags::Integer, max_excursion=250) = autocorrelation(data, nlags, autocorr_seqlen_heuristic(nlags), max_excursion)

"`calculate_filter(pulse_model, noise_autocorr, f_3db, dt)`
Return filter anf expected variance on pulse height. Filter is rolled off at `f_3db` (hertz) given that the point spacing is `dt` (seconds)"
function calculate_filter(pulse_model, noise_autocorr, f_3db, dt)
    M = hcat(pulse_model, ones(Float64, length(pulse_model)))
    # M will have 2 columns, or more. 1st column is pulse model. All other columns are
    # the things you want to be insensitive to (like exponentials).
    R = toeplitz(noise_autocorr)
    RinvM = (R \ M)
    A = M' * RinvM # "design matrix"
    variance = inv(A)[1,1] # (dV/V)^2, expected paramter covariance
    filter0 = A \ RinvM'
    filter = smoother(length(pulse_model), f_3db, dt)*vec(filter0[1,:])
    vec(filter), variance
end

"`smoother(n, f_3db, dt)`
Return an (n,n) array `s` such that `s*filter` gives a filter rolled off at `f_3db` (hertz). `dt` is the time (seconds) between samples in `filter`."
function smoother(n, f_3db, dt)
    s = Array(Float64, n, n)
    tau = 1 ./ (2*pi*f_3db*dt)
    for i = 1:n
        for j = 1:n
            s[i,j] = exp(-.5*((i-j)/tau)^2)
        end
        s[i,:] = s[i,:] ./ sum(s[i,:])
    end
    s
end

"Used to efficiently to polyfit in filter5lag."
const SG_filter = [[ -6 24  34 24 -6];
             [-14 -7   0  7 14];
             [ 10 -5 -10 -5 10]]./70.
"`filter5lag(filter, pulses)`
`pulses` can either be one pulse in a `Vector{T<:Number}`, or many pulses in a vector of pulses, or a 2d array where `a[1,:]` is the firse pulse.
Truncate the filter by 4 points, 2 on each end. Dot the trunacated filter with the pulse for 5 different lags. Fit a quadratic to the 5 lags.
Use the peak of the fit as the peak_value and the x center as the phase (aka subsample arrival time)."
function filter5lag{T<:Number}(filter, pulse::Vector{T})
    output = [dot(filter[3:end-2],pulse[i:end-5+i]) for i=1:5]
    a,b,c = SG_filter*output #quadratic fit to 5 lags
    filt_value = a-b^2/c/4 # peak of quadratic
    filt_phase = -b/c/2 # x location of peak
    filt_value = isfinite(filt_value) ? filt_value : -1
    filt_phase = isfinite(filt_phase) ? filt_phase : -10
    filt_value, filt_phase
end
function filter5lag{T<:Vector}(filter, pulses::Vector{T})
    filt_values = zeros(Float64, length(pulses))
    filt_phases = zeros(Float64, length(pulses))
    for i in 1:length(pulses)
        filt_values[i], filt_phases[i] = filter5lag(filter, pulses[i])
    end
    filt_values, filt_phases
end
function filter5lag{T<:Number}(filter, pulses::Array{T,2})
    filt_values = zeros(Float64, size(pulses, 2))
    filt_phases = zeros(Float64, size(pulses, 2))
    for j in 1:size(pulses,2)
        output = [dot(filter[3:end-2],pulses[i:end-5+i,j]) for i=1:5]
        pfit = SG_filter*output
        filt_value = pfit[1]-pfit[2]^2/4/pfit[3]
        filt_phase = -pfit[2]/pfit[3]/2
        filt_values[j] = isfinite(filt_value) ? filt_value : -1
        filt_phases[j] = isfinite(filt_phase) ? filt_phase : -10
    end
    filt_values, filt_phases
end

"Create and return a toeplitz matrix with c as the first column and
(optional) r as the first row. r[1] will be ignored in
favor of c[1] to determine the diagonal."
function toeplitz{T}(c::Vector{T}, r::Vector{T})
    nr = length(c)
    nc = length(r)
    m = Array(typeof(c[1]), nr, nc)
    for j=1:nc
        for i=1:j-1  # Strict upper triangle
            m[i,j] = r[j-i+1]
        end
        for i=j:nr
            m[i,j] = c[i-j+1]
        end
    end
    m
end
toeplitz{T<:Number}(c::Vector{T}) = toeplitz(c,c)

end # module


text:
/* Test boolean operators. */

int main() {
  printString("&&");
  printBool(test(-1) && test(0));
  printBool(test(-2) && test(1));
  printBool(test(3) && test(-5));
  printBool(test(234234) && test(21321));
  printString("||");
  printBool(test(-1) || test(0));
  printBool(test(-2) || test(1));
  printBool(test(3) || test(-5));
  printBool(test(234234) || test(21321));
  printString("!");
  printBool(true);
  printBool(false);
  return 0 ;

}

void printBool(boolean b) {
  if (!b) {
    printString("false");
  } else {
    printString("true");
 }
 return;
}

boolean test(int i) {
  printInt(i);
  return i > 0;
}

ext: jl

max_stars_repo_name: jwilk/jtc

max_issues_repo_name: jwilk/jtc

max_forks_repo_name: jwilk/jtc


text:
module ManifoldsBase

import Base:
    isapprox,
    exp,
    log,
    convert,
    copyto!,
    angle,
    eltype,
    isempty,
    length,
    similar,
    show,
    +,
    -,
    *
import LinearAlgebra: dot, norm, det, cross, I, UniformScaling, Diagonal

import Markdown: @doc_str
using LinearAlgebra

"""
    Manifold{F}

A manifold type. The `Manifold` is used to dispatch to different functions on a manifold,
usually as the first argument of the function. Examples are the [`exp`](@ref)onential and
[`log`](@ref)arithmic maps as well as more general functions that are built on them like the
[`geodesic`](@ref).

The manifold is parametrized by an [`AbstractNumbers`](@ref) to distinguish for example
real (ℝ) and complex (ℂ) manifolds.

For subtypes the preferred order of parameters is: size and simple value parameters,
followed by the [`AbstractNumbers`](@ref) `field`, followed by data type parameters,
which might depend on the abstract number field type.
"""
abstract type Manifold{𝔽} end

"""
    AbstractEstimationMethod

Abstract type for defining statistical estimation methods.
"""
abstract type AbstractEstimationMethod end

"""
    AbstractInverseRetractionMethod

Abstract type for methods for inverting a retraction (see [`inverse_retract`](@ref)).
"""
abstract type AbstractInverseRetractionMethod end

"""
    AbstractRetractionMethod

Abstract type for methods for [`retract`](@ref)ing a tangent vector to a manifold.
"""
abstract type AbstractRetractionMethod end

"""
    ExponentialRetraction

Retraction using the exponential map.
"""
struct ExponentialRetraction <: AbstractRetractionMethod end

"""
    PolarRetraction <: AbstractRetractionMethod

Retractions that are based on singular value decompositions of the matrix / matrices
for point and tangent vector on a [`Manifold`](@ref)
"""
struct PolarRetraction <: AbstractRetractionMethod end

"""
    ProjectionRetraction <: AbstractRetractionMethod

Retractions that are based on projection and usually addition in the embedding.
"""
struct ProjectionRetraction <: AbstractRetractionMethod end

"""
    QRRetraction <: AbstractRetractionMethod

Retractions that are based on a QR decomposition of the
matrix / matrices for point and tangent vector on a [`Manifold`](@ref)
"""
struct QRRetraction <: AbstractRetractionMethod end

"""
    LogarithmicInverseRetraction

Inverse retraction using the [`log`](@ref)arithmic map.
"""
struct LogarithmicInverseRetraction <: AbstractInverseRetractionMethod end

"""
    PolarInverseRetraction <: AbstractInverseRetractionMethod

Inverse retractions that are based on a singular value decomposition of the
matrix / matrices for point and tangent vector on a [`Manifold`](@ref)
"""
struct PolarInverseRetraction <: AbstractInverseRetractionMethod end

"""
    ProjectionInverseRetraction <: AbstractInverseRetractionMethod

Inverse retractions that are based on a projection (or its inversion).
"""
struct ProjectionInverseRetraction <: AbstractInverseRetractionMethod end

"""
    QRInverseRetraction <: AbstractInverseRetractionMethod

Inverse retractions that are based on a QR decomposition of the
matrix / matrices for point and tangent vector on a [`Manifold`](@ref)
"""
struct QRInverseRetraction <: AbstractInverseRetractionMethod end

"""
    MPoint

Type for a point on a manifold. While a [`Manifold`](@ref) does not necessarily require this
type, for example when it is implemented for `Vector`s or `Matrix` type elements, this type
can be used for more complicated representations, semantic verification, or even dispatch
for different representations of points on a manifold.
"""
abstract type MPoint end

"""
    OutOfInjectivityRadiusError

An error thrown when a function (for example [`log`](@ref)arithmic map or
[`inverse_retract`](@ref)) is given arguments outside of its [`injectivity_radius`](@ref).
"""
struct OutOfInjectivityRadiusError <: Exception end

"""
    TVector

Type for a tangent vector of a manifold. While a [`Manifold`](@ref) does not necessarily
require this type, for example when it is implemented for `Vector`s or `Matrix` type
elements, this type can be used for more complicated representations, semantic verification,
or even dispatch for different representations of tangent vectors and their types on a
manifold.
"""
abstract type TVector end

"""
    CoTVector

Type for a cotangent vector of a manifold. While a [`Manifold`](@ref) does not necessarily
require this type, for example when it is implemented for `Vector`s or `Matrix` type
elements, this type can be used for more complicated representations, semantic verification,
or even dispatch for different representations of cotangent vectors and their types on a
manifold.
"""
abstract type CoTVector end

"""
    allocate(a)
    allocate(a, dims::Integer...)
    allocate(a, dims::Tuple)
    allocate(a, T::Type)
    allocate(a, T::Type, dims::Integer...)
    allocate(a, T::Type, dims::Tuple)

Allocate an object similar to `a`. It is similar to function `similar`, although
instead of working only on the outermost layer of a nested structure, it maps recursively
through outer layers and calls `similar` on the innermost array-like object only.
Type `T` is the new number element type [`number_eltype`](@ref), if it is not given
the element type of `a` is retained. The `dims` argument can be given for non-nested
allocation and is forwarded to the function `similar`.
"""
allocate(a, args...)
allocate(a) = similar(a)
allocate(a, dims::Integer...) = similar(a, dims...)
allocate(a, dims::Tuple) = similar(a, dims)
allocate(a, T::Type) = similar(a, T)
allocate(a, T::Type, dims::Integer...) = similar(a, T, dims...)
allocate(a, T::Type, dims::Tuple) = similar(a, T, dims)
allocate(a::AbstractArray{<:AbstractArray}) = map(allocate, a)
allocate(a::AbstractArray{<:AbstractArray}, T::Type) = map(t -> allocate(t, T), a)
allocate(a::NTuple{N,AbstractArray} where {N}) = map(allocate, a)
allocate(a::NTuple{N,AbstractArray} where {N}, T::Type) = map(t -> allocate(t, T), a)

"""
    allocate_result(M::Manifold, f, x...)

Allocate an array for the result of function `f` on [`Manifold`](@ref) `M` and arguments
`x...` for implementing the non-modifying operation using the modifying operation.

Usefulness of passing a function is demonstrated by methods that allocate results of musical
isomorphisms.
"""
function allocate_result(M::Manifold, f, x...)
    T = allocate_result_type(M, f, x)
    return allocate(x[1], T)
end

"""
    allocate_result_type(M::Manifold, f, args::NTuple{N,Any}) where N

Return type of element of the array that will represent the result of function `f` and the
[`Manifold`](@ref) `M` on given arguments `args` (passed as a tuple).
"""
function allocate_result_type(M::Manifold, f, args::NTuple{N,Any}) where {N}
    return typeof(mapreduce(eti -> one(number_eltype(eti)), +, args))
end

"""
    angle(M::Manifold, p, X, Y)

Compute the angle between tangent vectors `X` and `Y` at point `p` from the
[`Manifold`](@ref) `M` with respect to the inner product from [`inner`](@ref).
"""
angle(M::Manifold, p, X, Y) = acos(real(inner(M, p, X, Y)) / norm(M, p, X) / norm(M, p, Y))

"""
    base_manifold(M::Manifold, depth = Val(-1))

Return the internally stored [`Manifold`](@ref) for decorated manifold `M` and the base
manifold for vector bundles or power manifolds. The optional parameter `depth` can be used
to remove only the first `depth` many decorators and return the [`Manifold`](@ref) from that
level, whether its decorated or not. Any negative value deactivates this depth limit.
"""
base_manifold(M::Manifold, depth = Val(-1)) = M

"""
    check_manifold_point(M::Manifold, p; kwargs...) -> Union{Nothing,String}

Return `nothing` when `p` is a point on the [`Manifold`](@ref) `M`. Otherwise, return an
error with description why the point does not belong to manifold `M`.

By default, `check_manifold_point` returns `nothing`, i.e. if no checks are implemented, the
assumption is to be optimistic for a point not deriving from the [`MPoint`](@ref) type.
"""
check_manifold_point(M::Manifold, p; kwargs...) = nothing

"""
    check_tangent_vector(M::Manifold, p, X; kwargs...) -> Union{Nothing,String}

Check whether `X` is a valid tangent vector in the tangent space of `p` on the
[`Manifold`](@ref) `M`. An implementation should first call [`check_manifold_point(M, p;
kwargs...)`](@ref) and then validate `X`. If it is not a tangent vector, an error string
should be returned.

By default, `check_tangent_vector` returns `nothing`, i.e. if no checks are implemented, the
assumption is to be optimistic for tangent vectors not deriving from the [`TVector`](@ref)
type.
"""
check_tangent_vector(M::Manifold, p, X; kwargs...) = nothing

"""
    check_size(M::Manifold, p)
    check_size(M::Manifold, p, X)

Check whether `p` has the right [`representation_size`](@ref) for a [`Manifold`](@ref) `M`.
Additionally if a tangent vector is given, both `p` and `X` are checked to be of
corresponding correct representation sizes for points and tangent vectors on `M`.

By default, `check_size` returns `nothing`, i.e. if no checks are implemented, the
assumption is to be optimistic.
"""
function check_size(M::Manifold, p)
    n = size(p)
    m = representation_size(M)
    if length(n) != length(m)
        return DomainError(
            length(n),
            "The point $(p) can not belong to the manifold $(M), since its size $(n) is not equal to the manifolds representation size ($(m)).",
        )
    end
    if n != m
        return DomainError(
            n,
            "The point $(p) can not belong to the manifold $(M), since its size $(n) is not equal to the manifolds representation size ($(m)).",
        )
    end
end
function check_size(M::Manifold, p, X)
    mse = check_size(M, p)
    mse === nothing || return mse
    n = size(X)
    m = representation_size(M)
    if length(n) != length(m)
        return DomainError(
            length(n),
            "The tangent vector $(X) can not belong to the manifold $(M), since its size $(n) is not equal to the manifolds representation size ($(m)).",
        )
    end
    if n != m
        return DomainError(
            n,
            "The tangent vector $(X) can not belong to the manifold $(M), since its size $(n) is not equal to the manifodls representation size ($(m)).",
        )
    end
end

@doc raw"""
    copyto!(M::Manifold, q, p)

Copy the value(s) from `p` to `q`, where both are points on the [`Manifold`](@ref) `M`.
This function defaults to calling `copyto!(q, p)`, but it might be useful to overwrite the
function at the level, where also information from `M` can be accessed.
"""
copyto!(::Manifold, q, p) = copyto!(q, p)

@doc raw"""
    copyto!(M::Manifold, Y, p, X)

Copy the value(s) from `X` to `Y`, where both are tangent vectors from the tangent space at
`p` on the [`Manifold`](@ref) `M`.
This function defaults to calling `copyto!(Y, X)`, but it might be useful to overwrite the
function at the level, where also information from `p` and `M` can be accessed.
"""
copyto!(::Manifold, Y, p, X) = copyto!(Y, X)

"""
    distance(M::Manifold, p, q)

Shortest distance between the points `p` and `q` on the [`Manifold`](@ref) `M`.
"""
distance(M::Manifold, p, q) = norm(M, p, log(M, p, q))

"""
    embed(M::Manifold, p)

Embed point `p` from the [`Manifold`](@ref) `M` into the ambient space.
This method is only available for manifolds where implicitly an embedding or ambient space
is given.
Additionally, `embed` includes changing data representation, if applicable, i.e.
if the points on `M` are not represented in the same way as points on the embedding,
the representation is changed accordingly.

See also: [`EmbeddedManifold`](@ref), [`project`](@ref project(M::Manifold,p))
"""
function embed(M::Manifold, p)
    q = allocate_result(M, embed, p)
    embed!(M, q, p)
    return q
end

"""
    embed!(M::Manifold, q, p)

Embed point `p` from the [`Manifold`](@ref) `M` into an ambient space.
This method is only available for manifolds where implicitly an embedding or ambient space
is given. Not implementing this function means, there is no proper embedding for your manifold.
Additionally, `embed` might include changing data representation, if applicable, i.e.
if points on `M` are not represented in the same way as their counterparts in the embedding,
the representation is changed accordingly.

If you have more than one embedding, see [`EmbeddedManifold`](@ref) for defining a second
embedding. If your point `p` is already represented in some embedding,
see [`AbstractEmbeddedManifold`](@ref) how you can avoid reimplementing code from the embedded manifold

See also: [`EmbeddedManifold`](@ref), [`project!`](@ref project!(M::Manifold, q, p))
"""
function embed!(M::Manifold, q, p)
    return error(manifold_function_not_implemented_message(M, embed!, q, p))
end

"""
    embed(M::Manifold, p, X)

Embed a tangent vector `X` at a point `p` on the [`Manifold`](@ref) `M` into an ambient space.
This method is only available for manifolds where implicitly an embedding or ambient space
is given. Not implementing this function means, there is no proper embedding for your tangent space(s).

Additionally, `embed` might include changing data representation, if applicable, i.e.
if tangent vectors on `M` are not represented in the same way as their counterparts in the
embedding, the representation is changed accordingly.

If you have more than one embedding, see [`EmbeddedManifold`](@ref) for defining a second
embedding. If your tangent vector `X` is already represented in some embedding,
see [`AbstractEmbeddedManifold`](@ref) how you can avoid reimplementing code from the embedded manifold

See also: [`EmbeddedManifold`](@ref), [`project`](@ref project(M::Manifold, p, X))
"""
function embed(M::Manifold, p, X)
    # Note that the order is switched,
    # since the allocation by default takes the type of the first.
    Y = allocate_result(M, embed, X, p)
    embed!(M, Y, p, X)
    return Y
end

"""
    embed!(M::Manifold, Y, p, X)

Embed a tangent vector `X` at a point `p` on the [`Manifold`](@ref) `M` into the ambient
space and return the result in `Y`.
This method is only available for manifolds where implicitly an embedding or ambient space
is given.
Additionally, `embed!` includes changing data representation, if applicable, i.e.
if the tangents on `M` are not represented in the same way as tangents on the embedding,
the representation is changed accordingly. This is the case for example for Lie groups,
when tangent vectors are represented in the Lie algebra. The embedded tangents are then in
the tangent spaces of the embedded base points.

See also: [`EmbeddedManifold`](@ref), [`project!`](@ref project!(M::Manifold, Y, p, X))
"""
function embed!(M::Manifold, Y, p, X)
    return error(manifold_function_not_implemented_message(M, embed!, Y, p, X))
end

"""
    exp(M::Manifold, p, X)
    exp(M::Manifold, p, X, t::Real = 1)

Compute the exponential map of tangent vector `X`, optionally scaled by `t`,  at point `p`
from manifold the [`Manifold`](@ref) `M`.
"""
function exp(M::Manifold, p, X)
    q = allocate_result(M, exp, p, X)
    exp!(M, q, p, X)
    return q
end
exp(M::Manifold, p, X, t::Real) = exp(M, p, t * X)

"""
    exp!(M::Manifold, q, p, X)
    exp!(M::Manifold, q, p, X, t::Real = 1)

Compute the exponential map of tangent vector `X`, optionally scaled by `t`,  at point `p`
from manifold the [`Manifold`](@ref) `M`.
The result is saved to `q`.
"""
function exp!(M::Manifold, q, p, X)
    return error(manifold_function_not_implemented_message(M, exp!, q, p, X))
end
exp!(M::Manifold, q, p, X, t::Real) = exp!(M, q, p, t * X)

"""
    geodesic(M::Manifold, p, X) -> Function

Get the geodesic with initial point `p` and velocity `X` on the [`Manifold`](@ref) `M`.
 The geodesic is the curve of constant velocity that is locally distance-minimizing.
 This function returns a function of (time) `t`.

    geodesic(M::Manifold, x, v, t::Real)
    geodesic(M::Manifold, x, v, T::AbstractVector) -> AbstractVector

Return the point at time `t` or points at times `t` in `T` along the geodesic.
"""
geodesic(M::Manifold, p, X) = t -> exp(M, p, X, t)
geodesic(M::Manifold, p, X, t::Real) = exp(M, p, X, t)
geodesic(M::Manifold, p, X, T::AbstractVector) = map(t -> exp(M, p, X, t), T)

@doc doc"""
    injectivity_radius(M::Manifold, p)

Return the distance $d$ such that [`exp(M, p, X)`](@ref exp(::Manifold, ::Any, ::Any)) is
injective for all tangent vectors shorter than $d$ (i.e. has an inverse).

    injectivity_radius(M::Manifold)

Infimum of the injectivity radius of all manifold points.

    injectivity_radius(M::Manifold[, x], method::AbstractRetractionMethod)
    injectivity_radius(M::Manifold, x, method::AbstractRetractionMethod)

Distance $d$ such that
[`retract(M, p, X, method)`](@ref retract(::Manifold, ::Any, ::Any, ::AbstractRetractionMethod))
is injective for all tangent vectors shorter than $d$ (i.e. has an inverse) for point `p`
if provided or all manifold points otherwise.
"""
function injectivity_radius(M::Manifold)
    return error(manifold_function_not_implemented_message(M, injectivity_radius))
end
injectivity_radius(M::Manifold, p) = injectivity_radius(M)
function injectivity_radius(M::Manifold, p, method::AbstractRetractionMethod)
    return injectivity_radius(M, method)
end
function injectivity_radius(M::Manifold, method::AbstractRetractionMethod)
    return error(manifold_function_not_implemented_message(M, injectivity_radius, method))
end
injectivity_radius(M::Manifold, p, ::ExponentialRetraction) = injectivity_radius(M, p)
injectivity_radius(M::Manifold, ::ExponentialRetraction) = injectivity_radius(M)

"""
    inner(M::Manifold, p, X, Y)

Compute the inner product of tangent vectors `X` and `Y` at point `p` from the
[`Manifold`](@ref) `M`.

See also: [`MetricManifold`](@ref Main.Manifolds.MetricManifold)
"""
function inner(M::Manifold, p, X, Y)
    return error(manifold_function_not_implemented_message(M, inner, p, X, Y))
end

"""
    inverse_retract!(M::Manifold, X, p, q[, method::AbstractInverseRetractionMethod])

Compute the inverse retraction, a cheaper, approximate version of the
[`log`](@ref)arithmic map), of points `p` and `q` on the [`Manifold`](@ref) `M`.
Result is saved to `X`.

Inverse retraction method can be specified by the last argument, defaulting to
[`LogarithmicInverseRetraction`](@ref). See the documentation of respective manifolds for
available methods.
"""
function inverse_retract!(M::Manifold, X, p, q)
    return inverse_retract!(M, X, p, q, LogarithmicInverseRetraction())
end
function inverse_retract!(M::Manifold, X, p, q, method::LogarithmicInverseRetraction)
    return log!(M, X, p, q)
end
function inverse_retract!(M::Manifold, X, p, q, method::AbstractInverseRetractionMethod)
    return error(
        manifold_function_not_implemented_message(M, inverse_retract!, X, p, q, method),
    )
end

"""
    inverse_retract(M::Manifold, x, y)
    inverse_retract(M::Manifold, x, y, method::AbstractInverseRetractionMethod

Compute the inverse retraction, a cheaper, approximate version of the
[`log`](@ref)arithmic map), of points `p` and `q` on the [`Manifold`](@ref) `M`.

Inverse retraction method can be specified by the last argument, defaulting to
[`LogarithmicInverseRetraction`](@ref). See the documentation of respective manifolds for
available methods.
"""
function inverse_retract(M::Manifold, p, q)
    X = allocate_result(M, inverse_retract, p, q)
    inverse_retract!(M, X, p, q)
    return X
end
function inverse_retract(M::Manifold, p, q, method::AbstractInverseRetractionMethod)
    X = allocate_result(M, inverse_retract, p, q)
    inverse_retract!(M, X, p, q, method)
    return X
end

"""
    isapprox(M::Manifold, p, q; kwargs...)

Check if points `p` and `q` from [`Manifold`](@ref) `M` are approximately equal.

Keyword arguments can be used to specify tolerances.
"""
isapprox(::Manifold, x, y; kwargs...) = isapprox(x, y; kwargs...)

"""
    isapprox(M::Manifold, p, X, Y; kwargs...)

Check if vectors `X` and `Y` tangent at `p` from [`Manifold`](@ref) `M` are approximately
equal.

Keyword arguments can be used to specify tolerances.
"""
isapprox(::Manifold, p, X, Y; kwargs...) = isapprox(X, Y; kwargs...)

"""
    is_manifold_point(M::Manifold, p, throw_error = false; kwargs...)

Return whether `p` is a valid point on the [`Manifold`](@ref) `M`.

If `throw_error` is `false`, the function returns either `true` or `false`. If `throw_error`
is `true`, the function either returns `true` or throws an error. By default the function
calls [`check_manifold_point(M, p; kwargs...)`](@ref) and checks whether the returned value
is `nothing` or an error.
"""
function is_manifold_point(M::Manifold, p, throw_error = false; kwargs...)
    mpe = check_manifold_point(M, p; kwargs...)
    mpe === nothing && return true
    return throw_error ? throw(mpe) : false
end

"""
    is_tangent_vector(M::Manifold, p, X, throw_error = false; kwargs...)

Return whether `X` is a valid tangent vector at point `p` on the [`Manifold`](@ref) `M`.
Returns either `true` or `false`.

If `throw_error` is `false`, the function returns either `true` or `false`. If `throw_error`
is `true`, the function either returns `true` or throws an error. By default the function
calls [`check_tangent_vector(M, p, X; kwargs...)`](@ref) and checks whether the returned
value is `nothing` or an error.
"""
function is_tangent_vector(M::Manifold, p, X, throw_error = false; kwargs...)
    mtve = check_tangent_vector(M, p, X; kwargs...)
    mtve === nothing && return true
    return throw_error ? throw(mtve) : false
end

"""
    log(M::Manifold, p, q)

Compute the logarithmic map of point `q` at base point `p` on the [`Manifold`](@ref) `M`.
"""
function log(M::Manifold, p, q)
    X = allocate_result(M, log, p, q)
    log!(M, X, p, q)
    return X
end

"""
    log!(M::Manifold, X, p, q)

Compute the logarithmic map of point `q` at base point `p` on the [`Manifold`](@ref) `M`.
The result is saved to `X`.
"""
function log!(M::Manifold, X, p, q)
    return error(manifold_function_not_implemented_message(M, log!, X, p, q))
end

@doc doc"""
    manifold_dimension(M::Manifold)

The dimension $n=\dim_{\mathcal M}$ of real space $\mathbb R^n$ to which the neighborhood of
each point of the [`Manifold`](@ref) `M` is homeomorphic.
"""
function manifold_dimension(M::Manifold)
    return error(manifold_function_not_implemented_message(M, manifold_dimension))
end

function manifold_function_not_implemented_message(M::Manifold, f, x...)
    s = join(map(string, map(typeof, x)), ", ", " and ")
    a = length(x) > 1 ? "arguments" : "argument"
    m = length(x) > 0 ? " for $(a) $(s)." : "."
    return "$(f) not implemented on $(M)$(m)"
end

"""
    mid_point(M::Manifold, p1, p2)

Calculate the middle between the two point `p1` and `p2` from manifold `M`.
By default uses [`log`](@ref), divides the vector by 2 and uses [`exp`](@ref).
"""
function mid_point(M::Manifold, p1, p2)
    q = allocate(p1)
    return mid_point!(M, q, p1, p2)
end

"""
    mid_point!(M::Manifold, q, p1, p2)

Calculate the middle between the two point `p1` and `p2` from manifold `M`.
By default uses [`log`](@ref), divides the vector by 2 and uses [`exp!`](@ref).
Saves the result in `q`.
"""
function mid_point!(M::Manifold, q, p1, p2)
    X = log(M, p1, p2)
    return exp!(M, q, p1, X / 2)
end

@static if VERSION <= v"1.1"
    function mid_point!(
        M::Manifold,
        q::AbstractArray{T1,0},
        p1::AbstractArray{T2,0},
        p2::AbstractArray{T3,0},
    ) where {T1,T2,T3}
        X = log(M, p1, p2)
        return exp!(M, q, p1, fill(X / 2))
    end
end

"""
    norm(M::Manifold, p, X)

Compute the norm of tangent vector `X` at point `p` from a [`Manifold`](@ref) `M`.
By default this is computed using [`inner`](@ref).
"""
norm(M::Manifold, p, X) = sqrt(real(inner(M, p, X, X)))

"""
    number_eltype(x)

Numeric element type of the a nested representation of a point or a vector.
To be used in conjuntion with [`allocate`](@ref) or [`allocate_result`](@ref).
"""
number_eltype(x) = eltype(x)
function number_eltype(x::AbstractArray)
    return typeof(mapreduce(eti -> one(number_eltype(eti)), +, x))
end
function number_eltype(x::Tuple)
    return typeof(mapreduce(eti -> one(number_eltype(eti)), +, x))
end

"""
    project(M::Manifold, p)

Project point `p` from the ambient space of the [`Manifold`](@ref) `M` to `M`.
This method is only available for manifolds where implicitly an embedding or ambient space
is given. Additionally, the projection includes changing data representation, if applicable,
i.e. if the points on `M` are not represented in the same array data, the data is changed
accordingly.

See also: [`EmbeddedManifold`](@ref), [`embed`](@ref embed(M::Manifold, p))
"""
function project(M::Manifold, p)
    q = allocate_result(M, project, p)
    project!(M, q, p)
    return q
end

"""
    project!(M::Manifold, q, p)

Project point `p` from the ambient space onto the [`Manifold`](@ref) `M`. The result is
storedin `q`.
This method is only available for manifolds where implicitly an embedding or ambient space
is given. Additionally, the projection includes changing data representation, if applicable,
i.e. if the points on `M` are not represented in the same array data, the data is changed
accordingly.

See also: [`EmbeddedManifold`](@ref), [`embed!`](@ref embed!(M::Manifold, q, p))
"""
function project!(M::Manifold, q, p)
    return error(manifold_function_not_implemented_message(M, project!, q, p))
end

"""
    project(M::Manifold, p, X)

Project ambient space representation of a vector `X` to a tangent vector at point `p` on
the [`Manifold`](@ref) `M`.
This method is only available for manifolds where implicitly an embedding or ambient space
is given.
Additionally, `project` includes changing data representation, if applicable, i.e.
if the tangents on `M` are not represented in the same way as points on the embedding,
the representation is changed accordingly. This is the case for example for Lie groups,
when tangent vectors are represented in the Lie algebra. after projection the change to the
Lie algebra is perfomed, too.

See also: [`EmbeddedManifold`](@ref), [`embed`](@ref embed(M::Manifold, p, X))
"""
function project(M::Manifold, p, X)
    # Note that the order is switched,
    # since the allocation by default takes the type of the first.
    Y = allocate_result(M, project, X, p)
    project!(M, Y, p, X)
    return Y
end

"""
    project!(M::Manifold, Y, p, X)

Project ambient space representation of a vector `X` to a tangent vector at point `p` on
the [`Manifold`](@ref) `M`. The result is saved in vector `Y`.
This method is only available for manifolds where implicitly an embedding or ambient space
is given.
Additionally, `project!` includes changing data representation, if applicable, i.e.
if the tangents on `M` are not represented in the same way as points on the embedding,
the representation is changed accordingly. This is the case for example for Lie groups,
when tangent vectors are represented in the Lie algebra. after projection the change to the
Lie algebra is perfomed, too.

See also: [`EmbeddedManifold`](@ref), [`embed!`](@ref embed!(M::Manifold, Y, p, X))
"""
function project!(M::Manifold, Y, p, X)
    return error(manifold_function_not_implemented_message(M, project!, Y, p, X))
end

@doc doc"""
    representation_size(M::Manifold)

The size of an array representing a point on [`Manifold`](@ref) `M`.
Returns `nothing` by default indicating that points are not represented using an
`AbstractArray`.
"""
function representation_size(M::Manifold)
    return nothing
end

"""
    retract(M::Manifold, p, X)
    retract(M::Manifold, p, X, t::Real=1)
    retract(M::Manifold, p, X, method::AbstractRetractionMethod)
    retract(M::Manifold, p, X, t::Real=1, method::AbstractRetractionMethod)

Compute a retraction, a cheaper, approximate version of the [`exp`](@ref)onential map,
from `p` into direction `X`, scaled by `t`, on the [`Manifold`](@ref) `M`.

Retraction method can be specified by the last argument, defaulting to
[`ExponentialRetraction`](@ref). See the documentation of respective manifolds for available
methods.
"""
function retract(M::Manifold, p, X)
    q = allocate_result(M, retract, p, X)
    retract!(M, q, p, X)
    return q
end
retract(M::Manifold, p, X, t::Real) = retract(M, p, t * X)
function retract(M::Manifold, p, X, method::AbstractRetractionMethod)
    q = allocate_result(M, retract, p, X)
    retract!(M, q, p, X, method)
    return q
end
function retract(M::Manifold, p, X, t::Real, method::AbstractRetractionMethod)
    return retract(M, p, t * X, method)
end

"""
    retract!(M::Manifold, q, p, X)
    retract!(M::Manifold, q, p, X, t::Real=1)
    retract!(M::Manifold, q, p, X, method::AbstractRetractionMethod)
    retract!(M::Manifold, q, p, X, t::Real=1, method::AbstractRetractionMethod)

Compute a retraction, a cheaper, approximate version of the [`exp`](@ref)onential map,
from `p` into direction `X`, scaled by `t`, on the [`Manifold`](@ref) manifold `M`.
Result is saved to `q`.

Retraction method can be specified by the last argument, defaulting to
[`ExponentialRetraction`](@ref). See the documentation of respective manifolds for available
methods.
"""
retract!(M::Manifold, q, p, X) = retract!(M, q, p, X, ExponentialRetraction())
retract!(M::Manifold, q, p, X, t::Real) = retract!(M, q, p, t * X)
retract!(M::Manifold, q, p, X, method::ExponentialRetraction) = exp!(M, q, p, X)
function retract!(M::Manifold, q, p, X, t::Real, method::AbstractRetractionMethod)
    return retract!(M, q, p, t * X, method)
end
function retract!(M::Manifold, q, p, X, method::AbstractRetractionMethod)
    return error(manifold_function_not_implemented_message(M, retract!, q, p, method))
end

@doc doc"""
    shortest_geodesic(M::Manifold, p, q) -> Function

Get a [`geodesic`](@ref) $\gamma_{p,q}(t)$ whose length is the shortest path between the
points `p`and `q`, where $\gamma_{p,q}(0)=p$ and $\gamma_{p,q}(1)=q$. When there are
multiple shortest geodesics, there is no guarantee which will be returned.

This function returns a function of time, which may be a `Real` or an `AbstractVector`.

    shortest_geodesic(M::Manifold, p, q, t::Real)
    shortest_geodesic(M::Manifold, p, q, T::AbstractVector) -> AbstractVector

Return the point at time `t` or points at times `t` in `T` along the shortest geodesic.
"""
shortest_geodesic(M::Manifold, p, q) = geodesic(M, p, log(M, p, q))
shortest_geodesic(M::Manifold, p, q, t::Real) = geodesic(M, p, log(M, p, q), t)
shortest_geodesic(M::Manifold, p, q, T::AbstractVector) = geodesic(M, p, log(M, p, q), T)

"""
    size_to_tuple(::Type{S}) where S<:Tuple

Converts a size given by `Tuple{N, M, ...}` into a tuple `(N, M, ...)`.
"""
Base.@pure size_to_tuple(::Type{S}) where {S<:Tuple} = tuple(S.parameters...)

"""
    zero_tangent_vector!(M::Manifold, X, p)

Save to `X` a vector such that retracting `X` to the [`Manifold`](@ref) `M` at `p`
produces `p`.
"""
zero_tangent_vector!(M::Manifold, X, p) = log!(M, X, p, p)

"""
    zero_tangent_vector(M::Manifold, p)

Return the tangent vector from the tangent space at `p` on the [`Manifold`](@ref) `M`, that
represents the zero vector, i.e. such that a retraction at `p` produces `p`.
"""
function zero_tangent_vector(M::Manifold, p)
    X = allocate_result(M, zero_tangent_vector, p)
    zero_tangent_vector!(M, X, p)
    return X
end
include("errors.jl")
include("numbers.jl")
include("vector_transport.jl")
include("DecoratorManifold.jl")
include("bases.jl")
include("ValidationManifold.jl")
include("EmbeddedManifold.jl")
include("DefaultManifold.jl")
include("PowerManifold.jl")

export Manifold, MPoint, TVector, CoTVector
export AbstractDecoratorManifold
export ValidationManifold, ValidationMPoint, ValidationTVector, ValidationCoTVector
export AbstractEmbeddingType,
    TransparentIsometricEmbedding, DefaultIsometricEmbeddingType, DefaultEmbeddingType
export AbstractEmbeddedManifold, EmbeddedManifold, TransparentIsometricEmbedding
export AbstractPowerManifold, PowerManifold
export AbstractPowerRepresentation, NestedPowerRepresentation

export OutOfInjectivityRadiusError

export AbstractRetractionMethod,
    ExponentialRetraction,
    QRRetraction,
    PolarRetraction,
    ProjectionRetraction,
    PowerRetraction,
    InversePowerRetraction

export AbstractInverseRetractionMethod,
    LogarithmicInverseRetraction,
    QRInverseRetraction,
    PolarInverseRetraction,
    ProjectionInverseRetraction

export AbstractVectorTransportMethod,
    DifferentiatedRetractionVectorTransport,
    ParallelTransport,
    PoleLadderTransport,
    PowerVectorTransport,
    ProjectionTransport,
    ScaledVectorTransport,
    SchildsLadderTransport

export CachedBasis,
    DefaultBasis,
    DefaultOrthogonalBasis,
    DefaultOrthonormalBasis,
    DiagonalizingOrthonormalBasis,
    DefaultOrthonormalBasis,
    GramSchmidtOrthonormalBasis,
    ProjectedOrthonormalBasis

export CompositeManifoldError, ComponentManifoldError

export allocate,
    base_manifold,
    check_manifold_point,
    check_tangent_vector,
    check_size,
    distance,
    exp,
    exp!,
    embed,
    embed!,
    geodesic,
    get_basis,
    get_component,
    get_component!,
    get_coordinates,
    get_coordinates!,
    get_embedding,
    get_vector,
    get_vector!,
    get_vectors,
    hat,
    hat!,
    shortest_geodesic,
    injectivity_radius,
    inner,
    inverse_retract,
    inverse_retract!,
    isapprox,
    is_manifold_point,
    is_tangent_vector,
    isempty,
    length,
    log,
    log!,
    manifold_dimension,
    mid_point,
    mid_point!,
    norm,
    number_eltype,
    number_of_coordinates,
    number_system,
    power_dimensions,
    project,
    project!,
    real_dimension,
    representation_size,
    set_component!,
    show,
    retract,
    retract!,
    vector_transport_along,
    vector_transport_along!,
    vector_transport_direction,
    vector_transport_direction!,
    vector_transport_to,
    vector_transport_to!,
    vee,
    vee!,
    zero_tangent_vector,
    zero_tangent_vector!

end # module


text:
# Example 2.5 from A. Law (2013): Simulation Modeling and Analysis, 5ed

using EventSimulation
using Distributions   # Exponential distribution
using StatsBase       # weighted mean
using DataFrames      # nicer ouptut of results

mutable struct Job
    tin::Float64  # time when job entered the system
    tout::Float64 # time when job left the system
    ttrue::Float64 # time needed without cycling
    remaining::Float64 # time remaining to finish the job
end

mutable struct SharedComputer <: AbstractState
    a::Exponential{Float64}    # arrival time distribution
    p::Exponential{Float64}    # processing time distribution
    τ::Float64                 # switching time
    q::Float64                 # processing quantum time
    n::Int                     # number of terminals
    reps::Int                  # number of finished jobs to collect
    job_queue::Vector{Job}     # queue of jobs waiting for processing
    finished_jobs::Vector{Job} # list of finished jobs
    job::Union{Job, Nothing}  # current job on computer
    # this structure holds information for reporting
    # for each perior between events it holds the length of queue,
    # if the computer were busy and the length of time period
    state_chunks::Vector{Tuple{Int, Bool, Float64}}
    SharedComputer(a, p, τ, q, n, reps) =
        new(Exponential(a), Exponential(p), τ, q, n, reps,
            Job[], Job[], nothing, Tuple{Int, Bool, Float64}[])
end

function arrival(s)
    ttrue = rand(s.state.p)
    j = Job(s.now, NaN, ttrue, ttrue)
    isnothing(s.state.job) && register!(s, start_compute, 0.0)
    push!(s.state.job_queue, j)
end

function start_compute(s)
    j = popfirst!(s.state.job_queue)
    s.state.job = j
    processed = j.remaining < s.state.q ? j.remaining : s.state.q
    j.remaining -= processed
    register!(s, stop_compute, processed + s.state.τ)
end

function stop_compute(s)
    j = s.state.job
    s.state.job = nothing
    if j.remaining == 0
        j.tout = s.now
        push!(s.state.finished_jobs, j)
        register!(s, arrival, rand(s.state.a))
    else
        push!(s.state.job_queue, j)
    end
    isempty(s.state.job_queue) || register!(s, start_compute, 0.0)
    length(s.state.finished_jobs) == s.state.reps && terminate!(s)
end

monitor(s, Δ) =
    push!(s.state.state_chunks,
          (length(s.state.job_queue), !isnothing(s.state.job), Δ))

function run(a, p, τ, q, n, reps)
    sc = SharedComputer(a, p, τ, q, n, reps)
    s = Scheduler(sc, Float64, monitor)
    for _ in 1:n
        register!(s, arrival, rand(sc.a))
    end
    go!(s)
    chunks = s.state.state_chunks
    (n=n,
     response=mean(j.tout - j.tin for j in s.state.finished_jobs),
     queue=mean(getindex.(chunks, 1), Weights(getindex.(chunks, 3))),
     CPU=mean(getindex.(chunks, 2), Weights(getindex.(chunks, 3))))
end

DataFrame([run(25.0, 0.8, 0.015, 0.1, n, 1000) for n in 10:10:80])


text:
# Create worker processes
nproc = 4 # The desired number of CPUs to run with #Sys.CPU_THREADS
# Clear any existing worker processes if present
Distributed.nprocs() > 1 && Distributed.rmprocs(Distributed.workers())
# Create worker processes
Distributed.addprocs(nproc - 1; exeflags="--project")
# Import functions used on all worker processes
Distributed.@everywhere using OPFLearn
#Distributed.@everywhere import OPFLearn: sample_producer, sample_processor

results = OPFLearn.dist_create_samples("test//data//pglib_opf_case14_ieee.m", 10, nproc=4)


text:
isdefined(Base, :__precompile__) && __precompile__()

module FITSIO

using Compat

export FITS,
       HDU,
       ImageHDU,
       TableHDU,
       ASCIITableHDU,
       FITSHeader,
       read_key,
       read_header,
       get_comment,
       set_comment!,
       copy_section

import Base: getindex,
             setindex!,
             length,
             show,
             read,
             write,
             close,
             ndims,
             size,
             endof,
             haskey,
             keys,
             values,
             start,
             next,
             done

# Libcfitsio submodule
include("libcfitsio.jl")

using .Libcfitsio

# There are a few direct `ccall`s to libcfitsio in this module. For this, we
# need a few non-exported things from Libcfitsio: the shared library handle,
# and a helper function for raising errors. TYPE_FROM_BITPIX is awkwardly
# defined in Libcfitsio, even though it is not used there.
import .Libcfitsio: libcfitsio,
                    fits_assert_ok,
                    TYPE_FROM_BITPIX

# HDU Types
abstract HDU

type ImageHDU <: HDU
    fitsfile::FITSFile
    ext::Int
end

type TableHDU <: HDU
    fitsfile::FITSFile
    ext::Int
end

type ASCIITableHDU <: HDU
    fitsfile::FITSFile
    ext::Int
end

# FITS
#
# The FITS type represents a FITS file. It holds a reference to a
# FITSFile object (basically the low-level CFITSIO pointer). It also
# holds a reference to all of the previously accessed HDU
# objects. This is so that only a single HDU object is created for
# each extension in the file. It also allows a FITS object to tell
# previously created HDUs about events that happen to the file, such
# as deleting extensions. This could be done by, e.g., setting ext=-1
# in the HDU object.
type FITS
    fitsfile::FITSFile
    filename::AbstractString
    mode::AbstractString
    hdus::Dict{Int, HDU}

    function FITS(filename::AbstractString, mode::AbstractString="r")
        f = (mode == "r"                     ? fits_open_file(filename, 0)   :
             mode == "r+" && isfile(filename)? fits_open_file(filename, 1)   :
             mode == "r+"                    ? fits_create_file(filename)    :
             mode == "w"                     ? fits_create_file("!"*filename):
             error("invalid open mode: $mode"))

        new(f, filename, mode, Dict{Int, HDU}())
    end
end

# FITSHeader
#
# An in-memory representation of the header of an HDU. It stores the
# (key, value, comment) information for each card in a header. We
# could almost just use an OrderedDict for this, but we need to store
# comments.
type FITSHeader
    keys::Vector{Compat.ASCIIString}
    values::Vector{Any}
    comments::Vector{Compat.ASCIIString}
    map::Dict{Compat.ASCIIString, Int}

    function FITSHeader(keys::Vector{Compat.ASCIIString}, values::Vector,
                        comments::Vector{Compat.ASCIIString})
        if ((length(keys) != length(values)) ||
            (length(keys) != length(comments)))
            error("keys, values, comments must be same length")
        end
        map = Dict{Compat.ASCIIString, Int64}()
        for i in 1:length(keys)
          map[keys[i]] = i
        end
        new(keys, convert(Vector{Any}, values), comments, map)
    end
end

include("fits.jl")  # FITS methods
include("header.jl")  # FITSHeader methods
include("image.jl")  # ImageHDU methods
include("table.jl")  # TableHDU & ASCIITableHDU methods

function libcfitsio_version()
    # fits_get_version returns a float. e.g., 3.341f0. We parse that
    # into a proper version number. E.g., 3.341 -> v"3.34.1"
    v = convert(Int, round(1000 * fits_get_version()))
    x = div(v, 1000)
    y = div(rem(v, 1000), 10)
    z = rem(v, 10)
    VersionNumber(x, y, z)
end

include("deprecations.jl")

end # module


text:
export message, ctype

Base.eltype(::Type{Type{T}}) where T = T

"""
    ccal(f, R, ts, vals...)
Calls a C-function `f` with return type `R`, argument types `ts` and argument
values `vals`.
"""
function ccal(f, R, ts, vals...)
  AS = Expr(:tuple, ts...)
  @eval ccall($f, $R, $AS, $(vals...))
end

# Doesn't make a lot of sense, but we use this to convert
# Julia UTF-8 strings to objects in Foundation
toobject(o::Object) = o
toobject(c::Class) = c
toobject(p::Ptr) = p

"""
    ctype(t::Type) -> Type

Returns what C type should be used for our Julia ObjC types such as
`Object` and `Selector`.
"""
ctype(x) = x
ctype(o::Type{Object}) = Ptr{Cvoid}
ctype(s::Type{Selector}) = Ptr{Cvoid}
ctype(a::AbstractArray) = map(ctype, a)

const cmsgsend = cglobal(:objc_msgSend)

"""
    message(obj, sel, args...)
Send a message to an objective-c object `obj` using selector `sel`.
"""
function message(obj::Union{Class, Object}, sel::Selector, args...)
  obj = toobject(obj)
  clarse = class(obj)
  m = method(clarse, sel)
  m == C_NULL && error("$clarse does not respond to $sel")
  types = signature(m)
  ctypes = ctype(types)

  args = Any[args...]
  for i = 1:length(args)
    types[i+3] == Object && (args[i] = toobject(args[i]))
  end

  result = ccal(cmsgsend, ctypes[1], tuple(ctypes[2:end]...),
                obj, sel, args...)
  types[1] in (Object, Selector) && return types[1](result)
  return result
end


text:
# writer.jl
# =========
#
# A writer for VCF formatted files.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/GeneticVariation.jl/blob/master/LICENSE

mutable struct Writer{T<:IO} <: BioCore.IO.AbstractWriter
    stream::T
end

"""
    VCF.Writer(output::IO, header::VCF.Header)
Create a data writer of the VCF file format.
# Arguments
* `output`: data sink
* `header`: VCF header object
"""
function Writer(output::IO, header::Header)
    writer = Writer(output)
    write(writer, header)
    return writer
end

function BioCore.IO.stream(writer::Writer)
    return writer.stream
end

function Base.write(writer::Writer, header::Header)
    n = 0
    for metainfo in header.metainfo
        n += write(writer.stream, metainfo, '\n')
    end
    n += write(writer.stream, "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO")
    if !isempty(header.sampleID)
        n += write(writer.stream, "\tFORMAT")
    end
    for id in header.sampleID
        n += write(writer.stream, '\t', id)
    end
    n += write(writer.stream, '\n')
    return n
end

function Base.write(writer::Writer, record::Record)
    return write(writer.stream, record, '\n')
end


text:
using Test

include("luhn.jl")

@testset "invalid" begin
    @testset "single digit strings can not be valid" begin
        @test !luhn("1")
    end

    @testset "A single zero is invalid" begin
        @test !luhn("0")
    end

    @testset "invalid Canadian SIN" begin
        @test !luhn("055 444 286")
    end

    @testset "invalid credit card" begin
        @test !luhn("8273 1232 7352 0569")
    end

    @testset "valid strings with a non-digit added become invalid" begin
        @test !luhn("046a 454 286")
    end

    @testset "invalid long number with an even remainder" begin
        @test !luhn("1 2345 6789 1234 5678 9012")
    end

    @testset "valid strings with a non-digit added at the end become invalid" begin
        @test !luhn("059a")
    end

    @testset "valid strings with punctuation included become invalid" begin
        @test !luhn("055-444-285")
    end

    @testset "valid strings with symbols included become invalid" begin
        @test !luhn("055# 444\$ 285")
    end

    @testset "single zero with space is invalid" begin
        @test !luhn(" 0")
    end

    @testset "using ascii value for non-doubled non-digit isn't allowed" begin
        @test !luhn("055b 444 285")
    end

    @testset "using ascii value for doubled non-digit isn't allowed" begin
        @test !luhn(":9")
    end
end

@testset "valid" begin
    @testset "valid Canadian SIN" begin
        @test luhn("055 444 285")
    end

    @testset "a simple valid SIN that remains valid if reversed" begin
        @test luhn("059")
    end

    @testset "a simple valid SIN that becomes invalid if reversed" begin
        @test luhn("59")
    end

    @testset "valid number with an even number of digits" begin
        @test luhn("095 245 88")
    end

    @testset "valid number with an odd number of spaces" begin
        @test luhn("234 567 891 234")
    end

    @testset "more than a single zero is valid" begin
        @test luhn("0000 0")
    end

    @testset "input digit 9 is correctly converted to output digit 9" begin
        @test luhn("091")
    end
end


text:
# 60x59 med rectangular (uses QR Lanczos)
using Pseudospectra, Test

@testset "Medium Rectangular" begin
    G = Pseudospectra.grcar(60)
    # drop a column to test rectangular case
    A = G[:,1:end-1]
    # for a rectangular mtx, must provide axes
    opts = Dict{Symbol,Any}(:ax => [-1.0,3.0,-3.0,3.0],:npts=>60)

    ps_data = new_matrix(A,opts)
    driver!(ps_data,opts,gs)
    @test iscomputed(ps_data)
end


text:
# Run a single op over XRTArrays
const AnyXLA = Union{XRTArray, XLAScalar, HloToken}

function build_computation(op::HloOp, args::AnyXLA...)
    comp = HloComputationProto(
        name = "comp",
        instructions = HloInstructionProto[ ],
        id = 0
    )
    root = HloInstructionProto(comp, op, map(i->(
        HloInstructionProto(comp, HloParameter(eltype(args[i]), size(args[i]), i-1))
    ), 1:length(args))...)
    comp.root_id = root.id
    pshape = ProgramShape(
        parameters = collect(map(shape, args)),
        result = root.shape
    )
    config = XLAComputationConfig(
        program_shape = pshape
    )
    hlo_module = HloModuleProto(
        name = "test",
        computations = [ comp ],
        entry_computation_name = "op",
        entry_computation_id = 0,
        id = 0,
        host_program_shape = pshape,
    )
    hlo = HloProto(
        hlo_module = hlo_module
    )
    hlo_snap = HloSnapshot(
        hlo = hlo
    )
    xlac = XLAComputation(
        config=config,
        hlo_snapshot = hlo_snap
    )
end

function _execute(op::HloOp, args...)
    sess = nothing
    for x in args
        if x.storage.remotestorage !== nothing
            sess = x.storage.remotestorage.sess
            break
        end
    end
    if sess === nothing
        ret = emulate(op, args...)
    else
        xrt = XLA.compile(sess, build_computation(op, args...))
        ret = run(xrt, args...)
    end
    ret
end

function execute(op::HloOp, args...)
    # This acts as a temporary inference barrier. We don't benefit from looking
    # through this (since we use infer_rt to figure out the return type) and when
    # aggressive optimizations are on, compiling the execute code can be *very*
    # expensive.
    Base.invokelatest(_execute, op, args...)::infer_rt(op, map(typeof, args)...)
end

@noinline (op::GenericHloOp)(args::AnyXLA...) = execute(op, args...)
@noinline (op::HloDot)(args::AnyXLA...) = execute(op, args...)
@noinline (op::HloReshape)(args::AnyXLA...) = execute(op, args...)
@noinline (op::HloBroadcast)(args::AnyXLA...) = execute(op, args...)
@noinline (op::HloConv)(args::AnyXLA...) = execute(op, args...)
@noinline (op::HloSlice)(args::AnyXLA...) = execute(op, args...)
@noinline (op::HloRng)(args::AnyXLA...) = execute(op, args...)
@noinline (op::HloTranspose)(args::AnyXLA...) = execute(op, args...)
@noinline (op::HloRev)(args::AnyXLA...) = execute(op, args...)
@noinline (op::HloDynamicSlice)(args::AnyXLA...) = execute(op, args...)
@noinline (op::HloDynamicUpdateSlice)(args::AnyXLA...) = execute(op, args...)
@noinline (op::HloConcatenate)(args::AnyXLA...) = execute(op, args...)
@noinline (op::HloInfeed)(args::AnyXLA...) = execute(op, args...)
@noinline (op::HloOutfeed)(args...) = execute(op, args...)
@noinline (op::HloAfterAll)(args::AnyXLA...) = execute(op, args...)
@noinline (op::HloIota)(args::AnyXLA...) = execute(op, args...)
@noinline (op::HloPad)(args::AnyXLA...) = execute(op, args...)
@noinline (op::HloGather)(args::AnyXLA...) = execute(op, args...)

# This function is invoked via invokelatest which acts as an inference barrier.
# Thus statically, we get the type given by `infer_rt`, while dynamically we get
# this nice and informative error. TODO: Actually implement dynamic semantics for these
function dynamic_not_implemented(op)
    error("Dynamic semantics not implemented for HLO operation of type $(typeof(op))\n"*
          "Try compiled mode or implement dynamic semantics.")
end

@noinline function (m::HloReduceWindow{fT})(f::fT, arg::XRTArray, init::XRTArray) where {fT}
    Base.invokelatest(dynamic_not_implemented, m)::infer_rt(m, typeof(f), typeof(arg))
end

@noinline function (m::HloReduce{fT})(f::fT, arg::XRTArray, init::XRTArray) where {fT}
    Base.invokelatest(dynamic_not_implemented, m)::infer_rt(m, typeof(f), typeof(arg), typeof(init))
end

@noinline function (m::HloSort{fT})(f::fT, args...) where {fT}
    rT = infer_rt(m, fT, map(typeof, args)...)
    Base.invokelatest(dynamic_not_implemented, m)::rT
end

@noinline function (m::HloSelectAndScatter{T,S})(select::T, scatter::S, op::XRTArray, source::XRTArray, init::XRTArray) where {T,S}
    Base.invokelatest(dynamic_not_implemented, m)::infer_rt(m, T, S, typeof(op), typeof(source), typeof(init))
end

@noinline function (m::HloScatter{T})(update::T, args...) where {T}
    Base.invokelatest(dynamic_not_implemented, m)::infer_rt(m, T, map(typeof, args)...)
end

@noinline function (m::HloMap{fT})(f::fT, args::XRTArray...) where {fT}
    Base.invokelatest(dynamic_not_implemented, m)::infer_rt(m, f, map(typeof, args)...)
end

@noinline function (m::HloCrossReplicaSum)(f, args::AnyXLA...)
    xargs = args
    Base.invokelatest(dynamic_not_implemented, m, args)::infer_rt(m, Core.Typeof(f), typeof(args).parameters...)
end


text:
"""
### Component <: Servable
name::String
f::Function
properties::Dict
------------------
- name::String - The name field is the way that a component is denoted in code.
- f::Function - The function that gets called with the Connection as an
argument.
- properties::Dict - A dictionary of symbols and values.
------------------
##### constructors
Component(name::String, tag::String, properties::Dict)
"""
mutable struct Component <: Servable
    name::String
    f::Function
    properties::Dict{Any, Any}
    function Component(name::String = "", tag::String = "",
         properties::Dict = Dict{Any, Any}())
         properties[:children] = Vector{Any}()
         f(c::Connection) = begin
             open_tag::String = "<$tag id = $name "
             text::String = ""
             write!(c, open_tag)
             for property in keys(properties)
                 special_keys = [:text, :children]
                 if ~(property in special_keys)
                     prop::String = string(properties[property])
                     propkey::String = string(property)
                     write!(c, " $propkey = $prop")
                 else
                     if property == :text
                         text = properties[property]
                     end
                 end
             end
             write!(c, ">")
             if :children in keys(properties)
                 [write!(c, s) for s in properties[:children]]
            end
            write!(c, "$text</$tag>")
         end
         new(name, f, properties)::Component
    end
    Component(name::String, tag::String, props::Base.Pairs) = begin
        Component(name, tag, Dict{Any, Any}(props))
    end
end
#==
Base
    Components
==#
img(name::String = ""; args ...) = Component(name, "img", args)::Component
link(name::String = ""; args ...) = Component(name, "link", args)::Component
meta(name::String = ""; args ...) = Component(name, "meta", args)::Component
input(name::String = ""; args ...) = Component(name, "input", args)::Component
a(name::String = ""; args ...) = Component(name, "a", args)::Component
p(name::String = ""; args ...) = Component(name, "p", args)::Component
h(name::String = "", n::Int64 = 1; args ...) = Component(name, "h$n", args)::Component
ul(name::String = ""; args ...) = Component(name, "ul", args)::Component
li(name::String = ""; args ...) = Component(name, "li", args)::Component
divider(name::String = ""; args ...) = Component(name, "div", args)::Component
br(name::String = ""; args ...) = Component(name, "/br", args)::Component
i(name::String = ""; args ...) = Component(name, "i", args)::Component
title(name::String = ""; args ...) = Component(name, "title", args)::Component
span(name::String = ""; args ...) = Component(name, "span", args)::Component
iframe(name::String = ""; args ...) = Component(name, "iframe", args)::Component
svg(name::String = ""; args ...) = Component(name, "svg", args)::Component
element(name::String = ""; args ...) = Component(name, "element", args)::Component
label(name::String = ""; args ...) = Component(name, "label", args)::Component
script(name::String = ""; args ...) = Component(name, "script", args)::Component
nav(name::String = ""; args ...) = Component(name, "nav", args)::Component
button(name::String = ""; args ...) = Component(name, "button", args)::Component
form(name::String = ""; args ...) = Component(name, "form", args)::Component
#==
Style
    Components
    ==#
mutable struct Animation <: StyleComponent
    name::String
    keyframes::Dict
    f::Function
    delay::Float64
    length::Float64
    iterations::Integer
    function Animation(name::String = "animation"; delay::Float64 = 0.0,
        length::Float64 = 5.2, iterations::Integer = 1)
        f(c::Connection) = begin
            s::String = "<style> @keyframes $name {"
            for anim in keys(keyframes)
                vals = keyframes[anim]
                s = s * "$anim {" * vals * "}"
            end
            write!(c, string(s * "}</style>"))
        end
        f() = begin
            s::String = "<style> @keyframes $name {"
            for anim in keys(keyframes)
                vals = keyframes[anim]
                s = s * "$anim {" * vals * "}"
            end
            string(s * "}</style>")::String
        end
        keyframes::Dict = Dict()
        new(name, keyframes, f, delay, length, iterations)
    end
end

mutable struct Style <: StyleComponent
    name::String
    f::Function
    properties::Dict{Any, Any}
    extras::String
    function Style(name::String; props ...)
        properties::Dict = Dict{Any, Any}(props)
        extras::String = ""
        f(c::Connection) = begin
            css = "<style>.$name { "
            for rule in keys(properties)
                property = string(rule)
                value = string(properties[rule])
                css = css * "$property: $value; "
            end
            css = css * "}</style>" * extras
            write!(c, css)
        end
        new(name::String, f::Function, properties::Dict, extras)::Style
    end
end


text:
push!(LOAD_PATH, joinpath(homedir(),"GitHub/KenyaCoV/src"))
using Plots,Parameters,Distributions,DifferentialEquations,JLD2,DataFrames,StatsPlots,FileIO,MAT,RecursiveArrayTools
using Revise
import KenyaCoV
using LinearAlgebra:eigen
using Statistics: median, quantile

gr()#Plotting frontend
"""
IMPERFECT social distancing
All households reduce contact outside household or workplace by 75%.
School contact rates set to zero, workplace contact rates reduced by 25%.
Household contact rates assumed to increase by 25%.
"""

"""
Load age structured data, define initial state and declare the KenyaCoV problem
"""

u0,P,P_dest = KenyaCoV.model_ingredients_from_data("data/data_for_age_structuredmodel.jld2",
                                            "data/flight_numbers.csv",
                                            "data/projected_global_prevelance.csv")

@load "data/detection_rates_for_different_taus.jld2" d_0 d_01 d_025 d_05 d_1
@load "data/susceptibility_rates.jld2" σ

#Initial infecteds

"""
Declare the treatment/isolation rates considered, and a callback that limits intervention
"""
treatment_rates = [(0.,1),(0.,0.5),(1/7.,1.),(1/7,0.5),(1/3.5,1.),(1/3.5,0.5)]

function isolation_limit(u,t,integrator) # Isolation can only continue until 1000 detected cases
  integrator.p.isolating_detecteds && sum(u[:,:,8]) > 1e3
end
function affect_isolation_limit!(integrator)
  integrator.p.τ = 0
  integrator.p.isolating_detecteds = false
end
cb_iso_limit = DiscreteCallback(isolation_limit,affect_isolation_limit!)

"""
Load age mixing matrices (these are all in to (row) from (col) format)
"""

@load "data/agemixingmatrix_Kenya_all_types.jld2" M_Kenya M_Kenya_ho M_Kenya_other M_Kenya_school M_Kenya_work
reduced_treatment_rates = [(0.,1),(1/3.5,0.5)]

"""
SCENARIO 1 --- Imperfect social distancing for 3 months:
* No age-specific susceptibilties --- its disease difference
* Asymptomatics are 0% as infectious as symptomatics
"""

u0,P,P_dest = KenyaCoV.model_ingredients_from_data("data/data_for_age_structuredmodel.jld2",
                                            "data/flight_numbers.csv",
                                            "data/projected_global_prevelance.csv")
#Redistribute susceptibility JUST to rescale infectiousness so we get the correct R₀/r
P.χ = ones(KenyaCoV.n_a)
P.rel_detection_rate = d_0
P.dt = 0.25;
P.ext_inf_rate = 0.;
P.ϵ = 0.
P.γ = 1/2.5
R₀_scale = KenyaCoV.calculate_R₀_scale(P)
P.χ = ones(KenyaCoV.n_a)/R₀_scale

P.σ = 1/rand(KenyaCoV.d_incubation)
P.β = 2.5*P.γ

#Choosing initial conditions - scenario 2
rand_chosen_R,eq_case_distribution,fifth_gen_case_distribution = KenyaCoV.calculate_R₀(P)
num_initial_detected = ceil(Int64,sum(fifth_gen_case_distribution.*d_01))
detected_age_group = argmax(fifth_gen_case_distribution.*d_01)
undetected_age_profile = round.(Int64,fifth_gen_case_distribution.*(1 .- d_01))

#Implement social distancing
KenyaCoV.calculate_R₀(P)
P.M = M_Kenya_ho*(1.25) .+ M_Kenya_work*(0.75) .+ M_Kenya_other*(0.25)
KenyaCoV.calculate_R₀(P)

u0[KenyaCoV.ind_nairobi_as,detected_age_group,4] = num_initial_detected
u0[KenyaCoV.ind_nairobi_as,:,3] = undetected_age_profile
prob = KenyaCoV.create_KenyaCoV_non_neg_prob(u0,(0.,90.),P)

results_1SI_short = KenyaCoV.run_scenario(P,prob,1000,reduced_treatment_rates,cb_iso_limit)
@save joinpath(homedir(),"Github/KenyaCoVOutputs/results_1SI_short.jld2") results_1SI_short

println("Finished 1 imperfect social distancing for 3 months")

"""
SCENARIO 2 --- Imperfect social distancing:
* No age-specific susceptibilties --- its disease difference
* Asymptomatics are 10% as infectious as symptomatics
"""

u0,P,P_dest = KenyaCoV.model_ingredients_from_data("data/data_for_age_structuredmodel.jld2",
                                            "data/flight_numbers.csv",
                                            "data/projected_global_prevelance.csv")
#Redistribute susceptibility JUST to rescale infectiousness so we get the correct R₀/r
P.χ = ones(KenyaCoV.n_a)
P.rel_detection_rate = d_01
P.dt = 0.25;
P.ext_inf_rate = 0.;
P.ϵ = 0.1
P.γ = 1/2.5
R₀_scale = KenyaCoV.calculate_R₀_scale(P)
P.χ = ones(KenyaCoV.n_a)/R₀_scale

P.σ = 1/rand(KenyaCoV.d_incubation)
P.β = 2.5*P.γ

#Choosing initial conditions - scenario 2
rand_chosen_R,eq_case_distribution,fifth_gen_case_distribution = KenyaCoV.calculate_R₀(P)
num_initial_detected = ceil(Int64,sum(fifth_gen_case_distribution.*d_01))
detected_age_group = argmax(fifth_gen_case_distribution.*d_01)
undetected_age_profile = round.(Int64,fifth_gen_case_distribution.*(1 .- d_01))

#Implement social distancing
KenyaCoV.calculate_R₀(P)
P.M = M_Kenya_ho*(1.25) .+ M_Kenya_work*(0.75) .+ M_Kenya_other*(0.25)
KenyaCoV.calculate_R₀(P)

u0[KenyaCoV.ind_nairobi_as,detected_age_group,4] = num_initial_detected
u0[KenyaCoV.ind_nairobi_as,:,3] = undetected_age_profile
prob = KenyaCoV.create_KenyaCoV_non_neg_prob(u0,(0.,90.),P)

results_2SI_short = KenyaCoV.run_scenario(P,prob,1000,reduced_treatment_rates,cb_iso_limit)
@save joinpath(homedir(),"Github/KenyaCoVOutputs/results_2SI_short.jld2") results_2SI_short

println("Finished 2 imperfect social distancing for 3 months")
"""
SCENARIO 3 --- Social distancing:
* No age-specific susceptibilties --- its disease difference
* Asymptomatics are 25% as infectious as symptomatics
"""
u0,P,P_dest = KenyaCoV.model_ingredients_from_data("data/data_for_age_structuredmodel.jld2",
                                            "data/flight_numbers.csv",
                                            "data/projected_global_prevelance.csv")
#Redistribute susceptibility JUST to rescale infectiousness so we get the correct R₀/r
P.χ = ones(KenyaCoV.n_a)
P.rel_detection_rate = d_025
P.dt = 0.25;
P.ext_inf_rate = 0.;
P.ϵ = 0.25
P.γ = 1/2.5
R₀_scale = KenyaCoV.calculate_R₀_scale(P)
P.χ = ones(KenyaCoV.n_a)/R₀_scale

P.σ = 1/rand(KenyaCoV.d_incubation)
P.β = 2.5*P.γ

#Choosing initial conditions - scenario 3
rand_chosen_R,eq_case_distribution,fifth_gen_case_distribution = KenyaCoV.calculate_R₀(P)
num_initial_detected = ceil(Int64,sum(fifth_gen_case_distribution.*d_025))
detected_age_group = argmax(fifth_gen_case_distribution.*d_025)
undetected_age_profile = round.(Int64,fifth_gen_case_distribution.*(1 .- d_025))

#Implement social distancing
KenyaCoV.calculate_R₀(P)
P.M = M_Kenya_ho*(1.25) .+ M_Kenya_work*(0.75) .+ M_Kenya_other*(0.25)
KenyaCoV.calculate_R₀(P)

u0[KenyaCoV.ind_nairobi_as,detected_age_group,4] = num_initial_detected
u0[KenyaCoV.ind_nairobi_as,:,3] = undetected_age_profile
prob = KenyaCoV.create_KenyaCoV_non_neg_prob(u0,(0.,90.),P)

results_3SI_short = KenyaCoV.run_scenario(P,prob,1000,reduced_treatment_rates,cb_iso_limit)
@save joinpath(homedir(),"Github/KenyaCoVOutputs/results_3SI_short.jld2") results_3SI_short

println("Finished 3 imperfect social distancing for 3 months")

"""
SCENARIO 4 --- Social distancing for 3 months:
* No age-specific susceptibilties --- its disease difference
* Asymptomatics are 50% as infectious as symptomatics
"""
u0,P,P_dest = KenyaCoV.model_ingredients_from_data("data/data_for_age_structuredmodel.jld2",
                                            "data/flight_numbers.csv",
                                            "data/projected_global_prevelance.csv")
#Redistribute susceptibility JUST to rescale infectiousness so we get the correct R₀/r
P.χ = ones(KenyaCoV.n_a)
P.rel_detection_rate = d_05
P.dt = 0.25;
P.ext_inf_rate = 0.;
P.ϵ = 0.5
P.γ = 1/2.5
R₀_scale = KenyaCoV.calculate_R₀_scale(P)
P.χ = ones(KenyaCoV.n_a)/R₀_scale

P.σ = 1/rand(KenyaCoV.d_incubation)
P.β = 2.5*P.γ

#Choosing initial conditions - scenario 4
rand_chosen_R,eq_case_distribution,fifth_gen_case_distribution = KenyaCoV.calculate_R₀(P)
num_initial_detected = ceil(Int64,sum(fifth_gen_case_distribution.*d_05))
detected_age_group = argmax(fifth_gen_case_distribution.*d_05)
undetected_age_profile = round.(Int64,fifth_gen_case_distribution.*(1 .- d_05))

#Implement social distancing
KenyaCoV.calculate_R₀(P)
P.M = M_Kenya_ho*(1.25) .+ M_Kenya_work*(0.75) .+ M_Kenya_other*(0.25)
KenyaCoV.calculate_R₀(P)

u0[KenyaCoV.ind_nairobi_as,detected_age_group,4] = num_initial_detected
u0[KenyaCoV.ind_nairobi_as,:,3] = undetected_age_profile
prob = KenyaCoV.create_KenyaCoV_non_neg_prob(u0,(0.,90.),P)

results_4SI_short = KenyaCoV.run_scenario(P,prob,1000,reduced_treatment_rates,cb_iso_limit)
@save joinpath(homedir(),"Github/KenyaCoVOutputs/results_4SI_short.jld2") results_4SI_short

println("Finished 4 imperfect social distancing for 3 months")

"""
SCENARIO 5 --- Social distancing for 3 months:
* No age-specific susceptibilties --- its disease difference
* Asymptomatics are 100% as infectious as symptomatics
"""
u0,P,P_dest = KenyaCoV.model_ingredients_from_data("data/data_for_age_structuredmodel.jld2",
                                            "data/flight_numbers.csv",
                                            "data/projected_global_prevelance.csv")
#Redistribute susceptibility JUST to rescale infectiousness so we get the correct R₀/r
P.χ = ones(KenyaCoV.n_a)
P.rel_detection_rate = d_1
P.dt = 0.25;
P.ext_inf_rate = 0.;
P.ϵ = 1.
P.γ = 1/2.5
R₀_scale = KenyaCoV.calculate_R₀_scale(P)
P.χ = ones(KenyaCoV.n_a)/R₀_scale

P.σ = 1/rand(KenyaCoV.d_incubation)
P.β = 2.5*P.γ

#Choosing initial conditions - scenario 4
rand_chosen_R,eq_case_distribution,fifth_gen_case_distribution = KenyaCoV.calculate_R₀(P)
num_initial_detected = ceil(Int64,sum(fifth_gen_case_distribution.*d_1))
detected_age_group = argmax(fifth_gen_case_distribution.*d_1)
undetected_age_profile = round.(Int64,fifth_gen_case_distribution.*(1 .- d_1))

#Implement social distancing
KenyaCoV.calculate_R₀(P)
P.M = M_Kenya_ho*(1.25) .+ M_Kenya_work*(0.75) .+ M_Kenya_other*(0.25)
KenyaCoV.calculate_R₀(P)

u0[KenyaCoV.ind_nairobi_as,detected_age_group,4] = num_initial_detected
u0[KenyaCoV.ind_nairobi_as,:,3] = undetected_age_profile
prob = KenyaCoV.create_KenyaCoV_non_neg_prob(u0,(0.,90.),P)

results_5SI_short = KenyaCoV.run_scenario(P,prob,1000,reduced_treatment_rates,cb_iso_limit)
@save joinpath(homedir(),"Github/KenyaCoVOutputs/results_5SI_short.jld2") results_5SI_short

println("Finished 5 imperfect social distancing for 3 months")


text:
export local_coordinate_coding

import ..LocalCoordinateCoding

using mlpack._Internal.io

import mlpack_jll
const local_coordinate_codingLibrary = mlpack_jll.libmlpack_julia_local_coordinate_coding

# Call the C binding of the mlpack local_coordinate_coding binding.
function local_coordinate_coding_mlpackMain()
  success = ccall((:local_coordinate_coding, local_coordinate_codingLibrary), Bool, ())
  if !success
    # Throw an exception---false means there was a C++ exception.
    throw(ErrorException("mlpack binding error; see output"))
  end
end

" Internal module to hold utility functions. "
module local_coordinate_coding_internal
  import ..local_coordinate_codingLibrary

import ...LocalCoordinateCoding

# Get the value of a model pointer parameter of type LocalCoordinateCoding.
function IOGetParamLocalCoordinateCoding(paramName::String)::LocalCoordinateCoding
  LocalCoordinateCoding(ccall((:IO_GetParamLocalCoordinateCodingPtr, local_coordinate_codingLibrary), Ptr{Nothing}, (Cstring,), paramName))
end

# Set the value of a model pointer parameter of type LocalCoordinateCoding.
function IOSetParamLocalCoordinateCoding(paramName::String, model::LocalCoordinateCoding)
  ccall((:IO_SetParamLocalCoordinateCodingPtr, local_coordinate_codingLibrary), Nothing, (Cstring, Ptr{Nothing}), paramName, model.ptr)
end

# Serialize a model to the given stream.
function serializeLocalCoordinateCoding(stream::IO, model::LocalCoordinateCoding)
  buf_len = UInt[0]
  buf_ptr = ccall((:SerializeLocalCoordinateCodingPtr, local_coordinate_codingLibrary), Ptr{UInt8}, (Ptr{Nothing}, Ptr{UInt}), model.ptr, Base.pointer(buf_len))
  buf = Base.unsafe_wrap(Vector{UInt8}, buf_ptr, buf_len[1]; own=true)
  write(stream, buf)
end
# Deserialize a model from the given stream.
function deserializeLocalCoordinateCoding(stream::IO)::LocalCoordinateCoding
  buffer = read(stream)
  LocalCoordinateCoding(ccall((:DeserializeLocalCoordinateCodingPtr, local_coordinate_codingLibrary), Ptr{Nothing}, (Ptr{UInt8}, UInt), Base.pointer(buffer), length(buffer)))
end
end # module

"""
    local_coordinate_coding(; [atoms, initial_dictionary, input_model, lambda, max_iterations, normalize, seed, test, tolerance, training, verbose])

An implementation of Local Coordinate Coding (LCC), which codes data that
approximately lives on a manifold using a variation of l1-norm regularized
sparse coding.  Given a dense data matrix X with n points and d dimensions, LCC
seeks to find a dense dictionary matrix D with k atoms in d dimensions, and a
coding matrix Z with n points in k dimensions.  Because of the regularization
method used, the atoms in D should lie close to the manifold on which the data
points lie.

The original data matrix X can then be reconstructed as D * Z.  Therefore, this
program finds a representation of each point in X as a sparse linear combination
of atoms in the dictionary D.

The coding is found with an algorithm which alternates between a dictionary
step, which updates the dictionary D, and a coding step, which updates the
coding matrix Z.

To run this program, the input matrix X must be specified (with -i), along with
the number of atoms in the dictionary (-k).  An initial dictionary may also be
specified with the `initial_dictionary` parameter.  The l1-norm regularization
parameter is specified with the `lambda` parameter.

For example, to run LCC on the dataset `data` using 200 atoms and an
l1-regularization parameter of 0.1, saving the dictionary `dictionary` and the
codes into `codes`, use

```julia
julia> using CSV
julia> data = CSV.read("data.csv")
julia> codes, dict, _ = local_coordinate_coding(atoms=200,
            lambda=0.1, training=data)
```

The maximum number of iterations may be specified with the `max_iterations`
parameter. Optionally, the input data matrix X can be normalized before coding
with the `normalize` parameter.

An LCC model may be saved using the `output_model` output parameter.  Then, to
encode new points from the dataset `points` with the previously saved model
`lcc_model`, saving the new codes to `new_codes`, the following command can be
used:

```julia
julia> using CSV
julia> points = CSV.read("points.csv")
julia> new_codes, _, _ =
            local_coordinate_coding(input_model=lcc_model, test=points)
```

# Arguments

 - `atoms::Int`: Number of atoms in the dictionary.  Default value `0`.

 - `initial_dictionary::Array{Float64, 2}`: Optional initial dictionary.
 - `input_model::LocalCoordinateCoding`: Input LCC model.
 - `lambda::Float64`: Weighted l1-norm regularization parameter.  Default
      value `0`.

 - `max_iterations::Int`: Maximum number of iterations for LCC (0
      indicates no limit).  Default value `0`.

 - `normalize::Bool`: If set, the input data matrix will be normalized
      before coding.  Default value `false`.

 - `seed::Int`: Random seed.  If 0, 'std::time(NULL)' is used.  Default
      value `0`.

 - `test::Array{Float64, 2}`: Test points to encode.
 - `tolerance::Float64`: Tolerance for objective function.  Default value
      `0.01`.

 - `training::Array{Float64, 2}`: Matrix of training data (X).
 - `verbose::Bool`: Display informational messages and the full list of
      parameters and timers at the end of execution.  Default value `false`.

# Return values

 - `codes::Array{Float64, 2}`: Output codes matrix.
 - `dictionary::Array{Float64, 2}`: Output dictionary matrix.
 - `output_model::LocalCoordinateCoding`: Output for trained LCC model.

"""
function local_coordinate_coding(;
                                 atoms::Union{Int, Missing} = missing,
                                 initial_dictionary = missing,
                                 input_model::Union{LocalCoordinateCoding, Missing} = missing,
                                 lambda::Union{Float64, Missing} = missing,
                                 max_iterations::Union{Int, Missing} = missing,
                                 normalize::Union{Bool, Missing} = missing,
                                 seed::Union{Int, Missing} = missing,
                                 test = missing,
                                 tolerance::Union{Float64, Missing} = missing,
                                 training = missing,
                                 verbose::Union{Bool, Missing} = missing,
                                 points_are_rows::Bool = true)
  # Force the symbols to load.
  ccall((:loadSymbols, local_coordinate_codingLibrary), Nothing, ());

  IORestoreSettings("Local Coordinate Coding")

  # Process each input argument before calling mlpackMain().
  if !ismissing(atoms)
    IOSetParam("atoms", convert(Int, atoms))
  end
  if !ismissing(initial_dictionary)
    IOSetParamMat("initial_dictionary", initial_dictionary, points_are_rows)
  end
  if !ismissing(input_model)
    local_coordinate_coding_internal.IOSetParamLocalCoordinateCoding("input_model", convert(LocalCoordinateCoding, input_model))
  end
  if !ismissing(lambda)
    IOSetParam("lambda", convert(Float64, lambda))
  end
  if !ismissing(max_iterations)
    IOSetParam("max_iterations", convert(Int, max_iterations))
  end
  if !ismissing(normalize)
    IOSetParam("normalize", convert(Bool, normalize))
  end
  if !ismissing(seed)
    IOSetParam("seed", convert(Int, seed))
  end
  if !ismissing(test)
    IOSetParamMat("test", test, points_are_rows)
  end
  if !ismissing(tolerance)
    IOSetParam("tolerance", convert(Float64, tolerance))
  end
  if !ismissing(training)
    IOSetParamMat("training", training, points_are_rows)
  end
  if verbose !== nothing && verbose === true
    IOEnableVerbose()
  else
    IODisableVerbose()
  end

  IOSetPassed("codes")
  IOSetPassed("dictionary")
  IOSetPassed("output_model")
  # Call the program.
  local_coordinate_coding_mlpackMain()

  return IOGetParamMat("codes", points_are_rows),
         IOGetParamMat("dictionary", points_are_rows),
         local_coordinate_coding_internal.IOGetParamLocalCoordinateCoding("output_model")
end


text:
"""
    gradient(f, p::AbstractParticles)

Calculate the gradient of `f` in `p`. This corresponds to a smoothed finite-difference approximation where the smoothing kernel is given by the distribution of `p`.
Return mean and std.
"""
function gradient(f,p::MonteCarloMeasurements.AbstractParticles)
    r = 2(p\f(p))
    pmean(r), pstd(r)
end
function gradient(f,p::Union{Integer, AbstractFloat})
    p = p  ± 0.000001
    r = 2(p\f(p))
    pmean(r)
end

function gradient(f::Function,p::MonteCarloMeasurements.MvParticles)
    r = (p-pmean.(p))\(f(p) - f(pmean.(p)))
end

"""
    jacobian(f::Function, p::Vector{<:AbstractParticles})

Calculate the Jacobian of `f` in `p`. This corresponds to a smoothed finite-difference approximation where the smoothing kernel is given by the distribution of `p`.
"""
function jacobian(f::Function,p::MonteCarloMeasurements.MvParticles)
    r = (p-pmean.(p))\(f(p) - f(pmean.(p)))
end


text:
module StanRun

using ArgCheck
using StanDump

import StanDump: standump

"""
    find_executable_dir(name)

Searching the directories in `ENV["PATH"]`, return the directory of
the first executable file with name `name`, or `nothing` if not found.

NOTE: Testing for executables is heuristic.
"""
function find_executable_dir(name)
    for dir in split(get(ENV, "PATH", ""), @static is_windows() ? ";" : ":")
        if isdir(dir) && name ∈ readdir(dir)
            path = joinpath(dir, name)
            if ((uperm(path) | gperm(path) | operm(path)) & 1) ≠ 0
                return dir
            end
        end
    end
    nothing
end

"""
   find_cmdstan_home()

Try to find the home directory of `CmdStan`, either from
`ENV[CMDSTAN_HOME]` or falling back to searching for the `stanc`
executable. Error when not found.
"""
function find_cmdstan_home()
    envkey = "CMDSTAN_HOME"
    if haskey(ENV, envkey)
        ENV[envkey]
    else
        dir = find_executable_dir("stanc")
        if dir == nothing
            error("stanc not found in path, could not determine CMDSTAN_HOME")
        end
        realpath(joinpath(dir, ".."))
    end
end

"""
A Stan program, represented by its source file. Also contains
information on the CmdStan home directory.
"""
struct Program
    cmdstan_home::String
    program_file::String
    function Program(program_file; cmdstan_home = find_cmdstan_home())
        @argcheck isabspath(cmdstan_home) "Need an absolute path for CMDSTAN_HOME."
        @argcheck isabspath(program_file) "Need an absolute path for the program file."
        ext = splitext(basename(program_file))[2]
        if !isempty(ext)
            warn("Program path $(program_file) contains an extension $(ext).")
        end
        new(cmdstan_home, program_file)
    end
end

# constants for various parts
"The `stanc` executable."
const STANC = Val{:stanc}

"The `stansummary` executable."
const STANSUMMARY = Val{:stansummary}

"Stan source code."
const SOURCE = Val{:source}

"Data file."
const DATA = Val{:data}

"Compiled executable."
const EXECUTABLE = Val{:executable}

"Regular expression for sample files (only the filename part)."
const SAMPLEBASERX = Val{:samplebaserx}

"Directory for all the program-specific files."
const DIR = Val{:dir}

"""
    Samples (MCMC draws). The chain ID is appended to the end.
"""
struct Samples
    id::Int
end

"""
    getpath(sp, part)

Return the path for the given `part`. Path does not necessarily point
to an existing file.
"""
getpath(sp::Program, ::Type{STANC}) =
    joinpath(sp.cmdstan_home, "bin/stanc")
getpath(sp::Program, ::Type{STANSUMMARY}) =
    joinpath(sp.cmdstan_home, "bin/stansummary")
getpath(sp::Program, ::Type{DIR}) = realpath(dirname(sp.program_file))
getpath(sp::Program, ::Type{SOURCE}) = sp.program_file * ".stan"
getpath(sp::Program, ::Type{DATA}) = sp.program_file * ".data.R"
getpath(sp::Program, ::Type{EXECUTABLE}) = sp.program_file
getpath(sp::Program, samples::Samples) =
    sp.program_file * "-samples-" * string(samples.id) * ".csv"
getpath(sp::Program, ::Type{SAMPLEBASERX}) =
    Regex(basename(sp.program_file) * "-samples-([[:digit:]]+).csv")

"""
    getparents(part)

Return the part identifier for the other parts that are required to make it.
"""
getparents(::Any) = ()
getparents(::Type{EXECUTABLE}) = (STANC, SOURCE)
getparents(::Samples) = (EXECUTABLE, DATA)

"""
    make(sp, part; force = false)

Make `part` of a Stan program, conditional on parent parts being more
recent. Traverses dependencies recursively, using `mtime`s, or
unconditionally when `force`. Return the `mtime` of the part.

This function implements the traversal algorithm, `_make` does the
actual work.
"""
function make(sp::Program, part; force = false)
    parents = getparents(part)
    parent_timestamps = make.(sp, parents; force = force)
    path = getpath(sp, part)
    if isempty(parent_timestamps)
        @assert isfile(path) "Could not find $(part) at $(getpath(sp, part))."
    else
        if force || !isfile(path) || (mtime(path) < maximum(parent_timestamps))
            _make(sp, part)
        end
    end
    mtime(path)
end

"""
    _make(sp, part)

Make `part` of a Stan program. Assumes that parent parts are available.
"""
function _make(sp::Program, ::Type{EXECUTABLE})
    executable_path = getpath(sp, EXECUTABLE)
    run(Cmd(`make $(executable_path)`; dir = sp.cmdstan_home))
end

function _make(sp::Program, samples::Samples)
    executable_path = getpath(sp, EXECUTABLE)
    data_path = getpath(sp, DATA)
    samples_path = getpath(sp, samples)
    run(`$(executable_path) sample data file=$(data_path) output file=$(samples_path)`)
end

function standump(sp::Program, data...)
    open(io -> standump(io, data...), getpath(sp, DATA), "w+")
end

"""
    sample(sp, id; force = false)

Sample samples with `id` from a Stan program `sp`. When `force`,
regenerate samples.
"""
sample(sp::Program, id; force = false) = make(sp, Samples(id); force = force)

"""
    sample_ids(sp)

Return the samples found associated with Stan program `sp`. Samples
are not checked to be more recent than either the data or the source,
use `sample` to enforce that.
"""
function sample_ids(sp::Program)
    dir = getpath(sp, DIR)
    files = filter(f -> isfile(joinpath(dir, f)), readdir(dir))
    matches = match.(getpath(sp, SAMPLEBASERX), files)
    [parse(Int, m.captures[1]) for m in matches if m != nothing]
end

end # module


text:
@testset "Utility" begin

   duex = similar(tmpv.dubV_σ)
   fduex = similar(tmpv.ubVfull)
   roy = 0.2
   geoid = 2
   itype=(geoid, roy,)
   p = prim
   h = peturb(σv)
   σ1 = σv - h
   σ2 = σv + h
   hh = σ2 - σ1
   st = state(2,2,0,0)
   stidx = state_idx(p.wp,2,2,0,0)

   fillflows!(duex,  flow, p, θt, σ1, stidx, itype...)
   fillflows!(fduex, flow, p, θt, σ2, stidx, itype...)
   fduex .-= duex
   fduex ./= hh
   fillflows!(duex, flowdσ, p, θt, σv, stidx, itype...)

   @views maxv, idx = findmax(abs.(duex .- fduex))
   sub = CartesianIndices(duex)[idx]
   @show "worst value is $maxv at $sub for duσ_add"
   @test duex ≈ fduex
end


text:
"""

$(TYPEDEF)

Structures to contain the details of a solute or solvent to
store in the results of the MDDF calculation.

$(TYPEDFIELDS)

"""
struct SolSummary
    natoms::Int
    nmols::Int
    natomspermol::Int
end
SolSummary(s::Selection) = SolSummary(s.natoms, s.nmols, s.natomspermol)


text:
#=

Implementation of Gated Recurrent Neural network as described in
J. Chung, Ç. Gülçehre, K. Cho, and Y. Bengio, “Empirical Evaluation of Gated Recurrent Neural Networks on Sequence Modeling.,”
arXiv, vol. cs.NE, 2014.

=#

export GRU

type GRUParameters{D<:Device, F<:Float} <: Parameters
  device::Device

  # Update gate
  updateweight_h::matrixParameters(D,F)
  updateweight_x::matrixParameters(D,F)

  # Reset gate
  resetweight_h::matrixParameters(D,F)
  resetweight_x::matrixParameters(D,F)
end

type GRU{D<:Device, F<:Float} <: Layer
  parameters::GRUParameters{D,F}

  # State
  output::denseRealMatrix(D,F)
  input_gradient::denseRealMatrix(D,F)

  function GRU(p)
    new(device, weight, bias, array(device, F, 0, 0), array(device, F, 0, 0))
  end
end

@doc doc"Initialize a linear module with the input and output size" ->
function GRU{D<:Device, F<:Float}(d::D, ::Type{F}, inputSize::Int64, outputSize::Int64)
  GRU{D,F}(d, matrixParameters(d, F, inputSize, outputSize), matrixParameters(d, F, 1, outputSize))
end

function forward!{D<:Device, F<:Float}(linear::GRU{D,F}, input::DenseMatrix)
  @ensuresize linear.output, (size(input,1), size(linear.weight.values, 2))

  gemm!('N', 'N', 1., input, linear.weight.values, 0., linear.output)
  broadcast!(+, linear.output, linear.output, linear.bias.values)
end

function compute_inputgradient!{D<:Device, F<:Float}(linear::GRU{D,F}, input::DenseMatrix, gradOutput::DenseMatrix)
  @ensuresize linear.input_gradient, (size(gradOutput,1), size(linear.weight.values, 1))

  gemm!('N', 'T', 1., gradOutput, linear.weight.values, 0., linear.input_gradient)
end

function update_gradient!{D<:Device, F<:Float}(linear::GRU{D,F}, input::DenseMatrix, gradOutput::DenseMatrix, scale::F=1.)
  if !isnull(linear.weight.gradient)
    gemm!('T', 'N', scale, input, gradOutput, 0., get(linear.weight.gradient))
  end

  if !isnull(linear.bias.gradient)
    axpy!(scale, sum(gradOutput, 1), get(linear.bias.gradient))
  end
end


text:
function makesignals1(dtype, nmax)
    x0 = rand(dtype, nmax)
    x_exp = copy(x0)
    return x_exp, x0
end
function makesignals2(dtype, nmax)
    x0 = rand(dtype, nmax)
    y0 = rand(dtype, nmax)
    x_exp = copy(x0)
    return x_exp, x0, y0
end
function makesignals3(dtype, nmax)
    x0 = rand(dtype, nmax)
    y0 = rand(dtype, nmax)
    z0 = rand(dtype, nmax)
    x_exp = copy(x0)
    return x_exp, x0, y0, z0
end

macro printtest1(dtype, nmax, ix, incx, n)
    quote
        nmax=$(esc(nmax))
        ix=$(esc(ix))
        incx=$(esc(incx))
        n=$(esc(n))
        $(esc(dtype))<:Float64 && println("  nmax: $nmax, ix: $ix, incx: $incx, n: $n")
    end
end
macro printtest2(dtype, nmax, ix, incx, iy, incy, n)
    quote
        nmax=$(esc(nmax))
        ix=$(esc(ix))
        incx=$(esc(incx))
        iy=$(esc(iy))
        incy=$(esc(incy))
        n=$(esc(n))
        $(esc(dtype))<:Float64 && println("  nmax: $nmax, ix: $ix, incx: $incx, iy: $iy, incy: $incy, n: $n")
    end
end
macro printtest3(dtype, nmax, ix, incx, iy, incy, iz, incz, n)
    quote
        nmax=$(esc(nmax))
        ix=$(esc(ix))
        incx=$(esc(incx))
        iy=$(esc(iy))
        incy=$(esc(incy))
        iz=$(esc(iz))
        incz=$(esc(incz))
        n=$(esc(n))
        $(esc(dtype))<:Float64 && println("  nmax: $nmax, ix: $ix, incx: $incx, iy: $iy, incy: $incy, iz: $iz, incz: $incz, n: $n")
    end
end

# 1 Array, scalar
macro testfunc_arr1a(func, x_exp, x0, ix, incx, a, n)
    quote
        func = $(esc(func))
        x_exp = $(esc(x_exp))
        x0 = $(esc(x0))
        ix = $(esc(ix))
        incx = $(esc(incx))
        a = $(esc(a))
        n = $(esc(n))

        x = copy(x0)
        func(x,ix,incx,a,n)
        @test_approx_eq x x_exp
        # compare kinds
        # inc1
        if incx == 1
            x = copy(x0)
            func(x,ix,a,n)
            @test_approx_eq x x_exp
        end
    end
end
# 2 Arrays
macro testfunc_arr2(func, x_exp, x0, ix, incx, y0, iy, incy, n)
    quote
        func = $(esc(func))
        x_exp = $(esc(x_exp))
        x0 = $(esc(x0))
        ix = $(esc(ix))
        incx = $(esc(incx))
        y0 = $(esc(y0))
        iy = $(esc(iy))
        incy = $(esc(incy))
        n = $(esc(n))

        x = copy(x0)
        y = copy(y0)
        func(x,ix,incx,y,iy,incy,n)
        @test_approx_eq x x_exp
        @test_approx_eq y y0

        # compare kinds
        # inceq
        if incx == incy
            x = copy(x0)
            y = copy(y0)
            func(x,ix,incx,y,iy,n)
            @test_approx_eq x x_exp
            @test_approx_eq y y0
            # inc1
            if incx == 1
                x = copy(x0)
                y = copy(y0)
                func(x,ix,y,iy,n)
                @test_approx_eq x x_exp
                @test_approx_eq y y0
                # inc1ieq
                if ix == iy
                    x = copy(x0)
                    y = copy(y0)
                    func(x,ix,y,n)
                    @test_approx_eq x x_exp
                    @test_approx_eq y y0
                end
            end
        end
    end
end
# 2 Arrays, scalar
macro testfunc_arr2a(func, x_exp, x0, ix, incx, y0, iy, incy, a, n)
    quote
        func = $(esc(func))
        x_exp = $(esc(x_exp))
        x0 = $(esc(x0))
        ix = $(esc(ix))
        incx = $(esc(incx))
        y0 = $(esc(y0))
        iy = $(esc(iy))
        incy = $(esc(incy))
        a = $(esc(a))
        n = $(esc(n))

        x = copy(x0)
        y = copy(y0)
        func(x,ix,incx,y,iy,incy,a,n)
        @test_approx_eq x x_exp
        @test_approx_eq y y0

        # compare kinds
        # inceq
        if incx == incy
            x = copy(x0)
            y = copy(y0)
            func(x,ix,incx,y,iy,a,n)
            @test_approx_eq x x_exp
            @test_approx_eq y y0
            # inc1
            if incx == 1
                x = copy(x0)
                y = copy(y0)
                func(x,ix,y,iy,a,n)
                @test_approx_eq x x_exp
                @test_approx_eq y y0
                # inc1ieq
                if ix == iy
                    x = copy(x0)
                    y = copy(y0)
                    func(x,ix,y,a,n)
                    @test_approx_eq x x_exp
                    @test_approx_eq y y0
                end
            end
        end
    end
end
# 3 Arrays
macro testfunc_arr3(func, x_exp, x0, ix, incx, y0, iy, incy, z0, iz, incz, n)
    quote
        func = $(esc(func))
        x_exp = $(esc(x_exp))
        x0 = $(esc(x0))
        ix = $(esc(ix))
        incx = $(esc(incx))
        y0 = $(esc(y0))
        iy = $(esc(iy))
        incy = $(esc(incy))
        z0 = $(esc(z0))
        iz = $(esc(iz))
        incz = $(esc(incz))
        n = $(esc(n))

        x = copy(x0)
        y = copy(y0)
        z = copy(z0)
        func(x,ix,incx,y,iy,incy,z,iz,incz,n)
        @test_approx_eq x x_exp
        @test_approx_eq y y0
        @test_approx_eq z z0

        # compare kinds
        # inceq
        if incx == incy
            x = copy(x0)
            y = copy(y0)
            func(x,ix,incx,y,iy,z,iz,n)
            @test_approx_eq x x_exp
            @test_approx_eq y y0
            @test_approx_eq z z0
            # inc1
            if incx == 1
                x = copy(x0)
                y = copy(y0)
                func(x,ix,y,iy,z,iz,n)
                @test_approx_eq x x_exp
                @test_approx_eq y y0
                @test_approx_eq z z0
                # inc1ieq
                if ix == iy
                    x = copy(x0)
                    y = copy(y0)
                    func(x,ix,y,z,n)
                    @test_approx_eq x x_exp
                    @test_approx_eq y y0
                    @test_approx_eq z z0
                end
            end
        end
    end
end
# 3 Arrays, scalar
macro testfunc_arr3a(func, x_exp, x0, ix, incx, y0, iy, incy, z0, iz, incz, a, n)
    quote
        func = $(esc(func))
        x_exp = $(esc(x_exp))
        x0 = $(esc(x0))
        ix = $(esc(ix))
        incx = $(esc(incx))
        y0 = $(esc(y0))
        iy = $(esc(iy))
        incy = $(esc(incy))
        z0 = $(esc(z0))
        iz = $(esc(iz))
        incz = $(esc(incz))
        a = $(esc(a))
        n = $(esc(n))

        x = copy(x0)
        y = copy(y0)
        z = copy(z0)
        func(x,ix,incx,y,iy,incy,z,iz,incz,a,n)
        @test_approx_eq x x_exp
        @test_approx_eq y y0
        @test_approx_eq z z0

        # compare kinds
        # inceq
        if incx == incy
            x = copy(x0)
            y = copy(y0)
            func(x,ix,incx,y,iy,z,iz,a,n)
            @test_approx_eq x x_exp
            @test_approx_eq y y0
            @test_approx_eq z z0
            # inc1
            if incx == 1
                x = copy(x0)
                y = copy(y0)
                func(x,ix,y,iy,z,iz,a,n)
                @test_approx_eq x x_exp
                @test_approx_eq y y0
                @test_approx_eq z z0
                # inc1ieq
                if ix == iy
                    x = copy(x0)
                    y = copy(y0)
                    func(x,ix,y,z,a,n)
                    @test_approx_eq x x_exp
                    @test_approx_eq y y0
                    @test_approx_eq z z0
                end
            end
        end
    end
end


text:
# Benchmark the calculation of a distance in a PDB file
# The distance is the closest distance between any atoms of residues 50 and 60
#   of chain A in 1AKE

using MIToS.PDB

pdb_filepath = "pdbs/1AKE.pdb"
struc = @residuesdict read(pdb_filepath, PDBFile) model "1" chain "A" group "ATOM" residue "*"

# Run to JIT compile
distance(struc["50"], struc["60"])

elapsed = @elapsed distance(struc["50"], struc["60"])

println(elapsed)

max_stars_repo_name: chemlove/pdb-benchmarks

max_issues_repo_name: chemlove/pdb-benchmarks

max_forks_repo_name: chemlove/pdb-benchmarks