text:
struct Quaternion{T<:Real} <: Number
    s::T
    v1::T
    v2::T
    v3::T
    norm::Bool
end

const QuaternionF16 = Quaternion{Float16}
const QuaternionF32 = Quaternion{Float32}
const QuaternionF64 = Quaternion{Float64}

Quaternion{T}(x::Real) where {T<:Real} = Quaternion(convert(T, x))
Quaternion{T}(x::Complex) where {T<:Real} = Quaternion(convert(Complex{T}, x))
Quaternion{T}(q::Quaternion) where {T<:Real} = Quaternion{T}(q.s, q.v1, q.v2, q.v3, q.norm)
Quaternion(s::Real, v1::Real, v2::Real, v3::Real, n::Bool = false) =
    Quaternion(promote(s, v1, v2, v3)..., n)
Quaternion(x::Real) = Quaternion(x, zero(x), zero(x), zero(x), abs(x) == one(x))
Quaternion(z::Complex) = Quaternion(z.re, z.im, zero(z.re), zero(z.re), abs(z) == one(z.re))
Quaternion(s::Real, a::AbstractVector) = Quaternion(s, a[1], a[2], a[3])
Quaternion(a::AbstractVector) = Quaternion(0, a[1], a[2], a[3])

promote_rule(::Type{Quaternion{T}}, ::Type{S}) where {T <: Real, S <: Real} = Quaternion{promote_type(T, S)}
promote_rule(::Type{Quaternion{T}}, ::Type{Complex{S}}) where {T <: Real, S <: Real} = Quaternion{promote_type(T, S)}
promote_rule(::Type{Quaternion{T}}, ::Type{Quaternion{S}}) where {T <: Real, S <: Real} = Quaternion{promote_type(T, S)}

quat(p, v1, v2, v3) = Quaternion(p, v1, v2, v3)
quat(p, v1, v2, v3, n) = Quaternion(p, v1, v2, v3, n)
quat(x) = Quaternion(x)
quat(s, a) = Quaternion(s, a)

function show(io::IO, q::Quaternion)
    pm(x) = x < 0 ? " - $(-x)" : " + $x"
    print(io, q.s, pm(q.v1), "im", pm(q.v2), "jm", pm(q.v3), "km")
end

real(::Type{Quaternion{T}}) where {T} = T
real(q::Quaternion) = q.s
imag_part(q::Quaternion) = (q.v1, q.v2, q.v3)
@deprecate imag(q::Quaternion) collect(imag_part(q)) false

(/)(q::Quaternion, x::Real) = Quaternion(q.s / x, q.v1 / x, q.v2 / x, q.v3 / x)

conj(q::Quaternion) = Quaternion(q.s, -q.v1, -q.v2, -q.v3, q.norm)
abs(q::Quaternion) = sqrt(q.s * q.s + q.v1 * q.v1 + q.v2 * q.v2 + q.v3 * q.v3)
float(q::Quaternion{T}) where T = convert(Quaternion{float(T)}, q)
abs_imag(q::Quaternion) = sqrt(q.v1 * q.v1 + q.v2 * q.v2 + q.v3 * q.v3)
abs2(q::Quaternion) = q.s * q.s + q.v1 * q.v1 + q.v2 * q.v2 + q.v3 * q.v3
inv(q::Quaternion) = q.norm ? conj(q) : conj(q) / abs2(q)

isreal(q::Quaternion) = iszero(q.v1) & iszero(q.v2) & iszero(q.v3)
isfinite(q::Quaternion) = q.norm | (isfinite(q.s) & isfinite(q.v1) & isfinite(q.v2) & isfinite(q.v3))
iszero(q::Quaternion) = ~q.norm & iszero(real(q)) & iszero(q.v1) & iszero(q.v2) & iszero(q.v3)
isnan(q::Quaternion) = isnan(real(q)) | isnan(q.v1) | isnan(q.v2) | isnan(q.v3)
isinf(q::Quaternion) = ~q.norm & (isinf(q.s) | isinf(q.v1) | isinf(q.v2) | isinf(q.v3))

function normalize(q::Quaternion)
    if (q.norm)
        return q
    end
    q = q / abs(q)
    Quaternion(q.s, q.v1, q.v2, q.v3, true)
end

function normalizea(q::Quaternion)
    if (q.norm)
        return (q, one(q.s))
    end
    a = abs(q)
    q = q / a
    (Quaternion(q.s, q.v1, q.v2, q.v3, true), a)
end

function normalizeq(q::Quaternion)
    a = abs(q)
    if a > 0
        q = q / a
        Quaternion(q.s, q.v1, q.v2, q.v3, true)
    else
        Quaternion(0.0, 1.0, 0.0, 0.0, true)
    end
end

(-)(q::Quaternion) = Quaternion(-q.s, -q.v1, -q.v2, -q.v3, q.norm)

(+)(q::Quaternion, w::Quaternion) =
    Quaternion(q.s + w.s, q.v1 + w.v1, q.v2 + w.v2, q.v3 + w.v3)

(-)(q::Quaternion, w::Quaternion) =
    Quaternion(q.s - w.s, q.v1 - w.v1, q.v2 - w.v2, q.v3 - w.v3)

(*)(q::Quaternion, w::Quaternion) = Quaternion(q.s * w.s - q.v1 * w.v1 - q.v2 * w.v2 - q.v3 * w.v3,
                                               q.s * w.v1 + q.v1 * w.s + q.v2 * w.v3 - q.v3 * w.v2,
                                               q.s * w.v2 - q.v1 * w.v3 + q.v2 * w.s + q.v3 * w.v1,
                                               q.s * w.v3 + q.v1 * w.v2 - q.v2 * w.v1 + q.v3 * w.s,
                                               q.norm && w.norm)
(/)(q::Quaternion, w::Quaternion) = q * inv(w)

(==)(q::Quaternion, w::Quaternion) = (q.s == w.s) & (q.v1 == w.v1) & (q.v2 == w.v2) & (q.v3 == w.v3) # ignore .norm field

angleaxis(q::Quaternion) = angle(q), axis(q)

angle(q::Quaternion) = 2 * atan(√(q.v1^2 + q.v2^2 + q.v3^2), q.s)

function axis(q::Quaternion)
    q = normalize(q)
    s = sin(angle(q) / 2)
    abs(s) > 0 ?
        [q.v1, q.v2, q.v3] / s :
        [1.0, 0.0, 0.0]
end

argq(q::Quaternion) = normalizeq(Quaternion(0, q.v1, q.v2, q.v3))

"""
    extend_analytic(f, q::Quaternion)

Evaluate the extension of the complex analytic function `f` to the quaternions at `q`.

Given ``q = s + a u``, where ``s`` is the real part, ``u`` is a pure unit quaternion,
and ``a \\ge 0`` is the magnitude of the imaginary part of ``q``,

```math
f(q) = \\Re(f(z)) + \\Im(f(z)) u,
```
is the extension of `f` to the quaternions, where ``z = a + s i`` is a complex analog to
``q``.

See Theorem 5 of [^Sudbery1970] for details.

[^Sudbery1970]
    Sudbery (1979). Quaternionic analysis. Mathematical Proceedings of the Cambridge
    Philosophical Society,85, pp 199225
    doi:[10.1017/S030500410005563](https://doi.org/10.1017/S0305004100055638)
"""
function extend_analytic(f, q::Quaternion)
    a = abs_imag(q)
    s = q.s
    z = complex(s, a)
    w = f(z)
    wr, wi = reim(w)
    scale = wi / a
    norm = _isexpfun(f) && iszero(s)
    if a > 0
        return Quaternion(wr, scale * q.v1, scale * q.v2, scale * q.v3, norm)
    else
        # q == real(q), so f(real(q)) may be real or complex, i.e. wi may be nonzero.
        # we choose to embed complex numbers in the quaternions by identifying the first
        # imaginary quaternion basis with the complex imaginary basis.
        return Quaternion(wr, oftype(scale, wi), zero(scale), zero(scale), norm)
    end
end

_isexpfun(::Union{typeof(exp),typeof(exp2),typeof(exp10)}) = true
_isexpfun(::Any) = false

"""
    cis(q::Quaternion)

Return ``\\exp(u * q)``, where ``u`` is the normalized imaginary part of `q`.

Let ``q = s + a u``, where ``s`` is the real part, ``u`` is a pure unit quaternion,
and ``a \\ge 0`` is the magnitude of the imaginary part of ``q``.

!!! Note
    This is the extension of `cis(z)` for complex `z` to the quaternions and is not
    equivalent to `exp(im * q)`. As a result, `cis(Quaternion(z)) ≠ cis(z)` when
    `imag(z) < 0`.
"""
cis(q::Quaternion)

if VERSION ≥ v"1.6"
    """
        cispi(q::Quaternion)

    Compute `cis(π * q)` more accurately.

    !!! Note
        This is not equivalent to `exp(π*im*q)`. See [cis(::Quaternion)](@ref) for details.
    """
    Base.cispi(q::Quaternion) = extend_analytic(cispi, q)
end

for f in (
    :sqrt, :exp, :exp2, :exp10, :expm1, :log2, :log10, :log1p, :cis,
    :sin, :cos, :tan, :asin, :acos, :atan, :sinh, :cosh, :tanh, :asinh, :acosh, :atanh,
    :csc, :sec, :cot, :acsc, :asec, :acot, :csch, :sech, :coth, :acsch, :asech, :acoth,
    :sinpi, :cospi,
)
    @eval Base.$f(q::Quaternion) = extend_analytic($f, q)
end

for f in (@static(VERSION ≥ v"1.6" ? (:sincos, :sincospi) : (:sincos,)))
    @eval begin
        function Base.$f(q::Quaternion)
            a = abs_imag(q)
            z = complex(q.s, a)
            s, c = $f(z)
            sr, si = reim(s)
            cr, ci = reim(c)
            sscale = si / a
            cscale = ci / a
            if a > 0
                return (
                    Quaternion(sr, sscale * q.v1, sscale * q.v2, sscale * q.v3),
                    Quaternion(cr, cscale * q.v1, cscale * q.v2, cscale * q.v3),
                )
            else
                return (
                    Quaternion(sr, oftype(sscale, si), zero(sscale), zero(sscale)),
                    Quaternion(cr, oftype(cscale, ci), zero(cscale), zero(cscale)),
                )
            end
        end
    end
end

function log(q::Quaternion)
    a = abs(q)
    M = abs_imag(q)
    theta = atan(M, q.s)
    scale = theta / ifelse(iszero(M), oneunit(M), M)
    return Quaternion(log(a), q.v1 * scale, q.v2 * scale, q.v3 * scale)
end

(^)(q::Quaternion, w::Quaternion) = exp(w * log(q))

function linpol(p::Quaternion, q::Quaternion, t::Real)
    p = normalize(p)
    q = normalize(q)
    qm = -q
    if abs(p - q) > abs(p - qm)
        q = qm
    end
    c = p.s * q.s + p.v1 * q.v1 + p.v2 * q.v2 + p.v3 * q.v3
    if c < 1.0
        o = acos(c)
        s = sin(o)
        sp = sin((1 - t) * o) / s
        sq = sin(t * o) / s
    else
        sp = 1 - t
        sq = t
    end
    Quaternion(sp * p.s  + sq * q.s,
                sp * p.v1 + sq * q.v1,
                sp * p.v2 + sq * q.v2,
                sp * p.v3 + sq * q.v3, true)
end

quatrand(rng = Random.GLOBAL_RNG)  = quat(randn(rng), randn(rng), randn(rng), randn(rng))
nquatrand(rng = Random.GLOBAL_RNG) = normalize(quatrand(rng))

function rand(rng::AbstractRNG, ::Random.SamplerType{Quaternion{T}}) where {T<:Real}
    Quaternion{T}(rand(rng, T), rand(rng, T), rand(rng, T), rand(rng, T), false)
end

function randn(rng::AbstractRNG, ::Type{Quaternion{T}}) where {T<:AbstractFloat}
    Quaternion{T}(
        randn(rng, T) * 1//2,
        randn(rng, T) * 1//2,
        randn(rng, T) * 1//2,
        randn(rng, T) * 1//2,
        false,
    )
end

## Rotations

function qrotation(axis::AbstractVector{T}, theta) where {T <: Real}
    if length(axis) != 3
        error("Must be a 3-vector")
    end
    normaxis = norm(axis)
    if iszero(normaxis)
        normaxis = oneunit(normaxis)
        theta = zero(theta)
    end
    s,c = sincos(theta / 2)
    scaleby = s / normaxis
    Quaternion(c, scaleby * axis[1], scaleby * axis[2], scaleby * axis[3], true)
end

# Variant of the above where norm(rotvec) encodes theta
function qrotation(rotvec::AbstractVector{T}) where {T <: Real}
    if length(rotvec) != 3
        error("Must be a 3-vector")
    end
    theta = norm(rotvec)
    s,c = sincos(theta / 2)
    scaleby = s / (iszero(theta) ? one(theta) : theta)
    Quaternion(c, scaleby * rotvec[1], scaleby * rotvec[2], scaleby * rotvec[3], true)
end

function qrotation(dcm::AbstractMatrix{T}) where {T<:Real}
    # See https://arxiv.org/pdf/math/0701759.pdf
    a2 = 1 + dcm[1,1] + dcm[2,2] + dcm[3,3]
    b2 = 1 + dcm[1,1] - dcm[2,2] - dcm[3,3]
    c2 = 1 - dcm[1,1] + dcm[2,2] - dcm[3,3]
    d2 = 1 - dcm[1,1] - dcm[2,2] + dcm[3,3]

    if a2 >= max(b2, c2, d2)
        a = sqrt(a2)/2
        return Quaternion(a, (dcm[3,2]-dcm[2,3])/4a, (dcm[1,3]-dcm[3,1])/4a, (dcm[2,1]-dcm[1,2])/4a)
    elseif b2 >= max(c2, d2)
        b = sqrt(b2)/2
        return Quaternion((dcm[3,2]-dcm[2,3])/4b, b, (dcm[2,1]+dcm[1,2])/4b, (dcm[1,3]+dcm[3,1])/4b)
    elseif c2 >= d2
        c = sqrt(c2)/2
        return Quaternion((dcm[1,3]-dcm[3,1])/4c, (dcm[2,1]+dcm[1,2])/4c, c, (dcm[3,2]+dcm[2,3])/4c)
    else
        d = sqrt(d2)/2
        return Quaternion((dcm[2,1]-dcm[1,2])/4d, (dcm[1,3]+dcm[3,1])/4d, (dcm[3,2]+dcm[2,3])/4d, d)
    end
end

function qrotation(dcm::AbstractMatrix{T}, qa::Quaternion) where {T<:Real}
    q = qrotation(dcm)
    abs(q-qa) < abs(q+qa) ? q : -q
end

rotationmatrix(q::Quaternion) = rotationmatrix_normalized(normalize(q))

function rotationmatrix_normalized(q::Quaternion)
    sx, sy, sz = 2q.s * q.v1, 2q.s * q.v2, 2q.s * q.v3
    xx, xy, xz = 2q.v1^2, 2q.v1 * q.v2, 2q.v1 * q.v3
    yy, yz, zz = 2q.v2^2, 2q.v2 * q.v3, 2q.v3^2
    [1 - (yy + zz)     xy - sz     xz + sy;
        xy + sz   1 - (xx + zz)    yz - sx;
        xz - sy      yz + sx  1 - (xx + yy)]
end

function slerp(qa::Quaternion{T}, qb::Quaternion{T}, t::T) where {T}
    # http://www.euclideanspace.com/maths/algebra/realNormedAlgebra/quaternions/slerp/
    coshalftheta = qa.s * qb.s + qa.v1 * qb.v1 + qa.v2 * qb.v2 + qa.v3 * qb.v3;

    if coshalftheta < 0
        qm = -qb
        coshalftheta = -coshalftheta
    else
        qm = qb
    end
    abs(coshalftheta) >= 1.0 && return qa

    halftheta    = acos(coshalftheta)
    sinhalftheta = sqrt(one(T) - coshalftheta * coshalftheta)

    if abs(sinhalftheta) < 0.001
        return Quaternion(
            T(0.5) * (qa.s  + qb.s),
            T(0.5) * (qa.v1 + qb.v1),
            T(0.5) * (qa.v2 + qb.v2),
            T(0.5) * (qa.v3 + qb.v3),
        )
    end

    ratio_a = sin((one(T) - t) * halftheta) / sinhalftheta
    ratio_b = sin(t * halftheta) / sinhalftheta

    Quaternion(
        qa.s  * ratio_a + qm.s  * ratio_b,
        qa.v1 * ratio_a + qm.v1 * ratio_b,
        qa.v2 * ratio_a + qm.v2 * ratio_b,
        qa.v3 * ratio_a + qm.v3 * ratio_b,
    )
end

function sylvester(a::Quaternion{T}, b::Quaternion{T}, c::Quaternion{T}) where {T<:Real}
    isreal(a) && return sylvester(real(a), b, c)
    isreal(b) && return sylvester(a, real(b), c)
    abs2a = abs2(a)
    abs2b = abs2(b)
    if abs2a > abs2b
        inva = conj(a) / abs2a
        d1 = -2real(b) - a - inva * abs2b
        x = d1 \ (c + inva * c * conj(b))
    else
        invb = conj(b) / abs2b
        d2 = -2real(a) - b - invb * abs2a
        x = (c + conj(a) * c * invb) / d2
    end
    return x
end
sylvester(a::Quaternion, b::Quaternion, c::Quaternion) = sylvester(promote(a, b, c)...)
sylvester(a::Quaternion, b::Quaternion, c::Real) = sylvester(promote(a, b, c)...)
# if either a or b commute with x, use a simpler expression
sylvester(a::Real, b::Real, c::Quaternion) = c / -(a + b)
sylvester(a::Real, b::Quaternion, c::Quaternion) = c / -(a + b)
sylvester(a::Quaternion, b::Real, c::Quaternion) = -(a + b) \ c
sylvester(a::Real, b::Quaternion, c::Real) = -c / (a + b)
sylvester(a::Quaternion, b::Real, c::Real) = (a + b) \ -c

function lyap(a::Quaternion{T}, c::Quaternion{T}) where {T<:Real}
    # if a commutes with c, use a simpler expression
    (isreal(a) || isreal(c)) && return c / -2real(a)
    return (c + a \ c * a) / -4real(a)
end
lyap(a::Quaternion, c::Quaternion) = lyap(promote(a, c)...)
# if a commutes with c, use a simpler expression
lyap(a::Real, c::Quaternion) = c / -2a
lyap(a::Quaternion, c::Real) = c / -2real(a)


text:
# This file is part of the TaylorIntegration.jl package; MIT licensed

using TaylorIntegration
using Test
using LinearAlgebra: norm
using Logging
import Logging: Warn

@testset "Testing `one_ode.jl`" begin

    local _order = 28
    local _abstol = 1.0E-20
    local tT = Taylor1(_order)

    max_iters_reached() = "Maximum number of integration steps reached; exiting.\n"

    @testset "Tests: dot{x}=x^2, x(0) = 1" begin
        eqs_mov(x, p, t) = x^2
        t0 = 0.0
        x0 = 1.0
        x0T = Taylor1(x0, _order)
        tT[1] = t0
        TaylorIntegration.jetcoeffs!(eqs_mov, tT, x0T, nothing)
        @test x0T.coeffs[end] == 1.0
        δt = _abstol^inv(_order-1)
        @test TaylorIntegration.stepsize(x0T, _abstol) == δt

        tv, xv = (@test_logs (Warn, max_iters_reached()) taylorinteg(
            eqs_mov, 1, 0.0, 1.0, _order, _abstol))
        @test length(tv) == 501
        @test length(xv) == 501
        @test xv[1] == x0
        @test tv[end] < 1.0

        tv, xv = (@test_logs (Warn, max_iters_reached()) taylorinteg(
            eqs_mov, x0, 0.0, 1.0, _order, _abstol, nothing))
        @test length(tv) == 501
        @test length(xv) == 501
        @test xv[1] == x0
        @test tv[end] < 1.0

        trange = 0.0:1/8:1.0
        xv = (@test_logs (Warn, max_iters_reached()) taylorinteg(
            eqs_mov, 1, trange, _order, _abstol))
        @test length(xv) == length(trange)
        @test typeof(xv) == Array{typeof(x0),1}
        @test xv[1] == x0
        @test isnan(xv[end])
        @test abs(xv[5] - 2.0) ≤ eps(2.0)
        tvr, xvr = (@test_logs min_level=Logging.Warn taylorinteg(
            eqs_mov, x0, trange[1], trange[end-1], _order, _abstol))
        @test tvr[1] == trange[1]
        @test tvr[end] == trange[end-1]
        @test xvr[1] == xv[1]
        @test xvr[end] == xv[end-1]

        trange = 0.0:1/8:1.0
        xv = (@test_logs (Warn, max_iters_reached()) taylorinteg(
            eqs_mov, x0, trange, _order, _abstol, nothing))
        @test length(xv) == length(trange)
        @test typeof(xv) == Array{typeof(x0),1}
        @test xv[1] == x0
        @test isnan(xv[end])
        @test abs(xv[5] - 2.0) ≤ eps(2.0)
        tvr, xvr = (@test_logs min_level=Logging.Warn taylorinteg(
            eqs_mov, x0, trange[1], trange[end-1], _order, _abstol, nothing))
        @test tvr[1] == trange[1]
        @test tvr[end] == trange[end-1]
        @test xvr[1] == xv[1]
        @test xvr[end] == xv[end-1]

        tarray = vec(trange)
        xv2 = (@test_logs (Warn, max_iters_reached()) taylorinteg(
            eqs_mov, x0, tarray, _order, _abstol))
        @test xv[1:end-1] == xv2[1:end-1]
        @test xv2[1:end-1] ≈ xv[1:end-1] atol=eps() rtol=0.0
        @test length(xv2) == length(tarray)
        @test typeof(xv2) == Array{typeof(x0),1}
        @test xv2[1] == x0
        @test isnan(xv2[end])
        @test abs(xv2[5] - 2.0) ≤ eps(2.0)

        # Output includes Taylor polynomial solution
        tv, xv, polynV = (@test_logs (Warn, max_iters_reached()) taylorinteg(
            eqs_mov, x0, 0, 0.5, _order, _abstol, Val(true), maxsteps=2))
        @test length(polynV) == 3
        @test xv[1] == x0
        @test polynV[1] == Taylor1(x0, _order)
        @test xv[2] == evaluate(polynV[2], tv[2]-tv[1])
        @test polynV[2] == Taylor1(ones(_order+1))
    end

    @testset "Tests: dot{x}=x^2, x(0) = 3; nsteps <= maxsteps" begin
        eqs_mov(x, p, t) = x.^2 #the ODE (i.e., the equations of motion)
        exactsol(t, x0) = x0/(1.0-x0*t) #the analytical solution
        t0 = 0.0
        tmax = 0.3
        x0 = 3.0
        x0T = Taylor1(x0, _order)
        tT[1] = t0
        TaylorIntegration.jetcoeffs!(eqs_mov, tT, x0T, nothing)
        @test x0T.coeffs[end] == 3.0^(_order+1)
        δt = (_abstol/x0T.coeffs[end-1])^inv(_order-1)
        @test TaylorIntegration.stepsize(x0T, _abstol) == δt

        tv, xv = (@test_logs min_level=Logging.Warn taylorinteg(
            eqs_mov, x0, 0, tmax, _order, _abstol))
        @test length(tv) < 501
        @test length(xv) < 501
        @test length(tv) == 14
        @test length(xv) == 14
        @test xv[1] == x0
        @test tv[end] < 1/3
        @test tv[end] == tmax
        @test abs(xv[end]-exactsol(tv[end], xv[1])) < 5e-14

        tmax = 0.33
        tv, xv = (@test_logs min_level=Logging.Warn taylorinteg(
            eqs_mov, x0, t0, tmax, _order, _abstol))
        @test length(tv) < 501
        @test length(xv) < 501
        @test length(tv) == 28
        @test length(xv) == 28
        @test xv[1] == x0
        @test tv[end] < 1/3
        @test tv[end] == tmax
        @test abs(xv[end]-exactsol(tv[end], xv[1])) < 1e-11
    end

    @testset "Test non-autonomous ODE (1): dot{x}=cos(t)" begin
        let f(x, p, t) = cos(t)
            t0 = 0//1
            tmax = 10.25*(2pi)
            abstol = 1e-20
            order = 25
            x0 = 0.0 #initial conditions such that x(t)=sin(t)

            tv, xv = (@test_logs min_level=Logging.Warn taylorinteg(
                f, x0, t0, tmax, order, abstol))
            @test length(tv) < 501
            @test length(xv) < 501
            @test xv[1] == x0
            @test tv[1] == t0
            @test abs(sin(tmax)-xv[end]) < 1e-14

            # Backward integration
            tb, xb = (@test_logs min_level=Logging.Warn taylorinteg(
                f, sin(tmax), tmax, t0, order, abstol))
            @test length(tb) < 501
            @test length(xb) < 501
            @test xb[1] == sin(tmax)
            @test tb[1] > tb[end]
            @test abs(sin(t0)-xb[end]) < 5e-14

            # Tests with a range, for comparison with backward integration
            tmax = 15*(2pi)
            Δt = (tmax-t0)/1024
            tspan = t0:Δt:tmax
            xv = (@test_logs min_level=Logging.Warn taylorinteg(
                f, x0, tspan, order, abstol))
            @test xv[1] == x0
            @test abs(sin(tmax)-xv[end]) < 1e-14

            # Backward integration
            xback = (@test_logs min_level=Logging.Warn taylorinteg(
                f, xv[end], reverse(tspan), order, abstol))
            @test xback[1] == xv[end]
            @test abs(sin(t0)-xback[end]) < 5e-14
            @test norm(xv[:]-xback[end:-1:1], Inf) < 5.0e-14

            # Tests if trange is properly sorted
            @test_throws AssertionError taylorinteg(f, x0, rand(t0:Δt:tmax, 100), order, abstol)
        end
    end

end

ext: jl


text:
# diffl.jl

using MIRT: diffl, diffl!, diffl_adj, diffl_adj!, diffl_map

using LinearMapsAA: LinearMapAM, LinearMapAO
using Test: @test, @testset, @test_throws, @inferred

@testset "1D" begin
	x = rand(4)
	@test diffl(x)[2:end] == diff(x)
end

@testset "2D" begin
	x = [2 4; 6 16]
	g = @inferred diffl(x, 1)
	@test all(g[1,:] .== 0)
	@test g[2:end,:] == diff(x, dims=1)
	g = @inferred diffl(x, 2)
	@test all(g[:,1] .== 0)
	@test g[:,2:end] == diff(x, dims=2)

	g = @inferred diffl(x, 1 ; edge=:none)
	@test g[2:end,:] == diff(x, dims=1)

	g = @inferred diffl(x, 1 ; add=true)
	@test g[2:end,:][:] == x[1,:] + x[2,:]

	x = rand(3,4)
	g1 = diffl(x)
	g2 = similar(g1)
	@inferred diffl!(g2, x) # test default 1
	@test g2[2:end,:] == g1[2:end,:]
end

@testset "stack" begin
	x = reshape((1:(2*3*4)).^2, 2, 3, 4)
	g = diffl(x, [3, 1])
	@test all(g[:,:,1,1] .== 0)
	@test all(g[1,:,:,2] .== 0)
	@test g[:,:,2:end,1] == diff(x, dims=3)
	@test g[2:end,:,:,2] == diff(x, dims=1)
end

@testset "adj" begin
	x = rand(3)
	@inferred diffl_adj(rand(3))
	@test_throws ArgumentError diffl_adj(rand(3) ; edge=:test)
#	@inferred diffl_adj(rand(3,4,2), 1:2) # fails
	@test size(diffl_adj(rand(3,4,2), 1:2)) == (3,4)
end

N = (2,3); d = 2

@testset "basics" begin
	O = diffl_map(N, d ; T=Int32, edge=:zero, add=false)
	@test O isa LinearMapAO
	@test Matrix(O)' == Matrix(O')
	@test O.name == "diffl_map"
	@test_throws String diffl_map(N ; edge=:none) # unsupported
end

@testset "adjoint" begin
	for N in [(3,), (3,4), (2,3,4)]
		dlist = [1, [1,]]
		length(N) > 1 && push!(dlist, 2:-1:1, 1:length(N))
		length(N) > 2 && push!(dlist, [length(N), 1])
		isinteractive() && @show dlist
		for d in dlist
			for edge in (:zero, :circ)
				for add in (false, true)
					for op in (false, true)
						O = diffl_map(N, d ;
							T=Int32, edge=edge, add=add, operator = op)
						@test Matrix(O)' == Matrix(O')
						@test O isa (op ? LinearMapAO : LinearMapAM)
					end
				end
			end
		end
	end
end


text:
using BigCombinatorics
using Test

@test Fibonacci(10)==Fibonacci(9)+Fibonacci(8)

@test Factorial(10)==big(factorial(10))

@test FallingFactorial(10,3) == 10*9*8
@test FallingFactorial(10,10) == Factorial(10)
@test FallingFactorial(10,12) == 0

@test RisingFactorial(10,3) == 10*11*12

@test DoubleFactorial(9) == 9*7*5*3
@test DoubleFactorial(10) == 10*8*6*4*2

@test Catalan(12) == 208012

@test Derangements(1) == 0
@test Derangements(12) == 176214841

@test MultiChoose(10,1) == 10
@test MultiChoose(10,0) == 1
@test MultiChoose(10,10) == 92378

@test Multinomial(5,5,5) == div(Factorial(15), Factorial(5)^3)
@test Multinomial([6,6,6,6]) == Multinomial(6,6,6,6)

@test Bell(10)==115975
@test sum(Stirling2(10,k) for k=0:10) == Bell(10)
@test Stirling1(10,10) == 1
@test Stirling1(10,0) == 0
@test sum(Stirling1(10,k) for k=0:10) == 0

@test IntPartitions(10) == 42
@test IntPartitionsDistinct(10) == 10

@test Euler(12) == 2702765

@test PowerSum(10,3) == sum(k^3 for k=1:10)

@test BigCombinatorics.cache_clear(Fibonacci)

@test sum(Eulerian(10,k) for k=1:10) == Factorial(10)


text:
using PtFEM, Test

data = Dict(
  # Beam(nels, nn, ndim, nst, nip, finite_element(nod, nodof))
  :struc_el => Beam(2, 1, 4, 1, :x, Line(2, 2), false),
  :properties => [4.0e4; 2.0e4],
  :etype => [1, 1, 2, 2],
  :x_coords => [0.0, 2.5, 5.0, 8.0, 10.0],
  :loaded_nodes => [
  (2, [-20.0 0.0]),
  (3, [-6.0 -3.0]),
  (4, [-8.8 2.2]),
  (5, [-1.2 0.5333])],
  :fixed_freedoms => [
    (1, 1, 0.0),
    (1, 2, -0.001),
    (3, 1, -0.005),
    (4, 1, 0.0)],
  :penalty => 1e19
)

@time fem, dis_df, fm_df = p43(data)

@test round.(fem.displacements[:, 5], digits=7) == [0.0047132, 0.0023433]


text:
## generic `Tables.rows` and `Tables.columns` fallbacks
## if a table provides Tables.rows or Tables.columns,
## we'll provide a default implementation of the other

# Turn any AbstractColumns into an AbstractRow iterator

# get the number of rows in the incoming table
function rowcount(cols)
    names = columnnames(cols)
    isempty(names) && return 0
    return length(getcolumn(cols, names[1]))
end

# a lazy row view into a AbstractColumns object
struct ColumnsRow{T} <: AbstractRow
    columns::T # an `AbstractColumns`-compatible object
    row::Int # row number
end

getcolumns(c::ColumnsRow) = getfield(c, :columns)
getrow(c::ColumnsRow) = getfield(c, :row)

# AbstractRow interface
Base.@propagate_inbounds getcolumn(c::ColumnsRow, ::Type{T}, col::Int, nm::Symbol) where {T} = getcolumn(getcolumns(c), T, col, nm)[getrow(c)]
Base.@propagate_inbounds getcolumn(c::ColumnsRow, i::Int) = getcolumn(getcolumns(c), i)[getrow(c)]
Base.@propagate_inbounds getcolumn(c::ColumnsRow, nm::Symbol) = getcolumn(getcolumns(c), nm)[getrow(c)]
columnnames(c::ColumnsRow) = columnnames(getcolumns(c))

@generated function Base.isless(c::ColumnsRow{T}, d::ColumnsRow{T}) where {T <: NamedTuple{names}} where names
    exprs = Expr[]
    for n in names
        var1 = Expr(:., :c, QuoteNode(n))
        var2 = Expr(:., :d, QuoteNode(n))
        bl = quote
            a, b = $var1, $var2
            isless(a, b) && return true
            isequal(a, b) || return false
        end
        push!(exprs, bl)
    end
    push!(exprs, :(return false))
    Expr(:block, exprs...)
end

@generated function Base.isequal(c::ColumnsRow{T}, d::ColumnsRow{T}) where {T <: NamedTuple{names}} where names
    exprs = Expr[]
    for n in names
        var1 = Expr(:., :c, QuoteNode(n))
        var2 = Expr(:., :d, QuoteNode(n))
        push!(exprs, :(isequal($var1, $var2) || return false))
    end
    push!(exprs, :(return true))
    Expr(:block, exprs...)
end

# RowIterator wraps an AbstractColumns object and provides row iteration via lazy row views
struct RowIterator{T}
    columns::T
    len::Int
end

Base.eltype(x::RowIterator{T}) where {T} = ColumnsRow{T}
Base.length(x::RowIterator) = getfield(x, :len)
Base.getproperty(x::RowIterator, nm::Symbol) = getcolumn(x, nm)
Base.getproperty(x::RowIterator, i::Int) = getcolumn(x, i)
Base.propertynames(x::RowIterator) = columnnames(x)
isrowtable(::Type{<:RowIterator}) = true

columnaccess(::Type{<:RowIterator}) = true
columns(x::RowIterator) = getfield(x, :columns)
columnnames(x::RowIterator) = columnnames(columns(x))
getcolumn(x::RowIterator, nm::Symbol) = getcolumn(columns(x), nm)
getcolumn(x::RowIterator, i::Int) = getcolumn(columns(x), i)
materializer(x::RowIterator) = materializer(columns(x))
schema(x::RowIterator) = schema(columns(x))

@inline function Base.iterate(rows::RowIterator, st=1)
    st > length(rows) && return nothing
    return ColumnsRow(columns(rows), st), st + 1
end

# this is our generic Tables.rows fallback definition
function rows(x::T) where {T}
    isrowtable(x) && return x
    # because this method is being called, we know `x` didn't define it's own Tables.rows
    # first check if it supports column access, and if so, wrap it in a RowIterator
    if columnaccess(T)
        cols = columns(x)
        return RowIterator(cols, Int(rowcount(cols)))
    # otherwise, if the input is at least iterable, we'll wrap it in an IteratorWrapper
    # which will iterate the input, validating that elements support the AbstractRow interface
    # and unwrapping any DataValues that are encountered
    elseif IteratorInterfaceExtensions.isiterable(x)
        return nondatavaluerows(x)
    end
    throw(ArgumentError("no default `Tables.rows` implementation for type: $T"))
end

# for AbstractRow iterators, we define a "collect"-like routine to build up columns from iterated rows

"""
    Tables.allocatecolumn(::Type{T}, len) => returns a column type (usually `AbstractVector`) with size to hold `len` elements

Custom column types can override with an appropriate "scalar" element type that should dispatch to their column allocator.
Alternatively, and more generally, custom scalars can overload `DataAPI.defaultarray` to signal the default array type.
"""
allocatecolumn(T, len) = DataAPI.defaultarray(T, 1)(undef, len)

@inline function allocatecolumns(::Schema{names, types}, len) where {names, types}
    if @generated
        vals = Tuple(:(allocatecolumn($(fieldtype(types, i)), len)) for i = 1:fieldcount(types))
        return :(NamedTuple{Base.map(Symbol, names)}(($(vals...),)))
    else
        return NamedTuple{Base.map(Symbol, names)}(Tuple(allocatecolumn(fieldtype(types, i), len) for i = 1:fieldcount(types)))
    end
end

# add! will push! or setindex! a value depending on if the row-iterator HasLength or not
@inline add!(val, col::Int, nm, dest::AbstractArray, ::Union{Base.HasLength, Base.HasShape}, row) = setindex!(dest, val, row)
@inline add!(val, col::Int, nm, dest::AbstractArray, L, row) = push!(dest, val)

@inline function buildcolumns(schema, rowitr::T) where {T}
    L = Base.IteratorSize(T)
    len = Base.haslength(T) ? length(rowitr) : 0
    nt = allocatecolumns(schema, len)
    for (i, row) in enumerate(rowitr)
        eachcolumns(add!, schema, row, nt, L, i)
    end
    return nt
end

@inline add!(dest::AbstractArray, val, ::Union{Base.HasLength, Base.HasShape}, row) = setindex!(dest, val, row)
@inline add!(dest::AbstractArray, val, T, row) = push!(dest, val)

replacex(t, col::Int, x) = ntuple(i->i == col ? x : t[i], length(t))

@inline function add_or_widen!(val, col::Int, nm, dest::AbstractArray{T}, row, updated, L) where {T}
    if val isa T
        add!(dest, val, L, row)
        return
    else
        new = allocatecolumn(promote_type(T, typeof(val)), length(dest))
        row > 1 && copyto!(new, 1, dest, 1, row - 1)
        add!(new, val, L, row)
        updated[] = replacex(updated[], col, new)
        return
    end
end

function __buildcolumns(rowitr, st, sch, columns, rownbr, updated)
    while true
        state = iterate(rowitr, st)
        state === nothing && break
        row, st = state
        rownbr += 1
        eachcolumns(add_or_widen!, sch, row, columns, rownbr, updated, Base.IteratorSize(rowitr))
        # little explanation here: we just called add_or_widen! for each column value of our row
        # note that when a column's type is widened, `updated` is set w/ the new set of columns
        # we then check if our current `columns` isn't the same object as our `updated` ref
        # if it isn't, we're going to call __buildcolumns again, passing our new updated ref as
        # columns, which allows __buildcolumns to specialize (i.e. recompile) based on the new types
        # of updated. So a new __buildcolumns will be compiled for each widening event.
        columns !== updated[] && return __buildcolumns(rowitr, st, sch, updated[], rownbr, updated)
    end
    return updated
end

# for the schema-less case, we do one extra step of initializing each column as an `EmptyVector`
# and doing an initial widening for each column in _buildcolumns, before passing the widened
# set of columns on to __buildcolumns
struct EmptyVector <: AbstractVector{Union{}}
    len::Int
end
Base.IndexStyle(::Type{EmptyVector}) = Base.IndexLinear()
Base.size(x::EmptyVector) = (x.len,)
Base.getindex(x::EmptyVector, i::Int) = throw(UndefRefError())

function _buildcolumns(rowitr, row, st, sch, columns, updated)
    eachcolumns(add_or_widen!, sch, row, columns, 1, updated, Base.IteratorSize(rowitr))
    return __buildcolumns(rowitr, st, sch, updated[], 1, updated)
end

# when Tables.schema(x) === nothing
@inline function buildcolumns(::Nothing, rowitr::T) where {T}
    state = iterate(rowitr)
    state === nothing && return NamedTuple()
    row, st = state
    names = Tuple(columnnames(row))
    len = Base.haslength(T) ? length(rowitr) : 0
    sch = Schema(names, nothing)
    columns = Tuple(EmptyVector(len) for _ = 1:length(names))
    return NamedTuple{Base.map(Symbol, names)}(_buildcolumns(rowitr, row, st, sch, columns, Ref{Any}(columns))[])
end

"""
    Tables.CopiedColumns

For some sinks, there's a concern about whether they can safely "own" columns from the input.
If mutation will be allowed, to be safe, they should always copy input columns, to avoid unintended mutation
to the original source.
When we've called `buildcolumns`, however, Tables.jl essentially built/owns the columns,
and it's happy to pass ownership to the sink. Thus, any built columns will be wrapped
in a `CopiedColumns` struct to signal to the sink that essentially "a copy has already been made"
and they're safe to assume ownership.
"""
struct CopiedColumns{T} <: AbstractColumns
    x::T
end

source(x::CopiedColumns) = getfield(x, :x)
istable(::Type{<:CopiedColumns}) = true
columnaccess(::Type{<:CopiedColumns}) = true
columns(x::CopiedColumns) = x
schema(x::CopiedColumns) = schema(source(x))
materializer(x::CopiedColumns) = materializer(source(x))

getcolumn(x::CopiedColumns, ::Type{T}, col::Int, nm::Symbol) where {T} = getcolumn(source(x), T, col, nm)
getcolumn(x::CopiedColumns, i::Int) = getcolumn(source(x), i)
getcolumn(x::CopiedColumns, nm::Symbol) = getcolumn(source(x), nm)
columnnames(x::CopiedColumns) = columnnames(source(x))

# here's our generic fallback Tables.columns definition
@inline function columns(x::T) where {T}
    # because this method is being called, we know `x` didn't define it's own Tables.columns method
    # first check if it explicitly supports row access, and if so, build up the desired columns
    if rowaccess(T)
        r = rows(x)
        return CopiedColumns(buildcolumns(schema(r), r))
    # though not widely supported, if a source supports the TableTraits column interface, use it
    elseif TableTraits.supports_get_columns_copy_using_missing(x)
        return CopiedColumns(TableTraits.get_columns_copy_using_missing(x))
    # otherwise, if the source is at least iterable, we'll wrap it in an IteratorWrapper and
    # build columns from that, which will check if the source correctly iterates valid AbstractRow objects
    # and unwraps DataValues for us
    elseif IteratorInterfaceExtensions.isiterable(x)
        iw = nondatavaluerows(x)
        return CopiedColumns(buildcolumns(schema(iw), iw))
    end
    throw(ArgumentError("no default `Tables.columns` implementation for type: $T"))
end


text:
function variable_dcgrid_voltage_magnitude(pm::AbstractWModels; kwargs...)
    variable_dcgrid_voltage_magnitude_sqr(pm; kwargs...)
end
"""
```
sum(p[a] for a in bus_arcs) + sum(pconv_grid_ac[c] for c in bus_convs_ac) == sum(pg[g] for g in bus_gens)  - pd - gs*w
sum(q[a] for a in bus_arcs) + sum(qconv_grid_ac[c] for c in bus_convs_ac) == sum(qg[g] for g in bus_gens)  - qd + bs*w
```
"""
function constraint_kcl_shunt(pm::AbstractWModels, n::Int, cnd::Int, i::Int, bus_arcs, bus_arcs_dc, bus_gens, bus_convs_ac, bus_loads, bus_shunts, pd, qd, gs, bs)
    w = PowerModels.var(pm, n, cnd, :w, i)
    p = PowerModels.var(pm, n, cnd, :p)
    q = PowerModels.var(pm, n, cnd, :q)
    pg = PowerModels.var(pm, n, cnd, :pg)
    qg = PowerModels.var(pm, n, cnd, :qg)
    pconv_grid_ac = PowerModels.var(pm, n, cnd, :pconv_tf_fr)
    qconv_grid_ac = PowerModels.var(pm, n, cnd, :qconv_tf_fr)

    PowerModels.con(pm, n, cnd, :kcl_p)[i] = @constraint(pm.model, sum(p[a] for a in bus_arcs) + sum(pconv_grid_ac[c] for c in bus_convs_ac)  == sum(pg[g] for g in bus_gens)  - sum(pd[d] for d in bus_loads) - sum(gs[s] for s in bus_shunts)*w)
    PowerModels.con(pm, n, cnd, :kcl_q)[i] = @constraint(pm.model, sum(q[a] for a in bus_arcs) + sum(qconv_grid_ac[c] for c in bus_convs_ac)  == sum(qg[g] for g in bus_gens)  - sum(qd[d] for d in bus_loads) + sum(bs[s] for s in bus_shunts)*w)
end
"""
Creates Ohms constraints for DC branches

```
p[f_idx] + p[t_idx] == p * g[l] * (wdc[f_bus] - wdcr[f_bus,t_bus])
```
"""
function constraint_ohms_dc_branch(pm::AbstractWRModels, n::Int, cnd::Int, f_bus, t_bus, f_idx, t_idx, r, p)
    p_dc_fr = PowerModels.var(pm, n, cnd, :p_dcgrid, f_idx)
    p_dc_to = PowerModels.var(pm, n, cnd, :p_dcgrid, t_idx)
    wdc_fr = PowerModels.var(pm, n, cnd, :wdc, f_bus)
    wdc_to = PowerModels.var(pm, n, cnd, :wdc, t_bus)
    wdc_frto = PowerModels.var(pm, n, cnd, :wdcr, (f_bus, t_bus))

    if r == 0
        @constraint(pm.model, p_dc_fr + p_dc_to == 0)
    else
        g = 1 / r
        @constraint(pm.model, p_dc_fr == p * g *  (wdc_fr - wdc_frto))
        @constraint(pm.model, p_dc_to == p * g *  (wdc_to - wdc_frto))
    end
end
"`wdc[i] == vdcm^2`"
function constraint_dc_voltage_magnitude_setpoint(pm::AbstractWModels, n::Int, cnd::Int, i, vdcm)
    wdc = PowerModels.var(pm, n, cnd, :wdc, i)

    PowerModels.con(pm, n, cnd, :v_dc)[i] = @constraint(pm.model, wdc == vdcm^2)
end

function add_dc_bus_voltage_setpoint(sol, pm::AbstractWModels)
    PowerModels.add_setpoint!(sol, pm, "busdc", "vm", :wdc, status_name="Vdc", inactive_status_value = 4; scale = (x,item,cnd) -> sqrt(x))
end

"""
Limits dc branch current

```
p[f_idx] <= wdc[f_bus] * Imax
```
"""
function constraint_dc_branch_current(pm::AbstractWModels, n::Int, cnd::Int, f_bus, f_idx, ccm_max, p)
    p_dc_fr = PowerModels.var(pm, n, cnd, :p_dcgrid, f_idx)
    wdc_fr = PowerModels.var(pm, n, cnd, :wdc, f_bus)

    @constraint(pm.model, p_dc_fr <= wdc_fr * ccm_max * p^2)
end


text:
#=
Version = 0.03
Author = "William Herrera"
Copyright = "Copyright: 2018 William Herrera"
Created = "24 Jan 2018"
Purpose = "EEG file OpenBCI Ganglion view data example"
=#

using OpenBCI
using EDFPlus
using Plots
pyplot()

boardIP = "192.168.1.23"
myIP = "192.168.1.1"
idfile = "../test/patientdata.json"

const PLOTINTERVAL = 10

linspace(start, stop, len) = LinRange{Float64}(start, stop, len)

function plottwobipolars(bdfh, pcount, maxpackets)
    if pcount > PLOTINTERVAL && pcount % PLOTINTERVAL == 1
        c1data = bdfh.BDFsignals[pcount-PLOTINTERVAL:pcount, 1:250][:]
        c2data = bdfh.BDFsignals[pcount-PLOTINTERVAL:pcount, 251:500][:]
        c3data = bdfh.BDFsignals[pcount-PLOTINTERVAL:pcount, 501:750][:]
        c4data = bdfh.BDFsignals[pcount-PLOTINTERVAL:pcount, 751:1000][:]
        Fp1T3 = c2data .- c1data
        Fp2T4 = c4data .- c3data
        Fp1T3 = EDFPlus.lowpassfilter(Fp1T3, 250, 40.0)
        Fp2T4 = EDFPlus.lowpassfilter(Fp2T4, 250, 40.0)
        Fp1T3 = EDFPlus.highpassfilter(Fp1T3, 250, 0.5)
        Fp2T4 = EDFPlus.highpassfilter(Fp2T4, 250, 0.5)
        ydata = [Fp1T3, Fp2T4]
        timepoints = linspace(0.0, PLOTINTERVAL, length(Fp1T3))

        # detach plotting, return quickly now so as to avoid dropped packets
        @async(begin
        plt = Plots.plot(timepoints, ydata, layout=(2,1),
                   xticks=collect(timepoints[1]:1:timepoints[end]),
                   yticks=false, legend=false, title="Interval from $(pcount-PLOTINTERVAL) to $pcount")
        Plots.plot!(yaxis=true, xaxis=false, tight_layout=true, ylabel = "Fp1-T3", subplot=1)
        Plots.plot!(yaxis=true, tight_layout=true, ylabel = "Fp2-T4", subplot=2)
        PyPlot.display(plt)
        end)
    end
end

makeganglionbdfplus("examplefile.bdf", boardIP, myIP, 120,
                             idfile=idfile, inspector=plottwobipolars)


text:
using Compat: @compat
using QuickTypes
using QuickTypes: construct, roottypeof, fieldsof, type_parameters, roottype,
      tuple_parameters
using Compat.Test
using ConstructionBase: setproperties

@compat abstract type Vehicle end

@qstruct Car{T<:Number, U}(size::T, nwheels::Int=4; manufacturer::U=nothing,
                              brand::String="off-brand") <: Vehicle

c = Car(10; manufacturer=("Danone", "Hershey"))

@test c.size==10
@test c.nwheels==4
@test c.manufacturer==("Danone", "Hershey")
@test c.brand=="off-brand"
c2 = @inferred setproperties(c, (size=42, nwheels=8))
@test c2.nwheels == 8
@test c2.size == 42
@test c2.brand == c.brand
# Check that the fields are in the right order
@test collect(fieldnames(Car)) == [:size, :nwheels, :manufacturer, :brand]
# This is essentially the definition of these functions.
@test construct(roottypeof(c), fieldsof(c)...) == c
@test type_parameters(Vector{Int}) == Base.Core.svec(Int64, 1)
@test tuple_parameters(Tuple{Int, Float64}) == Base.Core.svec(Int64, Float64)
@inferred roottypeof(1=>2) == Pair

################################################################################

@qstruct Empty()
Empty()
@test setproperties(Empty(), NamedTuple()) === Empty()

# Used to yield:
#     WARNING: static parameter T does not occur in signature for Type.
#     The method will not be callable.
@qstruct Blah{T}()

################################################################################

@qstruct Boring(x::Int)
@inferred Boring(10)
@test Boring(10).x == 10
@test Boring(10.0).x == 10   # check that convert is called correctly
@qstruct ParametricBoring{X}(x::X; _concise_show=true)
@inferred ParametricBoring(10)
@test ParametricBoring(10).x === 10
o = ParametricBoring(1)
@test setproperties(o, x=:one).x === :one

@qstruct Kwaroo(x; y=10)
@test Kwaroo(5) == Kwaroo(5; y=10)
o = Kwaroo(5, y=10)
o2 = @inferred setproperties(o, (x=:five, y=100.0))
@test o2 isa Kwaroo
@test o2.x === :five
@test o2.y === 100.0

################################################################################
# Slurping

@qstruct Slurp(x, y=1, args...; kwargs...)
s = Slurp(1,2,3,4,5,6,7; x=1, y=10+2)
@test s.args == (3,4,5,6,7)
@test s.kwargs == [(:x => 1), (:y => 12)]
s2 = @inferred setproperties(s, x=:hello)
@test s2 isa Slurp
@test s2.x == :hello
@test s2.y == s.y

################################################################################

@qmutable Foo{T}(x::T; y=2) do
    @assert x < 10
end

@test_throws AssertionError Foo(11; y=10.0)
@test_throws AssertionError construct(Foo, 11, 10.0)

################################################################################
# Fully-parametric

@qstruct_fp Plane(nwheels::Number; brand=:zoomba) do
    @assert nwheels < 100
end <: Vehicle

@test_throws MethodError Plane{Int, Symbol}(2; brand=12)
@test Plane{Int, Symbol}(2; brand=:zoomba).brand == :zoomba
@test supertype(Plane) == Vehicle
# This used to be a MethodError, but since we moved the outer constructor inside
# the type, it has become a TypeError. Not sure why!
@test_throws TypeError Plane("happy")

@qstruct_fp NoFields()   # was an error before it was special-cased

o = Plane(4)
o2 = @inferred setproperties(o, brand=10, nwheels=o.nwheels)
@test o2 isa Plane
@test o2.brand === 10
@test o2.nwheels === o.nwheels

################################################################################
# Narrowly-parametric

@qstruct_fp Foo_fp(a, b)
@qstruct_np Foo_np(a, b)
convert_f(foo) = convert(foo.a, 10)
@test_throws(Exception, @inferred convert_f(Foo_fp(Int, 2)))
@inferred convert_f(Foo_np(Int, 2))
@test fieldtype(typeof(Foo_np(Int, 2)), :a) == Type{Int64}

@qstruct Issue11(;no_default_value)
@test_throws UndefKeywordError Issue11()


text:
module QuantileRegression

    import StatsModels: TableRegressionModel, FormulaTerm, coef, @formula
    export qreg, coef, vcov, stderr, quantiles, IP, IRLS, @formula

    using DataFrames, Distributions, LinearAlgebra, LinearAlgebra.BLAS, StatsModels, StatsBase

    mutable struct QRegModel
        beta::Vector{Float64}
        vcov::Matrix{Float64}
        stderr::Vector{Float64}
        q
    end

    abstract type Solver end

    struct IP <: Solver
    end

    struct IRLS <: Solver
        tol
        maxIter
        threshold
    end
    IRLS(;tol::Real = 1e-12, maxIter::Integer = 1_000, threshold::Real = 1e-5) = IRLS(tol, maxIter, threshold)

    include("InteriorPoint.jl")
    include("IRLS.jl")
    include("Covariance.jl")

    function qreg(f::FormulaTerm, df::AbstractDataFrame, q, s::Solver = IP();beta=nothing)
        mf = ModelFrame(f, df)
        mm = ModelMatrix(mf)
        mr = response(mf)
        if beta === nothing
            beta = qreg_coef(mr, mm.m, q, s)
        end
        vcov = qreg_vcov(mr, mm.m, beta, q)
        stderror = sqrt.(diag(vcov))
        return TableRegressionModel(QRegModel(beta, vcov, stderror, q), mf, mm)
    end
    qreg(f::FormulaTerm, df::AbstractDataFrame, s::Solver) = qreg(f, df, 0.5, s)

    StatsBase.coef(x::QRegModel) = x.beta
    StatsBase.coef(rm::TableRegressionModel) = coef(rm.model)

    StatsBase.vcov(x::QRegModel) = x.vcov

    StatsBase.stderror(x::QRegModel) = x.stderr

    quantiles(x::QRegModel) = x.q
    quantiles(rm::TableRegressionModel) = quantiles(rm.model)

    function StatsBase.coeftable(mm::QRegModel)
        cc = coef(mm)
        se = stderror(mm)
        tt = cc./se
        CoefTable(hcat(kron(mm.q,ones(length(se))), cc,se,tt),
                 ["Quantile",  "Estimate","Std.Error","t value"],
                 ["x$i" for i = 1:length(cc)])
    end
end


text:
"""
    vantablack()

Vantablack spectrum from 264.6 to 16310.0 nm
"""
function vantablack()
    RSpec(vantablack[:,1],vantablack[:,2])
end

#= Vantablack hemispherical reflectance spectrum from 264.6 to 16310.0 nm

    sources:

    https://vantablack.co.jp/wp-content/uploads/2020/04/webdata-full-s-vis-v2-2020.pdf

    https://www.surreynanosystems.com/assets/media/vantablack-vb-a4-data-brochure-2016.pdf

    https://www.surreynanosystems.com/assets/media/vantablack-vbx-2-a4-data-sheet-008-downloadpdf.pdf

=#

const vantablack =
       [264.6    0.00591;
        273.8    0.00594;
        284.6    0.00594;
        294.9    0.00593;
        300.1    0.00591;
        311.3    0.00585;
        322.6    0.00573;
        336.4    0.00552;
        350.1    0.00529;
        362.7    0.00507;
        378.0    0.00487;
        392.8    0.00473;
        400.0    0.00466;
        407.7    0.00460;
        420.0    0.00448;
        433.9    0.00433;
        450.1    0.00414;
        470.8    0.00388;
        486.9    0.00371;
        500.0    0.00359;
        526.7    0.00328;
        560.2    0.00288;
        576.9    0.00269;
        590.1    0.00258;
        599.9    0.00252;
        628.0    0.00238;
        668.4    0.00222;
        699.5    0.00211;
        750.9    0.00197;
        799.8    0.00185;
        822.3    0.00183;
        841.5    0.00185;
        853.5    0.00190;
        867.8    0.00197;
        882.1    0.00200;
        900.1    0.00195;
        930.8    0.00186;
        966.9    0.00177;
        1000.5   0.00171;
        1031.4   0.00168;
        1070.5   0.00166;
        1100.8   0.00166;
        1200.4   0.00166;
        1228.9   0.00166;
        1266.7   0.00169;
        1300.9   0.00168;
        1339.0   0.00171;
        1400.0   0.00180;
        4350.0   0.00191;
        4510.0   0.00199;
        4640.0   0.00206;
        4890.0   0.00233;
        5030.0   0.00238;
        5520.0   0.00271;
        5810.0   0.00286;
        6010.0   0.00308;
        6220.0   0.00330;
        6380.0   0.00337;
        6670.0   0.00363;
        7020.0   0.00396;
        7240.0   0.00399;
        7480.0   0.00405;
        7700.0   0.00418;
        8040.0   0.00452;
        8300.0   0.00475;
        8560.0   0.00512;
        8800.0   0.00545;
        9030.0   0.00578;
        9250.0   0.00600;
        9470.0   0.00609;
        9790.0   0.00626;
        10120.0  0.00643;
        10390.0  0.00652;
        10880.0  0.00672;
        11570.0  0.00676;
        12130.0  0.00702;
        12530.0  0.00724;
        13010.0  0.00754;
        13610.0  0.00787;
        14440.0  0.00839;
        15010.0  0.00876;
        15670.0  0.00915;
        15970.0  0.00930;
        16310.0  0.00944]


text:
module TestConstructors
    using Compat, Compat.Test, DataFrames
    using DataFrames: Index
    const ≅ = isequal

    #
    # DataFrame
    #

    df = DataFrame()
    @test df.columns == Any[]
    @test df.colindex == Index()

    df = DataFrame(Any[CategoricalVector{Union{Float64, Missing}}(zeros(3)),
                       CategoricalVector{Union{Float64, Missing}}(ones(3))],
                   Index([:x1, :x2]))
    @test size(df, 1) == 3
    @test size(df, 2) == 2

    @test df == DataFrame(Any[CategoricalVector{Union{Float64, Missing}}(zeros(3)),
                              CategoricalVector{Union{Float64, Missing}}(ones(3))])
    @test df == DataFrame(x1 = Union{Int, Missing}[0.0, 0.0, 0.0],
                          x2 = Union{Int, Missing}[1.0, 1.0, 1.0])

    df2 = convert(DataFrame, Union{Float64, Missing}[0.0 1.0;
                                                  0.0 1.0;
                                                  0.0 1.0])
    names!(df2, [:x1, :x2])
    @test df[:x1] == df2[:x1]
    @test df[:x2] == df2[:x2]

    df2 = DataFrame([0.0 1.0;
                     0.0 1.0;
                     0.0 1.0])
    names!(df2, [:x1, :x2])
    @test df[:x1] == df2[:x1]
    @test df[:x2] == df2[:x2]

    df2 = DataFrame([0.0 1.0;
                     0.0 1.0;
                     0.0 1.0], [:a, :b])
    names!(df2, [:a, :b])
    @test df[:x1] == df2[:a]
    @test df[:x2] == df2[:b]

    @test df == DataFrame(x1 = Union{Float64, Missing}[0.0, 0.0, 0.0],
                          x2 = Union{Float64, Missing}[1.0, 1.0, 1.0])
    @test df == DataFrame(x1 = Union{Float64, Missing}[0.0, 0.0, 0.0],
                          x2 = Union{Float64, Missing}[1.0, 1.0, 1.0],
                          x3 = Union{Float64, Missing}[2.0, 2.0, 2.0])[[:x1, :x2]]

    df = DataFrame(Union{Int, Missing}, 2, 2)
    @test size(df) == (2, 2)
    @test eltypes(df) == [Union{Int, Missing}, Union{Int, Missing}]

    df = DataFrame([Union{Int, Missing}, Union{Float64, Missing}], [:x1, :x2], 2)
    @test size(df) == (2, 2)
    @test eltypes(df) == [Union{Int, Missing}, Union{Float64, Missing}]

    @test df ≅ DataFrame([Union{Int, Missing}, Union{Float64, Missing}], 2)

    @test_throws BoundsError SubDataFrame(DataFrame(A=1), 0)
    @test_throws BoundsError SubDataFrame(DataFrame(A=1), 0)
    @test SubDataFrame(DataFrame(A=1), 1) == DataFrame(A=1)
    @test SubDataFrame(DataFrame(A=1:10), 1:4) == DataFrame(A=1:4)
    @test view(SubDataFrame(DataFrame(A=1:10), 1:4), 2) == DataFrame(A=2)
    @test view(SubDataFrame(DataFrame(A=1:10), 1:4), [true, true, false, false]) == DataFrame(A=1:2)

    @test DataFrame(a=1, b=1:2) == DataFrame(a=[1,1], b=[1,2])

    @testset "pair constructor" begin
        df = DataFrame(:x1 => zeros(3), :x2 => ones(3))
        @test size(df, 1) == 3
        @test size(df, 2) == 2
        @test isequal(df, DataFrame(x1 = [0.0, 0.0, 0.0], x2 = [1.0, 1.0, 1.0]))

        df = DataFrame(:type => [], :begin => [])
        @test names(df) == [:type, :begin]
    end

    @testset "associative" begin
        df = DataFrame(Dict(:A => 1:3, :B => 4:6))
        @test df == DataFrame(A = 1:3, B = 4:6)
        @test eltypes(df) == [Int, Int]
    end

    @testset "recyclers" begin
        @test DataFrame(a = 1:5, b = 1) == DataFrame(a = collect(1:5), b = fill(1, 5))
        @test DataFrame(a = 1, b = 1:5) == DataFrame(a = fill(1, 5), b = collect(1:5))
    end

    @testset "constructor errors" begin
        @test_throws DimensionMismatch DataFrame(a=1, b=[])
        @test_throws DimensionMismatch DataFrame(Any[collect(1:10)], DataFrames.Index([:A, :B]))
        @test_throws DimensionMismatch DataFrame(A = rand(2,2))
        @test_throws DimensionMismatch DataFrame(A = rand(2,1))
    end

    @testset "column types" begin
        df = DataFrame(A = 1:3, B = 2:4, C = 3:5)
        answer = [Array{Int,1}, Array{Int,1}, Array{Int,1}]
        @test map(typeof, df.columns) == answer
        df[:D] = [4, 5, missing]
        push!(answer, Vector{Union{Int, Missing}})
        @test map(typeof, df.columns) == answer
        df[:E] = 'c'
        push!(answer, Vector{Char})
        @test map(typeof, df.columns) == answer
    end

    @testset "categorical constructor" begin
        df = DataFrame([Int, String], [:a, :b], [false, true], 3)
        @test !(df[:a] isa CategoricalVector)
        @test df[:b] isa CategoricalVector
        @test_throws DimensionMismatch DataFrame([Int, String], [:a, :b], [true], 3)
    end
end


text:
"""
    S_score_perms(data, nPerms; permFun, isVarFloor)

Obtains permutation p-values for S scores as described by Collins et al.

# Arguments

- data = RawData object, where the slot for the `X` predictor matrix (row
  covariates) is assumed to be the over-parameterized treatment contrasts for
  the conditions and the slot for the `Z` predictor matrix (column covariates)
  is assumed to be the over-parameterized treatment contrasts for the clones
- nPerms = number of permutations. Defaults to `1000`.

# Keyword arguments

- permFun = function used to permute `Y`. Defaults to `shuffle_rows`
  (shuffles rows of `Y`).
- isVarFloor = boolean indicating whether to include variance flooring and
  adjustments. Defaults to `true`.

# Value

Tuple
- Sscores: 2d array of floats; S scores
- pvals: 2d array of floats; permutation p-values

# Some notes

Permutations are computed in parallel when possible.

# Reference

Collins, S. R., Schuldiner, M., Krogan, N. J., & Weissman, J. S. (2006). A
    strategy for extracting and analyzing large-scale quantitative epistatic
    interaction data. Genome biology, 7(7), R63.

"""
function S_score_perms(data::RawData, nPerms::Int64=1000;
                       permFun::Function=shuffle_rows, isVarFloor::Bool=true)

    # Get S scores and permutation p-values
    return perm_pvals(S_score, data, nPerms;
                      permFun=permFun, isVarFloor=isVarFloor)
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export dashbio_speck

"""
    dashbio_speck(;kwargs...)

A Speck component.
The Speck component is a WebGL-based 3D molecule renderer.
Read more about the component here:
https://github.com/wwwtyro/speck
Keyword arguments:
- `id` (String; optional): The ID used to identify this component in Dash callbacks.
- `data` (optional): The xyz file data; a list of atoms such that each atom
has a dictionary defining the x, y, and z coordinates
along with the atom's symbol.. data has the following type: Array of lists containing elements 'symbol', 'x', 'y', 'z'.
Those elements have the following types:
  - `symbol` (String; optional)
  - `x` (Real; optional)
  - `y` (Real; optional)
  - `z` (Real; optional)s
- `presetView` (a value equal to: 'default', 'stickball', 'toon', 'licorice'; optional): One of several pre-loaded views: default, stick-ball, toon,
and licorice
- `scrollZoom` (Bool; optional): The option of whether or not to allow scrolling to control
the zoom.
- `view` (optional): An object that determines and controls various parameters
related to how the molecule is displayed.. view has the following type: lists containing elements 'aspect', 'zoom', 'translation', 'atomScale', 'relativeAtomScale', 'bondScale', 'rotation', 'ao', 'aoRes', 'brightness', 'outline', 'spf', 'bonds', 'bondThreshold', 'bondShade', 'atomShade', 'resolution', 'dofStrength', 'dofPosition', 'fxaa'.
Those elements have the following types:
  - `aspect` (Real; optional)
  - `zoom` (Real; optional)
  - `translation` (optional): . translation has the following type: lists containing elements 'x', 'y'.
Those elements have the following types:
  - `x` (Real; optional)
  - `y` (Real; optional)
  - `atomScale` (Real; optional)
  - `relativeAtomScale` (Real; optional)
  - `bondScale` (Real; optional)
  - `rotation` (optional): . rotation has the following type: lists containing elements .
Those elements have the following types:

  - `ao` (Real; optional)
  - `aoRes` (Real; optional)
  - `brightness` (Real; optional)
  - `outline` (Real; optional)
  - `spf` (Real; optional)
  - `bonds` (Bool; optional)
  - `bondThreshold` (Real; optional)
  - `bondShade` (Real; optional)
  - `atomShade` (Real; optional)
  - `resolution` (Real; optional)
  - `dofStrength` (Real; optional)
  - `dofPosition` (Real; optional)
  - `fxaa` (Real; optional)
"""
function dashbio_speck(; kwargs...)
        available_props = Symbol[:id, :data, :presetView, :scrollZoom, :view]
        wild_props = Symbol[]
        return Component("dashbio_speck", "Speck", "dash_bio", available_props, wild_props; kwargs...)
end


text:
#=
This file is auto-generated. Do not edit.
=#

#! format: off

"""
    mutable struct TwoPartCost <: OperationalCost
        variable::VariableCost
        fixed::Float64
    end

Data Structure Operational Cost Data in two parts: fixed and variable cost.

# Arguments
- `variable::VariableCost`: variable cost
- `fixed::Float64`: fixed cost, validation range: `(0, nothing)`, action if invalid: `warn`
"""
mutable struct TwoPartCost <: OperationalCost
    "variable cost"
    variable::VariableCost
    "fixed cost"
    fixed::Float64
end

function TwoPartCost(; variable, fixed, )
    TwoPartCost(variable, fixed, )
end

# Constructor for demo purposes; non-functional.
function TwoPartCost(::Nothing)
    TwoPartCost(;
        variable=VariableCost((0.0, 0.0)),
        fixed=0.0,
    )
end

"""Get [`TwoPartCost`](@ref) `variable`."""
get_variable(value::TwoPartCost) = value.variable
"""Get [`TwoPartCost`](@ref) `fixed`."""
get_fixed(value::TwoPartCost) = value.fixed

"""Set [`TwoPartCost`](@ref) `variable`."""
set_variable!(value::TwoPartCost, val) = value.variable = val
"""Set [`TwoPartCost`](@ref) `fixed`."""
set_fixed!(value::TwoPartCost, val) = value.fixed = val


text:
require("Pagerank")
using Pagerank

run(`test/bin/setup`)

tests = [
  "test/io/fast.jl",
  "test/io/dangling.jl",
  "test/rank/dangling.jl",
  "test/rank/gaps.jl",
  "test/rank/simple.jl"
]

println("Running suite.")

for test in tests
  println(" * $(test)")
  include(test)
end


text:
using Gadfly

plot(layer(x=[1,2,3], y=[4,5,6], Theme(default_color=colorant"red"), Geom.point),
     layer(x=[4,5,6], y=[1,2,3], Theme(default_color=colorant"blue"), Geom.point),
     Guide.XLabel("XLabel"),
     Guide.YLabel("YLabel"),
     Guide.Title("Title")
     )


text:
#####################
## Sequence OBJECT ##
#####################
"""
    Sequence(GR)

Sequence object.

# Argument
 - `GR::Array{Grad,2}` := Gradient array, fist axis is the dimension (x,y,z) and the second one is time.

# Examples
```julia-repl
julia> A = Sequence([Grad(1,1) for i=1:2,j=1:2])
Sequence(Grad[Grad(1, 1) Grad(1, 1); Grad(1, 1) Grad(1, 1)])

julia> B = Sequence([Grad(2,2) for i=1:2,j=1:2])
Sequence(Grad[Grad(2, 2) Grad(2, 2); Grad(2, 2) Grad(2, 2)])

julia> C = A + B
Sequence(Grad[Grad(1, 1) Grad(1, 1) Grad(2, 2) Grad(2, 2); Grad(1, 1) Grad(1, 1) Grad(2, 2) Grad(2, 2)])
```
"""
mutable struct Sequence
	GR::Array{Grad,2}		#Sequence in (X, Y and Z) and time
	RF::Array{RF,2}			#RF pulses in coil and time
	DAC::Array{DAC,1}		#DAC in time
	function Sequence(GR)	#If no RF is defined, just use a zero amplitude pulse
		M,N = size(GR)
		new([i <= M ? GR[i,j] : Grad(0,GR[1,j].T) for i=1:3, j=1:N],
			reshape(
				[RF(0, GR[1,i].T) for i = 1:N]
				,1,N),
			[DAC(0,GR[1,i].T) for i = 1:N])
	end
	function Sequence(GR,RF)	#If no DAC is defined, just use a zero amplitude pulse
		M,N = size(GR)
		new([i <= M ? GR[i,j] : Grad(0,GR[1,j].T) for i=1:3, j=1:N],
			RF,
			[DAC(0,GR[1,i].T) for i = 1:N])
	end
	Sequence(GR,RF,DAC) = begin
		if any(GR.T .!= RF.T .!= DAC.T)
			error("Gradient, RF, and DAC objects must have the same duration. ")
		end
		M,N = size(GR)
		if M != 3
			new([i <= M ? GR[i,j] : Grad(0,GR[1,j].T) for i=1:3, j=1:N],RF,DAC)
		else
			new(GR,RF,DAC)
		end
	end
end
#TODO: Add trapezoidal grads MACRO
Sequence(GR::Array{Grad,1}) = Sequence(reshape(GR,(length(GR),1)))
Sequence(GR::Array{Grad,1}, RF::Array{RF,1}) = Sequence(
																reshape(GR,(length(GR),1)),
																reshape(RF,(length(RF),1)),
																[DAC(0,GR[1,i].T) for i = 1:size(GR,2)]
																)
Sequence(GR::Array{Grad,1}, RF::Array{RF,1}, D::DAC) = Sequence(
																reshape(GR,(length(GR),1)),
																reshape(RF,(length(RF),1)),
																[D]
																)
Sequence() = Sequence(reshape([Grad(0,0); Grad(0,0)],2,1))
#Sequence operations. TODO: CHANGE RF to :,i whrn including coil sensitivities
Base.length(x::Sequence) = size(x.GR,2)
Base.iterate(x::Sequence) = (Sequence(x.GR[:,1], x.RF[:,1], x.DAC[1]), 2)
Base.iterate(x::Sequence, i::Integer) = (i <= length(x)) ? (Sequence(x.GR[:,i], x.RF[:,i], x.DAC[i]), i+1) : nothing
Base.getindex(x::Sequence, i::UnitRange{Int}) = Sequence(x.GR[:,i], x.RF[:,i], x.DAC[i])
Base.getindex(x::Sequence, i::Int) = Sequence(x.GR[:,i], x.RF[:,i], x.DAC[i])
Base.getindex(x::Sequence, i::BitArray{1}) = any(i) ? Sequence(x.GR[:,i], x.RF[:,i], x.DAC[i]) : nothing
Base.getindex(x::Sequence, i::Array{Bool,1}) = any(i) ? Sequence(x.GR[:,i], x.RF[:,i], x.DAC[i]) : nothing
Base.lastindex(x::Sequence) = size(x.GR,2)
Base.copy(x::Sequence) where Sequence = Sequence([deepcopy(getfield(x, k)) for k ∈ fieldnames(Sequence)]...)
Base.show(io::IO,s::Sequence) = begin
	print(io, "Sequence[ τ = $(round(dur(s)*1e3)) ms | DAC: $(is_DAC_on(s)) | GR: $(size(s.GR)) | RF: $(size(s.RF)) ]")
end
#Arithmetic operations
+(x::Sequence, y::Sequence) = Sequence([x.GR y.GR], [x.RF y.RF], [x.DAC; y.DAC])
-(x::Sequence, y::Sequence) = Sequence([x.GR -y.GR], [x.RF y.RF], [x.DAC; y.DAC])
-(x::Sequence) = Sequence(-x.GR, x.RF, x.DAC)
*(x::Sequence, α::Real) = Sequence(α*x.GR, x.RF, x.DAC)
*(α::Real, x::Sequence) = Sequence(α*x.GR, x.RF, x.DAC)
*(x::Sequence, A::Matrix{Float64}) = Sequence(A*x.GR, x.RF, x.DAC)
*(A::Matrix{Float64}, x::Sequence) = Sequence(A*x.GR, x.RF, x.DAC)
/(x::Sequence, α::Real) = Sequence(x.GR/α, x.RF, x.DAC)
#Sequence object functions
is_DAC_on(x::Sequence) = any(x.DAC.N .> 0)
"Tells if the sequence has elements with DAC on during t."
is_DAC_on(x::Sequence, t::Union{Array{Float64,1},Array{Float64,2}}) = begin
	N = length(x.DAC)
	T = [i==0 ? 0 : x.DAC[i].T for i=0:N]
	ts = [sum(T[1:i])  for i=1:N+1]
	activeGrads = [any(ts[i] .<= t .< ts[i+1]) for i=1:N]
	is_DAC_on(x[activeGrads]) #Get elements of sequence active during time period t
end
"Tells if the sequence has elements with RF on."
is_RF_on(x::Sequence) = any(abs.(getproperty.(x.RF,:A)) .> 0)
"Tells if the sequence has elements with RF on during t."
is_RF_on(x::Sequence, t::Vector{Float64}) = begin
	t = reshape(t,1,length(t));
	N = size(x.RF,2)
	T = [i==0 ? 0 : x.RF[i].T for i=0:N]
	ts = [sum(T[1:i])  for i=1:N+1]
 	activeRFs = [any(ts[i] .<= t .< ts[i+1]) for i=1:N]
	is_RF_on(x[activeRFs]) #Get elements of sequence active during time period t
end
"Duration `T` [s] of the gradient array."
dur(x::Array{Grad,1}) = sum(x[i].T for i=1:size(x,1))
"Duration `T` [s] of the gradient array."
dur(x::Array{Grad,2}) = maximum(sum([x[i,j].T for i=1:size(x,1),j=1:size(x,2)],dims=2))
"Duration `T` [s] of the sequence."
dur(x::Sequence) = dur(x.GR)
"""
	get_grad(seq,dim,t)

Get gradient array from sequence `seq` evaluated in time points `t` for the dimension `dim`.

# Arguments
 - `seq::Sequence`:= Sequence object.
 - `dim::Integer` := Dimension to retrieve, i.e. (1,2,3) = (Gx,Gy,Gz).
 - `t`            := Time vector to evaluate.

# Examples
```julia-repl
julia> seq = Sequence([Grad(1,1) Grad(-1,1); delay(1) delay(1)])
Sequence(Grad[Grad(1, 1) Grad(-1, 1); Grad(0, 1) Grad(0, 1)])

julia> t = 0:.5:2
0.0:0.5:2.0

julia> get_grad(seq,1,t)
5-element Array{Float64,1}:
  1.0
  1.0
  0.0
 -1.0
 -1.0
```
"""
get_grads(seq::Sequence,t) = begin
	M, N = size(seq.GR);
	A = [i <= M ? seq.GR[i,j].A : 0 for i=1:3, j=1:N];
	T = [seq.GR[1,j].T for j=1:N]; T = [0; T];
	⊓(t) = (abs.(t.-1/2) .<= 1/2)*1.;
	(sum([A[j,i]*⊓((t.-sum(T[1:i]))/T[i+1]) for i=1:N]) for j=1:3)
end

get_rfs(seq::Sequence,t) = begin
	M, N = size(seq.RF);
	A = [seq.RF[i,j].A for i=1:M, j=1:N];
	T = [seq.RF[1,j].T for j=1:N]; T = [0; T];
	⊓(t) = (abs.(t.-1/2) .<= 1/2)*1.;
	[sum([A[1,i]*⊓((t.-sum(T[1:i]))/T[i+1]) for i=1:N])]
end

"""
Linear interpolation of the gradients, only works for gradients with uniform duration.
"""
get_grads_linear(seq::Sequence,t) = begin
	M, N = size(seq.GR);
	A = [seq.GR[i,j].A for i=1:M, j=1:N];
	T = [seq.GR[1,j].T for j=1:N]; T = [0; T];
	∧(t) = max.(1 .- abs.(t),0);
	(sum([A[j,i]*∧((t.-sum(T[1:i]))/T[i+1]) for i=1:N]) for j=1:M)
end

get_DAC_on(seq::Sequence,t::Array{Float64,1}) = begin
	Δt = t[2]-t[1]
	M, N = size(seq.DAC)
	DAC = seq.DAC.N .> 0
	T = [seq.DAC[i].T for i=1:N]; T = [0; T]
	⊓(t) = (abs.(t.-1/2.) .<= 1/2)*1.
	DAC_bool = sum([DAC[i]*⊓((t.-sum(T[1:i]))/(T[i+1]+Δt)) for i=1:N])
	DAC_bool = (DAC_bool.>0)
	circshift(convert.(Bool,	[:]),-1)
end
"Calculates the `b`-value, in the PGSE sequnce is b = (2πγGδ)²(Δ-δ/3) and the normalized diffusion vector `n` from a *diffusion sequence*."
get_bvalue(DIF::Sequence) = begin
	M, N = size(DIF.GR)
	G = getproperty.(DIF.GR,:A) #[DIF.GR[i,j].A for i=1:M,j=1:N] #Strength of pulse
	δ = getproperty.(DIF.GR,:T)[1,:] #[DIF.GR[1,j].T for j=1:N]; #Duration of pulse
	T = [sum(δ[1:j]) for j=1:N]; T = [0; T] #Position of pulse
	τ = dur(DIF) #End of sequence
	# B-value
	b = 0
	for k=1:M,i=1:N,j=1:N
		ij = max(i,j)
		α = (i==j) ? 2/3 : 1/2
		b += G[k,i]*G[k,j]*δ[i]*δ[j]*(τ-T[ij]-α*δ[ij])
	end
	b_value = (2π*γ)^2*b # Trace of B tensor
end
"Calculates the `b`-matrix, such as `b`-value = g' B g [s/mm2] with g [T/m]."
get_Bmatrix(DIF::Sequence) = begin
	M, N = size(DIF.GR)
	δ = getproperty.(DIF.GR,:T)[1,:] #[DIF.GR[1,j].T for j=1:N]; #Duration of pulse
	T = [sum(δ[1:j]) for j=1:N]; T = [0; T] #Position of pulse
	τ = dur(DIF) #End of sequence
	# B-value, slower way of doing it
	# b = zeros(M,N,N)
	# for k=1,i=1:N,j=1:N
	# 	ij = max(i,j)
	# 	α = (i==j) ? 2/3 : 1/2
	# 	b[k,i,j] = δ[i]*δ[j]*(τ-T[ij]-α*δ[ij])
	# end
	ij = [max(i,j) for i=1:N, j=1:N]
	α = [(i==j) ? 2/3 : 1/2 for i=1:N, j=1:N]
	b = (δ' .* δ) .* (τ .- T[ij] .- α .* δ[ij])
	b_value = (2π*γ)^2*1e-6*b # Trace of B tensor
	b_value
end

"Calculates the `b` value for PGSE and a moment nulled sequence."
get_bvalue(G,Δ,δ;null::Bool=false) = (null ? 1/9 : 1)*(2π*γ*G*δ)^2*(Δ-δ/3)
"Calculates `q` = γδG diffusion vector from a *diffusion sequence*."
get_qvector(DIF::Sequence;type::String="val") = begin
	M, N = size(DIF.GR)
	G = getproperty.(DIF.GR,:A) #[DIF.GR[i,j].A for i=1:M,j=1:N] #Strength of pulse
	δ = getproperty.(DIF.GR,:T)[1,:] #[DIF.GR[1,j].T for j=1:N]; #Duration of pulse
	T = [sum(δ[1:j]) for j=1:N]; T = [0; T] #Position of pulse
	τ = dur(DIF) #End of sequence
	qτ = γ*sum([G[i,j]*δ[j] for i=1:M,j=1:N],dims=2) #Final q-value
	norm(qτ) ≥ 1e-10 ? error("This is not a diffusion sequence, M0{G(t)}=∫G(t)dt!=0.") : 0
	# q-trajectory
	int_rect(t,A,δ) = (t.≥0) ? A*(t.-(t.≥δ).*(t.-δ)) : 0
	q(t) = γ*sum([int_rect(t.-T[j],G[i,j],δ[j]) for i=1:M,j=1:N],dims=2)
 	# Output q-vector and q-trajectory
	if type=="val"
		q(τ/2) #1/m
	elseif type=="traj"
		t = range(0,stop=τ,length=10^3)
		qt = [q(t)[i] for i=1:M, t=t]
		l = PlotlyJS.Layout(;title="q(t) trajectory", yaxis_title="q [1/mm]",
	    xaxis_title="t [ms]",height=300)
		p = [PlotlyJS.scatter() for j=1:M]
		idx = ["qx" "qy" "qz"]
		for j=1:M
			p[j] = PlotlyJS.scatter(x=t*1e3, y=qt[j,:]*1e-3,name=idx[j])
		end
		PlotlyJS.plot(p, l)
	end
end
get_M0_M1_M2(SEQ::Sequence) = begin
	M, N = size(SEQ.GR)
	G = [SEQ.GR[i,j].A for i=1:M,j=1:N] #Strength of pulse
	δ = [SEQ.GR[1,j].T for j=1:N]; #Duration of pulse
	T = [sum(δ[1:j]) for j=1:N]; T = [0; T] #Position of pulse
	τ = dur(SEQ) #End of sequence
	#M0
	M0i(t,T,A,δ) = (t.≥T) ? A*(t.-T.-(t.≥δ+T).*(t.-δ.-T)) : 0
	M0(t) = sum([M0i(t,T[j],G[i,j],δ[j]) for i=1:M,j=1:N],dims=2)
	#M1
	M1i(t,T,A,δ) = (t.≥T) ? A/2*(t.^2 .-T^2 .-(t.≥δ+T).*(t.^2 .-(δ+T).^2))./2 : 0
	M1(t) = sum([M1i(t,T[j],G[i,j],δ[j]) for i=1:M,j=1:N],dims=2)
	#M1
	M2i(t,T,A,δ) = (t.≥T) ? A/3*(t.^3 .-T^3 .-(t.≥δ+T).*(t.^3 .-(δ+T).^3))./3 : 0
	M2(t) = sum([M2i(t,T[j],G[i,j],δ[j]) for i=1:M,j=1:N],dims=2)
	M0, M1, M2
end
get_max_grad(x::Sequence) = begin
	G = [x.GR[i,j].A for i=1:size(x.GR,1),j=1:size(x.GR,2)]
	Gmax = maximum(sqrt.(sum(G.^2,dims=1)))
end
"Outputs designed k-space trajectory from sequence object."
get_designed_kspace(x::Sequence) = begin
	M, N = size(x.GR)
	F = [x.GR[i,j].A*x.GR[i,j].T for j=1:N, i=1:M]
	F = [zeros(1,M); F]
	if M < 3
		F = [F zeros(N+1,3-M)]
	end
	k = γ*cumsum(F, dims=1)
end
"Outputs actual k-space trajectory for time vector `t` from sequence object."
get_actual_kspace(x::Sequence,t) = begin
	M, N = size(x.GR)
	Nt = length(t)
	Δt = t[2]-t[1]
	G = [get_grad(x,i,t[:]) for i=1:M]
	F = [G[j][i]*Δt for i=1:Nt, j=1:M]
	F = [zeros(1,M); F]
	k = γ*cumsum(F, dims=1)
end

""""Calculates the normalized moments at the end of the sequence. """
function get_Mmatrix(seq::Sequence;order=0)
    M, N = size(seq.GR)
    τ = dur(seq)
    δ = getproperty.(seq.GR,:T)[1,:] #[DIF.GR[1,j].T for j=1:N]; #Duration of pulse
	T = [sum(δ[1:j]) for j=1:N]; T = [0; T[1:end-1]] #Position of pulse
	if order == 0
		M0 = δ/τ
		M1 = δ.*(T .+ δ/2)/τ^2
		M2 = δ.*(T.^2 .+ T.*δ .+ δ.^2/3)/τ^3
		M3 = δ.*(T.^3 .+ 3/2 * T.^2 .*δ .+ T.*δ.^2 .+ δ.^3/4)/τ^4
	elseif order == 1
		M0 = δ/τ
		M1 = δ.*T/τ^2
		M2 = δ.*(T.^2 .+ δ.^2/6)/τ^3
		M3 = δ.*(T.^3 .+ T.*δ.^2/2)/τ^4
	end
    [M0'; M1'; M2'; M3']
end

using LinearAlgebra: I, Bidiagonal, norm
"""Slew rate matrix: |SR*g| ≤ Smax."""
function get_SRmatrix(seq::Sequence)
    _, N = size(seq.GR)
    T = getproperty.(seq.GR[1,:],:T)
    Δt = [T[i] for i=1:N-1]
    dv = [Δt; T[end]]
    ev = Δt;
	SR = Array(Bidiagonal(-1 ./ dv, 1 ./ ev, :U))
	SR = [SR[1,:]'*0 ; SR]; SR[1,1] = T[1]
	SR
end

## TO SCANNER (Philips)
"""Duration in [s] => samples, with dwell-time of Δt = 6.4 μs."""
δ2N(δ) = floor(Int64, δ * 156250) + 2

function write_diff_fwf(DIF, idx180, Gmax, bmax; filename="./qte_vectors_input.txt", name="Maxwell2")
    open(filename, "w") do io
        G1 = DIF[1:idx180-1]
        G2 = DIF[idx180 + 1:end]
        t1 = range(0,dur(G1),length=δ2N(dur(G1)))
		t2 = range(0,dur(G2),length=δ2N(dur(G2)))
		Gx2, Gy2, Gz2 = get_grads_linear(G2,t2)
		Gx1, Gy1, Gz1 = get_grads_linear(G1,t1)
        write(io, "$name $(δ2N(dur(G1))) $(δ2N(dur(G2))) $(bmax)\n")
        for i = 1:length(t1)
            fx1, fy1, fz1 = (Gx1[i], Gy1[i], Gz1[i])./Gmax
            fx2, fy2, fz2 = i ≤ length(t2) ? (-Gx2[i], -Gy2[i], -Gz2[i])./Gmax : (0,0,0)
            line = @sprintf "% .4f % .4f % .4f % .4f % .4f % .4f\n" fx1 fy1 fz1 fx2 fy2 fz2
            write(io, line)
        end
    end
end

function read_diff_fwf(filename="./qte_vectors_input.txt")
    f = readlines(filename)
    G = zeros(size(f[2:end],1),6)
    n1, n2 = tryparse.(Int64, split(f[1]," ")[2:end])
    for (i, l) in enumerate(f[2:end])
        g = tryparse.(Float64, split(l," "))
        G[i,:] = g[g .≠ nothing]
    end
    G, n1, n2
end

## PulseDesigner
module PulseDesigner
using ..MRIsim
using ..MRIsim: γ, get_designed_kspace, get_bvalue, get_max_grad

###############
## RF Pulses ##
###############
RF_hard(B1,T; G=[0,0]) = begin
	EX = Sequence([	Grad(G[1],T);	#Gx
					Grad(G[2],T)],	#Gy
					[RF(B1,T)]		#RF
					)
end
##################
## Gradient SEQ ##
##################
# Without moment nulling
DIF_base(G0,Δ,δ;verbose=false) = begin
	DIF = Sequence([Grad(G0,δ) delay(Δ-δ) -Grad(G0,δ); #Gx
	         		Grad(0, δ) delay(Δ-δ)  Grad(0 ,δ)]) #Gy
	if verbose
	b_value = get_bvalue(DIF)[1]
	GDmax = get_max_grad(DIF)
	println("Sequence with diffusion gradient: "*string(round(GDmax*1e3,digits=2))*" mT/m") # s/mm2
	println("Sequence with b-value: "*string(round(b_value*1e-6,digits=2))*" s/mm2") # s/mm2
	end
	DIF
end
# With moment nulling: M0,M1,M2
DIF_null(G0,Δ,δ,verbose=false) = begin
	if abs(Δ-4δ/3)<1e-6 #Necessary condition
	#Stoecker
	DIF = Sequence([Grad(G0,δ/3) -Grad(G0,2δ/3) delay(Δ-δ) Grad(G0,2δ/3) -Grad(G0,δ/3); #Gx
             		delay(δ/3) delay(2δ/3) delay(Δ-δ) delay(2δ/3) delay(δ/3)]) #Gy
	else
	#Castillo
	DIF = Sequence([ Grad(G0,δ/4) -Grad(G0,δ/2) Grad(G0,δ/4) delay(Δ-δ)
					-Grad(G0,δ/4)  Grad(G0,δ/2)-Grad(G0,δ/4); #Gx
             		 delay(δ/4) delay(δ/2) delay(δ/4) delay(Δ-δ)
					 delay(δ/4) delay(δ/2) delay(δ/4)])
	end

	if verbose
	b_value = get_bvalue(DIF)[1]
	GDmax = get_max_grad(DIF)
	println("Sequence with diffusion gradient: "*string(round(GDmax*1e3,digits=2))*" mT/m") # s/mm2
	println("Sequence with b-value: "*string(round(b_value*1e-6,digits=2))*" s/mm2") # s/mm2
	end
	DIF
end
# EPI
EPI_base(FOV::Float64,N::Int,Δt::Float64,Gmax::Float64) = begin
    #TODO: consider when N is odd
	Nx = Ny = N #Square acquisition
	Δx = FOV/(Nx-1)
	Ta = Δt*(Nx-1) #4-8 us
    Δτ = Ta/(Ny-1)
	Ga = 1/(γ*Δt*FOV)
	Ga ≥ Gmax ? error("Error: Ga exceeds Gmax, increase Δt to at least Δt_min="
	*string(round(1/(γ*Gmax*FOV),digits=2))*" us.") : 0
	EPI = Sequence(vcat(
	    [mod(i,2)==0 ? Grad(Ga*(-1)^(i/2),Ta) : delay(Δτ) for i=0:2*Ny-2],#Gx
	 	[mod(i,2)==1 ? Grad(Ga,Δτ) : delay(Ta) for i=0:2*Ny-2])) #Gy
	DACs = [mod(i,2)==1 ? DAC(0,Δτ) : DAC(N,Ta) for i=0:2*Ny-2]
	EPI.DAC = DACs
	# Recon parameters
	Wx = maximum(abs.(get_designed_kspace(EPI)[:,1]))
	Wy = maximum(abs.(get_designed_kspace(EPI)[:,2]))
	Δx_pix, Δy_pix = 1/Wx, 1/Wy
	Δfx_pix = γ*Ga*Δx_pix
	Δt_phase = Ta*(Nx/(Nx-1)+(Ny-1)^-1) #Time between two echoes
	Δfx_pix_phase = 1/(Ny*Δt_phase)
	println("## EPI parameters ##")
	println("Pixel Δf in freq. direction "*string(round(Δfx_pix,digits=2))*" Hz")
	println("Pixel Δf in phase direction "*string(round(Δfx_pix_phase,digits=2))*" Hz")
    PHASE = Sequence(1/2*[Grad(-Ga, Ta); Grad(-Ga, Ta)])
    DEPHASE = Sequence(1/2*[Grad(-Ga, Ta); Grad(-Ga, Ta)])
	PHASE+EPI+DEPHASE, Δx, Ga, Ta
end

export DIF_base, DIF_null, EPI_base
end


text:
import DiscreteMarkovChains
using Test

@testset "Strongly Connected Components" begin
    V = [1, 2, 3]
    E = [(1, 2), (2, 1), (3, 1)]
    components = DiscreteMarkovChains.strongly_connected_components(V, E)
    @test components == [[1, 2], [3]]
end

@testset "Breadth First Search" begin
    T = [1][:,:]
    period = DiscreteMarkovChains.breadth_first_search(T)
    @test period == 1

    T = [
        0 1;
        1 0
    ]
    period = DiscreteMarkovChains.breadth_first_search(T)
    @test period == 2

    T = [
        0 0.5 0.5;
        0.5 0.5 0;
        0.5 0.5 0;
    ]
    period = DiscreteMarkovChains.breadth_first_search(T)
    @test period == 1

    # this was a crazy bug
    T = [
        0 4 0 0 0;
        1 0 3 0 0;
        0 2 0 2 0;
        0 0 3 0 1;
        0 0 0 4 0;
    ]
    period = DiscreteMarkovChains.breadth_first_search(T)
    @test period == 2
end


text:
#############################################################################
# rationalnorm.jl
#
# Handles the k-norms for k > 1 where k is rational. Reduces the
# k-norm constraint to at most 2d ceil(log2(n + m)) + O(d) second
# order cone constraints. Here d is the dimension of the problem and k
# is converted to a rational with value k = n / m. The procedure for
# this is based on the paper "Second-order cone programming" by
# F. Alizadeh and D. Goldfarb, Mathematical Programming, Series B,
# 95:3-51, 2001, and is documented in the pdf available at
# https://github.com/jump-dev/Convex.jl/raw/master/docs/supplementary/rational_to_socp.pdf
#############################################################################

### k-norm for rational k

struct RationalNormAtom <: AbstractExpr
    head::Symbol
    id_hash::UInt64
    children::Tuple{AbstractExpr}
    size::Tuple{Int,Int}
    k::Rational{Int64}

    function RationalNormAtom(x::AbstractExpr, k::Rational{Int})
        children = (x,)
        k >= 1 || error("p-norms not defined for p < 1")
        return new(:rationalnorm, hash(children), children, (1, 1), k)
    end
end

function sign(x::RationalNormAtom)
    return Positive()
end

# The monotonicity
function monotonicity(x::RationalNormAtom)
    return (sign(x.children[1]) * Nondecreasing(),)
end

function curvature(x::RationalNormAtom)
    return ConvexVexity()
end

function evaluate(x::RationalNormAtom)
    return sum(abs.(evaluate(x.children[1])) .^ x.k)^(1 / x.k)
end

# conic_form!(x::RationalNormAtom, unique_conic_forms)
#
# Formulate the conic constraint
#
# sum(abs(x.^k))^(1/k) <= t
#
# into constraints of the form
#
# v[j] >= abs(x[j]) for all j, s >= 0, sum(s) <= t,
# v[j]^n <= s[j]^m * t^(n-m)
#
# where k is the rational n // m. Then, we use the reduction
# delineated in Duchi & Namkoong (2014) to reduce the last constraint
# into SOC constraints.
function conic_form!(x::RationalNormAtom, unique_conic_forms)
    if !has_conic_form(unique_conic_forms, x)
        # Add extra variables and constraints
        d = length(x.children[1])
        v = Variable(d, 1)
        s = Variable(d, 1)
        t = Variable(1)

        # Adding non-negativity constraints for all variables to problem
        obj = conic_form!(t, unique_conic_forms)
        conic_form!(v >= abs(x.children[1]), unique_conic_forms)
        # conic_form!(t >= 0, unique_conic_forms)
        conic_form!(s >= 0, unique_conic_forms)
        conic_form!(sum(s) <= t, unique_conic_forms)

        # Reduce to SOC constraints (get powers for each element)
        num = Int(numerator(x.k))
        denom = Int(denominator(x.k))
        # Construct list of inequalities of form u^2 <= vw, where the list
        # is given by triples in ineq_list.
        (ineq_list, var_list) =
            psocp.ProductToSimpleInequalities(denom, num - denom)
        if (length(ineq_list) > 10)
            @warn string(
                "Rational norm generating ",
                length(ineq_list),
                " intermediate constraints.\n\tIncreasing ",
                ":max_iters or decreasing solver tolerance\n\tmay give ",
                "more accurate solutions",
            )
        end
        # u corresponds to "introduced" variables; make a matrix of them
        # and then add equality constraints for the first and second
        # variable in the matrix (i.e. first and second columns), as they
        # correspond to s and v.
        u = Variable(length(var_list), d)
        # v is the first variable
        conic_form!(v == u[var_list[1], :]', unique_conic_forms)
        # s[jj] is the second variable
        conic_form!(s == u[var_list[2], :]', unique_conic_forms)

        for jj in 1:d
            # t is the third variable
            conic_form!(t == u[var_list[3], jj], unique_conic_forms)
            for ii in 1:length(ineq_list)
                temp1 = u[ineq_list[ii].left_var_ind, jj]
                temp2 = u[ineq_list[ii].t_var_ind, jj]
                temp3 = u[ineq_list[ii].s_var_ind, jj]
                # Want: t1^2 <= t2 * t3, so add constraints of form
                # norm([2 * t1, t2 - t3], 2) <= t2 + t3.
                conic_form!(
                    SOCElemConstraint(temp2 + temp3, temp2 - temp3, 2 * temp1),
                    unique_conic_forms,
                )
            end
        end
        cache_conic_form!(unique_conic_forms, x, obj)
    end
    return get_conic_form(unique_conic_forms, x)
end
function rationalnorm(x::AbstractExpr, k::Rational{Int})
    if sign(x) == ComplexSign()
        row, col = size(x)
        if row == 1 || col == 1
            return RationalNormAtom(abs(x), k)
        end
    else
        return RationalNormAtom(x, k)
    end
end


text:
using HTTP
using Sockets
using UUIDs
using ContextTracking

# Start a web server that generates a UUID and delegate request to another function
function start_server(host, port)
    socket = Sockets.listen(Sockets.getaddrinfo(host), port)
    task = @async HTTP.listen(; server = socket) do http
        @memo correlation_id = uuid4()
        result = process_request(http)
        write(http, result)
        HTTP.setstatus(http, 200)
        return nothing
    end
    return (task = task, socket = socket)
end

@ctx function process_request(http)
    correlation_id = context()[:correlation_id]
    return "My correlation_id is $correlation_id"
end

function stop_server(server)
    try
        close(server.socket)
        Base.throwto(server.task, InterruptException())
    catch ex
        ex isa InterruptException && return
        rethrow(ex)
    end
end

server = start_server("localhost", 8000)

# Test from the shell prompt
#=
$ curl http://localhost:8000/
OK.  My correlation_id is 73da9213-47e1-4139-8232-7f5c8e43aa08
=#

stop_server(server)


text:
module BCF

import BioCore: BioCore, isfilled, header
import VariantCallFormat: VCF
import BGZFStreams
import BufferedStreams

include("record.jl")
include("reader.jl")
include("writer.jl")

end


text:
using Distributions
using OnlineStats

module VarianceBasedSensitivity

# Include the type definitions
include("types.jl")

# Include test functions
include("test_functions.jl")

export SobolDists, SobolSamples, SobolEvals
export SobolEvalsSingle, SobolEvalsJoint
export SobolEvalsSingleAndJoint
export makeSobolSamples
export makeNewSobolSamples!
export +
export makeC, makeC!, makeCij
export sobolSampler, sobolSamplerSij
export sobolSamplerSingle
export computeSensSi, computeSensSiT
export computeSensSij
export sensSiSTD, sensSiTSTD, sensSijSTD

function makeSobolSamples(S::SobolDists,N::Int)
    # Form the A and B matricies
    A = zeros(Float64,N,S.P)
    B = zeros(Float64,N,S.P)
    for i = 1:S.P
        A[:,i] = Base.rand(S.dist[i],N)
        B[:,i] = Base.rand(S.dist[i],N)
    end
    return SobolSamples(N,S.P,S.dist,A,B)
end

function makeNewSobolSamples!(S::SobolSamples,N::Int)
    # Form the A and B matricies
    A = zeros(N,S.P)
    B = zeros(N,S.P)
    for i = 1:S.P
        A[:,i] = rand(S.dist[i],N)
        B[:,i] = rand(S.dist[i],N)
    end
    S.A = vcat(S.A,A)
    S.B = vcat(S.B,B)
    S.N += N
end

import Base.+
function +(x::SobolSamples,y::SobolSamples)
    # If the distributions match, add together
    if x.dist == y.dist
	z = SobolSamples((x.N+y.N),x.P,x.dist,vcat(x.A,y.A),vcat(x.B,y.B))
    else
        error("Distributions don't match!")
    end
    z
end

function +(x::SobolEvalsSingle,y::SobolEvalsSingle)
    # If the distributions match, add together
    if (x.dist == y.dist)
	tempfC = Array{Array{Float64,2}}(x.P)
        for i = 1:x.P
            tempfC[i] = vcat(x.fC[i],y.fC[i])
        end
	z = SobolEvalsSingle((x.N+y.N),x.P,x.dist,vcat(x.fA,y.fA),vcat(x.fB,y.fB),tempfC)
    else
        error("Distributions don't match!")
    end
    z
end

function +(x::SobolEvalsSingleAndJoint,y::SobolEvalsSingleAndJoint)
    # If the distributions match, add together
    if (x.dist == y.dist)
	tempfC = Array{Array{Float64,2}}(x.P)
        for i = 1:x.P
            tempfC[i] = vcat(x.fC[i],y.fC[i])
        end
	tempfCij = Array{Array{Float64,2}}(length(x.Cpairs))
        for i = 1:length(x.Cpairs)
            tempfCij[i] = vcat(x.fCij[i],y.fCij[i])
        end
	z = SobolEvalsSingleAndJoint((x.N+y.N),x.P,x.dist,vcat(x.fA,y.fA),vcat(x.fB,y.fB),tempfC,x.Cpairs,tempfCij)
    else
        error("Distributions don't match!")
    end
    z
end

function makeC(S::SobolSamples, k::Int)
    C = deepcopy(S.B)
    for i = 1:S.N
        C[i,k] = S.A[i,k];
    end
    C
end

function makeC(S::SobolSamples, i::Int, k::Int)
    C = deepcopy(S.B[k,:])
    C[i] = S.A[k,i];
    C
end

function makeC!(C::Array{Float64,2},S::SobolSamples, k::Int)
    if size(C) == size(S.A)
        for i = 1:S.N, j = vcat(collect(1:k-1),collect(k+1:S.P))
            C[i,j] = S.B[i,j]
        end
        for i = 1:S.N
            C[i,k] = S.A[i,k];
        end
    else
        error("Given C is the wrong sie!")
    end
end

function makeCij(S::SobolSamples, i::Int, j::Int, k::Int)
    C = deepcopy(S.B[k,:])
    C[i] = S.A[k,i]
    C[j] = S.A[k,j]
    C
end

function sobolSamplerSingle(f::Function,sob::SobolSamples)

    fout = f(sob.A[1,:])

    #sample_results = zeros(Float64,sob.N*(sob.P + 2),length(fout))
    sample_results = SharedArray{Float64,2}((sob.N*(sob.P + 2),length(fout)), init = 0)

    @sync @distributed for i = 1:size(sample_results,1)
        if i < sob.N + 1
            sample_results[i,:] = f(sob.A[i,:])

        elseif  i < 2*sob.N + 1
            sample_results[i,:] = f(sob.B[i-sob.N,:])

        else
            n = Int(mod(i-2*sob.N,sob.N))
            n == 0 ? n = sob.N : n
            cval = Int(floor((i - 2*sob.N)/sob.N) + 1)
            n == sob.N ? cval = cval - 1 : cval
            #println("i = ",i,", n = ",n,", cval = ",cval)
            sample_results[i,:] = f(makeC(sob,cval,n))
        end
    end

    formC = Array{Array{Float64,2}}(sob.P)
    for i = 1:length(formC)
        formC[i] = sample_results[(2+i-1)*sob.N+1:(2+i)*sob.N,:]
    end

    #println(sample_results)
    println(sample_results[1:10,:])
    println("  ")

    SobolEvalsSingle(sob.N,sob.P,sob.dist,sample_results[1:sob.N,:],sample_results[sob.N+1:2*sob.N,:],formC)

end

function sobolSamplerJoint(f::Function,sob::SobolSamples)

    fout = f(sob.A[1,:])

    sample_results = zeros(binomial(sob.P,2)*sob.N,length(fout))
    formctemp = zeros(Int,binomial(sob.P,2),2)
    counter = 1
    for i = 1:sob.P-1, j = i+1:sob.P
        formctemp[counter,:] = [i,j]
        counter += 1
    end
    Cij = kron(formctemp, ones(Int,sob.N))

    for i = 1:size(sample_results,1)
            n = Int(mod(i,sob.N))
            n == 0 ? n = sob.N : n
            #println("i = ",i,", n = ",n,", cval = ",cval)
            sample_results[i,:] = f(makeCij(sob,Cij[i,1],Cij[i,2],n))
    end

    CijDict = Array{Array{Int64,1}}(binomial(sob.P,2))
    formCij = Array{Array{Float64,2}}(binomial(sob.P,2))
    for i = 1:length(formCij)
        CijDict[i] = formctemp[i,:]
        formCij[i] = sample_results[(i-1)*sob.N+1:(i)*sob.N,:]
    end

    SobolEvalsJoint(sob.N,sob.P,sob.dist,CijDict,formCij)

end

function singleSobolSampleEvaluator(i::Int,output::SharedArray{Float64,2},sob::SobolSamples,f::Function)
	if i < sob.N + 1
		output[i,:] = f(sob.A[i,:])

	elseif  i < 2*sob.N + 1
		output[i,:] = f(sob.B[i-sob.N,:])

	else
		n = Int(mod(i-2*sob.N,sob.N))
		n == 0 ? n = sob.N : n
		cval = Int(floor((i - 2*sob.N)/sob.N) + 1)
		n == sob.N ? cval = cval - 1 : cval
		output[i,:] = f(makeC(sob,cval,n))
	end
end

function jointSobolSampleEvaluator(i::Int,output::SharedArray{Float64,2},Cij::Array{Int64,2},sob::SobolSamples,f::Function)
	n = Int(mod(i,sob.N))
        n == 0 ? n = sob.N : n
        output[i,:] = f(makeCij(sob,Cij[i,1],Cij[i,2],n))
end

function sobolSampler(f::Function,sob::SobolSamples)

    fout = f(sob.A[1,:])

    #sample_results = zeros(sob.N*(sob.P + 2),length(fout))
    sample_results = SharedArray{Float64,2}((sob.N*(sob.P + 2),length(fout)), init = 0)

    # my pmap call is not working correctly and is FAR too slow
    #pmap((x)->singleSobolSampleEvaluator(x,sample_results,sob,f),1:sob.N*(sob.P+2); batch_size = 100)
    @sync @distributed for i = 1:size(sample_results,1)
        singleSobolSampleEvaluator(i,sample_results,sob,f)
    end

    formC = Array{Array{Float64,2}}(sob.P)
    for i = 1:length(formC)
        formC[i] = sample_results[(2+i-1)*sob.N+1:(2+i)*sob.N,:]
    end

    fA = sample_results[1:sob.N,:]
    fB = sample_results[sob.N+1:2*sob.N,:]

    #sample_results = zeros(binomial(sob.P,2)*sob.N,length(fout))
    sample_results = SharedArray{Float64,2}((binomial(sob.P,2)*sob.N,length(fout)), init = 0)

    formctemp = zeros(Int,binomial(sob.P,2),2)
    counter = 1
    for i = 1:sob.P-1, j = i+1:sob.P
        formctemp[counter,:] = [i,j]
        counter += 1
    end
    Cij = kron(formctemp, ones(Int,sob.N))

    # pmap((x)->jointSobolSampleEvaluator(x,sample_results,Cij,sob,f),1:(binomial(sob.P,2)*sob.N); batch_size = 100)
    @sync @distributed for i = 1:size(sample_results,1)
        jointSobolSampleEvaluator(i,sample_results,Cij,sob,f)
    end

    CijDict = Array{Array{Int64,1}}(binomial(sob.P,2))
    formCij = Array{Array{Float64,2}}(binomial(sob.P,2))
    for i = 1:length(formCij)
        CijDict[i] = formctemp[i,:]
        formCij[i] = sample_results[(i-1)*sob.N+1:(i)*sob.N,:]
    end

    SobolEvalsSingleAndJoint(sob.N,sob.P,sob.dist,fA,fB,formC,CijDict,formCij)

end

# Include the sensitivity metric calculations
include("sobol_inds.jl")

end # module

max_stars_repo_name: mstobb/VarianceBasedSensitivity


text:
# # Preprocessing of Data
#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/notebooks/preprocessing.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/notebooks/preprocessing.ipynb)

# RHEOS offers several functions for sampling and filtering data; this page is intended to be a brief tutorial of their use. For detailed descriptions of functions and their optional arguments, see the [API](@ref) section.

using RHEOS
using PyPlot

# ## Upsampling and Downsampling

# We generate a simple sinusoid with constant sampling rate (more details about data generation are discussed in the [Generating Data](@ref) section).

## Generate a sinusoidal data set
foo = timeline(t_start = 0, t_end = 10, step = 0.4)
foo = strainfunction(foo, t -> sin(t))
fig, ax = subplots(1, 1, figsize = (5, 5))
ax.plot(foo.t, foo.ϵ, "--", marker = "o", markersize = 8, color = "blue")
#!nb fig #hide

# To downsample the full data set , the function [`resample`](@ref) is used with a `scale` argument less than 1. Similarly, to increase the sample rate, the `scale` argument should be greater than 1.

#md # !!! note "Note"
#md #     When scaling, the argument can be an integer, floating point, or rational type. The function contains greater functionality and versatility than shown here. See the [`resample`](@ref) docstring for more information.

#nb # Note: when scaling, the argument can be an integer, floating point, or rational type. The function contains greater functionality and versatility than shown here. See the `resample` docstring for more information.

## Downsample
foo_dsamp = resample(foo, scale = 1//2)
## Upsample
foo_usamp = resample(foo, scale = 2)

## Plotting
fig, ax = subplots(1, 2, figsize = (10, 5))
ax[1].set_title("Downsampling")
ax[1].plot(foo.t, foo.ϵ, "--", marker = "o", markersize = 8, color = "blue")
ax[1].plot(foo_dsamp.t, foo_dsamp.ϵ, "--", marker = "x", markersize = 10, markeredgewidth = 2, color = "orange")
ax[2].set_title("Upsampling")
ax[2].plot(foo.t, foo.ϵ, "--", marker = "o", markersize = 8, color = "blue")
ax[2].plot(foo_usamp.t, foo_usamp.ϵ, "--", marker = "x", markersize = 10, markeredgewidth = 2, color = "orange")
#!nb fig #hide

# ## Cutting

# RHEOS provides a dedicated function, [`cutting`](@ref), to remove the data outside a specified time interval.

foo_cut = cutting(foo, 2.0, 8.0)

## Plotting
fig, ax = subplots(1, 1, figsize = (5, 5))
ax.plot(foo.t, foo.ϵ, "--", marker = "o", markersize = 8, color = "blue")
ax.plot(foo_cut.t, foo_cut.ϵ, "--", marker = "x", markersize = 10, markeredgewidth = 2, color = "orange")
#!nb fig #hide

# ## Smoothing

# Lastly, RHEOS provides a smoothing function, [`smooth`](@ref). The first argument is the data to smooth and the second argument is the (very) approximate time scale of smoothing. (It uses Gaussian smoothing and can be thought of as a low pass filter for information occuring on time scales shorter than the 2nd argument). The padding can be changed using a keyword argument if desired, see [API](@ref) and [ImageFiltering.jl Documentation](https://juliaimages.github.io/ImageFiltering.jl/stable/) for more details.

foo_s = timeline(t_start = 0, t_end = 10, step = 0.02)
foo_s = strainfunction(foo_s, t -> 3*sin(t))
noise = strainfunction(foo_s, t -> rand())

foo_noisy = foo_s + noise
foo_smooth = smooth(foo_noisy, 1)

fig, ax = subplots(1, 1, figsize = (5, 5))
ax.plot(foo_noisy.t, foo_noisy.ϵ, color = "blue")
ax.plot(foo_smooth.t, foo_smooth.ϵ, color = "orange")
#!nb fig #hide


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    Base.precompile(Tuple{typeof(PkgServer.get_registries),String})
    Base.precompile(Tuple{typeof(PkgServer.serve_json),HTTP.Streams.Stream{HTTP.Messages.Request,HTTP.ConnectionPool.Transaction{Sockets.TCPSocket}},Dict{String,Any}})
    Base.precompile(Tuple{typeof(PkgServer.stream_file),IOStream,Int64,Int64,Task,Base.PipeEndpoint})
    Base.precompile(Tuple{typeof(PkgServer.verify_registry_hash),String,String})
    Base.precompile(Tuple{typeof(PkgServer.wait_first),Task,Vararg{Task,N} where N})
    Base.precompile(Tuple{typeof(collect),Base.KeySet{String,Dict{String,PkgServer.RegistryMeta}}})
    Base.precompile(Tuple{typeof(getindex),Dict{String,PkgServer.RegistryMeta},String})
    Base.precompile(Tuple{typeof(in),SubString{String},Base.KeySet{String,Dict{String,PkgServer.RegistryMeta}}})
    isdefined(HTTP.Servers, Symbol("#13#14")) && Base.precompile(Tuple{getfield(HTTP.Servers, Symbol("#13#14"))})
    isdefined(HTTP.Servers, Symbol("#8#9")) && Base.precompile(Tuple{getfield(HTTP.Servers, Symbol("#8#9"))})
    isdefined(PkgServer, Symbol("#1#2")) && Base.precompile(Tuple{getfield(PkgServer, Symbol("#1#2"))})
    isdefined(PkgServer, Symbol("#22#26")) && Base.precompile(Tuple{getfield(PkgServer, Symbol("#22#26"))})
    isdefined(PkgServer, Symbol("#30#33")) && Base.precompile(Tuple{getfield(PkgServer, Symbol("#30#33"))})
    isdefined(PkgServer, Symbol("#38#46")) && Base.precompile(Tuple{getfield(PkgServer, Symbol("#38#46")),String})
    isdefined(PkgServer, Symbol("#39#47")) && Base.precompile(Tuple{getfield(PkgServer, Symbol("#39#47"))})
    isdefined(PkgServer, Symbol("#40#48")) && Base.precompile(Tuple{getfield(PkgServer, Symbol("#40#48"))})
    isdefined(PkgServer, Symbol("#41#49")) && Base.precompile(Tuple{getfield(PkgServer, Symbol("#41#49"))})
    isdefined(PkgServer, Symbol("#42#50")) && Base.precompile(Tuple{getfield(PkgServer, Symbol("#42#50"))})
    isdefined(PkgServer, Symbol("#65#68")) && Base.precompile(Tuple{getfield(PkgServer, Symbol("#65#68"))})
    isdefined(PkgServer, Symbol("#ServerConfig#62#63")) && Base.precompile(Tuple{getfield(PkgServer, Symbol("#ServerConfig#62#63")),Sockets.InetAddr{Sockets.IPv4},String,Dict{String,PkgServer.RegistryMeta},Array{String,1},Int64,Type{PkgServer.ServerConfig}})
end


text:
@testset "analysis" begin

Context() do ctx
Builder(ctx) do builder
LLVM.Module("SomeModule", ctx) do mod
    ft = LLVM.FunctionType(LLVM.Int32Type(ctx))
    fn = LLVM.Function(mod, "SomeFunction", ft)

    entry = BasicBlock(fn, "entry")
    position!(builder, entry)

    ret!(builder)

    @test_throws LLVMException verify(mod)
    @test_throws LLVMException verify(fn)
end
end
end

Context() do ctx
Builder(ctx) do builder
LLVM.Module("SomeModule", ctx) do mod
    ft = LLVM.FunctionType(LLVM.VoidType(ctx))
    fn = LLVM.Function(mod, "SomeFunction", ft)

    entry = BasicBlock(fn, "entry")
    position!(builder, entry)

    ret!(builder)

    verify(mod)
    verify(fn)
end
end
end

end


text:
using Raylib
using Raylib: RayCamera3D, rayvector, RayVector3, RayColor, raycolor,
    RAYWHITE, LIGHTGRAY, DARKGRAY, MAROON,
    RED, BLUE, LIME, GOLD, SKYBLUE, BLACK

function main()
    screenWidth = 800
    screenHeight = 450

    Raylib.InitWindow(
        screenWidth, screenHeight,
        "raylib [core] example - 3d camera first person"
    )

    camera = RayCamera3D(
        rayvector(4, 2, 4),
        rayvector(0, 1.8, 0),
        rayvector(0, 1, 0),
        60,
        0,
    )

    heights = Vector{Cfloat}(undef, 20)
    positions = Vector{RayVector3}(undef, 20)
    colors = Vector{RayColor}(undef, 20)

    for i = 1:20
        heights[i] = rand(1:12)
        positions[i] = rayvector(
            rand(-15:15),
            heights[i] / 2,
            rand(-15:15),
        )
        colors[i] = raycolor(
            rand(20:255),
            rand(10:55),
            30,
            255,
        )
    end

    Raylib.SetCameraMode(camera, Int(Raylib.CAMERA_FIRST_PERSON))
    Raylib.SetTargetFPS(60)

    while !Raylib.WindowShouldClose()
        Raylib.UpdateCamera!(camera)

        Raylib.BeginDrawing()
        Raylib.ClearBackground(RAYWHITE)

        Raylib.BeginMode3D(camera)

        Raylib.DrawPlane(
            rayvector(0,0,0),
            rayvector(32, 32),
            LIGHTGRAY)

        Raylib.DrawCube(
            rayvector(-16, 2.5, 0), 1f0, 5f0, 32f0,
            BLUE)
        Raylib.DrawCube(
            rayvector(16, 2.5, 0), 1f0, 5f0, 32f0,
            LIME)
        Raylib.DrawCube(
            rayvector(0, 2.5, 16), 32f0, 5f0, 1f0,
            GOLD)

        for i = 1:20
            Raylib.DrawCube(
                positions[i], 2f0, heights[i], 2f0,
                colors[i])
            Raylib.DrawCubeWires(positions[i], 2f0, heights[i], 2f0,
                                 MAROON)

        end

        Raylib.EndMode3D()

        Raylib.DrawRectangle(10, 10, 220, 70,
                             Raylib.Fade(SKYBLUE, 0.5) )
        Raylib.DrawRectangleLines(10, 10, 220, 70, BLUE)

        Raylib.DrawText("First person camera default controls:",
                        20, 20, 10, BLACK)
        Raylib.DrawText("- Move with keys: W, A, S, D",
                        40, 40, 10, DARKGRAY)
        Raylib.DrawText("- Mouse move to look around",
                        40, 60, 10, DARKGRAY)

        Raylib.EndDrawing()
    end

    Raylib.CloseWindow()
end


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

@testset "graph.cycle" begin

@testset "minimumcyclebasis" begin
    p5 = pathgraph(5)
    @test isempty(minimumcyclebasis(p5))
    c5 = cyclegraph(5)
    @test length(minimumcyclebasis(c5)[1]) == 5
    k44 = completebipartite(3, 3)
    @test sum(map(length, minimumcyclebasis(k44))) == 16
    k5 = completegraph(5)
    @test sum(map(length, minimumcyclebasis(k5))) == 18
    cycs = plaingraph(6, [(1, 2), (2, 3), (3, 1), (3, 4), (4, 5), (5, 6), (6, 4)])
    @test length(minimumcyclebasis(cycs)) == 2
end

end # graph.cycle


text:
module Counters

# stdlib
using Dates
using Mmap
using Random
using Serialization
using Sockets

# AutoTM dependencies
using AutoTM
using MaxLFSR
using nGraph
using PCM
using PersistentArrays
using Slabs
using SystemSnoop

# External Packages
using DataStructures
using PrettyTables
using SIMD
using StructArrays
using Tables
using PGFPlotsX

# We use a `NamedPipe` to communicate from the process that is actually running the
# benchmarks to the process that is recording the counters.
#
# The NamedPipe has the following commands:
#
# * start -- start sampling.
# * stop -- stop sampling and serialize data to the currently set filepath.
# * shutdown -- exit.
# * filepath <payload> -- Set the serialization path to <payload>
# * measurements -- Indicate that a collection of SystemSnoop measurements are ready
#                   to be deserialized from `TRANSFERPATH` (defined below)
# * params -- The parameters of this run to be recorded in the DataTable.
#             Transfer methodology is the same as `measurements`.
# * sampletime <time> -- Configure the time between samples.
const SRCDIR = @__DIR__
const PKGDIR = dirname(SRCDIR)
const DATADIR = joinpath(PKGDIR, "data")
const TEMPDIR = joinpath(PKGDIR, "temp")

const PIPEPATH = joinpath(TEMPDIR, "pipe")
const TRANSFERPATH = joinpath(TEMPDIR, "transfer.jls")
const PARAMPATH = joinpath(TEMPDIR, "params.jls")
const FIGDIR = joinpath(PKGDIR, "figures")

function __init__()
    # Create temporary directories if needed.
    for dir in (DATADIR, TEMPDIR, FIGDIR)
        !ispath(dir) && mkdir(dir)
    end
end

# General Utility Functions
include("util.jl")

# Bridge into PCM. Also contains code that allows PCM data structures to use the
# SystemSnoop API
include("pcm.jl")

# Setup a server on `PIPEPATH` that will listen and perform requested measurements.
include("server.jl")

# Code for running AutoTM workloads in 1LM or 2LM.
include("autotm.jl")

# Kernels for testing the DRAM cache system.
include("kernels.jl")
include("kernel_runner.jl")

# Plot results
include("plots.jl")
include("heapplot.jl")

# Top level entry point for running all benchmarks
include("measurements.jl")
include("benchmarks.jl")

end # module

max_stars_repo_name: darchr/AutoTM

max_issues_repo_name: darchr/AutoTM

max_forks_repo_name: darchr/AutoTM