text:
"""
$(TYPEDEF)

Struct to hold extension data of sparse matrix.

$(TYPEDFIELDS)
"""
mutable struct SparseMatrixExtension{Tv,Ti<:Integer} <: AbstractSparseMatrix{Tv,Ti}

    """
    Number of rows
    """
    m::Ti

    """
    Number of columns
    """
    n::Ti

    """
    Number of nonzeros
    """
    nnz::Ti

    """
    Linked list of column entries. Initial length is n,
    it grows with each new entry.

    colptr[index] contains the next
    index in the list or zero, terminating the list which
    starts at index 1<=j<=n for each column j.
    """
    colptr::Vector{Ti}

    """
    Row numbers. For each index it contains the zero (initial state)
    or the row numbers corresponding to the column entry list in colptr.

    Initial length is n,
    it grows with each new entry.
    """
    rowval::Vector{Ti}

    """
    Nonzero entry values correspondin to each pair
    (colptr[index],rowval[index])

    Initial length is n,  it grows with each new entry.
    """
    nzval::Vector{Tv}
end

"""
$(TYPEDSIGNATURES)

Constructor of empty extension
"""
SparseMatrixExtension{Tv,Ti}(m::Integer, n::Integer)  where {Tv,Ti<:Integer} =    SparseMatrixExtension{Tv,Ti}(m,n,0,zeros(Ti,n),zeros(Ti,n),zeros(Tv,n))

"""
$(SIGNATURES)

Return value stored for entry or zero if not found
"""
function Base.getindex(E::SparseMatrixExtension{Tv,Ti},i::Integer, j::Integer) where {Tv,Ti<:Integer}

    if !((1 <= i <= E.m) & (1 <= j <= E.n))
        throw(BoundsError(E, (i,j)))
    end

    k=j
    while k>0
        if E.rowval[k]==i
            return E.nzval[k]
        end
        k=E.colptr[k]
    end

    return zero(Tv)
end

"""
$(SIGNATURES)

Update value of existing entry, otherwise extend matrix.
"""
function Base.setindex!(E::SparseMatrixExtension{Tv,Ti}, _v, _i::Integer, _j::Integer) where {Tv,Ti<:Integer}
    v = convert(Tv, _v)
    i = convert(Ti, _i)
    j = convert(Ti, _j)

    if !((1 <= i <= E.m) & (1 <= j <= E.n))
        throw(BoundsError(E, (i,j)))
    end

    if iszero(v)
        return E
    end

    # Set the first  column entry if it was not yet set.
    if E.rowval[j]==0
        E.rowval[j]=i
        E.nzval[j]=v
        E.nnz+=1
        return E
    end

    # Traverse list for existing entry
    k=j
    k0=j
    while k>0
        # Update value and return if entry has been found
        if E.rowval[k]==i
            E.nzval[k]=v
            return E
        end
        k0=k
        # Next element in the list
        k=E.colptr[k]
    end

    # Append entry if not found
    push!(E.nzval,v)
    push!(E.rowval,i)

    # Shift the end of the list
    push!(E.colptr,0)
    E.colptr[k0]=length(E.nzval)

    # Update number of nonzero entries
    E.nnz+=1
    return E
end

"""
$(TYPEDSIGNATURES)

Return tuple containing size of the matrix.
"""
Base.size(E::SparseMatrixExtension) = (E.m, E.n)

"""
$(TYPEDSIGNATURES)

Return number of nonzero entries.
"""
SparseArrays.nnz(E::SparseMatrixExtension)=E.nnz

"""
$(TYPEDSIGNATURES)

Dummy flush! method for Sparse matrix extension. Just
used in thest methods
"""
function flush!(M::SparseMatrixExtension{Tv, Ti}) where{Tv, Ti}
    return M
end


text:
# tests for the various kinds of waits
using Base.Test
using MPI

MPI.Init()
myrank = MPI.Comm_rank(MPI.COMM_WORLD)
commsize = MPI.Comm_rank(MPI.COMM_WORLD)

nsends = 2
send_arr = Array(Array{Int, 1}, nsends)
recv_arr = Array(Array{Int, 1}, nsends)

for i=1:nsends
    send_arr[i] = [i]
    recv_arr[i] = Array(Int, 1)
end

send_reqs = Array(MPI.Request, nsends)
recv_reqs = Array(MPI.Request, nsends)

# send to self
for i=1:nsends
    send_reqs[i] = MPI.Isend(send_arr[i], myrank, i, MPI.COMM_WORLD)
    recv_reqs[i] = MPI.Irecv!(recv_arr[i], myrank, i, MPI.COMM_WORLD)
end

send_check = zeros(Int, nsends)
recv_check = zeros(Int, nsends)
for i=1:nsends
    idx, stat = MPI.Waitany!(send_reqs)
    send_check[idx] += 1
    send_reqs[idx] = MPI.REQUEST_NULL
end

for i=1:nsends
  @test send_check[i] == 1
end

for i=1:nsends
    idx, stat = MPI.Waitany!(recv_reqs)
    recv_check[idx] += 1
    recv_reqs[idx] = MPI.REQUEST_NULL
end

for i=1:nsends
  @test recv_check[i] == 1
end

MPI.Barrier(MPI.COMM_WORLD)
MPI.Finalize()

ext: jl


text:
using Test, BellScenario, LinearAlgebra

@testset "./src/certifying_signaling_dimension.jl" begin

using SignalingDimension

@testset "maximum_likelihood_lower_bound()" begin
    @test 3 == maximum_likelihood_lower_bound(Strategy([1 0 0;0 0.5 0;0 0 0.5;0 0.5 0.5]))
    @test 1 == maximum_likelihood_lower_bound(Strategy(1/5*ones(5,5)))
    @test 2 == maximum_likelihood_lower_bound(DeterministicStrategy([1 0 1 0;0 1 0 1;0 0 0 0;0 0 0 0]))

    @test 1 == maximum_likelihood_lower_bound(Strategy(ones(1000,1000)/1000))
end

@testset "ambiguous_lower_bound()" begin
    @testset "the number of ambiguous rows affects the lower bound" begin
        P = Strategy([
            0.35 0    0;
            0    0.35 0;
            0    0    0.35;
            0.35 0.35 0.35;
            0.3  0.3  0.3;
        ])

        @test 2 == ambiguous_lower_bound(P, 4)
        @test 3 == ambiguous_lower_bound(P, 3)

        # test over all numbers of ambiguous rows
        @test 3 == ambiguous_lower_bound(P)
    end

    @testset "witnesses d=2 signaling dimension" begin
        P = DeterministicStrategy([1 0 0 1;0 1 1 0;0 0 0 0;0 0 0 0])

        @test 2 == ambiguous_lower_bound(P)
    end

    @testset "pre-sorting edge cases" begin
        @testset "edge case: primary sort by row sum, secondary by max-min diff is not optimal" begin
            P1 = Strategy([
                0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1;
                0.9 0   0   0   0   0   0   0  ;
                0   0.9 0   0   0   0   0   0  ;
                0   0   0.9 0   0   0   0   0  ;
                0   0   0   0.9 0   0   0   0  ;
                0   0   0   0   0.9 0   0   0  ;
                0   0   0   0   0   0.9 0   0  ;
                0   0   0   0   0   0   0.9 0  ;
                0   0   0   0   0   0   0   0.9;
            ])

            P2 = Strategy([
                0.9 0   0   0   0   0   0   0  ;
                0   0.9 0   0   0   0   0   0  ;
                0   0   0.9 0   0   0   0   0  ;
                0   0   0   0.9 0   0   0   0  ;
                0   0   0   0   0.9 0   0   0  ;
                0   0   0   0   0   0.9 0   0  ;
                0   0   0   0   0   0   0.9 0  ;
                0   0   0   0   0   0   0   0.9;
                0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1;
            ])

            @test 8 == ambiguous_lower_bound(P1,8)
            @test 8 == ambiguous_lower_bound(P2,8)
        end

        @testset "edge case: primary sort by row sum, secondary by max-min diff is optimal" begin
            # primary sort by max - min diff, secondary by row sum
            P1a  = Strategy([
                0   0.1 0.5 0;
                0.6 0.6 0.5 1;
                0   0.3 0   0;
                0.4 0   0   0;
            ])
            P1b = Strategy([
                0.4 0.4 1 1;
                0.3 0   0 0;
                0.3 0   0 0;
                0   0.3 0 0;
                0   0.3 0 0;
            ])

            #  primary sort by row sum, secondary by max-min diff
            P2a = Strategy([
                0   0.3 0   0;
                0.4 0   0   0;
                0   0.1 0.5 0;
                0.6 0.6 0.5 1;
            ])
            P2b = Strategy([
                0.3 0   0 0;
                0.3 0   0 0;
                0   0.3 0 0;
                0   0.3 0 0;
                0.4 0.4 1 1;
            ])

            @test 3 == ambiguous_lower_bound(P1a,3)
            @test 3 == ambiguous_lower_bound(P2a,3)

            @test 3 == ambiguous_lower_bound(P1b,4)
            @test 3 == ambiguous_lower_bound(P2b,4)
        end
    end

    @testset "errors" begin
        P = Strategy([1 0 0;0 1 0;0 0 1])

        @test_throws DomainError ambiguous_lower_bound(P, 4)
        @test_throws DomainError ambiguous_lower_bound(P, 0)
    end
end

@testset "trivial_upper_bound()" begin
    @test 3 == trivial_upper_bound(Strategy([1 0 1;0 1 0;0 0 0;0 0 0]))
    @test 4 == trivial_upper_bound(Strategy([1 0 0 0 0 0;0 1 0 0 0 0;0 0 1 0 0 0;0 0 0 1 1 1]))
end

@testset "attains_trivial_upper_bound()"  begin
    @test attains_trivial_upper_bound(Strategy([1 0 0;0 1 0;0 0 1]))
    @test !attains_trivial_upper_bound(Strategy([1 0 1;0 1 0;0 0 0]))

    @testset "strategy where ambiguous game must be tested" begin
        P = Strategy(1/2*[1 0 0;0 1 0;0 0 1;1 1 1])
        @test 2 == maximum_likelihood_lower_bound(P)
        @test 3 == ambiguous_lower_bound(P,3)
        @test attains_trivial_upper_bound(P)
    end
end

@testset "upper_bound()" begin
    @test 4 == upper_bound(Strategy([1 0 0 0;0 1 0 0;0 0 1 0;0 0 0 1]))
    @test 3 == upper_bound(Strategy([1 0 1 0;0 1 0 1;0 0 0 0;0 0 0 0;0 0 0 0]))
end

end


text:
#
# This finds Normalized, Synodic, Cartesian states
# for Halo Orbits near L1 & L2 for a variety of
# normalized CR3BP systems.
#
# Currently, nondimensional mass parameters are
# used which are close to those in the Sun-Earth,
# Earth-Moon, Sun-Mars, and Sun-Jupiter systems.
#

using DrWatson; @quickactivate "CR3BP-Manifold-Research"

using CSV
using DataFrames
using GeneralAstrodynamics

massparameters = let
    planets = (Venus, Earth, Mars, Jupiter, Saturn, Uranus, Neptune)
    planet_params = map(
        planet -> let μ₁=massparameter(Sun), μ₂=massparameter(planet)
            min(μ₁, μ₂) / (μ₁+μ₂)
        end, planets
    )
    earth_moon  = let μ₁=massparameter(Earth), μ₂=massparameter(Moon)
        min(μ₁, μ₂) / (μ₁+μ₂)
    end
    vcat(planet_params[1:2]..., earth_moon, planet_params[3:end]...)
end

names = [
    "Sun-Venus", "Sun-Earth",
    "Earth-Moon", "Sun-Mars",
    "Sun-Jupiter", "Sun-Saturn",
    "Sun-Uranus", "Sun-Neptune"
]

amplitudes = 0.0:1e-6:0.01

@assert length(names) == length(massparameters) "Invalid configuration!"

Threads.@threads for i = 1:length(names)
    μ = massparameters[i]
    name = names[i]
    orbits = DataFrame(
        MassParameter=Float64[], LagrangePoint=Int[], ZAmplitude=Float64[], JacobiConstant=Float64[], Period=Float64[],
        Rx=Float64[], Ry=Float64[], Rz=Float64[], Vx=Float64[], Vy=Float64[], Vz=Float64[])
    for L ∈ (1,2)
        fail_count = 0
        for Az ∈ amplitudes
            fail_count < 10 || break
            r, v, T = halo(μ; Az = Az, L = L, max_iter=20, nan_on_fail=true, disable_warnings=true);
            if T ≥ one(T) && all(el->!isnan(el), r) && all(el->!isnan(el), v) && isperiodic(CR3BPOrbit(CartesianState(r,v),μ), T)
                push!(orbits, (μ, L, Az, jacobi_constant(r, v, μ), T, r..., v...))
            else
                fail_count += 1
            end
        end
    end
    @info "Finished processing $name Halos!"
    CSV.write(datadir("exp_pro", "halos", lowercase(string(name,"-Halos.csv"))), orbits)
end

open(datadir("processed", "halos", "README.md"), "w") do io
    write(io,
    """
    # How are the Halo CSV files formatted?
    * The columns from left to right are:
      1. Nondimensional mass parameter
      2. Lagrange point (1 or 2)
      3. Nondimensional Z amplitude
      4. Jacobi Constant
      5. Nondimensional orbital period
      6. Nondimensional X
      7. Nondimensional Y
      8. Nondimensional Z
      9. Nondimensional Ẋ
      10. Nondimensional Ẏ
      11. Nondimensional Ż
    """)
end

max_stars_repo_name: cadojo/Halo-Orbit-Explorations

max_issues_repo_name: cadojo/Halo-Orbit-Explorations

max_forks_repo_name: cadojo/Halo-Orbit-Explorations


text:
"""
$(TYPEDEF)
"""
abstract type AbstractDesign end

"""
$(TYPEDEF)
"""
abstract type AbstractScreeningDesign <: AbstractDesign end

"""
$(TYPEDEF)
"""
abstract type AbstractResponseSurfaceDesign <: AbstractDesign end

"""
$(TYPEDEF)
"""
abstract type AbstractFactorialDesign <: AbstractDesign end

"""
$(TYPEDEF)
"""
abstract type AbstractOptimalDesign <: AbstractDesign end

"""
$(TYPEDEF)
"""
abstract type AbstractRandomDesign <: AbstractDesign end

"""
$(TYPEDEF)

Encapsulates  a  Plackett-Burman  screening  design  constructed  using  Paley's
method.  Factor  levels are  encoded as  `:high` and  `:low` symbols,  and extra
dummy variables, possibly  aliasing interactions between factors,  will be added
to pad a design.

$(TYPEDFIELDS)
"""
struct PlackettBurman <: AbstractScreeningDesign
    matrix::DataFrame
    factors::Tuple
    dummy_factors::Tuple
    formula::FormulaTerm
end

"""
$(TYPEDSIGNATURES)

```jldoctest
julia> PlackettBurman(@formula(y ~ x1 + x2 + x3 + x4))
PlackettBurman
Dimension: (8, 7)
Factors: (:x1, :x2, :x3, :x4)
Dummy Factors: (:dummy1, :dummy2, :dummy3)
Formula: y ~ -1 + x1 + x2 + x3 + x4 + dummy1 + dummy2 + dummy3
Design Matrix:
8×7 DataFrame
 Row │ x1     x2     x3     x4     dummy1  dummy2  dummy3
     │ Int64  Int64  Int64  Int64  Int64   Int64   Int64
─────┼────────────────────────────────────────────────────
   1 │     1      1      1      1       1       1       1
   2 │    -1      1     -1      1       1      -1      -1
   3 │     1     -1      1      1      -1      -1      -1
   4 │    -1      1      1     -1      -1      -1       1
   5 │     1      1     -1     -1      -1       1      -1
   6 │     1     -1     -1     -1       1      -1       1
   7 │    -1     -1     -1      1      -1       1       1
   8 │    -1     -1      1     -1       1       1      -1

```
"""
function PlackettBurman(formula::FormulaTerm; symbol_encoding::Bool = false)
    symbol_factors = Tuple(r.sym for r in formula.rhs)

    initial_design = plackettburman(length(symbol_factors))
    categorical_design = similar(initial_design, Symbol)

    if symbol_encoding
        for i = 1:size(initial_design, 1), j = 1:size(initial_design, 2)
            categorical_design[i, j] = initial_design[i, j] == 1.0 ? :high : :low
        end

        design = DataFrame(categorical_design, :auto)
    else
        design = DataFrame(initial_design, :auto)
    end

    dummy_factors = Tuple(Symbol("dummy" * string(i)) for i = 1:(length(names(design)) -
                                                                 length(symbol_factors)))

    design_names = (symbol_factors..., dummy_factors...)
    rename!(design, collect(design_names))

    PlackettBurman(design,
                   symbol_factors,
                   dummy_factors,
                   term(formula.lhs) ~ term(-1) + sum(term.(design_names)))
end

"""
$(TYPEDSIGNATURES)

```jldoctest
julia> PlackettBurman(4)
PlackettBurman
Dimension: (8, 7)
Factors: (:factor1, :factor2, :factor3, :factor4)
Dummy Factors: (:dummy1, :dummy2, :dummy3)
Formula: 0 ~ -1 + factor1 + factor2 + factor3 + factor4 + dummy1 + dummy2 + dummy3
Design Matrix:
8×7 DataFrame
 Row │ factor1  factor2  factor3  factor4  dummy1  dummy2  dummy3
     │ Int64    Int64    Int64    Int64    Int64   Int64   Int64
─────┼────────────────────────────────────────────────────────────
   1 │       1        1        1        1       1       1       1
   2 │      -1        1       -1        1       1      -1      -1
   3 │       1       -1        1        1      -1      -1      -1
   4 │      -1        1        1       -1      -1      -1       1
   5 │       1        1       -1       -1      -1       1      -1
   6 │       1       -1       -1       -1       1      -1       1
   7 │      -1       -1       -1        1      -1       1       1
   8 │      -1       -1        1       -1       1       1      -1

```
"""
function PlackettBurman(factors::Int)
    PlackettBurman(ConstantTerm(0) ~ sum(term.(Symbol("factor" * string(i)) for i = 1:factors)))
end

"""
$(TYPEDEF)

Box-Behnken design for response surface methodology.

$(TYPEDFIELDS)
"""
struct BoxBehnken <: AbstractResponseSurfaceDesign
    matrix::DataFrame
    factors::Tuple
    formula::FormulaTerm
end

"""
$(TYPEDSIGNATURES)

```jldoctest
julia> BoxBehnken(@formula( y ~ -1 + factor1 & factor1 + factor2 & factor2 + factor3 & factor3 + factor1 & factor2 + factor1 & factor3 + factor2 & factor3 + factor1 + factor2 + factor3), center=3)
BoxBehnken
Dimension: (15, 3)
Factors: (:factor1, :factor2, :factor3)
Formula: y ~ -1 + factor1 + factor2 + factor3 + factor1 & factor1 + factor2 & factor2 + factor3 & factor3 + factor1 & factor2 + factor1 & factor3 + factor2 & factor3
Design Matrix:
15×3 DataFrame
 Row │ factor1  factor2  factor3
     │ Float64  Float64  Float64
─────┼───────────────────────────
   1 │    -1.0     -1.0      0.0
   2 │     1.0     -1.0      0.0
   3 │    -1.0      1.0      0.0
   4 │     1.0      1.0      0.0
   5 │    -1.0      0.0     -1.0
   6 │     1.0      0.0     -1.0
   7 │    -1.0      0.0      1.0
   8 │     1.0      0.0      1.0
   9 │     0.0     -1.0     -1.0
  10 │     0.0      1.0     -1.0
  11 │     0.0     -1.0      1.0
  12 │     0.0      1.0      1.0
  13 │     0.0      0.0      0.0
  14 │     0.0      0.0      0.0
  15 │     0.0      0.0      0.0

```
"""
function BoxBehnken(formula::FormulaTerm; center::Int=1)
    symbol_factors_array = []
    # TODO can only handle interaction terms not functional terms
    for r in formula.rhs
        if r isa InteractionTerm
            for s in r.terms
                append!(symbol_factors_array, [s.sym])
            end
        elseif r isa Term
            append!(symbol_factors_array, [r.sym])
        end
    end
    symbol_factors = Tuple(unique!(symbol_factors_array))

    initial_design = boxbehnken(length(symbol_factors_array), center)
    design = DataFrame(initial_design, :auto)

    rename!(design, collect(symbol_factors))

    BoxBehnken(design, symbol_factors, formula)
end

"""
$(TYPEDSIGNATURES)

```jldoctest
julia> BoxBehnken(3, center=3)
BoxBehnken
Dimension: (15, 3)
Factors: (:factor1, :factor2, :factor3)
Formula: y ~ -1 + factor1 & factor1 + factor2 & factor2 + factor3 & factor3 + factor1 & factor2 + factor1 & factor3 + factor2 & factor3 + factor1 + factor2 + factor3
Design Matrix:
15×3 DataFrame
 Row │ factor1  factor2  factor3
     │ Float64  Float64  Float64
─────┼───────────────────────────
   1 │    -1.0     -1.0      0.0
   2 │     1.0     -1.0      0.0
   3 │    -1.0      1.0      0.0
   4 │     1.0      1.0      0.0
   5 │    -1.0      0.0     -1.0
   6 │     1.0      0.0     -1.0
   7 │    -1.0      0.0      1.0
   8 │     1.0      0.0      1.0
   9 │     0.0     -1.0     -1.0
  10 │     0.0      1.0     -1.0
  11 │     0.0     -1.0      1.0
  12 │     0.0      1.0      1.0
  13 │     0.0      0.0      0.0
  14 │     0.0      0.0      0.0
  15 │     0.0      0.0      0.0

```
"""
function BoxBehnken(factors::Int; center::Int=1)
    factors_tuple = Tuple(Symbol("factor" * string(i)) for i = 1:factors)
    q_terms = sum(InteractionTerm(term.((only(x),only(x)))) for x ∈ combinations(factors_tuple, 1))
    i_terms = sum((&)(term.(x)...) for x ∈ combinations(factors_tuple, 2))
    l_terms = sum(term(only(x)) for x ∈ combinations(factors_tuple, 1))
    formula = term(:y) ~ term(-1) + q_terms + i_terms + l_terms
    BoxBehnken(formula, center=center)
end

"""
$(TYPEDEF)

Encapsulates a central composite design.

$(TYPEDFIELDS)
"""
struct CentralComposite <: AbstractResponseSurfaceDesign
    matrix::DataFrame
    factors::Tuple
    formula::FormulaTerm
    alpha::Symbol
    face::Symbol
end

"""
$(TYPEDSIGNATURES)

```jldoctest
julia> CentralComposite(@formula(y ~ -1 + factor1 & factor1 + factor2 & factor2 + factor3 & factor3 + factor1 & factor2 + factor1 & factor3 + factor2 & factor3 + factor1 + factor2 + factor3))
CentralComposite
Dimension: (22, 3)
Factors: (:factor1, :factor2, :factor3)
Formula: y ~ -1 + factor1 + factor2 + factor3 + factor1 & factor1 + factor2 & factor2 + factor3 & factor3 + factor1 & factor2 + factor1 & factor3 + factor2 & factor3
Alpha: orthogonal
Face: circumscribed
Design Matrix:
22×3 DataFrame
 Row │ factor1   factor2   factor3
     │ Float64   Float64   Float64
─────┼──────────────────────────────
   1 │ -1.0      -1.0      -1.0
   2 │  1.0      -1.0      -1.0
   3 │ -1.0       1.0      -1.0
   4 │  1.0       1.0      -1.0
   5 │ -1.0      -1.0       1.0
   6 │  1.0      -1.0       1.0
   7 │ -1.0       1.0       1.0
   8 │  1.0       1.0       1.0
  ⋮  │    ⋮         ⋮         ⋮
  16 │  0.0       1.82574   0.0
  17 │  0.0       0.0      -1.82574
  18 │  0.0       0.0       1.82574
  19 │  0.0       0.0       0.0
  20 │  0.0       0.0       0.0
  21 │  0.0       0.0       0.0
  22 │  0.0       0.0       0.0
                      7 rows omitted

```
"""
function CentralComposite(formula::FormulaTerm; center::Array{Int}=[4, 4], alpha::Symbol=:orthogonal, face::Symbol=:circumscribed)
    symbol_factors_array = []
    # TODO can only handle interaction terms not functional terms
    for r in formula.rhs
        if r isa InteractionTerm
            for s in r.terms
                append!(symbol_factors_array, [s.sym])
            end
        elseif r isa Term
            append!(symbol_factors_array, [r.sym])
        end
    end
    symbol_factors = Tuple(unique!(symbol_factors_array))

    initial_design = ccdesign(length(symbol_factors_array), center, alpha, face)
    design = DataFrame(initial_design, :auto)
    rename!(design, collect(symbol_factors))

    CentralComposite(design, symbol_factors, formula, alpha, face)
end

"""
$(TYPEDSIGNATURES)

```jldoctest
julia> CentralComposite(3)
CentralComposite
Dimension: (22, 3)
Factors: (:factor1, :factor2, :factor3)
Formula: y ~ -1 + factor1 & factor1 + factor2 & factor2 + factor3 & factor3 + factor1 & factor2 + factor1 & factor3 + factor2 & factor3 + factor1 + factor2 + factor3
Alpha: orthogonal
Face: circumscribed
Design Matrix:
22×3 DataFrame
 Row │ factor1   factor2   factor3
     │ Float64   Float64   Float64
─────┼──────────────────────────────
   1 │ -1.0      -1.0      -1.0
   2 │  1.0      -1.0      -1.0
   3 │ -1.0       1.0      -1.0
   4 │  1.0       1.0      -1.0
   5 │ -1.0      -1.0       1.0
   6 │  1.0      -1.0       1.0
   7 │ -1.0       1.0       1.0
   8 │  1.0       1.0       1.0
  ⋮  │    ⋮         ⋮         ⋮
  16 │  0.0       1.82574   0.0
  17 │  0.0       0.0      -1.82574
  18 │  0.0       0.0       1.82574
  19 │  0.0       0.0       0.0
  20 │  0.0       0.0       0.0
  21 │  0.0       0.0       0.0
  22 │  0.0       0.0       0.0
                      7 rows omitted

```
"""
function CentralComposite(factors::Int; center::Array{Int}=[4, 4], alpha::Symbol=:orthogonal, face::Symbol=:circumscribed)
    factors_tuple = Tuple(Symbol("factor" * string(i)) for i = 1:factors)
    q_terms = sum(InteractionTerm(term.((only(x),only(x)))) for x ∈ combinations(factors_tuple, 1))
    i_terms = sum((&)(term.(x)...) for x ∈ combinations(factors_tuple, 2))
    l_terms = sum(term(only(x)) for x ∈ combinations(factors_tuple, 1))
    formula = term(:y) ~ term(-1) + q_terms + i_terms + l_terms
    CentralComposite(formula, center=center,  alpha=alpha, face=face)
end

"""
$(TYPEDEF)

Encapsulates  an *explicit*  full factorial  design, where  all experiments  are
completely computed  and stored.  Use [`IterableFullFactorial`](@ref)  to obtain
an *implicit* factorial  design with arbitrary access to the  elements of a full
factorial design.

$(TYPEDFIELDS)
"""
struct FullFactorial <: AbstractFactorialDesign
    matrix::DataFrame
    factors::NamedTuple
    formula::FormulaTerm
end

"""
$(TYPEDEF)

$(TYPEDFIELDS)
"""
struct IterableFullFactorial <: AbstractFactorialDesign
    matrix::DataFrame
    factors::NamedTuple
    formula::FormulaTerm
end

"""
$(TYPEDSIGNATURES)

```jldoctest
julia> FullFactorial((A = [1, 2, 4], B = [:a, :b], C = [1.0, -1.0]), @formula(y ~ A + B +C))
FullFactorial
Dimension: (12, 3)
Factors: (A = [1, 2, 4], B = [:a, :b], C = [1.0, -1.0])
Formula: y ~ A + B + C
Design Matrix:
12×3 DataFrame
 Row │ A    B    C
     │ Any  Any  Any
─────┼────────────────
   1 │ 1    a    1.0
   2 │ 2    a    1.0
   3 │ 4    a    1.0
   4 │ 1    b    1.0
   5 │ 2    b    1.0
   6 │ 4    b    1.0
   7 │ 1    a    -1.0
   8 │ 2    a    -1.0
   9 │ 4    a    -1.0
  10 │ 1    b    -1.0
  11 │ 2    b    -1.0
  12 │ 4    b    -1.0

```
"""
function FullFactorial(factors::NamedTuple, formula::FormulaTerm)
    iterator = fullfactorial(values(factors))
    matrix = DataFrame(explicit_fullfactorial(iterator), :auto)
    rename!(matrix, [r.sym for r in formula.rhs])
    FullFactorial(matrix, factors, formula)
end

"""
$(TYPEDSIGNATURES)

```jldoctest
julia> FullFactorial((A = [1, 2, 4], B = [:a, :b], C = [1.0, -1.0]))
FullFactorial
Dimension: (12, 3)
Factors: (A = [1, 2, 4], B = [:a, :b], C = [1.0, -1.0])
Formula: 0 ~ A + B + C
Design Matrix:
12×3 DataFrame
 Row │ A    B    C
     │ Any  Any  Any
─────┼────────────────
   1 │ 1    a    1.0
   2 │ 2    a    1.0
   3 │ 4    a    1.0
   4 │ 1    b    1.0
   5 │ 2    b    1.0
   6 │ 4    b    1.0
   7 │ 1    a    -1.0
   8 │ 2    a    -1.0
   9 │ 4    a    -1.0
  10 │ 1    b    -1.0
  11 │ 2    b    -1.0
  12 │ 4    b    -1.0

```
"""
function FullFactorial(factors::NamedTuple)
    FullFactorial(factors, ConstantTerm(0) ~ sum(term.(keys(factors))))
end

"""
$(TYPEDSIGNATURES)

```jldoctest
julia> FullFactorial(([1, 2, 4], [:a, :b], [1.0, -1.0]))
FullFactorial
Dimension: (12, 3)
Factors: (factor1 = [1, 2, 4], factor2 = [:a, :b], factor3 = [1.0, -1.0])
Formula: 0 ~ factor1 + factor2 + factor3
Design Matrix:
12×3 DataFrame
 Row │ factor1  factor2  factor3
     │ Any      Any      Any
─────┼───────────────────────────
   1 │ 1        a        1.0
   2 │ 2        a        1.0
   3 │ 4        a        1.0
   4 │ 1        b        1.0
   5 │ 2        b        1.0
   6 │ 4        b        1.0
   7 │ 1        a        -1.0
   8 │ 2        a        -1.0
   9 │ 4        a        -1.0
  10 │ 1        b        -1.0
  11 │ 2        b        -1.0
  12 │ 4        b        -1.0

```
"""
function FullFactorial(factors::Tuple)
    FullFactorial(NamedTuple{Tuple(Symbol("factor" * string(i)) for i = 1:length(factors))}(factors))
end

"""
$(TYPEDSIGNATURES)

```jldoctest
julia> FullFactorial(fill([-1, 1], 3))
FullFactorial
Dimension: (8, 3)
Factors: (factor1 = [-1, 1], factor2 = [-1, 1], factor3 = [-1, 1])
Formula: 0 ~ factor1 + factor2 + factor3
Design Matrix:
8×3 DataFrame
 Row │ factor1  factor2  factor3
     │ Int64    Int64    Int64
─────┼───────────────────────────
   1 │      -1       -1       -1
   2 │       1       -1       -1
   3 │      -1        1       -1
   4 │       1        1       -1
   5 │      -1       -1        1
   6 │       1       -1        1
   7 │      -1        1        1
   8 │       1        1        1

```
"""
function FullFactorial(factors::Array)
    FullFactorial(tuple(factors...))
end

"""
$(TYPEDEF)

$(TYPEDFIELDS)
"""
struct FractionalFactorial2Level <: AbstractFactorialDesign
    matrix::DataFrame
    factors::NamedTuple
    formula::FormulaTerm
end

"""
$(TYPEDSIGNATURES)

```jldoctest
julia> FractionalFactorial2Level(@formula(y ~ a + b + a&b + c + a&c+ b&c + a&b&c))
FractionalFactorial2Level
Dimension: (8, 7)
Factors: (a = [-1, 1], b = [-1, 1], c = [-1, 1])
Formula: y ~ a + b + c + a & b + a & c + b & c + a & b & c
Design Matrix:
8×7 DataFrame
 Row │ a      b      c      a_b    a_c    b_c    a_b_c
     │ Int64  Int64  Int64  Int64  Int64  Int64  Int64
─────┼─────────────────────────────────────────────────
   1 │    -1     -1     -1      1      1      1     -1
   2 │     1     -1     -1     -1     -1      1      1
   3 │    -1      1     -1     -1      1     -1      1
   4 │     1      1     -1      1     -1     -1     -1
   5 │    -1     -1      1      1     -1     -1      1
   6 │     1     -1      1     -1      1     -1     -1
   7 │    -1      1      1     -1     -1      1     -1
   8 │     1      1      1      1      1      1      1

```
"""
function FractionalFactorial2Level(formula::FormulaTerm)
    symbol_main_factors_array = []
    symbol_interaction_factors_array = []
    # TODO can only handle simple interaction terms but not functional terms
    for r in formula.rhs
        if r isa InteractionTerm
            append!(symbol_interaction_factors_array, [r])
            for s in r.terms
                append!(symbol_main_factors_array, [s.sym])
            end
        elseif r isa Term
            append!(symbol_main_factors_array, [r.sym])
        end
    end
    symbol_factors = Tuple(unique!(symbol_main_factors_array))
    levels = [[-1,1] for i in symbol_factors]
    A_ff = FullFactorial(NamedTuple{symbol_factors}(levels))
    for (index, interactionterm) in enumerate(symbol_interaction_factors_array)
        interaction_factors = [term.sym for term in interactionterm.terms]
        switch = false
        fac = A_ff.matrix[:, interaction_factors[1]]
        for factor in interaction_factors
            if switch == false
                switch = true
            else
                fac =  fac .* A_ff.matrix[!, factor]
            end
        end
        A_ff.matrix[!, Symbol(join(interaction_factors, "_"))] = fac
    end
    FractionalFactorial2Level(A_ff.matrix, NamedTuple{symbol_factors}(levels), formula)
end

"""
$(TYPEDEF)

$(TYPEDFIELDS)
"""
struct RandomDesign <: AbstractRandomDesign
    matrix::DataFrame
    factors::NamedTuple
    formula::FormulaTerm
end

"""
$(TYPEDEF)

$(TYPEDFIELDS)

Encapsulates one or more `Distribution`s, generating random designs.  Receives a
`NamedTuple`  of  factor names  associated  with  `Distribution`s from  the  the
*Distributions* package.   After instantiating a `DesignDistribution`,  you must
request   samples   from   it    using   [`rand`](@ref),   which   generates   a
[`RandomDesign`](@ref)

"""
struct DesignDistribution
    factors::NamedTuple
    formula::FormulaTerm
end

"""
$(TYPEDSIGNATURES)

```jldoctest
julia> DesignDistribution((f1 = Uniform(2, 3), f2 = DiscreteUniform(-1, 5), f3 = Uniform(5, 10)))
DesignDistribution
Formula: 0 ~ f1 + f2 + f3
Factor Distributions:
f1: Distributions.Uniform{Float64}(a=2.0, b=3.0)
f2: Distributions.DiscreteUniform(a=-1, b=5)
f3: Distributions.Uniform{Float64}(a=5.0, b=10.0)

```
"""
function DesignDistribution(factors::NamedTuple)
    DesignDistribution(factors, ConstantTerm(0) ~ sum(term.(keys(factors))))
end

"""
$(TYPEDSIGNATURES)

```jldoctest
julia> DesignDistribution((Uniform(2, 3), DiscreteUniform(-1, 5), Uniform(5, 10)))
DesignDistribution
Formula: 0 ~ factor1 + factor2 + factor3
Factor Distributions:
factor1: Distributions.Uniform{Float64}(a=2.0, b=3.0)
factor2: Distributions.DiscreteUniform(a=-1, b=5)
factor3: Distributions.Uniform{Float64}(a=5.0, b=10.0)

```
"""
function DesignDistribution(factors::Tuple)
    DesignDistribution(NamedTuple{Tuple(Symbol("factor" * string(i)) for i = 1:length(factors))}(factors))
end

"""
$(TYPEDSIGNATURES)

```jldoctest
julia> DesignDistribution([Uniform(2, 3), DiscreteUniform(-1, 5), Uniform(5, 10)])
DesignDistribution
Formula: 0 ~ factor1 + factor2 + factor3
Factor Distributions:
factor1: Distributions.Uniform{Float64}(a=2.0, b=3.0)
factor2: Distributions.DiscreteUniform(a=-1, b=5)
factor3: Distributions.Uniform{Float64}(a=5.0, b=10.0)

```
"""
function DesignDistribution(factors::Array)
    DesignDistribution(tuple(factors...))
end

"""
$(TYPEDSIGNATURES)

```jldoctest
julia> DesignDistribution(DiscreteNonParametric([-1, 1], [0.5, 0.5]), 6)
DesignDistribution
Formula: 0 ~ factor1 + factor2 + factor3 + factor4 + factor5 + factor6
Factor Distributions:
factor1: Distributions.DiscreteNonParametric{Int64, Float64, Vector{Int64}, Vector{Float64}}(support=[-1, 1], p=[0.5, 0.5])
factor2: Distributions.DiscreteNonParametric{Int64, Float64, Vector{Int64}, Vector{Float64}}(support=[-1, 1], p=[0.5, 0.5])
factor3: Distributions.DiscreteNonParametric{Int64, Float64, Vector{Int64}, Vector{Float64}}(support=[-1, 1], p=[0.5, 0.5])
factor4: Distributions.DiscreteNonParametric{Int64, Float64, Vector{Int64}, Vector{Float64}}(support=[-1, 1], p=[0.5, 0.5])
factor5: Distributions.DiscreteNonParametric{Int64, Float64, Vector{Int64}, Vector{Float64}}(support=[-1, 1], p=[0.5, 0.5])
factor6: Distributions.DiscreteNonParametric{Int64, Float64, Vector{Int64}, Vector{Float64}}(support=[-1, 1], p=[0.5, 0.5])

```
"""
function DesignDistribution(distribution::D, n::Int) where D <: Distribution
    factors = tuple(fill(distribution, n)...)
    DesignDistribution(factors)
end

"""
$(TYPEDSIGNATURES)

```jldoctest
julia> rand(DesignDistribution((f1 = Uniform(2, 3), f2 = DiscreteUniform(-1, 5), f3 = Uniform(5, 10))), 12)
ExperimentalDesign.RandomDesign
Dimension: (12, 3)
Factors: (f1 = Distributions.Uniform{Float64}(a=2.0, b=3.0), f2 = Distributions.DiscreteUniform(a=-1, b=5), f3 = Distributions.Uniform{Float64}(a=5.0, b=10.0))
Formula: 0 ~ f1 + f2 + f3
Design Matrix:
12×3 DataFrame
 Row │ f1       f2       f3
     │ Float64  Float64  Float64
─────┼───────────────────────────
   1 │ 2.04922     -1.0  9.62061
   2 │ 2.59117      3.0  5.79254
   3 │ 2.77148     -1.0  6.22902
   4 │ 2.25659      4.0  6.00256
   5 │ 2.64968      2.0  9.21818
   6 │ 2.31523      5.0  8.80749
   7 │ 2.86526      2.0  8.47297
   8 │ 2.44753      0.0  6.11722
   9 │ 2.08284     -1.0  6.34626
  10 │ 2.81201     -1.0  6.44508
  11 │ 2.64232      3.0  5.57183
  12 │ 2.08189      1.0  8.72759

```
"""
function rand(distribution::DesignDistribution, n::Int = 1)
    RandomDesign(DataFrame(random_design!(values(distribution.factors), n),
                           collect(keys(distribution.factors))),
                 distribution.factors,
                 distribution.formula)
end

"""
$(TYPEDEF)

$(TYPEDFIELDS)
"""
struct RandomLHCDesign <: AbstractRandomDesign
    matrix::DataFrame
    factors::Tuple
end

"""
$(TYPEDSIGNATURES)

```jldoctest
julia> RandomLHCDesign(4,3)
RandomLHCDesign
Dimension: (4, 3)
Factors: (:factor1, :factor2, :factor3)
Design Matrix:
4×3 DataFrame
 Row │ factor1  factor2  factor3
     │ Int64    Int64    Int64
─────┼───────────────────────────
   1 │       1        4        4
   2 │       2        1        2
   3 │       3        3        3
   4 │       4        2        1

```
"""
function RandomLHCDesign(n::Int, d::Int)
    factors = Tuple(Symbol("factor" * string(i)) for i = 1:d)
    design = DataFrame(randomLHC(n,d), :auto)
    rename!(design, collect(factors))
    RandomLHCDesign(design, factors)
end

"""
$(TYPEDSIGNATURES)

```jldoctest
julia> RandomLHCDesign(4,Tuple((:A,:B,:C)))
RandomLHCDesign
Dimension: (4, 3)
Factors: (:A, :B, :C)
Design Matrix:
4×3 DataFrame
 Row │ A      B      C
     │ Int64  Int64  Int64
─────┼─────────────────────
   1 │     1      4      4
   2 │     2      1      2
   3 │     3      3      3
   4 │     4      2      1
```
"""
function RandomLHCDesign(n::Int, factors::Tuple)
    design = DataFrame(randomLHC(n,length(factors)), :auto)
    rename!(design, collect(factors))
    RandomLHCDesign(design, factors)
end

"""
$(TYPEDEF)

$(TYPEDFIELDS)
"""
struct OptimLHCDesign <: AbstractRandomDesign
    matrix::DataFrame
    factors::Tuple
    fitness::Vector{Float64}
end

"""
$(TYPEDSIGNATURES)

```jldoctest
julia> OptimLHCDesign(5,3,4)
OptimLHCDesign
Dimension: (5, 3)
Factors: (:factor1, :factor2, :factor3)
Fitness: [1.3760238272524201, 1.3760238272524201, 1.3760238272524201, 1.3760238272524201, 1.3760238272524201]
Design Matrix:
5×3 DataFrame
 Row │ factor1  factor2  factor3
     │ Int64    Int64    Int64
─────┼───────────────────────────
   1 │       3        5        2
   2 │       5        1        3
   3 │       4        4        5
   4 │       2        2        1
   5 │       1        3        4
```
"""
function OptimLHCDesign(n::Int, d::Int, gens)
    factors = Tuple(Symbol("factor" * string(i)) for i = 1:d)
    design_matrix, fitness = LHCoptim(n, d, gens)
    design = DataFrame(design_matrix, :auto)
    rename!(design, collect(factors))
    OptimLHCDesign(design, factors, fitness)
end

"""
$(TYPEDSIGNATURES)

```jldoctest
julia> OptimLHCDesign(5,Tuple((:A,:B,:C)),4)
OptimLHCDesign
Dimension: (5, 3)
Factors: (:A, :B, :C)
Fitness: [1.3760238272524201, 1.3760238272524201, 1.3760238272524201, 1.3760238272524201, 1.3760238272524201]
Design Matrix:
5×3 DataFrame
 Row │ A      B      C
     │ Int64  Int64  Int64
─────┼─────────────────────
   1 │     3      5      2
   2 │     5      1      3
   3 │     4      4      5
   4 │     2      2      1
   5 │     1      3      4

```
"""
function OptimLHCDesign(n::Int, factors::Tuple, gens)
    design_matrix, fitness = LHCoptim(n, length(factors), gens)
    design = DataFrame(design_matrix, :auto)
    rename!(design, collect(factors))
    OptimLHCDesign(design, factors, fitness)
end

"""
$(TYPEDEF)

$(TYPEDFIELDS)

Contains a  set of candidate experiments,  a target optimality criterion,  and a
model prior  formula.  A set  of experiments  that maximizes a  given optimality
criterion can selected from a candidate set using the [`kl_exchange`](@ref)
algorithm.

"""
struct OptimalDesign <: AbstractOptimalDesign
    matrix::DataFrame
    factors::NamedTuple
    formula::FormulaTerm
    selected_experiments::Array{Int}
    criteria::Dict{Symbol, Float64}
end

"""
$(TYPEDSIGNATURES)

```jldoctest
julia> design_distribution = DesignDistribution((f1 = Uniform(2, 3), f2 = DiscreteUniform(-1, 5), f3 = Uniform(5, 10)))
DesignDistribution
Formula: 0 ~ f1 + f2 + f3
Factor Distributions:
f1: Distributions.Uniform{Float64}(a=2.0, b=3.0)
f2: Distributions.DiscreteUniform(a=-1, b=5)
f3: Distributions.Uniform{Float64}(a=5.0, b=10.0)

julia> design = rand(design_distribution, 400);

julia> f = @formula 0 ~ f1 + f2 + f3 + f2 ^ 2;

julia> OptimalDesign(design, f, 10)
OptimalDesign
Dimension: (10, 3)
Factors: (f1 = Distributions.Uniform{Float64}(a=2.0, b=3.0), f2 = Distributions.DiscreteUniform(a=-1, b=5), f3 = Distributions.Uniform{Float64}(a=5.0, b=10.0))
Formula: 0 ~ f1 + f2 + f3 + :(f2 ^ 2)
Selected Candidate Rows: [244, 49, 375, 43, 369, 44, 16, 346, 175, 205]
Optimality Criteria: Dict(:D => 2.3940431912232483)
Design Matrix:
10×3 DataFrame
 Row │ f1       f2       f3
     │ Float64  Float64  Float64
─────┼───────────────────────────
   1 │ 2.99329      1.0  5.09246
   2 │ 2.96899      5.0  9.39802
   3 │ 2.96274     -1.0  5.31426
   4 │ 2.00285      2.0  5.40398
   5 │ 2.8491       1.0  9.90621
   6 │ 2.00309     -1.0  9.49394
   7 │ 2.20051      2.0  9.75605
   8 │ 2.06422      5.0  5.1759
   9 │ 2.037       -1.0  9.13114
  10 │ 2.82612      5.0  9.66349

```
"""
function OptimalDesign(candidates::AbstractDesign,
                       formula::FormulaTerm,
                       experiments::Int;
                       tolerance::Float64 = 1e-9,
                       seed_design_size::Int = 2,
                       max_iterations::Int = 1000,
                       design_k::Int = experiments,
                       candidates_l::Int = size(candidates.matrix, 1) - design_k)
    selected = kl_exchange(formula,
                           candidates.matrix,
                           experiments = experiments,
                           tolerance = tolerance,
                           seed_design_size = seed_design_size,
                           max_iterations = max_iterations,
                           design_k = design_k,
                           candidates_l = candidates_l)

    optimality = d_criterion(selected)

    OptimalDesign(candidates.matrix[selected.indices[1], :],
                  candidates.factors,
                  formula,
                  selected.indices[1],
                  Dict(:D => optimality))
end


text:
module FoldsPreview

export with_distributed_preview, with_threaded_preview, with_preview, IntervalThrottle

using Distributed
using Transducers
using Transducers:
    Executor,
    R_,
    Transducer,
    combine,
    complete,
    inner,
    next,
    start,
    unwrap,
    wrap,
    wrapping,
    xform

include("utils.jl")

include("throttles.jl")
using .Throttles: IntervalThrottle, Throttle

include("core.jl")

end # module FoldsPreview


text:
import Deno_jll
import HTTP
using Test
using Sockets

# Not all architectures have a Deno_jll
hasproperty(Deno_jll, :deno) && @testset "WebSocket server" begin

    # Set up the references needed for a clean shutdown
    server_task_ref = Ref{Any}(nothing)
    closing_ref = Ref{Bool}(false)

    server = listen(IPv4(0), 36984)
    function close_it()
        closing_ref[] = true
        close(server)
    end

    # Will contain a list of received messages
    server_received_messages = Any[]

    # Start the server async
    server_task_ref[] = @async try
        HTTP.WebSockets.listen("127.0.0.1", UInt16(36984); server=server) do ws
            while isopen(ws.io)
                data = readavailable(ws)
                msg = String(data)
                push!(server_received_messages, msg)
                if msg == "close"
                    close(ws.io)
                    close_it()
                else
                    response = "Hello, " * msg
                    write(ws, response)
                end
            end
        end
    catch e
        if closing_ref[] && (e isa Base.IOError)
            # this is "expected"
        else
            @error "WebSocket server error" exception=(e,catch_backtrace())
            rethrow(e)
        end
    end

    success = try
        # Run our client tests using Deno
        # this throws error if the Deno tests fail
        run(`$(Deno_jll.deno()) test --allow-net`)
        true
    catch e
        if e isa ProcessFailedException
            false
        else
            rethrow(e)
        end
    finally
        close_it()
        wait(server_task_ref[])
    end
    @test success

    @test server_received_messages == ["world", "close"]
end


text:
@doc raw"""
`TableauVPRK`: Tableau of a Variational Partitioned Runge-Kutta method
```math
\begin{align*}
P_{n,i} &= \dfrac{\partial L}{\partial v} (Q_{n,i}, V_{n,i}) , &
Q_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} a_{ij} \, V_{n,j} , &
q_{n+1} &= q_{n} + h \sum \limits_{i=1}^{s} b_{i} \, V_{n,i} , \\
F_{n,i} &= \dfrac{\partial L}{\partial q} (Q_{n,i}, V_{n,i}) , &
P_{n,i} &= p_{n} + h \sum \limits_{i=1}^{s} \bar{a}_{ij} \, F_{n,j} - d_i \lambda , &
p_{n+1} &= p_{n} + h \sum \limits_{i=1}^{s} \bar{b}_{i} \, F_{n,i} , \\
&&
0 &= \sum \limits_{i=1}^{s} d_i V_i , &&
\end{align*}
```
satisfying the symplecticity conditions
```math
\begin{align*}
b_{i} \bar{a}_{ij} + b_{j} a_{ji} &= b_{i} b_{j} , &
\bar{b}_i &= b_i .
\end{align*}
```
"""
struct TableauVPRK{T} <: AbstractTableauPRK{T}
    @HeaderTableau

    q::CoefficientsRK{T}
    p::CoefficientsRK{T}

    R∞::Int

    d::Vector{T}

    function TableauVPRK{T}(name, o, q, p, R∞, d) where {T}
        @assert q.s == p.s == length(d)
        new(name, o, q.s, q, p, R∞, d)
    end

    function TableauVPRK{T}(name, o, q, p, R∞) where {T}
        @assert q.s == p.s
        new(name, o, q.s, q, p, R∞)
    end
end

function TableauVPRK(name::Symbol, order::Int, q::CoefficientsRK{T}, p::CoefficientsRK{T}, R∞::Int, d::Vector{T}) where {T}
    TableauVPRK{T}(name, order, q, p, R∞, d)
end

function TableauVPRK(name::Symbol, order::Int, q::CoefficientsRK{T}, p::CoefficientsRK{T}, R∞::Int) where {T}
    TableauVPRK{T}(name, order, q, p, R∞)
end

function TableauVPRK(name::Symbol, order::Int, q::CoefficientsRK{T}, R∞::Int, d::Vector{T}) where {T}
    TableauVPRK{T}(name, order, q, get_symplectic_conjugate_coefficients(q), R∞, d)
end

function TableauVPRK(name::Symbol, order::Int, q::CoefficientsRK{T}, R∞::Int) where {T}
    TableauVPRK{T}(name, order, q, get_symplectic_conjugate_coefficients(q), R∞)
end

"Parameters for right-hand side function of variational partitioned Runge-Kutta methods."
mutable struct ParametersVPRK{DT, TT, ET <: IODE{DT,TT}, D, S} <: AbstractParametersVPRK{DT,TT,ET,D,S}
    equ::ET
    tab::TableauVPRK{TT}
    Δt::TT

    t̅::TT
    q̅::Vector{DT}
    p̅::Vector{DT}
end

function ParametersVPRK(equ::ET, tab::TableauVPRK{TT}, Δt::TT) where {DT, TT, ET <: IODE{DT,TT}}
    q = zeros(DT, equ.d)
    p = zeros(DT, equ.d)

    ParametersVPRK{DT, TT, ET, equ.d, tab.s}(equ, tab, Δt, 0, q, p)
end

"Compute stages of variational partitioned Runge-Kutta methods."
@generated function function_stages!(x::Vector{ST}, b::Vector{ST},
                params::ParametersVPRK{DT,TT,ET,D,S}
            ) where {ST,DT,TT,ET,D,S}

    cache = NonlinearFunctionCacheVPRK{ST}(D, S)

    quote
        @assert length(x) == length(b)

        compute_stages!(x, $cache.Q, $cache.V, $cache.P, $cache.F, params)
        compute_rhs_vprk!(b, $cache.P, $cache.F, params)
        compute_rhs_vprk_correction!(b, $cache.V, params)
    end
end

"Variational partitioned Runge-Kutta integrator."
struct IntegratorVPRK{DT, TT, PT <: ParametersVPRK{DT,TT},
                              ST <: NonlinearSolver{DT},
                              IT <: InitialGuessPODE{DT,TT}} <: AbstractIntegratorVPRK{DT,TT}
    params::PT
    solver::ST
    iguess::IT
end

function IntegratorVPRK(equation::ET, tableau::TableauVPRK{TT}, Δt::TT) where {DT, TT, ET <: IODE{DT,TT}}
    D = equation.d
    S = tableau.s

    # create params
    params = ParametersVPRK(equation, tableau, Δt)

    # create nonlinear solver
    solver = create_nonlinear_solver(DT, D*S, params)

    # create initial guess
    iguess = InitialGuessPODE(get_config(:ig_interpolation), equation, Δt)

    # create integrator
    IntegratorVPRK{DT, TT, typeof(params), typeof(solver), typeof(iguess)}(params, solver, iguess)
end

IntegratorVPRKpNone = IntegratorVPRK

equation(int::IntegratorVPRK) = int.params.equ
timestep(int::IntegratorVPRK) = int.params.Δt
tableau(integrator::IntegratorVPRK) = integrator.params.tab
nstages(integrator::IntegratorVPRK) = integrator.params.tab.s

function initial_guess!(int::IntegratorVPRK, cache::IntegratorCacheVPRK)
    for i in 1:nstages(int)
        evaluate!(int.iguess, cache.q, cache.p, cache.v, cache.f,
                              cache.q̅, cache.p̅, cache.v̅, cache.f̅,
                              cache.q̃, cache.ṽ,
                              tableau(int).q.c[i])

        for k in 1:ndims(int)
            int.solver.x[ndims(int)*(i-1)+k] = cache.ṽ[k]
        end
    end
end

"Integrate ODE with variational partitioned Runge-Kutta integrator."
function integrate_step!(int::IntegratorVPRK{DT,TT}, cache::IntegratorCacheVPRK{DT,TT}) where {DT,TT}
    # update nonlinear solver parameters from cache
    update_params!(int.params, cache)

    # compute initial guess
    initial_guess!(int, cache)

    # reset cache
    reset!(cache, timestep(int))

    # call nonlinear solver
    solve!(int.solver)

    # print solver status
    print_solver_status(int.solver.status, int.solver.params, cache.n)

    # check if solution contains NaNs or error bounds are violated
    check_solver_status(int.solver.status, int.solver.params, cache.n)

    # compute vector fields at internal stages
    compute_stages!(int.solver.x, cache.Q, cache.V, cache.P, cache.F, int.params)

    # compute final update
    update_solution!(int, cache)

    # copy solution to initial guess
    update!(int.iguess, cache.t, cache.q, cache.p, cache.v, cache.f)

    # take care of periodic solutions
    cut_periodic_solution!(cache, equation(int).periodicity)
end


text:
function stripcomments(input, markers)
    function popspace!(arr)
        while !isempty(arr) && isspace(arr[end])
            pop!(arr)
        end
    end
    markers = join(markers)
    output = Char[]
    incomment = false
    for c ∈ input
        if c == '\n'
            incomment = false
            popspace!(output)
            append!(output,c)
        elseif incomment
        elseif c ∈ markers
            incomment = true
        else
            append!(output,c)
        end
    end
    popspace!(output)
    String(output)
end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
#set omnifunc=syntaxcomplete#Complete
#let g:ycm_collect_identifiers_from_tags_files = 1
#set tags+=base.tags

module MyModule
using Lib

intersection

using BigLib: thing1, thing2

import Base.show

importall OtherLib

export MyType, foo

type MyType
    x
end

foobuzz = 5
tt=1
export foobuzz, if

macro sayhello()
	return :( println("Hello, world!") )
end

export dead,
	   theparot

function dead(x)
	2*x
end

function theparot(x)
	dead(x)
end

a,b = 4,3

bar(x) = 2x
foo(a::MyType) = bar(a.x) + 1

show(io, a::MyType) = print(io, "MyType $(a.x)")
end

max_stars_repo_name: oxinabox/julia-vim-omnicomplete

max_issues_repo_name: oxinabox/julia-vim-omnicomplete

max_forks_repo_name: oxinabox/julia-vim-omnicomplete


text:
module LOLToolsDataDragon

export DataDragon
include("DataDragon.jl")

end # module LOLToolsDataDragon


text:
######################################################################
# plotting.jl
# Handles all the calls to PyPlot to generate plots
#
# hold is deprecated as of Matplotlib 2.0
#  Figures now default to hold (true)
#  I think we should adhere to this for future compliance
######################################################################

"""
`plot(m::SearchDomain, f::AbstractFilter, x::Vehicle)`

Plots the belief, jammer, and vehicles.
"""
function plot(m::SearchDomain, f::AbstractFilter, x::Vehicle; show_mean::Bool=false, show_cov::Bool=false, alpha=0.2, color="k")
	plot(m, f, (x.x, x.y, x.heading); show_mean=show_mean, show_cov=show_cov, alpha=alpha, color=color)
end
# allow plotting of multiple vehicles with a single shared filter
function plot(m::SearchDomain, f::AbstractFilter, vx::Vector{Vehicle})
	plot(m, f, vx[1])
	num_vehicles = length(vx)
	for i = 2:num_vehicles
		plot_vehicle(m, vx[i])
	end
end
function plot(m::SearchDomain, f::AbstractFilter, p::Pose; show_mean::Bool=false, show_cov::Bool=false, alpha=1.0, color="b")
	plot_theta(m)
	#hold(true) # deprecated
	plot_vehicle(m, p; color=color)
	if show_mean
		plot_mean(centroid(f), color=color)
	end
	if show_cov
		plot(m, centroid(f), covariance(f), color=color)
	end
	plot(m, f, alpha=alpha)
	return # so it doesn't spit out result of axis
end

function plot(m::SearchDomain, b::Matrix{Float64}, x::Vehicle; show_mean::Bool=false, show_cov::Bool=false, alpha=1.0, color="b", obs=nothing)
	plot(m, b, (x.x, x.y, x.heading); show_mean=show_mean, show_cov=show_cov, alpha=alpha, color=color, obs=obs)
end

function plot(m::SearchDomain, b::Matrix{Float64}, p::Pose; show_mean::Bool=false, show_cov::Bool=false, alpha=1.0, color="b", obs=nothing)
	plot_theta(m)
	#hold(true)  # deprecated
	plot_vehicle(m, p; color=color)
	#title_string = "plot: "
	title_string = ""
	if obs != nothing
		# TODO: just temporary for rss submission
		#title_string = "$(title_string)o = $(obs), "
		title_string = "t = $(obs) s"
	end
	if show_mean
		plot_mean(centroid(b, m.length), color=color)
	end
	if show_cov
		cov = covariance(b,m.length)
		evals = eigvals(cov)
		e1 = round(evals[1],1)
		e2 = round(evals[2],1)
		plot(m, centroid(b,m.length), cov, color=color)
		title_string = "$(title_string)e1 = $e1, e2 = $e2"
		#title("e1 = $e1, e2 = $e2")
	end
	title(title_string, fontname="Times New Roman", fontsize=12)
	plot(m, b, alpha=alpha)
	return # so it doesn't spit out result of axis
end

function plot(m::SearchDomain, f::AbstractFilter; alpha=1.0)
	error(typeof(f), " has not implemented `plot`(::SearchDomain, ::AbstractFilter)")
end

function plot(m::SearchDomain, f::DF; alpha=1.0, cmap="Greys")
	plot(m, f.b, alpha=alpha, cmap=cmap)
end
function plot(m::SearchDomain, b::Matrix{Float64}; alpha=1.0, cmap="Greys")

	# alpha changed by LD for rss submission 2017
	#alpha = 0.5
	#alpha = 0.6
    # commented out for cdc 2018

	a = [0,m.length,0,m.length]
	imshow(b', interpolation="none",cmap=cmap,origin="lower",extent=a,vmin=0, alpha=alpha)
	#imshow(b', cmap=cmap,origin="lower",extent=a,vmin=0, alpha=alpha)
	labels()
	axis(a)
	tick_params(direction="in")
end

plot(m::SearchDomain, f::EKF; alpha=0.1) = plot(m, f.mu, f.Sigma)
plot(m::SearchDomain, f::UKF; alpha=0.1) = plot(m, f.mu, f.Sigma)

# Plots mean and 95% confidence ellipse
function plot(m::SearchDomain, mu::LocTuple, Sigma::Matrix{Float64}; color="b")
	plot(m, [mu[1],mu[2]], Sigma, color=color)
end
function plot(m::SearchDomain, mu::Vector{Float64}, Sigma::Matrix{Float64}; color="b")
	a = [0,m.length,0,m.length]

	Sigma_half = sqrtm(Sigma)
	Sh11 = Sigma_half[1,1]
	Sh12 = Sigma_half[1,2]
	Sh21 = Sigma_half[2,1]
	Sh22 = Sigma_half[2,2]
	m1 = mu[1]
	m2 = mu[2]

	thetas = 0:.02:2pi
	num_theta = length(thetas)
	xvals = zeros(num_theta)
	yvals = zeros(num_theta)
	c = sqrt(-2.0*log(0.05))
	for i = 1:num_theta
		theta = thetas[i]
		s_theta = sin(theta)
		c_theta = cos(theta)
		xvals[i] = c*(Sh11*s_theta + Sh12*c_theta) + m1
		yvals[i] = c*(Sh21*s_theta + Sh22*c_theta) + m2
	end
	plot(xvals, yvals, color)
	plot(m1, m2, "$(color)x", ms=10)

	labels()
	#axis("square")
	axis(a)
	tick_params(direction="in")
end

#function plot_ellipse()
#end

function plot(m::SearchDomain, f::PF; alpha=0.2, cmap="Greys", color="k")
	mark_size = 12
	a = [0,m.length,0,m.length]
	x = zeros(f.n)
	y = zeros(f.n)
	for i = 1:f.n
		x[i] = particle(f,i)[1]
		y[i] = particle(f,i)[2]
	end
	#scatter(x,y,c=f.W,cmap="Greys",vmin=0)
	#scatter(x,y,c=f.W,cmap="Blues",vmin=0, alpha=0.2)
	scatter(x, y, s=1, alpha=alpha, c=color)
	#xmean, ymean = centroid(f)
	#plot(xmean, ymean, "gx", ms=10, mew=2)
	labels()
	axis("scaled")
	axis(a)
	tick_params(direction="in")
end
function plot(m::SearchDomain, b::Vector{Tuple{Float64,Float64}}; alpha=0.2, cmap="Greys", color="k")
	plot_theta(m)
	mark_size = 12
	a = [0,m.length,0,m.length]
    x = zeros(length(b))
    y = zeros(length(b))
    for i = 1:length(b)
        x[i] = b[i][1]
        y[i] = b[i][2]
	end
	#scatter(x,y,c=f.W,cmap="Greys",vmin=0)
	#scatter(x,y,c=f.W,cmap="Blues",vmin=0, alpha=0.2)
	scatter(x, y, s=1, alpha=alpha, c=color)
	#xmean, ymean = centroid(f)
	#plot(xmean, ymean, "gx", ms=10, mew=2)
	labels()
	axis("scaled")
	axis(a)
	tick_params(direction="in")
end
function plot(m::SearchDomain, b::Vector{Tuple{Float64,Float64}}, p; alpha=0.2, cmap="Greys", color="k")
	plot_theta(m)
	plot_vehicle(m, p; color=color)
	mark_size = 12
	a = [0,m.length,0,m.length]
    x = zeros(length(b))
    y = zeros(length(b))
    for i = 1:length(b)
        x[i] = b[i][1]
        y[i] = b[i][2]
	end
	#scatter(x,y,c=f.W,cmap="Greys",vmin=0)
	#scatter(x,y,c=f.W,cmap="Blues",vmin=0, alpha=0.2)
	scatter(x, y, s=1, alpha=alpha, c=color)
	#xmean, ymean = centroid(f)
	#plot(xmean, ymean, "gx", ms=10, mew=2)
	labels()
	axis("scaled")
	axis(a)
	tick_params(direction="in")
end
function plot(m::SearchDomain, b::ParticleCollection; alpha=1.0)
	mark_size = 12
	a = [0,m.length,0,m.length]
    n = length(b.particles)
	x = zeros(n)
	y = zeros(n)
	for i = 1:n
		x[i] = particle(b,i)[1]
		y[i] = particle(b,i)[2]
	end
	#scatter(x,y,c=f.W,cmap="Greys",vmin=0)
	#scatter(x,y,c=f.W,cmap="Blues",vmin=0, alpha=0.2)
	scatter(x, y, s=1, alpha=0.2, c="k")
	#xmean, ymean = centroid(f)
	#plot(xmean, ymean, "gx", ms=10, mew=2)
	labels()
	axis("scaled")
	axis(a)
	tick_params(direction="in")
end

function plot_mean(mu::LocTuple; color="b")
	plot_mean(mu[1], mu[2], color=color)
end
function plot_mean(xmean::Float64, ymean::Float64; color="b")
	plot(xmean, ymean, "$(color)X", ms=10, mew=1.5, mec="k")
end

# This is needed for plot_sim
#function plot_b(m::SearchDomain, b::Belief, x::Vehicle)
#	plot_theta(m)
#	hold(true)
#	plot_vehicles(X)
#	imshow(b', interpolation="none", cmap="Greys", origin="lower")
#	labels()
#end

#"""
#`plot_eid(m::SearchDomain, X::VehicleSet)`
#
#Plots the eid using the Fisher information matrix.
#"""
#function plot_eid(m::SearchDomain, X::VehicleSet; contours=false)
#	eid = EID(m,m.b)
#	if contours
#		plot_contour(m, eid)
#	end
#	plot_theta(m)
#	plot_vehicles(X)
#	imshow(eid', interpolation="none", cmap="Greys", origin="lower")
#	labels()
#end

#"""
#`plot_mi(m::SearchDomain, X::VehicleSet)`
#
#Plots the mutual information.
#"""
#function plot_mi(m::SearchDomain, X::VehicleSet; contours=false)
#	mut_info = mutual_information(m)
#	plot_vehicles(X)
#	imshow(mut_info', interpolation="none", cmap="Greys", origin="lower")
#	labels()
#	if contours
#		plot_contour(m, mut_info)
#	end
#end

#"""
#`plot_sim(m::SearchDomain, s::Simulation)`
#
#Steps through a simulation.
#"""
#function plot_sim(m::SearchDomain, s::Simulation)
#	for t = 0:s.T
#		hold(false)
#		pause(1)
#		plot_b(m, s.belief_list[t+1], s.state_list[t+1])
#	end
#end

######################################################################
# Helper functions
# mec = markeredgecolor
# mfc = markerfacecolor
# lw = linewidth
######################################################################
# Plots locations of the vehicles
plot_vehicle(m::SearchDomain,x::Vehicle; color="b") = plot_vehicle(m,x.x,x.y,x.heading, color=color)
plot_vehicle(m::SearchDomain, p::Pose; color="b") = plot_vehicle(m, p[1], p[2], p[3], color=color)
function plot_vehicle(m::SearchDomain, x::Float64, y::Float64, h::Float64; color="k")
	mark_size = 10

	# Changed by LD for 01/30/2017 for rss submission
	#c = 0.01 * m.length  *sqrt(2)
	#c = 0.015 * m.length  *sqrt(2)
	#c = 0.012 * m.length  *sqrt(2)
	c = 0.013 * m.length  *sqrt(2)
	dx = c*sind(45 + h)
	dy = c*cosd(45 + h)

	# Plot rotors
	# LD change 01/30/2017 for rss submission
	#rotor_size = 5
	#rotor_size = 7.5
	rotor_size = 5.5
	rotor_size = 6.5
	theta1 = 45.0 + h
	theta2 = 135.0 + h
	theta3 = 225.0 +  h
	theta4 = 315.0 + h
	#cs = "$(color)o"
	cs = "wo"
	#color = color
	plot(x+c*sind(theta1), y+c*cosd(theta1), cs, ms=rotor_size, mew=2.5, mfc="none", mec=color)
	plot(x+c*sind(theta2), y+c*cosd(theta2), cs, ms=rotor_size, mew=2.5, mfc="none", mec=color)
	plot(x+c*sind(theta3), y+c*cosd(theta3), cs, ms=rotor_size, mew=2.5, mfc="none", mec=color)
	plot(x+c*sind(theta4), y+c*cosd(theta4), cs, ms=rotor_size, mew=2.5, mfc="none", mec=color)

	# plot heading direction
	xline = [x, x+2*c*sind(h)]
	yline = [y, y+2*c*cosd(h)]
	# Changed by LD 01/30/2017 for rss submission
	#plot(xline, yline, color, mew=2)
	#plot(xline, yline, color, lw=2)
	#plot(xline, yline, color, lw=2.5)
	plot(xline, yline, color, lw=2.)

	# Plot frame
	plot(x, y, marker=(2,0,-theta1), ms=mark_size, mew=1, mec=color)
	plot(x, y, marker=(2,0,-theta2), ms=mark_size, mew=1, mec=color)

end

# Plots jammer location
function plot_theta(m::SearchDomain)
	# changed by LD 01/30/2017 for RSS submission
	#mark_size = 11
	#mark_size = 13
	mark_size = 11
	#mark_size = 13
	plot(m.theta[1], m.theta[2], "k^", markersize=mark_size, mew=2.5,markerfacecolor="none")
	#mark_size = 11
	#plot(m.theta[1], m.theta[2], "r^", markersize=mark_size, markerfacecolor="none", markeredgecolor="r", mew=2)
end

# Plots contours of some distribution `d` (a matrix).
function plot_contour(m::SearchDomain, d::Matrix{Float64})
	X,Y = meshgrid(linspace(0.,m.length,100), linspace(0.,m.length,100))
	contour(X, Y, d')
end

# Sets the appropriate plot labels
function labels()
	#xlabel("East (m)", fontsize=12)
	#ylabel("North (m)", fontsize=12)
	xlabel("East (m)")
	ylabel("North (m)")
end

######################################################################
# Old, deprecated, or unused
######################################################################
function plot_eid2(m, b, x, theta)
	mark_size = 12
	eid = EID(m,theta[1],theta[2])
	plot(x[1], x[2], "b*", markersize=mark_size)
	plot(theta[1], theta[2], "r^", markersize=mark_size)
	imshow(eid', interpolation="none", cmap="Greys", origin="lower")
	xlabel("x")
	ylabel("y")
end

function plot_eid3(m, b, x, theta)
	mark_size = 12
	eid = EID(m,theta[1],theta[2])
	#plot(x[1], m.num_cells - 1 - x[2], "b*", markersize=mark_size)
	#plot(theta[1], m.num_cells - 1 - theta[2], "r^", markersize=mark_size)
	#eid_plot = (eid')[end:-1:1, :]

	#imshow(eid_plot, interpolation="none", cmap="Greys")
	plot_contour(m, eid)
	xlabel("x")
	ylabel("y")
end

"""
`meshgrid(x, y)`

Returns X and Y.
"""
function meshgrid(x, y)
	nx = length(x)
	ny = length(y)

	X = zeros(ny, nx)
	Y = zeros(ny, nx)
	for xi = 1:nx
		for xj = 1:ny
			X[xj, xi] = x[xi]
			Y[xj, xi] = y[xj]
		end
	end
	return X, Y
end

######################################################################
# plotting multiple plots...
######################################################################
# TODO: make this just call the other version
function plot{TP<:AbstractFilter}(m::SearchDomain, farr::Vector{TP}, xarr::Vector{Vehicle}; show_mean::Bool=false, show_cov::Bool=false)
	for i = 1:length(farr)
		f = farr[i]
		x = xarr[i]
		plot(m, f, (x.x, x.y, x.heading); show_mean=show_mean, show_cov=show_cov)
	end
end

function plot{TP<:AbstractFilter}(m::SearchDomain, farr::Vector{TP}, parr::Vector{Pose}; show_mean::Bool=false, show_cov::Bool=false, alpha=1.0)
	num_vehicles = length(farr)
	cmaps = ["Blues", "Reds"]
	colors = ["b", "r"]
	if num_vehicles == 1
		cmaps = ["Greys"]
		colors = ["b"]
	elseif num_vehicles == 2
		cmaps = ["Blues", "Reds"]
		colors = ["b", "r"]
	elseif num_vehicles == 3
		cmaps = ["Blues", "Reds", "Greens"]
		colors = ["b", "r", "g"]
	elseif num_vehicles == 4
		cmaps = ["Blues", "Reds", "Greens", "Greys"]
		colors = ["b", "r", "g", "k"]
	end
	plot_theta(m)
	#hold(true)  # deprecated
	#cmaps = ["Blues", "Reds"]
	#colors = ["b", "r"]
	for (fi,f) in enumerate(farr)
		p = parr[fi]
		plot_vehicle(m,p, color=colors[fi])
		if show_mean
			plot_mean(centroid(f), color=colors[fi])
		end
		if show_cov
			plot(m, centroid(f), covariance(f), color=colors[fi])
		end
		plot(m, f, alpha=alpha, cmap=cmaps[fi])
	end
	return # so it doesn't spit out result of axis
end

function plot(m::SearchDomain, vsu::Vector{SimUnit};
              show_mean::Bool=false,
              show_cov::Bool=false,
              alpha=1.0
             )
	num_vehicles = length(vsu)
	cmaps = ["Blues", "Reds"]
	colors = ["b", "r"]
	if num_vehicles == 1
		cmaps = ["Greys"]
		colors = ["b"]
	elseif num_vehicles == 2
		cmaps = ["Blues", "Reds"]
		colors = ["b", "r"]
	elseif num_vehicles == 3
		cmaps = ["Blues", "Reds", "Greens"]
		colors = ["b", "r", "g"]
	elseif num_vehicles == 4
		cmaps = ["Blues", "Reds", "Greens", "Greys"]
		colors = ["b", "r", "g", "k"]
	end
	plot_theta(m)
	#hold(true)  # deprecated
	#cmaps = ["Blues", "Reds"]
	#colors = ["b", "r"]
	for (fi,uav) in enumerate(vsu)
		#p = parr[fi]
		plot_vehicle(m, uav.x, color=colors[fi])
		if show_mean
			plot_mean(centroid(uav.f), color=colors[fi])
		end
		if show_cov
			plot(m, centroid(uav.f), covariance(uav.f), color=colors[fi])
		end
        plot(m, uav.f, alpha=alpha, cmap=cmaps[fi], color=colors[fi])
	end
	return # so it doesn't spit out result of axis
end


text:
function grid{T}(width::T, d::Int)
    x = Array{T}(d, d, d)
    y = Array{T}(d, d, d)
    z = Array{T}(d, d, d)

    g = collect(linspace(-width/2, width/2, d))

    @inbounds for ix in 1:d, iy in 1:d, iz in 1:d
        x[ix, iy, iz] = g[ix]
        y[ix, iy, iz] = g[iy]
        z[ix, iy, iz] = g[iz]
    end

    return x, y, z
end

function grid_coords{T}(width::T, d::Int)
    x, y, z = grid(width, d)
    hcat(vec(x), vec(y), vec(z))'
end


text:
using Documenter, DICE

makedocs()

deploydocs(
    deps   = Deps.pip("mkdocs", "python-markdown-math"),
    repo   = "github.com/Libbum/DICE.jl.git",
    julia  = "0.6"
)


text:
import BalancePlot

if length(ARGS) < 2
    println("Usage: balanceplot <config file> <transaction file>")
    exit(1)
end

config = BalancePlot.load_ing_config_from_file(ARGS[1])
trs = BalancePlot.load_ing_transactions_from_file(ARGS[2])
trs2 = BalancePlot.convert_ing_transactions(trs, config)

BalancePlot.plot_transactions(trs2)

max_stars_repo_name: CrazedProgrammer/BalancePlot

max_issues_repo_name: CrazedProgrammer/BalancePlot

max_forks_repo_name: CrazedProgrammer/BalancePlot


text:
using DelimitedFiles

words1 = readdlm("./examples/forceAlingned.jl/words.txt")

words2 = readdlm("./examples/forceAlingned.jl/aenesa.txt")

words = [words1[:,1] words2[:,3] words2[:,4]]

# set item info
class = "IntervalTier"
name = "Words"
size1 = size(words)[1]
xmin = words[1,:][2]
xmax = words[end,:][3]

# print(" "^4*"item [10]:\n")
# print(" "^8*"name = \""*name*"\" \n")
# print(" "^8*"xmin = "*string(xmin)*" \n")
# print(" "^8*"xmax = "*string(xmax)*" \n")
# print(" "^8*"intervals: size = "*string(size1)*" \n")
# for i in 1:size1
#     print(" "^8*"intervals ["*string(i)*"]:\n")
#     print(" "^12*"xmin = "*string(words[i,:][2])*" \n")
#     print(" "^12*"xmax = "*string(words[i,:][3])*" \n")
#     print(" "^12*"text = "*string(words[i,:][1])*" \n")
# end

open("./examples/forceAlingned.jl/aenesa_results.TextGrid", "a") do file
    write(file, " "^4*"item [10]:\n")
    write(file, " "^8*"class = \""*class*"\" \n")
    write(file, " "^8*"name = \""*name*"\" \n")
    write(file, " "^8*"xmin = "*string(xmin)*" \n")
    write(file, " "^8*"xmax = "*string(xmax)*" \n")
    write(file, " "^8*"intervals: size = "*string(size1)*" \n")
    for i in 1:size1
        write(file, " "^8*"intervals ["*string(i)*"]:\n")
        write(file, " "^12*"xmin = "*string(words[i,:][2])*" \n")
        write(file, " "^12*"xmax = "*string(words[i,:][3])*" \n")
        write(file, " "^12*"text = \""*string(words[i,:][1])*"\" \n")
    end
end


text:
# Univariate Gaussian Mixture Models

struct UnivariateGMM <: UnivariateMixture{Continuous,Normal}
    K::Int
    means::Vector{Float64}
    stds::Vector{Float64}
    prior::Categorical

    function UnivariateGMM(ms::Vector{Float64}, ss::Vector{Float64}, pri::Categorical)
        K = length(ms)
        length(ss) == K || throw(DimensionMismatch())
        ncategories(pri) == K ||
            error("The number of categories in pri should be equal to the number of components.")
        new(K, ms, ss, pri)
    end
end

@distr_support UnivariateGMM -Inf Inf

ncomponents(d::UnivariateGMM) = d.K

component(d::UnivariateGMM, k::Int) = Normal(d.means[k], d.stds[k])

probs(d::UnivariateGMM) = probs(d.prior)

mean(d::UnivariateGMM) = dot(d.means, probs(d))

rand(d::UnivariateGMM) = (k = rand(d.prior); d.means[k] + randn() * d.stds[k])

params(d::UnivariateGMM) = (d.means, d.stds, d.prior)

struct UnivariateGMMSampler <: Sampleable{Univariate,Continuous}
    means::Vector{Float64}
    stds::Vector{Float64}
    psampler::AliasTable
end

rand(s::UnivariateGMMSampler) = (k = rand(s.psampler); s.means[k] + randn() * s.stds[k])
sampler(d::UnivariateGMM) = UnivariateGMMSampler(d.means, d.stds, sampler(d.prior))


text:
#module Scenario
#using  ..Types

#import generic functions
#import .Ballistics: InternalBallistics!, ExternalBallistics!, Vulnerability!

#export createProjectile, createFragment, createTarget

function createTarget(size::Float64, position::Float64)
    Target1D(size,position)
end

"""
    createTargetRect(a,b, position)
Return a rectangular target by specifying the length 'a', the height 'b' and the 'position' (range)
"""
function createTargetRect(a::Float64,b::Float64, position::Float64)
    TargetRect(a,b,position)
end

"""
    createTargetCirc(size, position)
Returns a circular target by specifying the radius (size) and the position (range)
"""
function createTargetCirc(size::Float64, position::Float64)
    TargetCirc(size,position)
end

"""
    createTarget()
Retursns an object of type 'Target'.
Optional arguments are:
* ρ the size of the target
* position, the range to the target
"""
function createTarget(;ρ=nothing,position=nothing)
    Target(position,ρ)
end

"""
    createGun(u₀,lat,QE,AZ,tc; args ...)
Returns an object of type 'Gun' by specifying the muzzle velocity (u₀) its latitude (lat), the elevation angle (QE),
azimuth angle (AZ) and the twist (tc).
Optional arguments are
* lw, the gun length
* X2w, the height of the gun
"""
function createGun(u₀::Float64,lat::Float64,QE::Float64,AZ::Float64,tc::Float64; lw=nothing, X2w=nothing)
    Gun(u₀,lat,tc,lw,X2w,QE,AZ)
end

#end


text:
module DisjunctiveConstraints

import MathOptInterface
const MOI = MathOptInterface
const MOIU = MOI.Utilities

const VI = MOI.VariableIndex
const CI = MOI.ConstraintIndex
const SV = MOI.SingleVariable
const SAT = MOI.ScalarAffineTerm{Float64}
const SAF = MOI.ScalarAffineFunction{Float64}
const ET = MOI.EqualTo{Float64}
const LT = MOI.LessThan{Float64}
const GT = MOI.GreaterThan{Float64}
const ZO = MOI.ZeroOne

export DisjunctiveSet, CombinatorialDisjunctiveSet
export minimum_activity, maximum_activity

include("sets.jl")
include("activity.jl")
include("big_m.jl")

end # module


text:
import MathOptInterface
import MathOptInterface: Utilities

const MOI = MathOptInterface
const MOIU = MOI.Utilities

# indices
const VI = MOI.VariableIndex
const CI = MOI.ConstraintIndex

# function aliases
const SV = MOI.SingleVariable
const SAF = MOI.ScalarAffineFunction{Float64}
const SQF = MOI.ScalarQuadraticFunction{Float64}

# set aliases
const Bounds{T} = Union{
    MOI.EqualTo{T},
    MOI.GreaterThan{T},
    MOI.LessThan{T},
    MOI.Interval{T}
}

mutable struct Optimizer <: MOI.AbstractOptimizer
    inner::BaronModel
    nlp_block_data::Union{Nothing, MOI.NLPBlockData}
end

Optimizer(; options...) = Optimizer(BaronModel(; options...), nothing)

MOI.get(model::Optimizer, ::MOI.SolverName) = "BARON"

# empty
function MOI.is_empty(model::Optimizer)
    BARON.is_empty(model.inner) && model.nlp_block_data === nothing
end

function MOI.empty!(model::Optimizer)
    model.inner = BaronModel(; ((Symbol(key), val) for (key, val) in model.inner.options)...)
    model.nlp_block_data = nothing
    return
end

# copy
MOIU.supports_default_copy_to(model::Optimizer, copy_names::Bool) = !copy_names
function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike; kws...)
    return MOIU.automatic_copy_to(dest, src; kws...)
end

# optimize
function MOI.optimize!(model::Optimizer)
    if !IS_SOLVER_SET
        error(("""BARON.jl was not built correctly.
                 Set the environment variable `BARON_EXEC` and run `using Pkg; Pkg.build("BARON")`."""))
    end
    write_bar_file(model.inner)
    run(`$baron_exec $(model.inner.problem_file_name)`)
    read_results(model.inner)
end

# RawParameter
MOI.supports(::Optimizer, ::MOI.RawParameter) = true
function MOI.set(model::Optimizer, param::MOI.RawParameter, value)
    model.inner.options[param.name] = value
    return
end
function MOI.get(model::Optimizer, param::MOI.RawParameter)
    return get(model.inner.options, param.name) do
        throw(ErrorException("Requested parameter $(param.name) is not set."))
    end
end

# TimeLimitSec
MOI.supports(::Optimizer, ::MOI.TimeLimitSec) = true
function MOI.set(model::Optimizer, ::MOI.TimeLimitSec, val::Real)
    model.inner.options["MaxTime"] = Float64(val)
    return
end

# BARON's default time limit is 1000 seconds.
function MOI.get(model::Optimizer, ::MOI.TimeLimitSec)
    return get(model.inner.options, "MaxTime", 1000.0)
end

include(joinpath("moi", "util.jl"))
include(joinpath("moi", "variables.jl"))
include(joinpath("moi", "constraints.jl"))
include(joinpath("moi", "objective.jl"))
include(joinpath("moi", "nlp.jl"))
include(joinpath("moi", "results.jl"))


text:
# using CSV

#experiment number in filenameMTiXsens
print_bool=false

MTiXsens=Dict()
Xsensor=Dict()
IBS = Dict()

if size(filenameMTiXsens)!=size(filenameXsensor)
    @warn "MTiXsens and Xsensor do not have the same number of files"
end

if size(filenameMTiXsens)!=size(filenameIBS)
    @warn "MTiXsens and IBS do not have the same number of files"

end

for experiment=1:size(filenameMTiXsens,1)
    #Read CSV of Xsens
    MTiXsens[filenameMTiXsens[experiment,1]]=Dict()
    for n=1:5
        MTiXsens[filenameMTiXsens[experiment,1]][sensorXsensNames[n]] = CSV.read(string(MTiXsensFolderPath,filenameMTiXsens[experiment,1], filenameMTiXsensID[n],".txt");header=5)
    end

    if print_bool
        println(string("Xsens CSV read finished, experiment : "),experiment)
    end
end

for experiment=1:size(filenameXsensor,1)
    #Read CSV of Xsensor
    #Xsensor[filenameXsensor[experiment,1]]=Dict()
    Xsensor[filenameXsensor[experiment,1]]= CSV.read(string(XsensorFolderPath,filenameXsensor[experiment,1], ".csv");header=2)

    if print_bool
        println(string("Xsensor CSV read finished, experiment : "),experiment)
    end
end

for experiment=1:size(filenameIBS,1)
    #Read CSV of IBS
    #IBS[filenameIBS[experiment,1]]=Dict()
    IBS[filenameIBS[experiment,1]]= CSV.read(string(IBSFolderPath,filenameIBS[experiment,1], ".CSV");header=2)

    if print_bool
        println(string("Instrumented Backseat CSV read finished, experiment : "),experiment)
    end
end

println("CSV read finished")


text:
# This file is for functions that just need simple standard overloading.

## Helpers:

function nameddimsarray_result(original_nda, reduced_data, reduction_dims)
    L = names(original_nda)
    return NamedDimsArray{L}(reduced_data)
end

# if reducing over `:` then results is a scalar
function nameddimsarray_result(original_nda, reduced_data, reduction_dims::Colon)
    return reduced_data
end

###################################################################################
# Overloads

# 1 Arg
for (mod, funs) in (
    (:Base, (
        :sum, :prod, :count, :maximum, :minimum, :extrema, :cumsum, :cumprod,
        :sort, :sort!,)
    ),
    (:Statistics, (:mean, :std, :var, :median, :cov, :cor)),
)
    for fun in funs
        @eval function $mod.$fun(a::NamedDimsArray; dims=:, kwargs...)
            numerical_dims = dim(a, dims)
            data = $mod.$fun(parent(a); dims=numerical_dims, kwargs...)
            return nameddimsarray_result(a, data, numerical_dims)
        end
    end
end

# 1 arg before
for (mod, funs) in (
    (:Base, (:mapslices,)),
)
    for fun in funs
        @eval function $mod.$fun(f, a::NamedDimsArray; dims=:, kwargs...)
            numerical_dims = dim(a, dims)
            data = $mod.$fun(f, parent(a); dims=numerical_dims, kwargs...)
            return nameddimsarray_result(a, data, numerical_dims)
        end
    end
end

# 2 arg before
for (mod, funs) in (
    (:Base, (:mapreduce,)),
)
    for fun in funs
        @eval function $mod.$fun(f1, f2, a::NamedDimsArray; dims=:, kwargs...)
            numerical_dims = dim(a, dims)
            data = $mod.$fun(f1, f2, parent(a); dims=numerical_dims, kwargs...)
            return nameddimsarray_result(a, data, numerical_dims)
        end
    end
end

################################################
# Non-dim Overloads

# 1 Arg
for (mod, funs) in (
    (:Base, (:zero, :one, :copy,)),
)
    for fun in funs
        @eval function $mod.$fun(a::NamedDimsArray{L}) where L
            data = $mod.$fun(parent(a))
            return NamedDimsArray{L}(data)
        end
    end
end


text:
"""
    reorder!(s::SnpData, order::Vector{<:Integer}; famnm=s.srcfam * ".reordered")
reorders the SnpData in-place, including the saved files if the sources exist. overwrites the fam file if `writefam` is true.
SnpData should be opened with "r+" mode (e.g. `Snpdata("mouse", "r+")`).
WARNING: the data are directly modified on the disk.
"""
function reorder!(s::SnpData, order::Vector{<:Integer}; famnm=split(s.srcfam, ".fam")[1] * ".reordered.fam")
    s.snparray[:,:] .= s.snparray[order, :]
    s.person_info[:, :] .= s.person_info[order, :]
    writedlm(famnm, hcat([s.person_info[!, k]
                            for k in PERSON_INFO_KEYS]...))
end


text:
import Base: *, /, \, exp, ≈, ==

"""
    QuatRotation{T} <: Rotation

4-parameter attitute representation that is singularity-free. Quaternions with unit norm
represent a double-cover of SO(3). The `QuatRotation` does NOT strictly enforce the unit
norm constraint, but certain methods will assume you have a unit quaternion.
Follows the Hamilton convention for quaternions.

# Constructors
```julia
QuatRotation(w,x,y,z)
QuatRotation(q::AbstractVector)
```
where `w` is the scalar (real) part, `x`,`y`, and `z` are the vector (imaginary) part,
and `q = [w,x,y,z]`.
"""
struct QuatRotation{T} <: Rotation{3,T}
    q::Quaternion{T}

    @inline function QuatRotation{T}(w, x, y, z, normalize::Bool = true) where T
        if normalize
            inorm = inv(sqrt(w*w + x*x + y*y + z*z))
            new{T}(Quaternion(w*inorm, x*inorm, y*inorm, z*inorm, true))
        else
            new{T}(Quaternion(w, x, y, z, true))
        end
    end

    @inline function QuatRotation{T}(q::Quaternion) where T
        if q.norm
            new{T}(q)
        else
            throw(InexactError(nameof(T), T, q))
        end
    end
    QuatRotation{T}(q::QuatRotation) where T = new{T}(q.q)
end

# ~~~~~~~~~~~~~~~ Constructors ~~~~~~~~~~~~~~~ #
# Use default map
function QuatRotation(w,x,y,z, normalize::Bool = true)
    types = promote(w,x,y,z)
    QuatRotation{eltype(types)}(w,x,y,z, normalize)
end

function QuatRotation(q::T) where T<:Quaternion
    if q.norm
        return QuatRotation(q.s, q.v1, q.v2, q.v3, false)
    else
        throw(InexactError(nameof(T), T, q))
    end
end

function Quaternions.Quaternion(q::QuatRotation)
    return q.q
end

# Pass in Vectors
@inline function (::Type{Q})(q::AbstractVector, normalize::Bool = true) where Q <: QuatRotation
    check_length(q, 4)
    Q(q[1], q[2], q[3], q[4], normalize)
end
@inline (::Type{Q})(q::StaticVector{4}, normalize::Bool = true) where Q <: QuatRotation =
    Q(q[1], q[2], q[3], q[4], normalize)

# Copy constructors
QuatRotation(q::QuatRotation) = q

# QuatRotation <=> Quat
# (::Type{Q})(q::Quat) where Q <: QuatRotation = Q(q.w, q.x, q.y, q.z, false)
# (::Type{Q})(q::QuatRotation) where Q <: Quat = Q(q.w, q.x, q.y, q.z, false)
# const AllQuats{T} = Union{<:Quat{T}, <:QuatRotation{T}}

# ~~~~~~~~~~~~~~~ StaticArrays Interface ~~~~~~~~~~~~~~~ #
function (::Type{Q})(t::NTuple{9}) where Q<:QuatRotation
    #=
    This function solves the system of equations in Section 3.1
    of https://arxiv.org/pdf/math/0701759.pdf. This cheap method
    only works for matrices that are already orthonormal (orthogonal
    and unit length columns).
    Use `nearest_rotation` to get rotation matrix from the given matrix.
    =#

    a = 1 + t[1] + t[5] + t[9]
    b = 1 + t[1] - t[5] - t[9]
    c = 1 - t[1] + t[5] - t[9]
    d = 1 - t[1] - t[5] + t[9]
    max_abcd = max(a, b, c, d)
    if a == max_abcd
        b = t[6] - t[8]
        c = t[7] - t[3]
        d = t[2] - t[4]
    elseif b == max_abcd
        a = t[6] - t[8]
        c = t[2] + t[4]
        d = t[7] + t[3]
    elseif c == max_abcd
        a = t[7] - t[3]
        b = t[2] + t[4]
        d = t[6] + t[8]
    else
        a = t[2] - t[4]
        b = t[7] + t[3]
        c = t[6] + t[8]
    end
    return principal_value(Q(a, b, c, d))
end

function Base.getindex(q::QuatRotation, i::Int)
    w = q.q.s
    x = q.q.v1
    y = q.q.v2
    z = q.q.v3

    if i == 1
        ww = (w * w)
        xx = (x * x)
        yy = (y * y)
        zz = (z * z)

        ww + xx - yy - zz
    elseif i == 2
        xy = (x * y)
        zw = (w * z)

        2 * (xy + zw)
    elseif i == 3
        xz = (x * z)
        yw = (y * w)

        2 * (xz - yw)
    elseif i == 4
        xy = (x * y)
        zw = (w * z)

        2 * (xy - zw)
    elseif i == 5
        ww = (w * w)
        xx = (x * x)
        yy = (y * y)
        zz = (z * z)

        ww - xx + yy - zz
    elseif i == 6
        yz = (y * z)
        xw = (w * x)

        2 * (yz + xw)
    elseif i == 7
        xz = (x * z)
        yw = (y * w)

        2 * (xz + yw)
    elseif i == 8
        yz = (y * z)
        xw = (w * x)

        2 * (yz - xw)
    elseif i == 9
        ww = (w * w)
        xx = (x * x)
        yy = (y * y)
        zz = (z * z)

        ww - xx - yy + zz
    else
        throw(BoundsError(r,i))
    end
end

function Base.Tuple(q::QuatRotation)
    w = q.q.s
    x = q.q.v1
    y = q.q.v2
    z = q.q.v3

    ww = (w * w)
    xx = (x * x)
    yy = (y * y)
    zz = (z * z)
    xy = (x * y)
    zw = (w * z)
    xz = (x * z)
    yw = (y * w)
    yz = (y * z)
    xw = (w * x)

    # initialize rotation part
    return (ww + xx - yy - zz,
            2 * (xy + zw),
            2 * (xz - yw),
            2 * (xy - zw),
            ww - xx + yy - zz,
            2 * (yz + xw),
            2 * (xz + yw),
            2 * (yz - xw),
            ww - xx - yy + zz)
end

# ~~~~~~~~~~~~~~~ Getters ~~~~~~~~~~~~~~~ #
@inline scalar(q::QuatRotation) = real(q.q)
@inline vector(q::QuatRotation) = vector(q.q)
@inline vector(q::Quaternion) = SVector{3}(q.v1, q.v2, q.v3)

"""
    params(R::Rotation)

Return an `SVector` of the underlying parameters used by the rotation representation.

# Example
```julia
p = MRP(1.0, 2.0, 3.0)
Rotations.params(p) == @SVector [1.0, 2.0, 3.0]  # true
```
"""
@inline params(q::QuatRotation) = SVector{4}(q.q.s, q.q.v1, q.q.v2, q.q.v3)

# ~~~~~~~~~~~~~~~ Initializers ~~~~~~~~~~~~~~~ #
function Random.rand(rng::AbstractRNG, ::Random.SamplerType{<:QuatRotation{T}}) where T
    _normalize(QuatRotation{T}(randn(rng,T), randn(rng,T), randn(rng,T), randn(rng,T)))
end
@inline Base.one(::Type{Q}) where Q <: QuatRotation = Q(1.0, 0.0, 0.0, 0.0)

# ~~~~~~~~~~~~~~~ Math Operations ~~~~~~~~~~~~~~~ #

# Inverses
Base.inv(q::Q) where Q <: QuatRotation = Q(conj(q.q))

function _normalize(q::Q) where Q <: QuatRotation
    w = q.q.s
    x = q.q.v1
    y = q.q.v2
    z = q.q.v3

    n = norm(params(q))
    Q(w/n, x/n, y/n, z/n)
end

# Identity
(::Type{Q})(I::UniformScaling) where Q <: QuatRotation = one(Q)

# Exponentials and Logarithms
"""
    _pure_quaternion(v::AbstractVector)
    _pure_quaternion(x, y, z)

Create a `Quaternion` with zero scalar part (i.e. `q.q.s == 0`).
"""
function _pure_quaternion(v::AbstractVector)
    check_length(v, 3)
    Quaternion(zero(eltype(v)), v[1], v[2], v[3], false)
end

@inline _pure_quaternion(x::Real, y::Real, z::Real) =
    Quaternion(zero(x), x, y, z, false)

function expm(ϕ::AbstractVector)
    check_length(ϕ, 3)
    θ = norm(ϕ)
    sθ,cθ = sincos(θ/2)
    M = sinc(θ/π/2)/2
    QuatRotation(cθ, ϕ[1]*M, ϕ[2]*M, ϕ[3]*M, false)
end

function _log_as_quat(q::Q, eps=1e-6) where Q <: QuatRotation
    w = q.q.s
    x = q.q.v1
    y = q.q.v2
    z = q.q.v3

    # Assumes unit quaternion
    θ = sqrt(x^2 + y^2 + z^2)
    if θ > eps
        M = atan(θ, w)/θ
    else
        M = (1-(θ^2/(3w^2)))/w
    end
    _pure_quaternion(M*vector(q))
end

function logm(q::QuatRotation)
    # Assumes unit quaternion
    return 2*vector(_log_as_quat(q))
end

# Composition
"""
    (*)(q::QuatRotation, w::QuatRotation)

Quternion Composition

Equivalent to
```julia
lmult(q) * SVector(w)
rmult(w) * SVector(q)
```

Sets the output mapping equal to the mapping of `w`
"""
function (*)(q1::QuatRotation, q2::QuatRotation)
    return QuatRotation(q1.q*q2.q)
end

"""
    (*)(q::QuatRotation, r::StaticVector)

Rotate a vector

Equivalent to `hmat()' lmult(q) * rmult(q)' hmat() * r`
"""
function Base.:*(q::QuatRotation, r::StaticVector)  # must be StaticVector to avoid ambiguity
    check_length(r, 3)
    w = real(q.q)
    v = vector(q.q)
    (w^2 - v'v)*r + 2*v*(v'r) + 2*w*cross(v,r)
end

(\)(q1::QuatRotation, q2::QuatRotation) = inv(q1)*q2
(/)(q1::QuatRotation, q2::QuatRotation) = q1*inv(q2)

"""
    rotation_between(from, to)

Compute the quaternion that rotates vector `from` so that it aligns with vector
`to`, along the geodesic (shortest path).
"""
rotation_between(from::AbstractVector, to::AbstractVector) = rotation_between(SVector{3}(from), SVector{3}(to))
function rotation_between(from::SVector{3}, to::SVector{3})
    # Robustified version of implementation from https://www.gamedev.net/topic/429507-finding-the-quaternion-betwee-two-vectors/#entry3856228
    normprod = sqrt(dot(from, from) * dot(to, to))
    T = typeof(normprod)
    normprod < eps(T) && throw(ArgumentError("Input vectors must be nonzero."))
    w = normprod + dot(from, to)
    v = abs(w) < 100 * eps(T) ? perpendicular_vector(from) : cross(from, to)
    @inbounds return QuatRotation(w, v[1], v[2], v[3]) # relies on normalization in constructor
end

# ~~~~~~~~~~~~~~~ Kinematics ~~~~~~~~~~~~~~~ $
"""
    kinematics(R::Rotation{3}, ω::AbstractVector)

The time derivative of the rotation R, according to the definition

``Ṙ = \\lim_{Δt → 0} \\frac{R(t + Δt) - R(t)}{Δt}``

where `ω` is the angular velocity. This is equivalent to

``Ṙ = \\lim_{Δt → 0} \\frac{R δR - R}{Δt}``

where ``δR`` is some small rotation, parameterized by a small rotation ``δθ`` about
an axis ``r``, such that ``lim_{Δt → 0} \\frac{δθ r}{Δt} = ω``

The kinematics are extremely useful when computing the dynamics of rigid bodies, since
`Ṙ = kinematics(R,ω)` is the first-order ODE for the evolution of the attitude dynamics.

See "Fundamentals of Spacecraft Attitude Determination and Control" by Markley and Crassidis
Sections 3.1-3.2 for more details.
"""
function kinematics(q::Q, ω::AbstractVector) where Q <: QuatRotation
    params(q*Q(0.0, ω[1], ω[2], ω[3], false))/2
end

# ~~~~~~~~~~~~~~~ Linear Algebraic Conversions ~~~~~~~~~~~~~~~ #
"""
    lmult(q::QuatRotation)
    lmult(q::StaticVector{4})

`lmult(q2)*params(q1)` returns a vector equivalent to `q2*q1` (quaternion composition)
"""
function lmult(q::QuatRotation)
    w = q.q.s
    x = q.q.v1
    y = q.q.v2
    z = q.q.v3

    SA[
        w -x -y -z;
        x  w -z  y;
        y  z  w -x;
        z -y  x  w;
    ]
end
function lmult(q::Quaternion)
    SA[
        q.s  -q.v1 -q.v2 -q.v3;
        q.v1  q.s  -q.v3  q.v2;
        q.v2  q.v3  q.s  -q.v1;
        q.v3 -q.v2  q.v1  q.s;
    ]
end
lmult(q::StaticVector{4}) = lmult(QuatRotation(q, false))

"""
    rmult(q::QuatRotation)
    rmult(q::StaticVector{4})

`rmult(q1)*params(q2)` return a vector equivalent to `q2*q1` (quaternion composition)
"""
function rmult(q::QuatRotation)
    w = q.q.s
    x = q.q.v1
    y = q.q.v2
    z = q.q.v3

    SA[
        w -x -y -z;
        x  w  z -y;
        y -z  w  x;
        z  y -x  w;
    ]
end
function rmult(q::Quaternion)
    SA[
        q.s  -q.v1 -q.v2 -q.v3;
        q.v1  q.s   q.v3 -q.v2;
        q.v2 -q.v3  q.s   q.v1;
        q.v3  q.v2 -q.v1  q.s;
    ]
end
rmult(q::SVector{4}) = rmult(QuatRotation(q, false))

"""
    tmat()

`tmat()*params(q)`return a vector equivalent to `inv(q)`, where `q` is a `QuatRotation`
"""
function tmat(::Type{T}=Float64) where T
    SA{T}[
        1  0  0  0;
        0 -1  0  0;
        0  0 -1  0;
        0  0  0 -1;
    ]
end

"""
    vmat()

`vmat()*params(q)`` returns the imaginary
    (vector) part of the quaternion `q` (equivalent to `vector(q)``)
"""
function vmat(::Type{T}=Float64) where T
    SA{T}[
        0 1 0 0;
        0 0 1 0;
        0 0 0 1
    ]
end

"""
    hmat()
    hmat(r::AbstractVector)

`hmat()*r` or `hmat(r)` converts `r` into a pure quaternion, where `r` is 3-dimensional.

`hmat() == vmat()'`
"""
function hmat(::Type{T}=Float64) where T
    SA{T}[
        0 0 0;
        1 0 0;
        0 1 0;
        0 0 1.;
    ]
end

function hmat(r)
    @assert length(r) == 3
    SA[0, r[1], r[2], r[3]]
end

# ~~~~~~~~~~~~~~~ Useful Jacobians ~~~~~~~~~~~~~~~ #
"""
    ∇differential(q::QuatRotation)

Jacobian of `lmult(q) QuatMap(ϕ)`, when ϕ is near zero.

Useful for converting Jacobians from R⁴ to R³ and
    correctly account for unit norm constraint. Jacobians for different
    differential quaternion parameterization are the same up to a constant.
"""
function ∇differential(q::QuatRotation)
    w = q.q.s
    x = q.q.v1
    y = q.q.v2
    z = q.q.v3

    SA[
        -x -y -z;
         w -z  y;
         z  w -x;
        -y  x  w;
    ]
end

"""
    ∇²differential(q::QuatRotation, b::AbstractVector)

Jacobian of `(∂/∂ϕ lmult(q) QuatMap(ϕ))`b, evaluated at ϕ=0, and `b` has length 4.
"""
function ∇²differential(q::QuatRotation, b::AbstractVector)
    check_length(b, 4)
    b1 = -params(q)'b
    Diagonal(@SVector fill(b1,3))
end

"""
    ∇rotate(R::Rotation{3}, r::AbstractVector)

Jacobian of `R*r` with respect to the rotation
"""
function ∇rotate(q::QuatRotation, r::AbstractVector)
    check_length(r, 3)
    rhat = QuatRotation(zero(eltype(r)), r[1], r[2], r[3], false)
    2vmat()*rmult(q)'rmult(rhat)
end

"""
    ∇composition1(R2::Rotation{3}, R1::Rotation{3})

Jacobian of `R2*R1` with respect to `R1`
"""
function ∇composition1(q2::QuatRotation, q1::QuatRotation)
    lmult(q2)
end

"""
    ∇composition2(R2::Rotation{3}, R1::Rotation{3})

Jacobian of `R2*R1` with respect to `R2`
"""
function ∇composition2(q2::QuatRotation, q1::QuatRotation)
    rmult(q1)
end


text:
module GeoArrays

using GeoFormatTypes
using ArchGDAL
using CoordinateTransformations
using StaticArrays
const GFT = GeoFormatTypes

include("geoarray.jl")
include("geoutils.jl")
include("utils.jl")
include("io.jl")
include("interpolate.jl")
include("crs.jl")
include("operations.jl")
include("plot.jl")

export GeoArray

export coords
export centercoords
export indices

export bbox
export bbox!
export bboxes

export compose!

export interpolate!

export epsg!
export crs!

export -,+,*,/

end


text:
export varsize

"""
    varsize(S)

The number of elements specified by the template type `S`.
"""
varsize(::Type{T}) where {T <: Real} = 1
varsize(::Type{Tuple{}}) = 0
varsize(::Type{NamedTuple{(), Tuple{}}}) = 0
varsize(::Type{SVector{N, T}}) where {N, T <: Real} = N

@generated function varsize(::Type{S}) where {S}
    types = fieldtypes(S)
    isempty(types) ? 0 : sum(varsize, types)
end

varsize(nt::Tuple) = sum([varsize(nt[i]) for i in 1:length(nt)])

function varsize(nt::NamedTuple)
    if Symbol(:type) in keys(nt)
        return sum(get_len(get_val(nt.type)))
    else
        v = values(nt)
        if v isa Tuple{}
            return 0
        else
            return sum(varsize.(v))
        end
    end
end


text:
#
# As of September 19 2021
# Tests against the Age UK Pensioner calculator.
# https://benefitscheck.ageuk.org.uk/
# ref: AC73A0470
#

using Test
using Dates
using ScottishTaxBenefitModel
using .ModelHousehold
using .STBParameters
using .STBIncomes
using .Definitions
using .GeneralTaxComponents
using .SingleHouseholdCalculations
using .RunSettings
using .Utils
using .ExampleHelpers

sys21_22 = load_file( "../params/sys_2021_22.jl" )
load_file!( sys21_22, "../params/sys_2021-uplift-removed.jl")
load_file!( sys21_22, "../params/budget_2021_uc_changes.jl")

println( "weeklyise start wpm=$PWPM wpy=52")
weeklyise!( sys21_22; wpy=52, wpm=PWPM  )
settings = Settings()

@testset "Single Pensioner from AGEUK" begin

    hh = make_hh(
        adults = 1,
        children = 0,
        earnings = 0,
        rent = 0,
        rooms = 0,
        age = 68,
        tenure = Private_Rented_Furnished )
    head = get_head( hh )
    settings.means_tested_routing = lmt_full

    head = get_head( hh )
    empty!( head.income )
    unemploy!( head )
    enable!( head )

    hres = do_one_calc( hh, sys21_22, settings )
    @test hres.ahc_net_income ≈ 179.60
    println( to_md_table(hres.bus[1].legacy_mtbens ))
    println( "Age 68; \n"*inctostr(  hres.bus[1].pers[head.pid].income ))

    head.age = 80

    hres = do_one_calc( hh, sys21_22, settings )
    @test hres.ahc_net_income ≈  (137.60*1.1) + 25.74
    println( "Age 80; old pension x 1.1 \n"*inctostr(  hres.bus[1].pers[head.pid].income ))

    head.benefit_ratios[state_pension] = 154.0/137.60 # just qualify for savings credit
    hres = do_one_calc( hh, sys21_22, settings )
    # pen / pen credit / savings credit
    @test hres.ahc_net_income ≈ 154 + 23.10 + 0.18
    println( "Age 80; old pension x 1.119 so qualify for savings credit \n"*inctostr(  hres.bus[1].pers[head.pid].income ))

    head.benefit_ratios[state_pension] = 1
    pid1 = add_child!( hh, 3, Female ) # 01/01/2018
    head.age=68 # 01/01/1953

    # £98.76 = 52.10PC + 10.00 SCP + 15.51 CTR + 21.15 CB + 179.60PENS
    hres = do_one_calc( hh, sys21_22, settings )
    @test hres.bhc_net_income ≈ (179.60 + 98.76) atol=0.04 # FIXME check weird rounding

    hh.gross_rent = 100
    hres = do_one_calc( hh, sys21_22, settings )
    @test hres.bhc_net_income ≈ (179.60 + 198.76) atol=0.04 # FIXME check weird rounding
    println( to_md_table(hres.bus[1].legacy_mtbens ))
    println( "Age 68; + 1 chld no rent \n"*inctostr(  hres.bus[1].income ))

end

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
"expected_arc(i,n) this polynomial indicates the expected number of NN-arcs that would pass over index `i` for a series of length`n` if the series was completely random. This is used to normalize the `segment_profile` to account for bias towards the edges."
expected_arc(i,n) = 2*(-i^2*n/(n-1)^2 + i*n/(n-1)) # polynomial coefficients solved for by using constraints e(0) = 0, e(n-1) = 0, e(n/2) = n/2

"""
    segment_profile(p::Profile)

Calculate the MP-index of a matrix profile. This index has the same length as the profile, and tells you how many nearest-neighbor arcs passes over index `i`. It's normalized by the expected number of arcs for a completely random time series, so that a value of 1 indicates a poor segmentation point, and a value close to 0 indicates a likely segmentation point.
"""
function segment_profile(p::Profile)
    I = p.I
    n = length(I)
    mpi = zeros(eltype(p.P), n)
    for i in 1:n
        Ii = I[i]
        if Ii > i
            for j = i:Ii-1
                mpi[j] += 1
            end
        else
            for j = Ii+1:i
                mpi[j] += 1
            end
        end
    end
    mpi .= min.(mpi ./ expected_arc.(0:n-1, n), 1)
    mpi
end

"""
    i = segment(p::Profile)

Returns an index `i` indicating the most likely segmentation point of profile `p`, i.e., the point which the fewst nearest-neighbor arcs passes over.

Ref: Matrix Profile VIII: Domain Agnostic Online Semantic
Segmentation at Superhuman Performance Levels
"""
function segment(p::Profile)
    argmin(segment_profile(p))
end


text:
#=
    The File contains types and interfaces definitions
=#

# All interfaces defined here must provide one-based indexing,
# yet internal indexing is implementation-defined

#------------------------------------------------------------------------------

# AbstractSparseObject<T>
#
# `T` is the type of the base_ring field
#
# Base interface for anything Sparse implemented in this project
"""
AbstractSparseObject{T}

Base interface for anything Sparse implemented in this project

`T` is the type of the ground field
"""
abstract type AbstractSparseObject{T} end

#------------------------------------------------------------------------------

# AbstractSparseVector<T>
#
# `T` is the type if the base_ring field
#
# Base interface for Sparse vectors family
# An object which implements `AbstractSparseVector` can be treated
# as an element of Rⁿ vectorspace. Thus, some vectorspace operations and
# left matrix multiplication are supported for such objects
"""
AbstractSparseObject{T}

Base interface for 1-dimensional Sparse vector. An object which implements `AbstractSparseVector` can be treated
as an element of Rⁿ vectorspace. Thus, some vectorspace operations and
left matvec are supported for such objects.

`T` is the type of the ground field
"""
abstract type AbstractSparseVector{T} <: AbstractSparseObject{T} end

#------------------------------------------------------------------------------

# AbstractSparseMatrix<T>
#
# `T` is the type if the base_ring field
#
# Base interface for Sparse matrices family
# An object which implements `AbstractSparseMatrix` can be treated
# as an element from a Matrix algebra
# An object which implements `AbstractSparseMatrix` can be treated
# as a linear operator in Rⁿ
#
#
# Implementing some functions for other sparse formats
# rather than DOK is tricky,
# probably, this interface it too strong to be applied to every class
"""
AbstractSparseObject{T}

Base interface for 2-dimensional Sparse vector, i.e Sparse Matrix. An object which implements `AbstractSparseMatrix` can be treated
as an element of Rⁿ⃰ⁿ vectorspace or as an element of a matrix algebra. Thus, some vectorspace operations and matrix prodducts are supported for such objects.

`T` is the type of the ground field
"""
abstract type AbstractSparseMatrix{T} <: AbstractSparseObject{T} end

#------------------------------------------------------------------------------

function Nemo.base_ring(::AbstractSparseMatrix) end

function first_nonzero(::T) where {T<:AbstractSparseMatrix} end
function getindex(::AbstractSparseMatrix, i::Int) end

function scale(::AbstractSparseMatrix, ::C) where {C} end
function scale!(::AbstractSparseMatrix, ::C) where {C} end

# To Be Changed sometime
# !!!
# A tender for a better name!
function Base.reduce(::T, ::T, ::C) where {T<:AbstractSparseMatrix{F}} where {F, C} end
function reduce!(::T, ::T, ::C) where {T<:AbstractSparseMatrix{F}} where {F, C} end

function inner(::AbstractSparseMatrix, ::AbstractSparseMatrix) end

# hmmm
function Base.prod(::T, ::T) where {T<:AbstractSparseMatrix{F}} where {F} end

function Base.zero(::AbstractSparseMatrix) end
function Base.isempty(::AbstractSparseMatrix) end
function Base.iszero(::AbstractSparseMatrix) end

function dim(::AbstractSparseMatrix) end
function Base.length(::AbstractSparseMatrix) end
function Base.size(::AbstractSparseMatrix) end
function Base.size(::AbstractSparseMatrix, i::Int) end

function density(::AbstractSparseMatrix) end

function ==(::T, ::T) where {T<:AbstractSparseMatrix{F}} where {F} end
function !=(::T, ::T) where {T<:AbstractSparseMatrix{F}} where {F} end

function Base.valtype(::AbstractSparseMatrix) end
function Base.eltype(::AbstractSparseMatrix) end

function Base.show(::AbstractSparseMatrix) end
function print_sparsik(::AbstractSparseMatrix) end

function modular_reduction(::AbstractSparseMatrix, ::T) where {T} end
function rational_reconstruction(::AbstractSparseMatrix) end

function Base.iterate(::AbstractSparseMatrix, state) end

function unit_vector(::Int, ::Int, ::T) where {T} end

function random_sparsik(::Tuple{Int}, ::T) where {T} end

#------------------------------------------------------------------------------