text:
# SimLynx/src/statistics.jl
# Licensed under the MIT License. See LICENSE.md file in the project root for
# full license information.

abstract type Stats{T<:Real} end

# Accumulated Statistics

"""
    AccumulatedStats{T<:Real} <: Stats{T}

    AccumulatedStats()
    AccumulatedStats(value::T)

Accumulates statistics of values over time.

# Fields
- `min::T` - the minimum value observed (or zero(T))
- `max::T` - the maximum value observed (or zero(T))
- `n::Float64` - the weighted number of values (i.e., total time)
- `sum::Float64` - the weighted sum of the values
- `sum_squares::Float64` - the weighted sum of squares of the values
"""
mutable struct AccumulatedStats{T<:Real} <: Stats{T}
    min::T
    max::T
    n::Float64
    sum::Float64
    sum_squares::Float64
    AccumulatedStats{T}() where {T<:Real} =
        new(zero(T), zero(T), 0.0, 0.0, 0.0)
    AccumulatedStats{T}(value::T) where {T<:Real} =
        new(value, value, 0.0, 0.0, 0.0)
end

"Accumulate running statistics over a time."
function update!(stats::AccumulatedStats{T}, x::T, t::Float64) where {T<:Real}
    if t < 0.0
        throw(ArgumentError("weight t cannot be negative, got $t"))
    end
    if stats.n > 0.0
        stats.min = min(x, stats.min)
        stats.max = max(x, stats.max)
    else
        stats.min = x
        stats.max = x
    end
    stats.n += t
    stats.sum += x * t
    stats.sum_squares += x^2 * t
    return nothing
end

# Tallied Statistics

"""
    TalliedStats{T<:Real} <: Stats{T}

    TalliedStats()
    TalliedStats(value::T)

Tally statistics of values.

# Fields
- `min::T` - the minimum value observed (or zero(T))
- `max::T` - the maximum value observed (or zero(T))
- `n::Float64` - the number of values
- `sum::Float64` - the sum of the values
- `sum_squares::Float64` - the sum of squares of the values
"""
mutable struct TalliedStats{T<:Real} <: Stats{T}
    min::T
    max::T
    n::Int64
    sum::Float64
    sum_squares::Float64
    TalliedStats{T}() where {T<:Real} =
        new(zero(T), zero(T), 0.0, 0.0, 0.0)
    TalliedStats{T}(value::T) where {T<:Real} =
        new(value, value, 1, value, value^2)
end

"Tally running statistics."
function update!(stats::TalliedStats{T}, x::T) where {T<:Real}
    if stats.n > 0
        stats.min = min(x, stats.min)
        stats.max = max(x, stats.max)
    else
        stats.min = x
        stats.max = x
    end
    stats.n += 1
    stats.sum += x
    stats.sum_squares += x^2
    return nothing
end

# Stats Operations

mean(stats::Stats) = stats.sum / stats.n
mean_square(stats::Stats) = stats.sum_squares / stats.n
variance(stats::Stats) = mean_square(stats) - mean(stats)^2
stddev(stats::Stats) = sqrt(variance(stats))

"Overloaded statistics properties for mean, variance, and stddev."
function Base.getproperty(stats::Stats, name::Symbol)
    if name === :mean
        return mean(stats)
    elseif name === :variance
        return variance(stats)
    elseif name === :stddev
        return stddev(stats)
   else # fallback to getfield
       return getfield(stats, name)
   end
end

"""
Print the accumulated statistics for a variable with optional format and
title."""
function print_stats(stats::Stats;
                     format="%f",
                     title="Statistics")
    println(title)
    println("     min = $(stats.min)")
    println("     max = $(stats.max)")
    println("       n = $(stats.n)")
    println("    mean = $(stats.mean)")
    println("variance = $(stats.variance)")
    println("  stddev = $(stats.stddev)")
    flush(stdout)
end


text:
include("../src/JOH.jl")
include("MDMKP.jl")

tight_problems = MDMKP.load_folder()
loose_problems =  MDMKP.loosen.(tight_problems)
all_problems = vcat(tight_problems, loose_problems)

decision = MDMKP.decide.(tight_problems)
[(i, count(==(i), decision)) for i in unique(decision)]

decision = MDMKP.decide.(loose_problems)
[(i, count(==(i), decision)) for i in unique(decision)]

decision = MDMKP.decide.(all_problems)
[(i, count(==(i), decision)) for i in unique(decision)]

println(decision)

ext: jl

max_stars_repo_name: Dibillilia/JuliaOptHeuristics

max_issues_repo_name: Dibillilia/JuliaOptHeuristics

max_forks_repo_name: Dibillilia/JuliaOptHeuristics


text:
using Base.Iterators
using Zygote

"""
    xlogx(x::Number)

Compute `x * log(x)`, returning zero if `x` is zero. (Copied from StatsFuns.jl package.)
"""
function xlogx(x::Number)
    result = x * log(x)
    ifelse(iszero(x), zero(result), result)
end

"""
Compute the von Neumann entropy of a density matrix `ρ`.
"""
neumann_entropy(ρ::AbstractMatrix) = -sum(xlogx.(real(eigvals(Matrix(ρ)))))

"""
    binary_entropy(p::Number)

Compute the binary entropy of `p` (natural logarithm).
"""
binary_entropy(p::Number) = -(xlogx(p) + xlogx(1 - p))

"""
Compute `-log(Z)`, with `Z` the partition function.
"""
thermal_logZ(βH::AbstractMatrix) = -log(tr(exp(-Matrix(βH))))

"""
Trace distance to target density matrix `σ`.
"""
trdistance(θlist, σ, cgc) = 0.5 * opnorm(real.(matrix(apply(latent_density(θlist), cgc))) - σ, 1)

"""
Use representation of density matrix in terms of Pauli matrices; cos(θ) is Bloch vector coefficient of Pauli-Z for each qubit
"""
latent_density(θlist) = DensityMatrix(kron([[1, 0, 0, cos(θ)] for θ in θlist]...), length(θlist))
# """
# Use representation of density matrix in terms of Pauli matrices; cos(θ) is Bloch vector coefficient of Pauli-Z for each qubit
# """
# function latent_density2(θlist)
#     ρ = kron([[θ1, 0, 0, θ2] for (θ1,θ2) in partition(θlist,2)]...)
#     return DensityMatrix(ρ ./ ρ[1], length(θlist)÷2)
# end

# Zygote.@adjoint function partition(θlist, step::Int)
#     return partition(θlist, step), θ->(collect(Base.Iterators.flatten(θ)), 0)
# end

function ngradient(f, xs::AbstractArray...)
    grads = zero.(xs)
    for (x, Δ) in zip(xs, grads), i in 1:length(x)
        δ = sqrt(eps())
        tmp = x[i]
        x[i] = tmp - δ/2
        y1 = f(xs...)
        x[i] = tmp + δ/2
        y2 = f(xs...)
        x[i] = tmp
        Δ[i] = (y2-y1)/δ
        if eltype(x) <: Complex
            # derivative with respect to imaginary part
            x[i] = tmp - im*δ/2
            y1 = f(xs...)
            x[i] = tmp + im*δ/2
            y2 = f(xs...)
            x[i] = tmp
            Δ[i] += im*(y2-y1)/δ
        end
    end
    return grads
end

max_stars_repo_name: keefehuang/vqt

max_issues_repo_name: keefehuang/vqt

max_forks_repo_name: keefehuang/vqt


text:
include("Graph_class.jl")
include("Tools.jl")

using Random
using SharedArrays
using Distributed

#addprocs(4)

# """
#     activated_nodes_of_this_round = getRRS_one_round(g::Array{Array{Tuple{Int32, Float64}, 1}, 1}, A, activated_nodes_of_last_round, round)

#     A: Activated nodes of all rounds
# """
# function getRRS_one_round(g :: Array{Array{Tuple{Int32, Float64}, 1}, 1}, A, activated_nodes_of_last_round, delta, round_seed)
#     activated_nodes_of_this_round = Set{Int32}()
#     tmp = 0
#     if round_seed == 1
#         for s in activated_nodes_of_last_round
#             for nb in g[s]
#                 if nb[1] in A
#                     continue
#                 end
#                 tmp = rand(1)[1]
#                 if tmp < nb[2]
#                     push!(activated_nodes_of_this_round, nb[1])
#                 end
#             end
#         end
#     else
#         for s in activated_nodes_of_last_round
#             for nb in g[s]
#                 if nb[1] in A
#                     continue
#                 end
#                 tmp = rand(1)[1]
#                 if tmp < delta * nb[2]
#                     push!(activated_nodes_of_this_round, nb[1])
#                 end
#             end
#         end
#     end

#     #activated_nodes_of_this_round = toarray(activated_nodes_of_this_round) # convert to vector
#     #A = vcat(A, activated_nodes_of_this_round) # add to original vector ``
#     return activated_nodes_of_this_round
# end

"""
    list_of_nodes = getRRS(G :: Graph, n, delta)
    A: Activated nodes of all rounds

    2021-07-18, also return the nodes get picked

    2021-10-17, combine with getRRS_one_round

"""
function getRRS(g:: Array{Array{Tuple{Int32, Float64}, 1}, 1}, n::Number, delta::Number)
    rr_seeds = rand(1:n, 1)
    tmp = rand(sum([length(g[i]) for i in 1:n])) #m

    A = Set{Int32}(rr_seeds[1])
    round_seed = 1
    activated_nodes_of_last_round = A
    activated_nodes_of_this_round = Set{Int32}()
    count = 1 # random number

    while length(activated_nodes_of_last_round) != 0
        for s in activated_nodes_of_last_round
            for nb in g[s]
                if nb[1] in A
                    continue
                end

                if round_seed == 1
                    if tmp[count] < nb[2]
                        push!(activated_nodes_of_this_round, nb[1])
                    end
                else
                    if tmp[count] < delta * nb[2]
                        push!(activated_nodes_of_this_round, nb[1])
                    end
                end
                count += 1
            end
        end

        union!(A, activated_nodes_of_this_round)
        activated_nodes_of_last_round = activated_nodes_of_this_round
        activated_nodes_of_this_round = Set{Int32}()
        round_seed += 1
    end

    return toarray(A), rr_seeds[1]
end

"""
    node_selection_influence_max(R, k)

    Node selection procedure

    Return Fᵣ(SEED), and SEED
"""
function node_selection_influence_max(R, k)
    copy_k = k
    copy_R = deepcopy(R)

    SEED = zeros(Int32, copy_k)
    size = length(copy_R) # initial size
    tmp_size,  cm, fc, ris_spread = 0, 0, 0, 0
    # The process of node-selection
    for i in 1:copy_k
        if length(copy_R) == 0
            copy_k = i - 1
            break
        end
        flat_list = vcat(copy_R...) # Flat the list
        (cm, fc) = most_common(flat_list)
        tmp_size += fc
        SEED[i] = cm
        filter!(e -> !(cm in e), copy_R)
    end

    ris_spread = tmp_size / size

    return ris_spread, SEED[1:copy_k]
end

"""
    ris_influence_max(G::Graph, k, delta)

    delta: network parameter (ack, spread, reject)
    ℓ: accuracy parameter

    ris only for inluence maximization algorithm.
    Return: influence_spread, seed nodes, sample_size

        Tang, Youze, Yanchen Shi, and Xiaokui Xiao.
        "Influence maximization in near-linear time: A martingale approach."
        Proceedings of the 2015 ACM SIGMOD International Conference on Management of Data. 2015.
"""
function ris_influence_max(G::Graph, k, delta = 0.5, ϵ = 0.3, ℓ = 1)
    # Create the "reversed" Link Table
    g = Array{Array{Tuple{Int32, Float64}, 1}, 1}(undef, G.n)
    for i = 1 : G.n
        g[i] = []
    end

    for (_, u, v, w) in G.E
        push!(g[v], (u, w))
    end

    T = time()

    n = G.n
    ϵ′ = √2 * ϵ
    λ′ = (2 + 2/3*ϵ′)*(logchoose(n, k) + ℓ * log(n) + log(log2(n)))*n / (ϵ′)^2
    λ⁺ = 4*n*(ϵ/3 + 2)*(ℓ * log(n) + log(2) + logchoose(n, k)) / (ϵ)^2
    LB = 1
    R, SEED = [], []
    ris_spread, x, θᵢ = 0, 0, 0

    for i in 1 : round(log2(n) - 1)
        x = n / (2^i)
        θᵢ = λ′ / x
        if length(R) <= θᵢ
            append!(R, [getRRS(g, n, delta)[1] for _ in 1 : (θᵢ - length(R) + 1)])
        end

        ris_spread, _ = node_selection_influence_max(R, k)

        if n*ris_spread >= (1 + ϵ′)*x
            LB = n * ris_spread / (1 + ϵ′)
            break
        end

    end

    θ = λ⁺ / LB
    if length(R) <= θ
        append!(R, [getRRS(g, n, delta)[1] for _ in 1 : (θ - length(R) + 1)])
    end

    ris_spread, SEED = node_selection_influence_max(R, k)

    T = time() - T
    return ris_spread * n, SEED, length(R), T
end

function choose_high_degree(G::Graph, k)
    T = time()
    g = zeros(G.n)
    for (_, u, _, _) in G.E
        g[u] += 1
    end

    seednodes = topk(g, k)
    T = time() - T

    return T, seednodes
end

function random_spread(G::Graph, k, delta, samplesize = 10)
    n = G.n
    tmp = 0

    T = time()

    for i = 1 : ceil(Int, log(n) * samplesize)
        rand_k_nodes = rand(1:n, k)
        (_, p4) = cascadeIC(G, rand_k_nodes, round(Int, 2*(10*10)*log(2*n)), delta)
        tmp = tmp .+ p4
    end

    T = time() - T

    return T / ceil(Int, log(n) * samplesize), tmp / ceil(Int, log(n) * samplesize)
end


text:
# A cache to preallocate everything for the GDTW distance
struct GDTWWorkspace{T1, T2, T3}
    τ::T2
    l::T1
    u::T1
    l_prev::T1
    u_prev::T1
    l₀::T1
    u₀::T1
    min_costs::T2
    costs::T3
end

"""
    GDTWWorkspace{T}(M, N)

Creates a cache of numeric type `T` for use in [`gdtw`](@ref).
"""
function GDTWWorkspace(::Type{T}, M, N) where {T}
    GDTWWorkspace{Vector{T}, Matrix{T}, Array{T, 3}}(
        zeros(T, M, N), zeros(T, N), zeros(T, N),
        zeros(T, N), zeros(T, N), zeros(T, N),
        zeros(T, N), zeros(T, M, N), zeros(T, M, M, N)
    )
end

GDTWWorkspace(M, N) = GDTWWorkspace(Float64, M, N)

# refine the bounds, as described in Section 4.2 of DB19
function refine!(l_current, u_current, l_prev, u_prev, l₀, u₀, warp; η)
    @avx for i in eachindex(l_current, u_current, l_prev, u_prev, l₀, u₀, warp)
        δ = η * (u_prev[i] - l_prev[i]) / 2
        l_current[i] = max(warp[i] - δ, l₀[i])
        u_current[i] = min(warp[i] + δ, u₀[i])
    end
    return nothing
end

# inital choices of `l` and `u`, modified from Eq (7) of DB19
function inital_bounds!(l, u, t, smin, smax, symmetric)
    # We need to loosen the bounds to account for floating point error
    # Otherwise these bounds can be too tight and disallow valid moves
    # which can lead to very wrong results.
    smin = .99*smin
    smax = 1.01*smax

    @inbounds for i in eachindex(t, l, u)
        # You must be able to get to `warp[i]` in time `t[i]`, so
        #   `smax >= warp[i] / t[i] >= smin`
        # This gives a lower and upper bound on `warp[i]`, i.e., on `τ[:, i]`.
        # You must be able to get to `1` from `warp[i]` in time `1-t[i]`, so
        #   `smin <= (1-warp[i])/(1-t[i]) <= smax`
        # this gives another lower and upper bound.
        # The resulting bounds:
        lower = max(smin * t[i], 1 - smax * (1 - t[i]))
        upper = min(smax * t[i], 1 - smin * (1 - t[i]))

        if symmetric
            # We need to apply the above bounds to `ψ(s) = 2s - ϕ(s) = 2t[i] - warp[i]`
            # as well Which leads to the following.
            l[i] = max(lower, 2*t[i] - upper)
            u[i] = min(upper, 2*t[i] - lower)
        else
            l[i] = lower
            u[i] = upper
        end
    end
    return nothing
end

# an inplace update to `τ` with the bounds `l` and `u`.
# produces `τ` as described in Section 4.1 of DB19.
function update_τ!(τ, t, M, l, u)
    N = length(t)
    @assert size(τ) == (M, N)
    @inbounds for t = 1:N, j = 1:M
        τ[j, t] = l[t] + ((j - 1) / (M - 1)) * (u[t] - l[t])
    end
    return nothing
end

"""
    gdtw(
        x,
        y,
        ::Type{T}  = Float64;
        symmetric::Bool = true,
        M::Int     = 100,
        N::Int     = 100,
        t          = range(T(0), stop = T(1), length = N),
        cache::GDTWWorkspace = GDTWWorkspace(T, M, length(t)),
        λcum       = T(0.01),
        λinst      = T(0.01),
        η          = T(1 / 8),
        max_iters  = 3,
        metric     = (x, y) -> norm(x - y),
        Rcum       = abs2,
        smin::Real = T(0.001),
        smax::Real = T(5.0),
        Rinst      = symmetric  ?
                        ϕ′ -> ( (smin <= ϕ′ <= smax)
                            && (smin <= 2 - ϕ′ <= smax) ) ? (ϕ′-1)^2 : typemax(T)
                                :
                        ϕ′ -> (smin <= ϕ′ <= smax) ? (ϕ′-1)^2 : typemax(T),
        verbose    = false,
        warp       = zeros(T, length(t)),
        callback   = nothing,
    ) where T -> cost, ϕ, ψ

Computes a general DTW distance following [DB19](https://arxiv.org/abs/1905.12893).

Aims to find `ϕ(s)` to minimize

    ∫ metric(x(ϕ(s)), y(ψ(s))) + λinst*Rinst(ϕ'(s) - 1) + λcum*Rcum(ϕ(s) - s) ds

over the interval `s ∈ [0,1]`, where `ψ(s) = 2s - ϕ(s)` (if `symmetric=true`) or `ψ(s) = s`
(if `symmetric = false`). The integral is discretized in time into `N` points (or according
to the times `t`, if `t` is specified). Additionally, the possible values obtained by `ϕ`
(and hence `ψ`) at each possible time `s` are discretized into `M` points.

If `max_iters > 1`, then after solving the doubly-discretized problem to obtain the optimal `ϕ`,
the problem is solved again by choosing a new discretization of `M` possible values
of `ϕ(s)` in an interval (whose width is governed by the parameter `η`) around the
previous optimal value. This is repeated until the problem has been solved `max_iters`
times in total. Setting `verbose=true` prints the cost at each iteration; a "high enough"
value of `max_iters` can be chosen by inspecting when the cost stabilizes sufficiently.

The parameters are:

* `x`: the continuous time signal to warp (see [`LinearInterpolation`](@ref) for generating such a signal from discrete data)
* `y`: the continuous-time signal to warp to
* `T`: the numeric type to be used in the problem
* `symmetric`: if true, `ψ(s) = 2s - ϕ(s)`, otherwise `ψ(s) = s`.
* `t`: the discretization of time on `[0,1]`; either `t` or `N` should be specified
* `M`: the discretization of the values obtained by the warping path
* `metric`:  a function `metric(u,v) -> ℝ` to compute differences between the signals at a time point (such as a Distances.jl distance)
* `Rcum`: penalty function on the cumulative warp
* `Rinst`: penalty function on the instantaenous warping. Should be infinite outside of `[smin, smax]`.
* `smin`, `smax`: minimum and maximum allowed instantaenous warping. Should have `smin > 0` and `smin < smax`.
* `λcum`, `λinst`: the regularization constants for `Rcum` and `Rinst`, respectively

The following may be pre-allocated and reused between distance computations with the same `M` and `N` (or `length(t)`).

* `cache`: a cache of matrices and vectors, generated by `GDTW.GDTWWorkspace{T}(N,M)`

"""
function gdtw(args...; kwargs...)
    data = prepare_gdtw(args...; kwargs...)
    cost = iterative_gdtw!(data)
    return cost, gdtw_warpings(data)...
end

"""
    prepare_gdtw(x, y; kwargs...)

Creates a NamedTuple of parameters, using the same keyword argments as `dist`.
A preprocessing step before calling `iterative_gdtw!`.
"""
function prepare_gdtw(
    x,
    y,
    ::Type{T}  = Float64;
    symmetric::Bool = true,
    M::Int     = 100,
    N::Int     = 100,
    t::AbstractVector{T} = range(T(0), stop = T(1), length = N),
    cache::GDTWWorkspace = GDTWWorkspace(T, M, length(t)),
    λcum::T    = T(0.01),
    λinst::T   = T(0.01),
    η::T       = T(1 / 8),
    max_iters  = 3,
    metric     = (x, y) -> norm(x - y),
    Rcum       = abs2,
    smin::T    = T(0.001),
    smax::T    = T(5.0),
    Rinst      = symmetric  ?
                    ϕ′ -> ( (smin <= ϕ′ <= smax)
                        && (smin <= 2 - ϕ′ <= smax) ) ? (ϕ′-1)^2 : typemax(T)
                            :
                    ϕ′ -> (smin <= ϕ′ <= smax) ? (ϕ′-1)^2 : typemax(T),
    verbose::Bool = false,
    warp       = zeros(T, length(t)),
    callback   = nothing,
) where T
    N = length(t)

    (M > N / smax) || @warn "`M <= N / smax`; problem may be infeasible" M N smax

    @unpack l₀, l_prev, l,  u₀, u_prev, u, τ = cache
    inital_bounds!(l₀, u₀, t, smin, smax, symmetric)
    l_prev .= l₀
    u_prev .= u₀
    u .= u₀
    l .= l₀
    update_τ!(τ, t, M, l, u)

    function node_weight(j, s)
        s == length(t) && return zero(T)
        Rval = Rcum(τ[j, s] - t[s])
        yval = symmetric ? 2*t[s] - τ[j, s] : t[s]
        (t[s+1] - t[s])*(metric(x(τ[j, s]), y(yval)) + λcum * Rval)
    end

    @inline function edge_weight((j, s), (k, s2))
        s + 1 ≠ s2 && return typemax(T)
        ϕ′ = (τ[k, s+1] - τ[j, s]) / (t[s+1] - t[s])
        (t[s+1] - t[s]) * (λinst * Rinst(ϕ′))
    end

    return (
        iter        = Ref(1),
        N           = N,
        M           = M,
        τ           = τ,
        node_weight = node_weight,
        edge_weight = edge_weight,
        η           = η,
        max_iters   = max_iters,
        t           = t,
        smin        = smin,
        smax        = smax,
        callback    = callback,
        verbose     = verbose,
        metric      = metric,
        cache       = cache,
        warp        = warp,
        symmetric   = symmetric,
    )
end

"""
    iterative_gdtw!(data; max_iters = data.max_iters, verbose = data.verbose) -> cost

Runs the GDTW algorithm with iterative refinement for `max_iters` iterations,
returning the resulting cost. Uses the [`GDTWWorkspace`](@ref) in `data.cache` and
updates the `data.warp` vector. Here, `data` is usually obtained by [`prepare_gdtw`](@ref).

This can be called multiple times on the same `data` with a higher value of `max_iters`
to refine a calculation without starting over, which can be useful for checking convergence.

## Example

```julia
data = prepare_gdtw(x,y; max_iters = 3)
cost3 = iterative_gdtw!(data) # performs 3 iterations, returns the cost
cost10 = iterative_gdtw!(data; max_iters = 10) # performs 7 more iterations, returns the cost
ϕ, ψ = gdtw_warpings(data)

# Equivalently,
cost10, ϕ, ψ = gdtw(x,y; max_iters = 10)
```
"""
function iterative_gdtw!(data; max_iters = data.max_iters, verbose = data.verbose)
    @unpack N, M, τ, η, iter, t, callback = data
    @unpack cache, warp, symmetric = data
    if iter[] > max_iters
        @warn "`iter[] > max_iters`; no iterations performed." iter[] max_iters
        return zero(eltype(warp))
    end
    local cost

    # First iteration is special cased because
    # we can't send the cost and warp to the callback yet,
    # and we can quit early if we don't need to do refinement.
    if iter[] == 1
        if callback !== nothing
            callback((iter=1, t=t, τ=τ))
        end

        cost = single_gdtw!(data)
        verbose && @info "Iteration" iter[] cost
        iter[] += 1
        max_iters == 1 && return cost
    end

    @unpack l_prev, u_prev, l, u, l₀, u₀ = cache

    while iter[] <= max_iters
        l_prev .= l
        u_prev .= u
        refine!(l, u, l_prev, u_prev, l₀, u₀, warp; η=η)
        update_τ!(τ, t, M, l, u)
        cost = single_gdtw!(data)
        if callback !== nothing
            callback((iter = iter, t = t, τ = τ, warp = warp, cost = cost))
        end
        verbose && @info "Iteration" iter[] cost

        iter[] += 1
    end

    return cost
end

"""
    gdtw_warpings(data) -> ϕ, ψ

Computes the interpolations from a `data` `NamedTuple`
with entries for the time points `t`, warping points `warp`,
and a boolean `symmetric`.
"""
function gdtw_warpings(data)
    @unpack t, warp, symmetric = data
    if symmetric
        ψ = LinearInterpolation(2*t - warp, t)
    else
        ψ = LinearInterpolation(t, t)
    end
    ϕ = LinearInterpolation(warp, t)
    return ϕ, ψ
end

## Dynamic programming to compute the distance

function single_gdtw!(data::T) where {T}
    @unpack N, M, node_weight, cache, edge_weight, τ, warp = data
    @unpack min_costs, costs = cache
    calc_costs!(min_costs, costs, N, M, node_weight, edge_weight)
    cost = min_costs[end, end]
    trackback!(warp, costs, τ)
    return cost
end

function calc_costs!(min_costs, costs, N, M, node_weight::F1, edge_weight::F2) where {F1,F2}
    @boundscheck checkbounds(min_costs, 1:M, 1:N)
    @boundscheck checkbounds(costs, 1:M, 1:M, 1:N)

    @inbounds begin
        min_costs .= node_weight.(1:M, permutedims(1:N))
        # t = 2 case
        for j = 1:M
            costs[1, j, 2] = min_costs[1, 1] + edge_weight((1, 1), (j, 2))
            min_costs[j, 2] += costs[1, j, 2]
        end
        for t = 3:N
            for j = 1:M
                mi = typemax(eltype(costs))
                for k = 1:M
                    c = min_costs[k, t-1] + edge_weight((k, t - 1), (j, t))
                    costs[k, j, t] = c
                    mi = ifelse(c < mi, c, mi)
                end
                min_costs[j, t] += mi
            end
        end
    end
    return nothing
end

function trackback!(warp, costs, τ)
    (M, _, N) = size(costs)
    @boundscheck checkbounds(costs, 1:M, 1:M, 1:N)
    c = M
    @inbounds for t = N:-1:3
        warp[t] = τ[c, t]
        c = argmin(@views costs[:, c, t])
    end
    warp[2] = τ[c, 2]
    warp[1] = τ[1, 1]

    return nothing
end

## Interpolations

"""
LinearInterpolation(x::AbstractVector) -> Function

Provides a linear interpolation of `x` on the interval `[0,1]`.
"""
struct LinearInterpolation{Tx,Tt} <: Function
    x::Tx
    t::Tt
    function LinearInterpolation(x::Tx, ts::Ts) where {Tx,Ts}
        issorted(ts) || throw(ArgumentError("Time parameter `ts` must be sorted in increasing order."))
        T = eltype(Tx)
        t = (ts .- T(first(ts))) ./ T( last(ts) - first(ts))
        Tt = typeof(t)
        new{Tx,Tt}(x, t)
    end
end

LinearInterpolation(x) = LinearInterpolation(x, axes(x, ndims(x)))

function (xt::LinearInterpolation)(s)
    x = xt.x
    0 <= s <= 1 || return zero(x[!, 1])
    t = xt.t
    i = searchsortedlast(t, s)
    (i == 0) && return x[!, 1]
    (i == lastlength(x)) && return x[!, lastlength(x)]
    (s == t[i]) && return x[!, i]
    weight = (s - t[i]) / (t[i+1] - t[i])
    omw = 1 - weight
    x[!, i] .* omw .+ x[!, i+1] .* weight
end


text:
"""
PUSH_TO_LAYER

A setinel to aid the creation of layers.
If set to true, all the objects are pushed to the current layer.
Can be accessed using PUSH_TO_LAYER[1]
"""
const PUSH_TO_LAYER = Array{Any,1}()
# by default push all the objects to the video
push!(PUSH_TO_LAYER, false)

"""
    @JLayer(frames, width, height, position, body)
Calls the [`to_layer_m`](@ref) method to create a [`Layer`](@ref) out of the arguments given.

# Arguments
- `frames`:a `UnitRange` that defines for which frames the layer is active
- `width`: defines the width of the layer
- `height`: defines the height of the layer
- `position`: location of the center of the layer on the main canvas
- `transparent` : Whether the layer should have a transparent background(:transparent or :opaque)
- `body`
    - It contains all the objects(and thier respective actions) definitions for a layer
    - A layer can have it's own separate background
    - Anything defined within the `begin...end` block stays in the layer
    - A Layer has it's own coordinate reference sysstem, i.e. it has it's own origin
    So eg : `Point(100, 100)` is different when defined in a layer and doesn't represent
    the location 100, 100 on the main canvas

`width`, `height`, `position` and `transparent` are optional and default to the video's width, height, origin and :opaque respectively.
Layer declaration should take place before pushing objects to it if one is not using the macro

# Example
```julia
function ground(args...)
    background("white")
    sethue("black")
end

video = Video(500, 500)
Background(1:100, ground)
object((args...)->circle(O, 50, :fill))

l1 = @JLayer 10:70 100 100 Point(150, 150) begin
    red_ball = Object(20:60, (args...)->object(O, "red"), Point(50,0))
    act!(red_ball, Action(anim_rotate_around(2π, O)))
end
render(video; pathname="test.gif")
```
"""
macro JLayer(
    frames::Expr,
    width::Int,
    height::Int,
    position,
    transparent::QuoteNode,
    body::Expr,
)
    # frame range 1:50 is an expr
    # transparent -> Symbol is a QuoteNode
    esc(
        to_layer_m(
            frames,
            body,
            width = width,
            height = height,
            position = position,
            transparent = transparent,
        ),
    )
end

macro JLayer(frames::Expr, width::Int, height::Int, position, body::Expr)
    esc(to_layer_m(frames, body, width = width, height = height, position = position))
end

macro JLayer(frames::Expr, body::Expr)
    esc(to_layer_m(frames, body))
end

macro JLayer(frames::Expr, position::Expr, body::Expr)
    esc(to_layer_m(frames, position = position, body))
end

macro JLayer(frames::Expr, position::Expr, transparent::QuoteNode, body::Expr)
    esc(to_layer_m(frames, position = position, transparent = transparent, body))
end

macro JLayer(frames::Expr, transparent::QuoteNode, body::Expr)
    esc(to_layer_m(frames, transparent = transparent, body))
end

macro JLayer(frames::Expr, width::Int, height::Int, body)
    esc(to_layer_m(frames, body, width = width, height = height))
end

macro JLayer(frames::Expr, width::Int, height::Int, transparent::QuoteNode, body::Expr)
    esc(to_layer_m(frames, body, width = width, height = height, transparent = transparent))
end

"""
    to_layer_m( frames, body; width, height, position)
Helper method for the [`JLayer`](@ref) macro
Returns an expression that creates a layer and pushes the objects defined withing the body to the layer
:transparent is the default while the other :opaque copies the video's background
"""
function to_layer_m(
    frames,
    body;
    width = nothing, # will be set to CURRENT_VIDEO[1].width inside Javis.Layer
    height = nothing,
    position = Point(0, 0),
    transparent = QuoteNode(:transparent),
)
    quote
        layer = Javis.Layer($frames, $width, $height, $position)

        if $transparent == :transparent
            push!(layer.opts, :transparent => true)
        end

        # by default fetch the video's background as a layer's background
        # this is overriden by passing another ground to the layer explicity in the begin end block
        # if no background is needed :transparent flag should be passed
        video_backgrounds =
            filter(x -> get(x.opts, :in_global_layer, false), $CURRENT_VIDEO[1].objects)
        push!(layer.layer_objects, video_backgrounds...)

        if isempty(Javis.CURRENT_LAYER)
            push!(Javis.CURRENT_LAYER, layer)
        else
            Javis.CURRENT_LAYER[1] = layer
        end
        Javis.PUSH_TO_LAYER[1] = true
        $body
        Javis.PUSH_TO_LAYER[1] = false

        # not a problem now but a todo for later
        # remove duplicate backgrounds
        # if layer's background is defined, delete the first object
        # which is the video's background
        layer
    end
end

"""
    show_layer_frame(frames::UnitRange, layer_frame::Union{UnitRange,Int}, layer::Layer)
Repeat a layer's frame/frames for a given frame range.

# Arguments
- `frames::UnitRange`: The frame range for which the layer should be repeated
- `layer_frame::Union{UnitRange,Int}`: The layer frame range to repeat
- `layer::Layer`:the layer to be repeated
"""
function show_layer_frame(
    frames::UnitRange,
    layer_frame::Union{UnitRange,Int},
    layer::Layer,
)
    lc = layer.layer_cache
    lc.frames = frames
    if layer_frame isa Int
        lc.layer_frames = layer_frame + first(get_frames(layer)) - 1
    else
        lc.layer_frames =
            (first(layer_frame) + first(get_frames(layer)) - 1):(last(
                layer_frame,
            ) + first(get_frames(layer)) - 1)
    end
end

"""
    get_position(obj::Object)

Get access to the position of a layer.

# Returns
- `Point`: the point stored by the layer.

# Throws
- If the function of Layer didn't return a Point or Transformation
"""
function get_position(l::Layer)
    return get_position(l.position)
end

"""
    to_layer!(l::Layer, object::Object)

Pushes an object into the layer and out of the list of independent objects.
This method is helpful in case one doesn't want to include an object in the
`begin...end` block of [`@JLayer`](@ref).
"""
function to_layer!(l::Layer, object::Object)
    remove_from_video(object)
    push!(l.layer_objects, object)
end

"""
    to_layer!(l::Layer, objects::Vector{Object})

Pushes a list of objects into the layer,
This method is helpful in case one doesn't want to include an object in the
`begin...end` block of [`@JLayer`](@ref).
"""
function to_layer!(l::Layer, objects::Vector{Object})
    remove_from_video(objects)
    push!(l.layer_objects, objects...)
end

"""
    remove_from_video(object::Object)

Removes an object or a list of objects from the main video.
This is a helper method for the [`@to_layer!`](@ref) method and is supposed to be used internally
"""
function remove_from_video(object::Object)
    filter!(x -> x != object, CURRENT_VIDEO[1].objects)
end

function remove_from_video(objects::Vector{Object})
    filter!(x -> x ∉ objects, CURRENT_VIDEO[1].objects)
end


text:
# This file is part of TaylorSeries.jl, MIT licensed
#

using TaylorSeries
using Test

@testset "Testing an identity proved by Euler (8 variables)" begin
    make_variable(name, index::Int) = string(name, TaylorSeries.subscriptify(index))

    variable_names = String[make_variable("α", i) for i in 1:4]
    append!(variable_names, [make_variable("β", i) for i in 1:4])

    a1, a2, a3, a4, b1, b2, b3, b4 = set_variables(variable_names, order=4)

    lhs1 = a1^2 + a2^2 + a3^2 + a4^2
    lhs2 = b1^2 + b2^2 + b3^2 + b4^2
    lhs = lhs1 * lhs2

    rhs1 = (a1*b1 - a2*b2 - a3*b3 - a4*b4)^2
    rhs2 = (a1*b2 + a2*b1 + a3*b4 - a4*b3)^2
    rhs3 = (a1*b3 - a2*b4 + a3*b1 + a4*b2)^2
    rhs4 = (a1*b4 + a2*b3 - a3*b2 + a4*b1)^2

    rhs = rhs1 + rhs2 + rhs3 + rhs4

    @test lhs == rhs
    v = randn(8)
    @test evaluate( rhs, v) == evaluate( lhs, v)
end


text:
export MongoClient

type MongoClient
    host::String
    port::Integer
    _mongo::Ptr{Void}

    function MongoClient(host::String, port::Integer)
        p = ccall((:mongo_alloc, MONGO_LIB), Ptr{Void},  ())
        if p == C_NULL
            error("Could not create a mongo object")
        end

        client = new(host, port, p)
        finalizer(client, destroy)

        connect(client)

        return client
    end
end

MongoClient() = MongoClient("localhost", 27017)

### Private methods ###

function connect(client::MongoClient)
    errno = ccall((:mongo_client, MONGO_LIB), Int32, (Ptr{Void},Ptr{Uint8},Int32), client._mongo, bytestring(client.host), client.port)
    if errno == MONGO_ERROR
        error("mongo_client error: $(last_error(client))")
    end
end

function destroy(client::MongoClient)
    ccall((:mongo_destroy, MONGO_LIB), Void, (Ptr{Void},), client._mongo)
end

function last_error(client::MongoClient)
    ccall((:mongo_get_err, MONGO_LIB), Int32, (Ptr{Void},), client._mongo)
end


text:
module ComponentArrays

using ArrayInterface
using LinearAlgebra
using Requires

const FlatIdx = Union{Int, CartesianIndex, AbstractArray}
const FlatOrColonIdx = Union{FlatIdx, Colon}

include("utils.jl")
export fastindices # Deprecated

include("lazyarray.jl")

include("axis.jl")
export AbstractAxis, Axis, PartitionedAxis, ShapedAxis, ViewAxis, FlatAxis

include("componentindex.jl")
export KeepIndex

include("componentarray.jl")
export ComponentArray, ComponentVector, ComponentMatrix, getaxes, getdata, valkeys

include("array_interface.jl")
# Base methods: parent, size, elsize, axes, reinterpret, hcat, vcat, permutedims, IndexStyle, to_indices, to_index, getindex, setindex!, view, pointer, unsafe_convert, strides, stride
# ArrayInterface methods: strides, size, lu_instance, parent_type

include("namedtuple_interface.jl")
# Base methods: hash, ==, keys, haskey, propertynames, getproperty, setproperty!

include("similar_convert_copy.jl")
# Base methods: similar, zero, copy, copyto!, deepcopy, convert (to Array and NamedTuple), promote

include("broadcasting.jl")
# Base methods: BroadcastStyle, convert(to Broadcasted{Nothing}), similar, map, dataids
# Broadcast methods: BroadcastStyle, broadcasted, broadcast_unalias

include("show.jl")

include("plot_utils.jl")
export labels, label2index

required(filename) = include(joinpath("if_required", filename))

function __init__()
    @require ConstructionBase="187b0558-2788-49d3-abe0-74a17ed4e7c9" required("constructionbase.jl")
    @require DiffEqBase="2b5f629d-d688-5b77-993f-72d75c75574e" begin
        @require RecursiveFactorization="f2c3362d-daeb-58d1-803e-2bc74f2840b4" required("diffeqbase.jl")
    end
    @require SciMLBase="0bca4576-84f4-4d90-8ffe-ffa030f20462" required("scimlbase.jl")
    @require RecursiveArrayTools="731186ca-8d62-57ce-b412-fbd966d074cd" required("recursivearraytools.jl")
    @require StaticArrays="90137ffa-7385-5640-81b9-e52037218182" required("staticarrays.jl")
    @require ChainRulesCore="d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" required("chainrulescore.jl")
    @require ReverseDiff="37e2e3b7-166d-5795-8a7a-e32c996b4267" required("reversediff.jl")
end

end


text:
begin
    using Markdown
    using Plots
end

md"""
Function that optimizes sinusoid parameters to fit $\delta^{18}O$ data.

The second core function of the ShellChron growth model. Loops
through all data windows and uses the `growth_model` function
to create d18O series that match the input data. This step is
iterated and optimized (minimizing the Sum of Squared Residuals)
through the SCEUA algorithm (by Duan et al., 1992) which finds
the optimal input parameters to the growth rate and Sea Surface
Temperature (SST) sinusoids to simulate d18O data.

[dat] Matrix containing the input data

[dynwindow] Information on the position and length of modeling
windows

[transfer_function] Transfer function used to convert d18Oc to temperature
data.

[$\delta^{18}O_w$] Either a single value (constant $\delta^{18}O_w$) or a vector of length
equal to the period in SST data (365 days by default) containing information
about seasonality in d18Ow. Defaults to constant d18Ow of 0 permille VSMOW
(the modern mean ocean value)

[T_per] Period of SST sinusoid (in days; default = 365)

[G_per] Period of growth rate sinusoid (in days; default = 365)

[t_int] Time interval (in days; default = 1)

[t_maxtemp] Timing of the warmest day of the year (in julian day;
default = 182.5, or May 26th halfway through the year)

[SCEUApar] Parameters for SCEUA optimization (iniflg, ngs, maxn, kstop
pcento, peps). For details, refer to Duan et al. (1992) in references

[sinfit] Apply sinusoidal fitting to guess initial parameters for SCEUA
optimization? \code{true/FALSE}

[MC] Number of Monte Carlo simulations to apply for error propagation
Default = 1000

[plot] Should results of modeling be plotted? \code{true/FALSE}
@return A list containing the \code{resultarray} which contains the full
result of all simulations on each data window and the \code{parmat} listing
all optimized growth rate and SST parameters used to model d18O in each data
window

seealso Duan, Qingyun, Soroosh Sorooshian, and Vijai Gupta. "Effective and
efficient global optimization for conceptual rainfall runoff models." Water
resources research 28.4 (1992): 1015-1031. https://doi.org/10.1029/91WR02985
@references package dependencies: ggplot2 3.2.1; rtop 0.5.14
Function dependencies: sinreg, d18O_model, growth_model

\doi{10.1029/91WR02985}

"""
function run_model(  # Core function to run the entire model
            dat, # data to run the model on
            dynwindow; # The window vetor resulting from reading in the data
            transfer_function = "KimONeil97",
            d18Ow = 0,
            T_per = 365, # Temperature sinusoid parameters
            G_per = 365, # Growth sinusoid parameters
            t_int = 1, # Default time interval = 1 day
            t_maxtemp = 182.5, # Default time (day) at which maximum temperature is reached is 182.5 (exactly halfway through the year, or 1st of June)
            SCEUApar = (1, 25, 10000, 5, 0.01, 0.01), # Set parameters for SCEUA optimization (iniflg, ngs, maxn, kstop, pcento, peps)
            sinfit = true, # Apply sinusoidal fitting to guess initial parameters for SCEUA optimization? (true/FALSE)
            MC = 1000, # If errors = true, give the number of iterations for Monte Carlo simulation used in error propagation (default = 1000, if MC = 0 no eror propagation is done)
            plot = FALSE # Should the progress of model fitting be plotted?
    )

    d18Oc = nothing
    d18Oc_err = nothing
    Omod = nothing # Predefine variables to circumvent global variable binding error

    # Prepare data arrays for storage of modeling results
    resultarray = [ # Create array to contain all modeling results of overlapping windows
        zeros(dat), # Replicate matrix of dat + extra columns to contain all variables
        dim = [length(dat.D), length(dynwindow.x) + length(dat[:, 1]), 9] # Six variables, being: Modeled d18O, residuals, Day of the Year, Growth between datapoints, Instantaneous growth rate at datapoint and Modeled temperature
    ]

    parmat = DataFrame(zeros(7, length(dynwindow.x), names=dynwindow.x) # DataFrame in which to store the modeling parameters
    parmat.index = ["T_amp", "T_pha", "T_av", "G_amp", "G_pha", "G_av", "G_skw"]

    # Prepare plot to show model progress
    if plot == true
        p = plot(dat.D, dat.d18Oc, m=:c, yerr=dat.d18Oc_err, xerr=dat.D_err)
        fitplot = ggplot2::ggplot(dat, ggplot2::aes(D, d18Oc)) + # Create a plot showing the fit of each window on the original data, start with a plot of the original data
        ggplot2::geom_point() +
        ggplot2::geom_line() +
        ggplot2::geom_errorbar(ggplot2::aes(ymin = d18Oc - d18Oc_err, # Add error bars on model result (1 SD)
            ymax = d18Oc + d18Oc_err),
            width = dat$D_err) +
        ggplot2::geom_errorbarh(ggplot2::aes(xmin = D - D_err,
            xmax = D + D_err),
            height = 0.05) +
        ggplot2::ggtitle("Plot of d18Oc values vs. depth. Black = data, Red = model, Errorbars = 1 SD")
        plot(fitplot)
    end

    # Estimate growth rate variability and round up to nearest higher magnitude of 10 for conservative boundary
    GRavest = max(diff(dat[dat$YEARMARKER == 1,1])) / 365 # Estimate maximum growth rate from yearmarkers
    GRavmax = 10 ^ (ceiling(log(GRavest, 10))) # Round up to nearest higher magnitude of 10

    # Find tailored range of temperatures from data
    d18Oc_range = range(dat$d18Oc) # Find d18Oc range in data
    if(transfer_function == "KimONeil97"){ # Find temperature range (to be superseded with inverse d18O_model function in later updates)
        T_range = 18.03 * 1000 / (log((d18Oc_range - (0.97002 * rev(range(d18Ow)) - 29.98)) / 1000 + 1) * 1000 + 32.42) - 273.15 # Use Kim and O'Neil (1997) with conversion between VSMOW and VPDB by Brand et al. (2014)
    }else if(transfer_function == "GrossmanKu86"){
        T_range =  20.6 - 4.34 * (d18Oc_range - rev(range(d18Ow)) - 0.2) # Use Grossmann and Ku (1986) modified by Dettmann et al. (1999)
    }else{
        print("ERROR: Supplied transfer function is not recognized")
    }
    T_max = max(T_range)
    T_amp_max = 2 * abs(diff(T_range))

    # Collate lower boundaries of parameters
    parl = c(
        T_amp = 0, # Minimum T amplitude in degrees C
        T_pha = 0, # Minimum phase in days
        T_av = -4, # Minimum average T in degrees C
        G_amp = 0, # Minimum seasonal GR range in um/d
        G_pha = 0, # Minimum GR phase in days
        G_av = -1 * GRavmax, # Minimum average GR in um/d.
        G_skw = 0 # Minimum skew factor
    )

    # Collate upper boundaries of parameters
    paru = c(
        T_amp = round(T_amp_max + 0.5, 0), # Maximum T amplitude in degrees C
        T_pha = 365, # Maximum phase in days
        T_av = round(T_max + 0.5, 0), # Maximum average T in degrees C
        G_amp = 2 * GRavmax, # Maximum seasonal GR range in um/d
        G_pha = 365, # Maximum GR phase in days
        G_av = GRavmax, # Maximum average GR in um/d (based on conservative boundaries of YEARMARKER indicators)
        G_skw = 100 # Maximum skew factor
    )

    # Set parameters for SCEUA optimization
    iniflg = SCEUApar[1] # Flag for initial parameter array (default = 1; included)
    ngs = SCEUApar[2] # Number of complexes (sub-populations, default = 25)
    maxn = SCEUApar[3] # Maximum number of function evaluations allowed during optimization (default = 10000)
    kstop = SCEUApar[4] # Maximum number of evolution loops before convergency (default = 5)
    pcento = SCEUApar[5] # Percentage change allowed in function value criterion before stop (default = 0.01)
    peps = SCEUApar[6] # Convergence level for parameter set (difference between parameters required for stop; default = 0.01)

    # Run the model on all windows

    for(i in 1:length(dynwindow$x)){ # Loop over shell record
        print(paste("Processing Datawindow:", i, "of", length(dynwindow$x))) # Keep track of progress

        # Isolate year of d18O data based on window data
        Dsam = dat[dynwindow$x[i]:(dynwindow$x[i] + dynwindow$y[i] - 1), 1]
        Osam = dat[dynwindow$x[i]:(dynwindow$x[i] + dynwindow$y[i] - 1), 2]
        if(MC > 0){
            D_err = dat[dynwindow$x[i]:(dynwindow$x[i] + dynwindow$y[i] - 1), 4] # Optional: include error on D
            O_err = dat[dynwindow$x[i]:(dynwindow$x[i] + dynwindow$y[i] - 1), 5] # Optional: include error on d18Oc
        }else{
            D_err = rep(0, dynwindow$y)
            O_err = rep(0, dynwindow$y)
            MC = 0
        }

        if(sinfit){ # If sinusoidal fitting is enabled
            sinlist = sinreg(Dsam, Osam) # Run sinusoidal regression to find initial parameter values
            # Estimate starting parameters from regression results
            O_av_start = sinlist[[1]][1] # Export starting value for annual d18O average
            O_amp_start = sinlist[[1]][2] # Export starting value for d18O amplitude
            O_pha_start = sinlist[[1]][4] %% sinlist[[1]][3] # Estimate position (in depth of the first peak in d18O)
            O_per_start = sinlist[[1]][3] # Export starting value for period in distance domain
        }else{
            O_av_start = mean(Osam) # Estimate starting value for annual d18O average by mean of d18O in record
            O_amp_start = diff(range(Osam)) / 2 # Estimate starting value for d18O amplitude by half the difference between minimum and maximum d18Oc
            O_per_start = diff(range(Dsam)) # Estimate starting period as thickness of isolated year
            O_pha_start = 0.25 * O_per_start # Set starting phase to one quarter of a cycle if it cannot be estimated from sinusoidal regression
        }

        if(transfer_function == "KimONeil97"){
            T_av_start = 18.03 * 1000 / (1000 * log((O_av_start - (0.97002 * mean(d18Ow) - 29.98)) / 1000 + 1) + 32.42) - 273.15  # Estimate mean temperature. Use Kim and O'Neil (1997) with conversion between VSMOW and VPDB by Brand et al. (2014)
            T_amp_start = 18.03 * 1000 / (1000 * log((O_av_start - O_amp_start - (0.97002 * mean(d18Ow) - 29.98)) / 1000 + 1) + 32.42) - 273.15 - T_av_start # Estimate temperature amplitude. Use Kim and O'Neil (1997) with conversion between VSMOW and VPDB by Brand et al. (2014)
        }else if(transfer_function == "GrossmanKu86"){
            T_av_start = 20.6 - 4.34 * (O_av_start - mean(d18Ow) - 0.2) # Estimate mean temperature. Use Grossmann and Ku (1986) modified by Dettmann et al. (1999)
            T_amp_start = 20.6 - 4.34 * (O_av_start - O_amp_start - mean(d18Ow) - 0.2) - T_av_start # Estimate mean temperature. Use Grossmann and Ku (1986) modified by Dettmann et al. (1999)
        }else{
            print("ERROR: Supplied transfer function is not recognized")
        }

        O_peak = O_pha_start + Dsam[1] # Find position of d18O peak in distance domain
        T_pha_start = ((O_pha_start - 0.5 * O_per_start) %% O_per_start) / O_per_start * T_per # Estimate position of first peak in temperature (low in d18O) relative to annual cycle (days)
        G_av_start = O_per_start / G_per # Estimate average growth rate in distance/day

        years = 3 # Set default number of years to 3

        # Collate starting parameters
        par0 = c(
            T_amp = T_amp_start,
            T_pha = T_pha_start,
            T_av = T_av_start,

            G_amp = G_av_start / 2, # Start by estimating growth rate changes by half the average
            G_pha = T_pha_start, # Start by estimating that the peak in growth rate coincides with the peak in temperature
            G_av = G_av_start,
            G_skw = 50 # Start with a no skew
        )

        invisible(capture.output( # Suppress the details on converging SCEUA
            sceua_list = rtop::sceua(growth_model,
                par0,
                T_per = T_per,
                G_per = G_per,
                years = years,
                t_int = t_int,
                transfer_function = transfer_function,
                d18Ow = d18Ow,
                Dsam = Dsam,
                Osam = Osam,
                t_maxtemp = t_maxtemp,
                parl,
                paru,
                maxn,
                kstop,
                pcento,
                ngs,
                iniflg = iniflg,
                peps = peps,
                implicit = function(pars){sum(pars[4]/2, pars[6]) < 1} # Make sure that the cumulative GR curve is not located below 0 (if G_av < - G_amp / 2)
                )
            ))

        par1 = sceua_list[[1]] # Eport parameters of final model
        names(par1) = names(par0)

        result = growth_model(par1, T_per, G_per, years, t_int, transfer_function, d18Ow, Dsam, Osam, t_maxtemp, plot = FALSE, MC, D_err, O_err, return = "result") # Calculate the end result of the best fit

        if(plot == true){
            fitplot = fitplot + # Add the new model fit to the plot to track progress of the model
                ggplot2::geom_point(data = as.data.frame(result), ggplot2::aes(Dsam, Omod), colour = "red") +
                ggplot2::geom_line(data = as.data.frame(result), ggplot2::aes(Dsam, Omod), colour = "red")
            print(fitplot)
        }

        resultarray[, i + length(dat[1, ]), ] = rbind(matrix(NA, nrow = i - 1, ncol = 9), result[, 3:11], matrix(NA, nrow = length(dat$D) - i - length(result[,1]) + 1, ncol = 9)) # Add results to result array
        parmat[, i] = par1 # Add parameters to parameter array
    }

    # Provide names for all dimensions in the result array
    dimnames(resultarray) = list(
        paste("sample", 1:length(resultarray[, 1, 3])),
        c(colnames(dat), paste("window", 1:length(dynwindow$x))),
        c("Modeled_d18O", "d18O_residuals", "Time_of_year", "Instantaneous_growth_rate", "Modeled temperature", "Modeled_d18O_SD", "Time_of_Year_SD", "Instantaneous_growth_rate_SD", "Modeled_temperature_SD")
    )

    colnames(parmat) = paste("window", 1:length(parmat[1,]))
    return(list(resultarray, parmat))
end


text:
module FlowNodeCouplingFlowsRadialFlowTest

using Test
using ReactiveMP
using ReactiveMP: getdim, getz0, getα, getβ, getall, setz0!, setα!, setβ!
using ReactiveMP: forward, forward!, jacobian, jacobian!, inv_jacobian, inv_jacobian!
using ReactiveMP: det_jacobian, absdet_jacobian, logdet_jacobian, logdet_jacobian, logabsdet_jacobian

@testset "Radial Flow" begin

    @testset "Constructor" begin

        # check for unspecified input
        f = RadialFlow()
        @test typeof(f) == ReactiveMP.RadialFlowPlaceholder

        # check when creating placeholder
        f = ReactiveMP.RadialFlowPlaceholder()
        @test typeof(f) == ReactiveMP.RadialFlowPlaceholder

        # check when creating empty planar flow
        f = ReactiveMP.RadialFlowEmpty(3)
        @test typeof(f) == ReactiveMP.RadialFlowEmpty{3}
        @test getdim(f) == 3

        # check for specified dimensionality
        f = RadialFlow(3)
        @test typeof(f.z0) == Array{Float64,1}
        @test typeof(f.α)  == Float64
        @test typeof(f.β)  == Float64
        @test length(f.z0) == 3

        f = RadialFlow(5)
        @test typeof(f.z0) == Array{Float64,1}
        @test typeof(f.α)  == Float64
        @test typeof(f.β)  == Float64
        @test length(f.z0) == 5

        # check for prespecified parameters (univariate)
        f = RadialFlow(1.0, 2.0, 3.0)
        @test typeof(f.z0) == Float64
        @test typeof(f.α)  == Float64
        @test typeof(f.β)  == Float64
        @test f.z0 == 1.0
        @test f.α  == 2.0
        @test f.β  == 3.0

        # check for prespecified parameters (multivariates)
        f = RadialFlow([2.0, 3.0], 5.0, 3.0)
        @test typeof(f.z0) == Array{Float64,1}
        @test typeof(f.α)  == Float64
        @test typeof(f.β)  == Float64
        @test f.z0 == [2.0, 3.0]
        @test f.α  == 5.0
        @test f.β  == 3.0

        # wrong constructors
        @test_throws MethodError RadialFlow(1.0, 2.0, [2.0, 3.0])
        @test_throws MethodError RadialFlow(0.5, 1.0, [5.0, 6.0])
        @test_throws MethodError RadialFlow(2.0, [2.0, 3.0], 1.0)
        @test_throws MethodError RadialFlow(1.0, [5.0, 6.0], 0.5)
        @test_throws AssertionError RadialFlow(1.0, -2.0, 1.0)
        @test_throws AssertionError RadialFlow(2.0, -1.0, -2.0)
        @test_throws AssertionError RadialFlow([2.0, 3.0], -1.0, 1.0)
        @test_throws AssertionError RadialFlow([5.0, 9.0], -5.0, -2.0)

    end

    @testset "Prepare-Compile" begin

        f = ReactiveMP.RadialFlowPlaceholder()
        @test ReactiveMP.prepare(2, f) == ReactiveMP.RadialFlowEmpty(2)
        @test ReactiveMP.prepare(5, f) == ReactiveMP.RadialFlowEmpty(5)
        @test ReactiveMP.prepare(9, f) == ReactiveMP.RadialFlowEmpty(9)

        f = ReactiveMP.RadialFlowEmpty(1)
        params = [1.0, 2.0, 3.0]
        fc = ReactiveMP.compile(f)
        fcp = ReactiveMP.compile(f, params)
        @test typeof(fc)     <: RadialFlow
        @test typeof(fc.z0)  == Float64
        @test typeof(fc.α)   == Float64
        @test typeof(fc.β)   == Float64
        @test typeof(fcp)    <: RadialFlow
        @test typeof(fcp.z0) == Float64
        @test typeof(fcp.α)  == Float64
        @test typeof(fcp.β)  == Float64
        @test fcp.z0         == 1.0
        @test fcp.α          == 2.0
        @test fcp.β          == 3.0

        f = ReactiveMP.RadialFlowEmpty(2)
        params = [1.0, 2.0, 3.0, 4.0]
        fc = ReactiveMP.compile(f)
        fcp = ReactiveMP.compile(f, params)
        @test typeof(fc)     <: RadialFlow
        @test typeof(fc.z0)  == Array{Float64,1}
        @test typeof(fc.α)   == Float64
        @test typeof(fc.β)   == Float64
        @test typeof(fcp)    <: RadialFlow
        @test typeof(fcp.z0) == Array{Float64,1}
        @test typeof(fcp.α)  == Float64
        @test typeof(fcp.β)  == Float64
        @test fcp.z0         == [1.0, 2.0]
        @test fcp.α          == 3.0
        @test fcp.β          == 4.0

    end

    @testset "nr_params" begin

        for k = 1:10
            f = ReactiveMP.RadialFlowEmpty(k)
            fc = compile(f)
            @test nr_params(f)  == k + 2
            @test nr_params(fc) == k + 2
        end

    end

    @testset "Get-Set" begin

        # check get functions for univariate RadialFlow
        f = RadialFlow(1.0, 2.0, 3.0)
        @test getz0(f)  == f.z0
        @test getα(f)   == f.α
        @test getβ(f)   == f.β
        @test getall(f) == (f.z0, f.α, f.β)
        @test getz0(f)  == 1.0
        @test getα(f)   == 2.0
        @test getβ(f)   == 3.0
        @test getall(f) == (1.0, 2.0, 3.0)

        # check get functions for multivariate RadialFlow
        f = RadialFlow([1.0, 2.0], 3.0, 4.0)
        @test getz0(f)  == f.z0
        @test getα(f)   == f.α
        @test getβ(f)   == f.β
        @test getall(f) == (f.z0, f.α, f.β)
        @test getz0(f)  == [1.0, 2.0]
        @test getα(f)   == 3.0
        @test getβ(f)   == 4.0
        @test getall(f) == ([1.0, 2.0], 3.0, 4.0)

        # check setz0! function (univariate)
        f = RadialFlow(1.0, 2.0, 3.0)
        setz0!(f, 4.0)
        @test getz0(f) == 4.0
        @test getα(f)  == 2.0
        @test getβ(f)  == 3.0

        # check setz0! function (multivariate)
        f = RadialFlow([1.0, 2.0], 3.0, 4.0)
        setz0!(f, [10.0, 11.0])
        @test getz0(f) == [10.0, 11.0]
        @test getα(f)  == 3.0
        @test getβ(f)  == 4.0

        # check setα! function (univariate)
        f = RadialFlow(1.0, 2.0, 3.0)
        setα!(f, 4.0)
        @test getz0(f) == 1.0
        @test getα(f)  == 4.0
        @test getβ(f)  == 3.0

        # check setα! function (multivariate)
        f = RadialFlow([1.0, 2.0], 3.0, 5.0)
        setα!(f, 4.0)
        @test getz0(f) == [1.0, 2.0]
        @test getα(f)  == 4.0
        @test getβ(f)  == 5.0

        # check setβ! function (univariate)
        f = RadialFlow(1.0, 2.0, 3.0)
        setβ!(f, 4.0)
        @test getz0(f) == 1.0
        @test getα(f)  == 2.0
        @test getβ(f)  == 4.0

        # check setβ! function (multivariate)
        f = RadialFlow([1.0, 2.0], 3.0, 5.0)
        setβ!(f, 6.0)
        @test getz0(f) == [1.0, 2.0]
        @test getα(f)  == 3.0
        @test getβ(f)  == 6.0

        # check getdim
        f = ReactiveMP.RadialFlowEmpty(2)
        @test ReactiveMP.getdim(f) == 2

        # check errors (univariate)
        f = RadialFlow(1.0, 2.0, 3.0)
        @test_throws MethodError setz0!(f, [1.0, 2.0])
        @test_throws MethodError setα!(f, [1.0, 2.0])
        @test_throws MethodError setβ!(f, [1.0, 2.0])

        # check errors (multivariate)
        f = RadialFlow([1.0, 2.0], 3.0, 5.0)
        @test_throws MethodError setz0!(f, 3.0)
        @test_throws MethodError setα!(f, [4.0, 3.0])
        @test_throws MethodError setβ!(f, [1.0, 2.0])
        @test_throws AssertionError setz0!(f, [1.0, 3.0, 6.0])

    end

    @testset "Base" begin

        # check base functions (univariate)
        f = RadialFlow(1.0, 2.0, 3.0)
        @test eltype(f) == Float64
        @test eltype(RadialFlow{Float64,Float64}) == Float64
        @test size(f) == 1
        @test length(f) == 1

        # check base functions (multivariate)
        f = RadialFlow([1.0, 2.0], 3.0, 5.0)
        @test eltype(f) == Float64
        @test eltype(RadialFlow{Array{Float64,1},Float64}) == Float64
        @test size(f) == 2
        @test length(f) == 2

        # check base functions empty object
        f = ReactiveMP.RadialFlowEmpty(2)
        @test size(f)   == 2
        @test length(f) == 2

    end

    @testset "Forward-Backward" begin

        # check forward function (univariate)
        f = RadialFlow(1.0, 2.0, -3.0)
        @test forward(f, 1.5) == 0.9
        @test forward(f, 2.5) == 1.2142857142857144
        @test forward.(f, [1.5, 2.5]) == [0.9, 1.2142857142857144]

        # check forward function (univariate mapping, multivariate input)
        f = RadialFlow(1.0, 2.0, -3.0)
        @test forward(f, [1.5]) == 0.9
        @test forward(f, [2.5]) == 1.2142857142857144
        @test forward.(f, [[1.5], [2.5]]) == [0.9, 1.2142857142857144]

        # check forward function (multivariate)
        f = RadialFlow([1.0, 2.0], 3.0, 1.0)
        @test forward(f, [-4.0, 3.0]) == [-4.617358680468466, 3.1234717360936934]
        @test forward(f, [-2.0, 1.5]) == [-2.4965751817893183, 1.4172374697017802]
        @test forward.(f, [[-4.0, 3.0], [-2.0, 1.5]]) == [[-4.617358680468466, 3.1234717360936934], [-2.4965751817893183, 1.4172374697017802]]

        # check forward! function (multivariate)
        f = RadialFlow([1.0, 2.0], 3.0, 1.0)
        output = zeros(2)
        forward!(output, f, [-4.0, 3.0])
        @test output == [-4.617358680468466, 3.1234717360936934]
        forward!(output, f, [-2.0, 1.5])
        @test output == [-2.4965751817893183, 1.4172374697017802]

    end

    @testset "Jacobian" begin

        # check jacobian function (univariate)
        f = RadialFlow(1.0, 2.0, -3.0)
        @test jacobian(f, 1.5) == 0.03999999999999987
        @test jacobian(f, 2.5) == 0.5102040816326531
        @test jacobian.(f, [1.5, 2.5]) == [0.03999999999999987, 0.5102040816326531]

        # check jacobian function (univariate mapping, multivariate input)
        f = RadialFlow(1.0, 2.0, -3.0)
        @test jacobian(f, [1.5]) == 0.03999999999999987
        @test jacobian(f, [2.5]) == 0.5102040816326531
        @test jacobian.(f, [[1.5], [2.5]]) == [0.03999999999999987, 0.5102040816326531]

        # check jacobian function (multivariate)
        f = RadialFlow([1.0, 2.0], 3.0, 1.0)
        @test jacobian(f, [-4.0, 3.0]) == [1.0487256521978074 0.014949216779177184; 0.014949216779177184 1.120481892737858]
        @test jacobian(f, [-2.0, 1.5]) == [1.084447783638529 -0.013512879492985073; -0.013512879492985073 1.1632729140142752]
        @test jacobian.(f, [[-4.0, 3.0], [-2.0, 1.5]]) == [[1.0487256521978074 0.014949216779177184; 0.014949216779177184 1.120481892737858], [1.084447783638529 -0.013512879492985073; -0.013512879492985073 1.1632729140142752]]

        # check jacobian! function (multivariate)
        f = RadialFlow([1.0, 2.0], 3.0, 1.0)
        output = zeros(2,2)
        jacobian!(output, f, [-4.0, 3.0])
        @test output == [1.0487256521978074 0.014949216779177184; 0.014949216779177184 1.120481892737858]
        jacobian!(output, f, [-2.0, 1.5])
        @test output == [1.084447783638529 -0.013512879492985073; -0.013512879492985073 1.1632729140142752]

    end

    @testset "Utility Jacobian" begin

        # check utility functions jacobian (univariate)
        f = RadialFlow(1.0, 2.0, -3.0)
        @test det_jacobian(f, 1.5) == 0.03999999999999987
        @test inv_jacobian(f, 1.5) ≈ 25
        @test absdet_jacobian(f, 1.5) == 0.03999999999999987
        @test logabsdet_jacobian(f, 1.5) == -3.218875824868204

        # check utility functions jacobian (multivariate)
        f = RadialFlow([1.0, 2.0], 3.0, 1.0)
        @test det_jacobian(f, [-4.0, 3.0]) == 1.1748546246550329
        @test inv_jacobian(f, [-4.0, 3.0]) ≈ [0.9537196085574075 -0.012724311983337218; -0.012724311983337218 0.8926429110373889]
        @test absdet_jacobian(f, [-4.0, 3.0]) == 1.1748546246550329
        @test logabsdet_jacobian(f, [-4.0, 3.0]) == (0.16114441624386985, 1.0)

    end

end

end


text:
include("preliminaries.jl")

@testset "smote" begin
    include("smote.jl")
end


text:
export iauGst00b
"""
Greenwich apparent sidereal time (consistent with IAU 2000
resolutions but using the truncated nutation model IAU 2000B).

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   uta,utb    double    UT1 as a 2-part Julian Date (Notes 1,2)

Returned (function value):
               double    Greenwich apparent sidereal time (radians)

Notes:

   1. The UT1 date uta+utb is a Julian Date, apportioned in any
   convenient way between the argument pair.  For example,
   JD=2450123.7 could be expressed in any of these ways, among
   others:

            uta            utb

      2450123.7           0.0       (JD method)
      2451545.0       -1421.3       (J2000 method)
      2400000.5       50123.2       (MJD method)
      2450123.5           0.2       (date & time method)

   The JD method is the most natural and convenient to use in cases
   where the loss of several decimal digits of resolution is
   acceptable.  The J2000 and MJD methods are good compromises
   between resolution and convenience.  For UT, the date & time
   method is best matched to the algorithm that is used by the Earth
   Rotation Angle function, called internally:  maximum precision is
   delivered when the uta argument is for 0hrs UT1 on the day in
   question and the utb argument lies in the range 0 to 1, or vice
   versa.

   2. The result is compatible with the IAU 2000 resolutions, except
   that accuracy has been compromised for the sake of speed and
   convenience in two respects:

   . UT is used instead of TDB (or TT) to compute the precession
      component of GMST and the equation of the equinoxes.  This
      results in errors of order 0.1 mas at present.

   . The IAU 2000B abridged nutation model (McCarthy & Luzum, 2001)
      is used, introducing errors of up to 1 mas.

   3. This GAST is compatible with the IAU 2000 resolutions and must be
   used only in conjunction with other IAU 2000 compatible
   components such as precession-nutation.

   4. The result is returned in the range 0 to 2pi.

   5. The algorithm is from Capitaine et al. (2003) and IERS
   Conventions 2003.

Called:
   iauGmst00    Greenwich mean sidereal time, IAU 2000
   iauEe00b     equation of the equinoxes, IAU 2000B
   iauAnp       normalize angle into range 0 to 2pi

References:

   Capitaine, N., Wallace, P.T. and McCarthy, D.D., "Expressions to
   implement the IAU 2000 definition of UT1", Astronomy &
   Astrophysics, 406, 1135-1149 (2003)

   McCarthy, D.D. & Luzum, B.J., "An abridged model of the
   precession-nutation of the celestial pole", Celestial Mechanics &
   Dynamical Astronomy, 85, 37-49 (2003)

   McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
   IERS Technical Note No. 32, BKG (2004)

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# double iauGst00b(double uta, double utb)

function iauGst00b(uta::Real, utb::Real)
   return ccall((:iauGst00b, libsofa_c), Cdouble,
               (Cdouble, Cdouble),
               convert(Float64, uta), convert(Float64, utb))
end


text:
using AntsArticle
using Documenter

makedocs(;
    modules=[AntsArticle],
    authors="yakir12 <12.yakir@gmail.com> and contributors",
    repo="https://github.com/yakir12/AntsArticle.jl/blob/{commit}{path}#L{line}",
    sitename="AntsArticle.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://yakir12.github.io/AntsArticle.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/yakir12/AntsArticle.jl",
)


text:
using LinearAlgebra

function matcond(n::Int, c::Float64)
# Function generates a random square matrix A of size n with
# a given condition number c.
# Inputs:
#	n: size of matrix A, n>1
#	c: condition of matrix A, c>= 1.0
#
# Usage: matcond(10, 100.0)
#
# Pawel Zielinski
        if n < 2
         error("size n should be > 1")
        end
        if c< 1.0
         error("condition number  c of a matrix  should be >= 1.0")
        end
        (U,S,V)=svd(rand(n,n))
        return U*diagm(0 =>[LinRange(1.0,c,n);])*V'
end

max_stars_repo_name: jerry-sky/academic-notebook

max_issues_repo_name: jerry-sky/academic-notebook

max_forks_repo_name: jerry-sky/academic-notebook


text:
using IndexedTables: collect_columns_flattened

@testset "collectnamedtuples" begin
    v = [@NT(a = 1, b = 2), @NT(a = 1, b = 3)]
    @test collect_columns(v) == Columns(@NT(a = Int[1, 1], b = Int[2, 3]))

    # test inferrability with constant eltype
    itr = [@NT(a = 1, b = 2), @NT(a = 1, b = 2), @NT(a = 1, b = 12)]
    st = start(itr)
    el, st = next(itr, st)
    dest = similar(IndexedTables.arrayof(typeof(el)), 3)
    dest[1] = el
    @inferred IndexedTables.collect_to_columns!(dest, itr, 2, st)

    v = [@NT(a = 1, b = 2), @NT(a = 1.2, b = 3)]
    @test collect_columns(v) == Columns(@NT(a = [1, 1.2], b = Int[2, 3]))
    @test typeof(collect_columns(v)) == typeof(Columns(@NT(a = [1, 1.2], b = Int[2, 3])))

    v = [@NT(a = 1, b = 2), @NT(a = 1.2, b = "3")]
    @test collect_columns(v) == Columns(@NT(a = [1, 1.2], b = Any[2, "3"]))
    @test typeof(collect_columns(v)) == typeof(Columns(@NT(a = [1, 1.2], b = Any[2, "3"])))

    v = [@NT(a = 1, b = 2), @NT(a = 1.2, b = 2), @NT(a = 1, b = "3")]
    @test collect_columns(v) == Columns(@NT(a = [1, 1.2, 1], b = Any[2, 2, "3"]))
    @test typeof(collect_columns(v)) == typeof(Columns(@NT(a = [1, 1.2, 1], b = Any[2, 2, "3"])))

    # length unknown
    itr = Iterators.filter(isodd, 1:8)
    tuple_itr = (@NT(a = i+1, b = i-1) for i in itr)
    @test collect_columns(tuple_itr) == Columns(@NT(a = [2, 4, 6, 8], b = [0, 2, 4, 6]))
    tuple_itr_real = (i == 1 ? @NT(a = 1.2, b =i-1) : @NT(a = i+1, b = i-1) for i in itr)
    @test collect_columns(tuple_itr_real) == Columns(@NT(a = Real[1.2, 4, 6, 8], b = [0, 2, 4, 6]))

    # empty
    itr = Iterators.filter(t -> t > 10, 1:8)
    tuple_itr = (@NT(a = i+1, b = i-1) for i in itr)
    @test collect_columns(tuple_itr) == Columns(@NT(a = Int[], b = Int[]))

    itr = (i for i in 0:-1)
    tuple_itr = (@NT(a = i+1, b = i-1) for i in itr)
    @test collect_columns(tuple_itr) == Columns(@NT(a = Int[], b = Int[]))
end

@testset "collecttuples" begin
    v = [(1, 2), (1, 3)]
    @test collect_columns(v) == Columns((Int[1, 1], Int[2, 3]))
    @inferred collect_columns(v)

    v = [(1, 2), (1.2, 3)]
    @test collect_columns(v) == Columns(([1, 1.2], Int[2, 3]))

    v = [(1, 2), (1.2, "3")]
    @test collect_columns(v) == Columns(([1, 1.2], Any[2, "3"]))
    @test typeof(collect_columns(v)) == typeof(Columns(([1, 1.2], Any[2, "3"])))

    v = [(1, 2), (1.2, 2), (1, "3")]
    @test collect_columns(v) == Columns(([1, 1.2, 1], Any[2, 2, "3"]))
    # length unknown
    itr = Iterators.filter(isodd, 1:8)
    tuple_itr = ((i+1, i-1) for i in itr)
    @test collect_columns(tuple_itr) == Columns(([2, 4, 6, 8], [0, 2, 4, 6]))
    tuple_itr_real = (i == 1 ? (1.2, i-1) : (i+1, i-1) for i in itr)
    @test collect_columns(tuple_itr_real) == Columns(([1.2, 4, 6, 8], [0, 2, 4, 6]))
    @test typeof(collect_columns(tuple_itr_real)) == typeof(Columns(([1.2, 4, 6, 8], [0, 2, 4, 6])))

    # empty
    itr = Iterators.filter(t -> t > 10, 1:8)
    tuple_itr = ((i+1, i-1) for i in itr)
    @test collect_columns(tuple_itr) == Columns(Int[], Int[])

    itr = (i for i in 0:-1)
    tuple_itr = ((i+1, i-1) for i in itr)
    @test collect_columns(tuple_itr) == Columns(Int[], Int[])
end

@testset "collectscalars" begin
    v = (i for i in 1:3)
    @test collect_columns(v) == [1,2,3]
    @inferred collect_columns(v)

    v = (i == 1 ? 1.2 : i for i in 1:3)
    @test collect_columns(v) == collect(v)

    itr = Iterators.filter(isodd, 1:100)
    @test collect_columns(itr) == collect(itr)
    real_itr = (i == 1 ? 1.5 : i for i in itr)
    @test collect_columns(real_itr) == collect(real_itr)
    @test eltype(collect_columns(real_itr)) == Float64

    #empty
    itr = Iterators.filter(t -> t > 10, 1:8)
    tuple_itr = (exp(i) for i in itr)
    @test collect_columns(tuple_itr) == Float64[]

    itr = (i for i in 0:-1)
    tuple_itr = (exp(i) for i in itr)
    @test collect_columns(tuple_itr) == Float64[]

    t = collect_columns(@NT(a = i) for i in (1, DataValue{Int}(), 3))
    @test columns(t, 1) isa DataValueArray
    @test isequal(columns(t, 1), DataValueArray([1, DataValue{Int}(), 3]))
end

@testset "collectpairs" begin
    v = (i=>i+1 for i in 1:3)
    @test collect_columns(v) == Columns([1,2,3]=>[2,3,4])
    @test eltype(collect_columns(v)) == Pair{Int, Int}

    v = (i == 1 ? (1.2 => i+1) : (i => i+1) for i in 1:3)
    @test collect_columns(v) == Columns([1.2,2,3]=>[2,3,4])
    @test eltype(collect_columns(v)) == Pair{Float64, Int}

    v = (@NT(a=i) => @NT(b="a$i") for i in 1:3)
    @test collect_columns(v) == Columns(Columns(@NT(a = [1,2,3]))=>Columns(@NT(b = ["a1","a2","a3"])))
    @test eltype(collect_columns(v)) == Pair{NamedTuples._NT_a{Int64}, NamedTuples._NT_b{String}}

    v = (i == 1 ? @NT(a="1") => @NT(b="a$i") : @NT(a=i) => @NT(b="a$i") for i in 1:3)
    @test collect_columns(v) == Columns(Columns(@NT(a = ["1",2,3]))=>Columns(@NT(b = ["a1","a2","a3"])))
    @test eltype(collect_columns(v)) == Pair{NamedTuples._NT_a{Any}, NamedTuples._NT_b{String}}

    # empty
    v = (@NT(a=i) => @NT(b="a$i") for i in 0:-1)
    @test collect_columns(v) == Columns(Columns(@NT(a = Int[]))=>Columns(@NT(b = String[])))
    @test eltype(collect_columns(v)) == Pair{NamedTuples._NT_a{Int}, NamedTuples._NT_b{String}}

    v = Iterators.filter(t -> t.first.a == 4, (@NT(a=i) => @NT(b="a$i") for i in 1:3))
    @test collect_columns(v) == Columns(Columns(@NT(a = Int[]))=>Columns(@NT(b = String[])))
    @test eltype(collect_columns(v)) == Pair{NamedTuples._NT_a{Int}, NamedTuples._NT_b{String}}

    t = table(@NT(b=1) => @NT(a = i) for i in (2, DataValue{Int}(), 3))
    @test t == table(@NT(b = [1,1,1], a = [2, DataValue{Int}(), 3]), pkey = :b)
end

@testset "issubtype" begin
    @test IndexedTables._is_subtype(Int, Int)
    @test IndexedTables._is_subtype(Int, DataValue{Int})
    @test !IndexedTables._is_subtype(DataValue{Int}, Int)
    @test IndexedTables._is_subtype(DataValue{Int}, DataValue{Int})
    @test !IndexedTables._is_subtype(DataValue{Int}, DataValue{String})
    @test !IndexedTables._is_subtype(Int, String)
end

@testset "collectflattened" begin
    t = [(:a => [1, 2]), (:b => [1, 3])]
    @test collect_columns_flattened(t) == Columns([:a, :a, :b, :b] => [1, 2, 1, 3])
    t = ([@NT(a = 1), @NT(a = 2)], [@NT(a = 1.1), @NT(a = 2.2)])
    @test collect_columns_flattened(t) == Columns(a = [1, 2, 1.1, 2.2])
    @test eltype(collect_columns_flattened(t)) == typeof(@NT(a=1.1))
    t = [(:a => table(1:2, ["a", "b"])), (:b => table(3:4, ["c", "d"]))]
    @test table(collect_columns_flattened(t)) == table([:a, :a, :b, :b], 1:4, ["a", "b", "c", "d"], pkey = 1)
end


text:
module Questionnaires

using Codex
using CSV
using DataFrames
using Dates

struct Items
    normal::Array{Int,1}
    reversed::Array{Int,1}
end

const NATO_NAMES = [
    "alfa",
    "bravo",
    "charlie",
    "delta",
    "echo",
    "foxtrot",
    "golf",
    "hotel",
    "india",
    "julliet",
    "kilo",
    "lima",
    "mike"
]

reverse(x::Int) = 6 - x
reverse(x::Missing) = missing

function get_scores(df::DataFrame, items::Items; average=true)::Array
    if length(items.normal) != 0
        v_normal = [Symbol("v$i") for i in items.normal]
        normal_scores = [sum(row[v_normal]) for row in eachrow(df[:, v_normal])]
    else
        normal_scores = repeat([0], nrow(df))
    end

    if length(items.reversed) != 0
        v_reversed = [Symbol("v$i") for i in items.reversed]
        reversed_scores = [sum(map(reverse, row)) for row in eachrow(df[:, v_reversed])]
    else
        reversed_scores = repeat([0], nrow(df))
    end

    totals = normal_scores .+ reversed_scores
    n_items = length(items.normal) + length(items.reversed)

    return average ? totals ./ n_items : totals
end

include("demographics.jl")
include("resilience.jl")
include("commitment.jl")
include("self-efficacy.jl")
include("personality.jl")
include("intelligence.jl")
include("toughness.jl")
include("optimism.jl")
include("coping.jl")
include("inspire.jl")

"""
    get_hnd(path::AbstractString)::DataFrame

Get HowNutsAreTheDutch data and select big five, age and more.
"""
function get_hnd(path::AbstractString)::DataFrame
    path = string(path)::String
    # Don't need the added complexity of weakrefstrings here.
    df = try
        CSV.read(path, DataFrame; stringtype=String)
    catch
        # CSV version 0.8.
        CSV.read(path, DataFrame)
    end
    # Assumes that all neo scores are missing if neo_neurot is 999.
    filter!(:neo_neurot => !=(999), df)
    df[!, :group] .= "civilians"
    domains = (:neo_neurot => :N, :neo_extraversion => :E, :neo_openness => :O,
        :neo_agreeable => :A, :neo_conscient => :C)
    select!(df, :id, :age, :start_educ => :education, domains...)
    @assert nrow(df) == 4_984
    df
end

transformation_map = Dict{String,Function}(
    "bravo" => unify_demographics,
    "charlie" => resilience2scores,
    "delta" => Commitment.delta2scores,
    "echo" => self_efficacy2scores,
    "foxtrot" => Intelligence.foxtrot2scores,
    "golf" => Intelligence.golf2scores,
    "india" => Toughness.india2scores,
    "julliet" => Coping.julliet2scores,
    "kilo" => Optimism.kilo2scores,
    "lima" => personality2scores,
    "mike" => Inspire.mike2scores
)

"""
    responses(data_dir::String, nato_name::String)::DataFrame

Responses for questionnaire `nato_name` as contained in directory `data_dir`.
Returns a DataFrame with rows `{ id, r...}` where `id` is a long identifier and not the one from the backend.
"""
function responses(data_dir::String, nato_name::String)::DataFrame
    responses_dir = joinpath(data_dir, "responses")
    responses_file = joinpath(responses_dir, "$nato_name.csv")
    responses_data = try
        CSV.read(responses_file, DataFrame; delim=';', stringtype=String)
    catch
        # CSV.jl 0.8.
        CSV.read(responses_file, DataFrame; delim=';')
    end
    # To speed up further processing.
    Codex.TransformExport.rm_timing!(responses_data)
    Codex.TransformExport.rm_boring_timestamps!(responses_data)
    Codex.TransformExport.rm_boring_foreign_keys!(responses_data)

    people_file = joinpath(data_dir, "people.csv")
    people_data = Codex.TransformExport.read_csv(people_file, delim=';')

    if nato_name != "mike" || !contains(data_dir, "2018")
        select!(people_data, :person_id => :backend_id, :first_name => :id)
        rename!(responses_data, Dict(:filled_out_by_id => "backend_id"))
        # If people_data is free from missings, then matchmissing equal should not introduce
        # missings, I think.
        dropmissing!(people_data)
        joined = nothing
        joined = innerjoin(people_data, responses_data, on=:backend_id, matchmissing=:equal)
        select!(joined, Not(:backend_id))
    else
        joined = responses_data
    end

    if "locale" in names(joined)
        select!(joined, Not(:locale))
    end

    if nato_name in keys(transformation_map)
        f = transformation_map[nato_name]
        joined = f(joined)
    end
    return joined
end

remove_auth_prefix(s) = ismissing(s) ? missing : s[7:end]

"""
    replace_usernames(df::DataFrame, id_col::Symbol, id_username::DataFrame)

For every username in `df[!, id_col]`, replace the name by the long identifier.
The long identifier is necessary, because usernames where not anonymized in 2018.
"""
function replace_usernames(df::DataFrame, id_col::Symbol, id_username::DataFrame)
    matchmissing = :notequal
    on = id_col => :username
    makeunique = true
    df = leftjoin(df, id_username; on, matchmissing, makeunique)
    new_col = Symbol("$(id_col)_1")
    df[!, new_col] = remove_auth_prefix.(df[!, new_col])
    df[!, :id] = [ismissing(new_id) ? id : new_id for (id, new_id) in zip(df[!, :id], df[!, new_col])]
    select!(df, Not(new_col))
    df
end

"""
    dropouts(raw_dir::String)

Returns dropout data where all IDs are in the long identifier format.
On 2021-07-19, returned a 197x7 dataset.
"""
function dropouts(raw_dir::String)
    dropouts_file = joinpath(raw_dir, "dropouts.csv")
    dropouts = Codex.TransformExport.read_csv(dropouts_file; delim=';')

    id_username_file = joinpath(raw_dir, "id-username.csv")
    id_username = Codex.TransformExport.read_csv(id_username_file; delim=',')

    dropouts = replace_usernames(dropouts, :id, id_username)
    transform!(dropouts, :dropout => ByRow(Bool) => :dropout)
    msg = "Unable to convert all usernames. Got:\n$dropouts"
    @assert all([length(id) == 24 for id in dropouts.id]) msg
    dropouts
end

function clean_graduates_dropouts(responses::DataFrame, dropouts::DataFrame, group, cohort)
    @assert group != "operators"

    responses = transform(responses, :id => ByRow(remove_auth_prefix) => :id)
    df = leftjoin(responses, dropouts; on=:id)

    function filter_group(dropout::Bool, dropout_reason::Union{Missing,String})
        if group == "graduates"
            !dropout
        else
            # Read as: if the query asks for "dropouts-medical", then get data with:
            medical_reason_match() = group == "dropouts-medical" ?
                dropout_reason == "B" || dropout_reason == "medisch" :
                dropout_reason == "A" || dropout_reason == "C" || dropout_reason == "ontheffing"
            dropout && medical_reason_match()
        end
    end
    function filter_group(dropout::Missing, dropout_reason)
        false
    end

    df = subset!(df, [:dropout, :dropout_reason] => ByRow(filter_group))
    df[!, :group] .= group
    df = select!(df, :group, names(responses)...)
end

"""
    validate_dropouts(df::DataFrame)

Validate "dropouts.csv".
"""
function validate_dropouts(df::DataFrame)
    disallowmissing!(df, :dropout)
    valid(dropout_reason::Missing, dropout::Bool) = !dropout
    valid(dropout_reason::String, dropout::Bool) = dropout
    invs = subset(df, [:dropout_reason, :dropout] => ByRow(!valid))
    @assert nrow(invs) == 0 "Found invalid rows:\n$invs"
end

"""
    responses(data_dir::String, nato_name::String, group::String; measurement=999)::DataFrame

Responses for group `group` and measurement `measurement`, where `group` is one of `graduates`, `operators`, `dropouts-medical` or `dropouts-non-medical`.
`measurement` is only used to split the 2018 data, for the later datasets it is ignored.

## Example in 2021-07:

```
julia> dir = joinpath(ysf_raw, "2020-first");

julia> df = Codex.Questionnaires.responses(dir, "kilo", "dropouts");

julia> nrow(df)
25

julia> select(first(df, 5), Not([:id, :completed_at]))
5×3 DataFrame
 Row │ group     optimism  pessimism
     │ String    Int64     Int64
─────┼───────────────────────────────
   1 │ dropouts        10         10
   2 │ dropouts        13          6
   3 │ dropouts        10          9
   4 │ dropouts        12          6
   5 │ dropouts        12          5
```
"""
function responses(data_dir::String, nato_name::String, group::String; measurement=999)
    responses_data = responses(data_dir, nato_name)

    cohort = parse(Int, match(r"[0-9]{4}", data_dir).match)
    dropouts_data = dropouts(dirname(data_dir))
    validate_dropouts(dropouts_data)

    if group == "operators"
        df = responses_data
        df[:, :group] .= "operators"
    else # Graduates and dropouts.
        df = clean_graduates_dropouts(responses_data, dropouts_data, group, cohort)
    end

    if cohort == 2018
        if !(measurement in [1, 2])
            throw(AssertionError("Measurement has to be specified for the 2018 data"))
        end
        threshold = Date("2019-01-01")
        function filter_date(x)::Bool
            date = Date(first(split(x, ' ')), DateFormat("dd-mm-yyyy"))
            measurement == 1 ? date < threshold : threshold < date
        end
        return nato_name != "mike" ? filter([:completed_at] => filter_date, df) : df
    else # After 2018, the datasets are already split on before and after selection.
        df
    end
end

function first_measurement(raw_dir::AbstractString, nato_name::AbstractString)
    parameters = [
        (raw_dir, "2018-first", "graduates"),
        (raw_dir, "2018-first", "dropouts-medical"),
        (raw_dir, "2018-first", "dropouts-non-medical"),
        (raw_dir, "2019-first", "graduates"),
        (raw_dir, "2019-first", "dropouts-medical"),
        (raw_dir, "2019-first", "dropouts-non-medical"),
        (raw_dir, "2020-operators", "operators"),
        (raw_dir, "2020-first", "graduates"),
        (raw_dir, "2020-first", "dropouts-medical"),
        (raw_dir, "2020-first", "dropouts-non-medical"),
        (raw_dir, "2021-march-first", "graduates"),
        (raw_dir, "2021-march-first", "dropouts-medical"),
        (raw_dir, "2021-march-first", "dropouts-non-medical"),
    ]
    measurement = 1
    function helper(dir, cohort_dir, group)
        data_dir = joinpath(dir, cohort_dir)
        responses(data_dir, nato_name, group; measurement)
    end
    vcat([helper(p...) for p in parameters]...)
end

function second_measurement(raw_dir::AbstractString, nato_name::AbstractString)
    parameters = [
        (raw_dir, "2018-second", "graduates", 2),
        (raw_dir, "2018-second", "dropouts-medical", 2),
        (raw_dir, "2018-second", "dropouts-non-medical", 2),
        (raw_dir, "2019-second", "graduates", 2),
        (raw_dir, "2019-second", "dropouts-medical", 2),
        (raw_dir, "2019-second", "dropouts-non-medical", 2),

#        (raw_dir, "2020-first", "dropouts-medical", 1),
#        (raw_dir, "2020-first", "dropouts-non-medical", 1),
    ]
    if nato_name == "mike"
        filter!(p -> p[2] != "2018-second", parameters)
    end
    helper(dir, cohort_dir, group, measurement::Int) =
        responses(joinpath(dir, cohort_dir), nato_name, group; measurement)
    vcat([helper(p...) for p in parameters]...)

end

"""
    join_questionnaires(raw_dir::String, questionnaires::Array{String,1}, groups::Array{String,1})::DataFrame

Joines multiple `questionnaires` for the members of `groups`, where `groups` ensures that the joins do not remove rows which shouldn't be removed.
For instance, given questionnaires `"A"` and `"B"` containing respectively columns `:a1`, `:a2` and `:b1`, this method will return `{ group, id, A_a1, A_a2, B_b1 }`.
"""
function join_questionnaires(raw_dir::String, questionnaires::Array{String,1}, groups::Array{String,1})::DataFrame
    function prepare_responses(q::String)::DataFrame
        df = first_measurement(raw_dir, q)
        filter!(:group => x -> x in groups, df)
        select!(df, Not(:completed_at))
    end

    function ysf_join(A::DataFrame, b::String)
        B = prepare_responses(b)
        # Group information is already provided by `A`.
        select!(B, Not(:group))
        rename!(s -> s == "id" ? s : "$(b)_$s", B)
        innerjoin(A, B, on = :id)
    end

    function ysf_join(a::String, b::String)
        A = prepare_responses(a)
        rename!(s -> s == "id" || s == "group" ? s : "$(a)_$s", A)
        ysf_join(A, b)
    end

    reduce(ysf_join, questionnaires)
end

"""
    join_vo_questionnaires(raw_dir::String)::DataFrame

Combine information from multiple questionnaires to allow model fitting.
"""
function join_vo_questionnaires(raw_dir::String)::DataFrame
    questionnaires = sort(collect(keys(transformation_map)))
    # Ignoring delta since only two operators participated in it.
    questionnaires = filter(!in(["bravo", "delta"]), questionnaires)
    df = Codex.Questionnaires.join_questionnaires(
        raw_dir,
        questionnaires,
        ["graduates", "operators", "dropouts-non-medical"]
    )
    df[:, :binary_group] = [x == "graduates" || x == "operators" ? 1 : 0 for x in df[:, :group]]
    df
end

function csv_files(dir::String)
    return filter(endswith(".csv"), readdir(dir))
end

"""
    unfinished_questionnaires(responses_dir::String, id::String)

Return questionnaire names in `responses_dir` for which `id` is not in `filled_out_by_id` column.
This assumes that unfinished questionnaires have a `missing` in the `filled_out_by_id` column.
"""
function unfinished_questionnaires(responses_dir::String, id::String)
    csvs = csv_files(responses_dir)
    unfinished_csv_files = filter(csvs) do csv_file
        csv_path = joinpath(responses_dir, csv_file)
        df = CSV.read(csv_path, DataFrame)
        if "filled_out_by_id" in names(df)
            return !(id in skipmissing(df.filled_out_by_id))
        else
            # We can ignore the dataset since it is not a valid u-can-act output.
            return true
        end
    end
    unfinished = first.(splitext.(unfinished_csv_files))
    return unfinished
end

"""
    all_ids(responses_dir)::Vector{String}

"""
function all_ids(responses_dir)::Vector{String}
    csvs = csv_files(responses_dir)
    ids = map(csvs) do csv_file
        csv_path = joinpath(responses_dir, csv_file)
        df = CSV.read(csv_path, DataFrame)
        if "filled_out_by_id" in names(df)
            return collect(skipmissing(df.filled_out_by_id))
        else
            return String[]
        end
    end
    return sort(unique(Iterators.flatten(ids)))
end

struct Unfinished
    id::String
    unfinished_questionnaires::Vector{String}
end

"""
    unfinished_info(responses_dir; required::Union{Nothing,Vector{String}}=missing)::Vector{Unfinished}

Return the ids for which not all questionnaires in `required` have been filled in.
When `ismissing(required)`, take all the questionnaires in `responses_dir`.

# Example
```
julia> responses_dir = joinpath(homedir(), "git", "ysf-raw", "2021-08", "responses");

julia> required = ["alfa", "bravo", "charlie", "delta", "echo", "foxtrot", "golf", "hotel", "india", "julliet", "kilo", "lima", "mike"];

julia> Codex.Questionnaires.unfinished_info(responses_dir; required)

"""
function unfinished_info(responses_dir; required::Union{Missing,Vector{String}}=missing)::Vector{Unfinished}
    if ismissing(required)
        required = first.(splitext.(csv_files(responses_dir)))
    end
    ids = all_ids(responses_dir)
    U = [unfinished_questionnaires(responses_dir, id) for id in ids]
    U = [filter(in(required), unfinished) for unfinished in U]
    Z = collect(zip(ids, U))
    Z = filter(t -> !isempty(last(t)), Z)
    return [Unfinished(id, unfinished) for (id, unfinished) in Z]
end

end # module


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    RefinementMethod

A method for refining meshes.
"""
abstract type RefinementMethod end

"""
    refine(mesh, method)

Refine `mesh` with refinement `method`.
"""
function refine end

# ----------------
# IMPLEMENTATIONS
# ----------------

include("refinement/quad.jl")
include("refinement/catmullclark.jl")


text:
struct stCoeffMX{d,mt<:CoefficientMatrix{d}} <: CoefficientMatrix{d}
    nID::Int64
    cMX::mt
end
(stcm::stCoeffMX)(t) = stcm.cMX(t)

struct stAdditive{d,T<:AdditiveVector{d}} <: AdditiveVector{d}
    nID::Int64
    V::T
end
(stV::stAdditive)(t) = stV.V(t)

# Linear Delay Differential Equation Problem
struct LDDEProblem{d,AT <: ProportionalMX{d,<:MatrixOrFunction},BT <: DelayMX{d,<:Any,<:Any}, cT <: Additive{d,<:VectorOrFunction}, αT<:stCoeffMX{d,<:ProportionalMX{d,<:Any}}, βT<:stCoeffMX{d,<:DelayMX{d,<:Any,<:Any}}, σT<:AdditiveVector{d}} <: AbstractLDDEProblem{d}
    A::AT # Coefficient of the proportional term (present /non-discretised term)
    Bs::Vector{BT} # Coefficient of the delayed terms
    αs::Vector{αT} # Noise from the delay matrix (only single delay!!!)
    βs::Vector{βT} # Noise from the delay matrix (only single delay!!!)
    c::cT # Additive vector
    σs::Vector{σT} # Additive vectors for the noises
    w::Int64 # Number of noise sources
end

function LDDEProblem(A::ProportionalMX, Bs::Vector{<:DelayMX}, αs::Vector{<:stCoeffMX}, βs::Vector{<:stCoeffMX}, c::Additive=Additive(size(A, 2)), σs::Vector{<:stAdditive}=Vector{AdditiveVector}(undef,0))
    w::Int64 = all(isempty.((αs,βs,σs))) ? 0 : maximum(αβ.nID for αβ in (αs...,βs...,σs...))
    LDDEProblem(A, Bs, αs, βs, c, σs, w)
end


text:
# if nworkers() > 1
#     rmprocs(workers())
# end
#
# if Base.JLOptions().code_coverage == 1
#     addprocs(Sys.CPU_CORES - 3 , exeflags = ["--code-coverage=user", "--inline=no", "--check-bounds=yes"])
# else
#     addprocs(Sys.CPU_CORES - 3, exeflags = "--check-bounds=yes")
# end
# @show nprocs()

using Simplices
using SimplexSplitting
using InvariantDistribution
using Distributions
using Base.Test

include("mm_dd_all.jl")
include("invdist.jl")


text:
module SumOfSquaresOptimization
import Base.Operators: +, *
using DataStructures # for SortedDict, for row-redundancy checks
import DataStructures: status # annoying that we have to overload this with a different meaning

# sos.jl
export Program
export constraint, partconstraint, objective, partobjective
export @constraint, @partconstraint, @objective, @partobjective
export perturb_objective, perturb_objective_sparse

# solve.jl
export sossolve

# solution.jl
export SoSSolution
export dumpsol, obj, primalobj, dualobj, status
export moment, @moment

# symmetrize.jl
export symmetrize!
export symmetrize_cyclic!, symmetrize_dihedral!, symmetrize_hypercube!, symmetrize_full!

include("monoms.jl")
include("program.jl")
include("sdp.jl")
include("solve.jl")
include("solution.jl")
include("symmetrize.jl")

end # module SumOfSquaresOptimization

max_stars_repo_name: willperry/SumOfSquaresOptimization


text:
""" Add from private repo folder once "copyright" situatoin is resolved"""
function maximal_RPI_set() end

max_stars_repo_name: mforets/InvariantSets

max_issues_repo_name: mforets/InvariantSets


text:
module BasicAuthRequest

using ..Base64

import ..Layer, ..request
using ..URIs
using ..Pairs: getkv, setkv
import ..@debug, ..DEBUG_LEVEL

"""
    request(BasicAuthLayer, method, ::URI, headers, body) -> HTTP.Response

Add `Authorization: Basic` header using credentials from url userinfo.
"""
abstract type BasicAuthLayer{Next <: Layer} <: Layer{Next} end
export BasicAuthLayer

function request(::Type{BasicAuthLayer{Next}},
                 method::String, url::URI, headers, body; kw...) where Next

    userinfo = unescapeuri(url.userinfo)

    if !isempty(userinfo) && getkv(headers, "Authorization", "") == ""
        @debug 1 "Adding Authorization: Basic header."
        setkv(headers, "Authorization", "Basic $(base64encode(userinfo))")
    end

    return request(Next, method, url, headers, body; kw...)
end

end # module BasicAuthRequest


text:
module GeoMapping

export distance, azimuth, reckon

"""
    d = distance(lat1,lon1,lat2,lon2)

Compute the great-circle distance between the points (`lat1`,`lon1`) and (`lat2`,`lon2`).
The units of all input and output parameters are degrees.
"""
function distance(lat1,lon1,lat2,lon2)
    #https://en.wikipedia.org/w/index.php?title=Great-circle_distance&oldid=749078136#Computational_formulas

    Δλ = π/180 * (lon2 - lon1)
    ϕ1 = π/180 * lat1
    ϕ2 = π/180 * lat2

    cosΔσ = sin(ϕ1)*sin(ϕ2) + cos(ϕ1)*cos(ϕ2)*cos(Δλ)

    eins = one(cosΔσ)
    cosΔσ = max(min(cosΔσ,eins),-eins)
    Δσ = acos(cosΔσ)
    return 180/π * Δσ
end

"""
    az = azimuth(lat1,lon1,lat2,lon2)

Compute azimuth, i.e. the angle between the line segment defined by the points (`lat1`,`lon1`) and (`lat2`,`lon2`)
and the North.
The units of all input and output parameters are degrees.
"""
function azimuth(lat1,lon1,lat2,lon2)
    # https://en.wikipedia.org/w/index.php?title=Azimuth&oldid=750059816#Calculating_azimuth

    Δλ = π/180 * (lon2 - lon1)
    ϕ1 = π/180 * lat1
    ϕ2 = π/180 * lat2

    α = atan(sin(Δλ), cos(ϕ1)*tan(ϕ2) - sin(ϕ1)*cos(Δλ))
    return 180/π * α
end

# Base on reckon from myself
# https://sourceforge.net/p/octave/mapping/ci/3f19801d4b93d3b3923df9fa62d268660e5cb4fa/tree/inst/reckon.m
# relicenced to LGPL-v3

"""
    lato,lono = reckon(lat,lon,range,azimuth)

Compute the coordinates of the end-point of a displacement on a
sphere. `lat`,`lon` are the coordinates of the starting point, `range`
is the covered distance of the displacements along a great circle and
`azimuth` is the direction of the displacement relative to the North.
The units of all input and output parameters are degrees.
This function can also be used to define a spherical coordinate system
with rotated poles.
"""
function reckon(lat,lon,range,azimuth)

    # convert to radian
    rad2deg = π/180

    lat = lat*rad2deg
    lon = lon*rad2deg
    range = range*rad2deg
    azimuth = azimuth*rad2deg

    tmp = sin.(lat).*cos.(range) + cos.(lat).*sin.(range).*cos.(azimuth)

    # clip tmp to -1 and 1
    eins = one(eltype(tmp))
    tmp = max.(min.(tmp,eins),-eins)

    lato = π/2 .- acos.(tmp)

    cos_γ = (cos.(range) - sin.(lato).*sin.(lat))./(cos.(lato).*cos.(lat))
    sin_γ = sin.(azimuth).*sin.(range)./cos.(lato)

    γ = atan.(sin_γ,cos_γ)

    lono = lon .+ γ

    # bring the lono in the interval [-π π[

    lono = mod.(lono .+ π,2*π) .- π

    # convert to degrees

    lono = lono/rad2deg
    lato = lato/rad2deg

    return lato,lono
end

end

max_stars_repo_name: rctngle/OpenWeatherDecoder

max_issues_repo_name: rctngle/OpenWeatherDecoder

max_forks_repo_name: rctngle/OpenWeatherDecoder