text:
"""
    AmoebaFiber2D(f::AbstractPolynomial, w=(0., 0.))

Construct the fiber `AmoebaFiber2D` of the bivariate polynomial `f` at `w`.
"""
function AmoebaFiber2D(p::MP.AbstractPolynomialLike{S}, w=(0.0, 0.0)) where S
    vars = MP.variables(p)
    @assert (length(vars) == 2) "Expected a bivariate polynomial, but got $(p)."

    T = realified_coefficient_type(S)
    f_re, f_im = SP.Polynomial.(1.0 .* realify(p))
    AmoebaFiber2D(f_re, f_im, w)
end

function AmoebaFiber2D(f_re::SP.Polynomial{T}, f_im::SP.Polynomial, w=(0.0, 0.0)) where T
    v = exp.(w)
    U = zeros(T, 2, 4)
    x = SVector{4, Float64}(0.0, 0.0, 0.0, 0.0)

    exponents = SP.exponents.((f_re, f_im))
    original_coefficients = copy.(SP.coefficients.((f_re, f_im)))

    AmoebaFiber2D(f_re, f_im, v, U, x, exponents, original_coefficients)
end

function update_fiber!(F::AmoebaFiber2D, w)
    w1, w2 = w
    F.v = (exp(w1), exp(w2))
    condition!(F)
    F
end

function condition!(fiber::AmoebaFiber2D)
    fiber.f_re.coefficients .= fiber.original_coefficients[1]
    fiber.f_im.coefficients .= fiber.original_coefficients[2]

    n1 = normalizer_amoeba_2d(fiber.exponents[1], fiber.f_re.coefficients, fiber.v)
    n2 = normalizer_amoeba_2d(fiber.exponents[2], fiber.f_im.coefficients, fiber.v)
    normalizer = inv(max(n1, n2))
    SP.scale_coefficients!(fiber.f_re, normalizer)
    SP.scale_coefficients!(fiber.f_im, normalizer)
    nothing
end

function normalizer_amoeba_2d(E, c, v)
    v1, v2 = v
    normalizer = sum(1:size(E, 2)) do j
        m1 = E[1,j]+E[3,j]
        m2 = E[2,j]+E[4,j]
        abs(c[j]) * Base.FastMath.pow_fast(v1, m1) *
        Base.FastMath.pow_fast(v2, m2)
    end
end

function phi!(fiber::AmoebaFiber2D, θ)
    v1, v2 = fiber.v
    x = SVector(v1 * cos(θ[1]), v2 * cos(θ[2]), v1 * sin(θ[1]), v2 * sin(θ[2]))
    fiber.x = x
    x
end

function evaluate(fiber::AmoebaFiber2D{T}, θ) where T
    x = phi!(fiber, θ)
    SVector(SP.evaluate(fiber.f_re, x), SP.evaluate(fiber.f_im, x))
end

function jacobian(fiber::AmoebaFiber2D, θ, precomputed=false)
    x = precomputed ? fiber.x : phi!(fiber, θ)
    U = fiber.U
    @inbounds begin
        ∇f_re = SP.gradient(fiber.f_re, x)
        for i=1:4
            U[1, i] = ∇f_re[i]
        end
        U[2, 3] = U[1, 1]
        U[2, 4] = U[1, 2]
        U[2, 1] = -U[1, 3]
        U[2, 2] = -U[1, 4]
    end

    @inbounds out = begin
        x3, x4 = -x[3], -x[4]
        ∂re∂θ1 = U[1, 1] * x3
        ∂re∂θ2 = U[1, 2] * x4
        ∂im∂θ1 = U[2, 1] * x3
        ∂im∂θ2 = U[2, 2] * x4

        ∂re∂θ1 = muladd(U[1, 3], x[1], ∂re∂θ1)
        ∂re∂θ2 = muladd(U[1, 4], x[2], ∂re∂θ2)
        ∂im∂θ1 = muladd(U[2, 3], x[1], ∂im∂θ1)
        ∂im∂θ2 = muladd(U[2, 4], x[2], ∂im∂θ2)

        # SMatrix has column-major storage
        SMatrix{2,2}(∂re∂θ1, ∂im∂θ1, ∂re∂θ2, ∂im∂θ2)
    end

    out
end


text:
# WAH Vector Cursors
# ==================
#
# A type that acts as a cursor through a WAHVector during some looping
# operations.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/WAHVectors.jl/blob/master/LICENSE.md

mutable struct WAHCursor
    wahvec::WAHVector
    len::Int
    word_i::Int
    fill::UInt32
    isfill::Bool
    fillbit::UInt32
    nwords::UInt32

    function WAHCursor(vec::WAHVector)
        x = new(vec, length(vec.data), 1, 0x00000000, false, 0x00000000, 0x00000000)
        decode!(x)
        return x
    end
end

@inline function decode!(x::WAHCursor)
    currentWord = x.wahvec.data[x.word_i]
    if isruns(currentWord)
        isones = is_ones_runs(currentWord)
        x.fill = ifelse(isones, 0x7FFFFFFF, 0x00000000)
        x.nwords = nwords(currentWord)
        x.isfill = true
        x.fillbit = ifelse(isones, 0x00000001, 0x00000000)
    else
        x.nwords = 0x00000001
        x.isfill = false
        x.fill = x.wahvec.data[x.word_i]
    end
end

@inline function check_to_move!(x::WAHCursor)
    if x.nwords == 0x00000000
        x.word_i += 1
        if x.word_i <= x.len
            decode!(x)
        end
    end
end

@inline function Base.:&(x::WAHCursor, y::WAHCursor)
    if x.isfill
        if y.isfill
            num_words = min(x.nwords, y.nwords)
            x.nwords -= num_words
            y.nwords -= num_words
            #answer = WAHElement(x.fillbit & y.fillbit, num_words)
            return WAHElement(x.fillbit & y.fillbit, num_words)
        else
            #answer = WAHElement(x.fill & y.wahvec.data[y.word_i])
            #answer = WAHElement(x.fill & y.fill)
            x.nwords -= 0x00000001
            y.nwords = 0x00000000
        end
    elseif y.isfill
        #answer = WAHElement(x.wahvec[x.word_i] & y.fill)
        #answer = WAHElement(x.fill & y.fill)
        y.nwords -= 0x00000001
        x.nwords = 0x00000000
    else
        #answer = WAHElement(x.wahvec[x.word_i] & y.wahvec[y.word_i])
        #answer = WAHElement(x.fill & y.fill)
        x.nwords = 0x00000000
        y.nwords = 0x00000000
    end
    #return answer
    return WAHElement(x.fill & y.fill)
end

@inline function _compute_and_vec(x::WAHCursor, y::WAHCursor)
    num_words = min(x.nwords, y.nwords)
    xyfill = x.isfill & y.isfill
    xfill = x.isfill & !y.isfill
    yfill = !x.isfill & x.isfill
    nofill = !x.isfill & y.isfill
    x.nwords -= ifelse(xyfill, num_words, 0)
    x.nwords -= ifelse(xfill, 1, 0)
    x.nwords = ifelse(yfill | nofill, 0, x.nwords)
    y.nwords -= ifelse(xyfill, num_words, 0)
    y.nwords -= ifelse(yfill, 1, 0)
    y.nwords = ifelse(xfill | nofill, 0, y.nwords)
    answer = ifelse(xyfill,
        WAHElement(x.fillbit & y.fillbit, num_words),
        WAHElement(x.fill & y.fill))
    return answer
end

ext: jl


text:
function FieldAngularCoordinates(arg0::Field, arg1::AngularCoordinates)
    return FieldAngularCoordinates((Field, AngularCoordinates), arg0, arg1)
end

function FieldAngularCoordinates(arg0::FieldPVCoordinates, arg1::FieldPVCoordinates, arg2::FieldPVCoordinates, arg3::FieldPVCoordinates, arg4::jdouble)
    return FieldAngularCoordinates((FieldPVCoordinates, FieldPVCoordinates, FieldPVCoordinates, FieldPVCoordinates, jdouble), arg0, arg1, arg2, arg3, arg4)
end

function FieldAngularCoordinates(arg0::FieldRotation)
    return FieldAngularCoordinates((FieldRotation,), arg0)
end

function FieldAngularCoordinates(arg0::FieldRotation, arg1::FieldVector3D)
    return FieldAngularCoordinates((FieldRotation, FieldVector3D), arg0, arg1)
end

function FieldAngularCoordinates(arg0::FieldRotation, arg1::FieldVector3D, arg2::FieldVector3D)
    return FieldAngularCoordinates((FieldRotation, FieldVector3D, FieldVector3D), arg0, arg1, arg2)
end

function add_offset(obj::FieldAngularCoordinates, arg0::FieldAngularCoordinates)
    return jcall(obj, "addOffset", FieldAngularCoordinates, (FieldAngularCoordinates,), arg0)
end

function apply_to(obj::FieldAngularCoordinates, arg0::FieldPVCoordinates)
    return jcall(obj, "applyTo", FieldPVCoordinates, (FieldPVCoordinates,), arg0)
end

function apply_to(obj::FieldAngularCoordinates, arg0::PVCoordinates)
    return jcall(obj, "applyTo", FieldPVCoordinates, (PVCoordinates,), arg0)
end

function apply_to(obj::FieldAngularCoordinates, arg0::TimeStampedFieldPVCoordinates)
    return jcall(obj, "applyTo", TimeStampedFieldPVCoordinates, (TimeStampedFieldPVCoordinates,), arg0)
end

function apply_to(obj::FieldAngularCoordinates, arg0::TimeStampedPVCoordinates)
    return jcall(obj, "applyTo", TimeStampedFieldPVCoordinates, (TimeStampedPVCoordinates,), arg0)
end

function create_from_modified_rodrigues(::Type{FieldAngularCoordinates}, arg0::Vector{Vector{RealFieldElement}})
    return jcall(FieldAngularCoordinates, "createFromModifiedRodrigues", FieldAngularCoordinates, (Vector{Vector{RealFieldElement}},), arg0)
end

function estimate_rate(::Type{FieldAngularCoordinates}, arg0::FieldRotation, arg1::FieldRotation, arg2::RealFieldElement)
    return jcall(FieldAngularCoordinates, "estimateRate", FieldVector3D, (FieldRotation, FieldRotation, RealFieldElement), arg0, arg1, arg2)
end

function estimate_rate(::Type{FieldAngularCoordinates}, arg0::FieldRotation, arg1::FieldRotation, arg2::jdouble)
    return jcall(FieldAngularCoordinates, "estimateRate", FieldVector3D, (FieldRotation, FieldRotation, jdouble), arg0, arg1, arg2)
end

function get_identity(::Type{FieldAngularCoordinates}, arg0::Field)
    return jcall(FieldAngularCoordinates, "getIdentity", FieldAngularCoordinates, (Field,), arg0)
end

function get_modified_rodrigues(obj::FieldAngularCoordinates, arg0::jdouble)
    return jcall(obj, "getModifiedRodrigues", Vector{Vector{RealFieldElement}}, (jdouble,), arg0)
end

function get_rotation(obj::FieldAngularCoordinates)
    return jcall(obj, "getRotation", FieldRotation, ())
end

function get_rotation_acceleration(obj::FieldAngularCoordinates)
    return jcall(obj, "getRotationAcceleration", FieldVector3D, ())
end

function get_rotation_rate(obj::FieldAngularCoordinates)
    return jcall(obj, "getRotationRate", FieldVector3D, ())
end

function revert(obj::FieldAngularCoordinates)
    return jcall(obj, "revert", FieldAngularCoordinates, ())
end

function shifted_by(obj::FieldAngularCoordinates, arg0::RealFieldElement)
    return jcall(obj, "shiftedBy", FieldAngularCoordinates, (RealFieldElement,), arg0)
end

function shifted_by(obj::FieldAngularCoordinates, arg0::jdouble)
    return jcall(obj, "shiftedBy", FieldAngularCoordinates, (jdouble,), arg0)
end

function subtract_offset(obj::FieldAngularCoordinates, arg0::FieldAngularCoordinates)
    return jcall(obj, "subtractOffset", FieldAngularCoordinates, (FieldAngularCoordinates,), arg0)
end

function to_angular_coordinates(obj::FieldAngularCoordinates)
    return jcall(obj, "toAngularCoordinates", AngularCoordinates, ())
end

function to_derivative_structure_rotation(obj::FieldAngularCoordinates, arg0::jint)
    return jcall(obj, "toDerivativeStructureRotation", FieldRotation, (jint,), arg0)
end

function to_univariate_derivative1_rotation(obj::FieldAngularCoordinates)
    return jcall(obj, "toUnivariateDerivative1Rotation", FieldRotation, ())
end


text:
module Billard_Balls16

using Modia3D
using Modia3D.StaticArrays
import Modia3D.ModiaMath

vmatBalls = Modia3D.Material(color="Red" , transparency=0.0)
vmatTable = Modia3D.Material(color="DarkGreen" , transparency=0.5)

cmatTable   = "BilliardTable"
cmatBall    = "BilliardBall"
cmatCushion = "BilliardCushion"

diameter = 0.06
radius = diameter/2
distance_balls = sqrt(3)/2*diameter
dist = 0.0001

TableX = 2.24
TableY = 1.12
TableZ = 0.05
LyBox = 0.05
LzBox = 0.05
hz  = 7*radius/5
hzz = hz/3

@assembly Table(world) begin
  box = Modia3D.Solid(Modia3D.SolidBox(TableX, TableY, TableZ) , "BilliardTable", vmatTable; contactMaterial = cmatTable)
  table = Modia3D.Object3D(world, box, r=[0.0, 0.0, -TableZ/2], fixed=true)
end

@assembly Cushion(world) begin
  cushionPart11 = Modia3D.Solid(Modia3D.SolidBox(LyBox, TableY, hz) , "BilliardCushion", vmatTable)
  cushionPart12 = Modia3D.Solid(Modia3D.SolidBox(1.5*LyBox, TableY, hzz,rsmall=0.0) , "BilliardCushion", vmatTable; contactMaterial = cmatCushion)
  cushionPart21 = Modia3D.Solid(Modia3D.SolidBox(TableX, LyBox , hz) , "BilliardCushion", vmatTable)
  cushionPart22 = Modia3D.Solid(Modia3D.SolidBox(TableX, 1.5*LyBox, hzz,rsmall=0.0) , "BilliardCushion", vmatTable; contactMaterial = cmatCushion)

  cushion11 = Modia3D.Object3D(world    , cushionPart11, r=[TableX/2 - LyBox/2, 0.0, hz/2], fixed=true)
  cushion12 = Modia3D.Object3D(cushion11, cushionPart12, r=[-LyBox/4     , 0.0, hz/2+hzz/2], fixed=true)
  cushion31 = Modia3D.Object3D(world    , cushionPart11, r=[-TableX/2 + LyBox/2, 0.0, hz/2], fixed=true)
  cushion32 = Modia3D.Object3D(cushion31, cushionPart12, r=[LyBox/4     , 0.0, hz/2+hzz/2], fixed=true)
  cushion21 = Modia3D.Object3D(world    , cushionPart21, r=[0.0, TableY/2 - LyBox/2 , hz/2], fixed=true)
  cushion22 = Modia3D.Object3D(cushion21, cushionPart22, r=[0.0, -LyBox/4 , hz/2+hzz/2], fixed=true)
  cushion41 = Modia3D.Object3D(world    , cushionPart21, r=[0.0, -TableY/2 + LyBox/2 , hz/2], fixed=true)
  cushion42 = Modia3D.Object3D(cushion41, cushionPart22, r=[0.0, LyBox/4 , hz/2+hzz/2], fixed=true)
end

@assembly BillardBall(world, xPos, yPos) begin
  ball = Modia3D.Object3D(world, Modia3D.Solid(Modia3D.SolidSphere(diameter), "BilliardBall", vmatBalls ; contactMaterial = cmatBall), fixed = false, r=[xPos, yPos, diameter/2])
end

@assembly Billard() begin
  world = Modia3D.Object3D(visualizeFrame=false)
  table = Table(world)
  cushion = Cushion(world)
  ballStart = Modia3D.Object3D(world, Modia3D.Solid(Modia3D.SolidSphere(diameter), "BilliardBall", vmatBalls ; contactMaterial = cmatBall), fixed = false, r=[-0.8, -0.1, diameter/2], v_start=[3.0, 0.1, 0.0])

  ball1 = BillardBall(world, TableX/6, 0.0)

  ball2  = BillardBall(world, TableX/6 + 1*distance_balls + dist,  1/2*(diameter+dist))
  ball3  = BillardBall(world, TableX/6 + 1*distance_balls + dist, -1/2*(diameter+dist))

  ball4  = BillardBall(world, TableX/6 + 2*(distance_balls+dist),  (diameter+dist))
  ball5  = BillardBall(world, TableX/6 + 2*(distance_balls+dist),  0.0)
  ball6  = BillardBall(world, TableX/6 + 2*(distance_balls+dist), -(diameter+dist))

  ball7  = BillardBall(world, TableX/6 + 3*(distance_balls+dist),  3/2*(diameter+dist))
  ball8  = BillardBall(world, TableX/6 + 3*(distance_balls+dist),  1/2*(diameter+dist))
  ball9  = BillardBall(world, TableX/6 + 3*(distance_balls+dist), -1/2*(diameter+dist))
  ball10 = BillardBall(world, TableX/6 + 3*(distance_balls+dist), -3/2*(diameter+dist))

  ball11 = BillardBall(world, TableX/6 + 4*(distance_balls+dist),  2*(diameter+dist))
  ball12 = BillardBall(world, TableX/6 + 4*(distance_balls+dist),  (diameter+dist))
  ball13 = BillardBall(world, TableX/6 + 4*(distance_balls+dist),  0.0)
  ball14 = BillardBall(world, TableX/6 + 4*(distance_balls+dist), -(diameter+dist))
  ball15 = BillardBall(world, TableX/6 + 4*(distance_balls+dist), -2*(diameter+dist))
end

gravField = Modia3D.UniformGravityField(g=9.81, n=[0,0,-1])
bill = Billard(sceneOptions=Modia3D.SceneOptions(gravityField=gravField,visualizeFrames=false, defaultFrameLength=0.2,nz_max = 100, enableContactDetection=true, visualizeContactPoints=false, visualizeSupportPoints=false))

model = Modia3D.SimulationModel( bill )
result = ModiaMath.simulate!(model; stopTime=5.0, tolerance=1e-8,interval=0.0001, log=false)

println("... success of Billard_Balls16.jl!")
end


text:
"""
    get_parameters_vals(items::Vector{<:AbstractItem})

# Description

Returns a matrix with item parameters displayed by row.
"""
function get_parameters_vals(items::Vector{<:AbstractItem})
ret = Vector{Vector{Float64}}(undef, size(items, 1))
max_length = 1
i_2 = 0
for i in items
    i_2 += 1
    local pars = get_parameters_vals(i)
    ret[i_2] = copy(pars)
    max_length = max_length < size(pars, 1) ? size(pars, 1) : max_length
end

for i_3 = 1:i_2
    local length_i = size(ret[i_3], 1)
    if length_i < max_length
        ret[i_3] = vcat(ret[i_3], zeros(Float64, max_length - length_i))
    end
end
return permutedims(reduce(hcat, ret))
end

"""
    get_parameters(items::Vector{<:AbstractItem})

# Description

Returns a vector of `<:AbstractParameters` objects.
"""
function get_parameters(items::Vector{<:AbstractItem})
return map( i -> i.parameters, items)
end


text:
using Documenter, Quaternionic

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

makedocs(
    sitename="Quaternionic.jl",
    modules = [Quaternionic],
    format = Documenter.HTML(
        prettyurls = !("local" in ARGS),  # Use clean URLs, unless built as a "local" build
        edit_link = "main",  # Link out to "main" branch on github
        canonical = "https://moble.github.io/Quaternionic.jl/stable/",
    ),
    pages = [
        "Introduction" => "index.md",
        "Basics" => "manual.md",
        "Functions of time" => "functions_of_time.md",
    ],
    # doctest = false
)

deploydocs(
    repo="github.com/moble/Quaternionic.jl.git",
    devbranch="main"
)


text:
# indexnotation/sum.jl
#
# A wrapper to store the sum of a set of indexed objects and evaluate lazily,
# i.e. evaluate upon calling `deindexify`.

immutable SumOfIndexedObjects{Os<:Tuple{Vararg{AbstractIndexedObject}}} <: AbstractIndexedObject
    objects::Os
end

@inline _conjtuple(objects::Tuple{Any}) = tuple(conj(objects[1]))
@inline _conjtuple(objects::Tuple) = tuple(conj(objects[1]), _conjtuple(Base.tail(objects))...)
@inline _eltypetuple(objects::Tuple{Any}) = eltype(objects[1])
@inline _eltypetuple(objects::Tuple) = promote_type(eltype(objects[1]), _eltypetuple(Base.tail(objects)))
@inline _multuple(a, objects::Tuple{Any}) = tuple(a*objects[1])
@inline _multuple(a, objects::Tuple) = tuple(a*objects[1], _multuple(a, Base.tail(objects))...)

Base.conj(A::SumOfIndexedObjects) = SumOfIndexedObjects(_conjtuple(A.objects))

Base.eltype(A::SumOfIndexedObjects) = _eltypetuple(A.objects)

*(A::SumOfIndexedObjects, β::Number) = SumOfIndexedObjects(_multuple(β, A.objects))
*(β::Number, A::SumOfIndexedObjects) = SumOfIndexedObjects(_multuple(β, A.objects))
-(A::SumOfIndexedObjects) = *(-1, A)

indices{Os}(A::SumOfIndexedObjects{Os}) = indices(Os.parameters[1])

function +(A::SumOfIndexedObjects, B::SumOfIndexedObjects)
    add_indices(indices(A), indices(B)) # performs index check
    SumOfIndexedObjects(tuple(A.objects..., B.objects...))
end
function +(A::SumOfIndexedObjects, b::AbstractIndexedObject)
    add_indices(indices(A), indices(b)) # performs index check
    SumOfIndexedObjects(tuple(A.objects..., b))
end
function +(a::AbstractIndexedObject, B::SumOfIndexedObjects)
    add_indices(indices(a), indices(B)) # performs index check
    SumOfIndexedObjects(tuple(a, B.objects...))
end
function +(a::AbstractIndexedObject, b::AbstractIndexedObject)
    add_indices(indices(a), indices(b)) # performs index check
    SumOfIndexedObjects(tuple(a, b))
end

-(a::AbstractIndexedObject, b::AbstractIndexedObject) = +(a, -b)

@generated function deindexify{Os}(A::SumOfIndexedObjects{Os}, I::Indices)
    addex = Expr(:block, [:(deindexify!(dst, A.objects[$j], I, +1)) for j=2:length(Os.parameters)]...)
    quote
        dst = deindexify(A.objects[1], I, eltype(A))
        $addex
        return dst
    end
end

@generated function deindexify!{Os}(dst, A::SumOfIndexedObjects{Os}, I::Indices, β=0)
    addex = Expr(:block, [:(deindexify!(dst, A.objects[$j], I, +1)) for j=2:length(Os.parameters)]...)
    quote
        deindexify!(dst, A.objects[1], I, β)
        $addex
        return dst
    end
end


text:
immutable FidMtx
  MVxxBVxx   :: Array{Float64,1}
  MVxxBTxx   :: SparseMatrixCSC{Float64,Int64}
  MVxzBVxz   :: Array{Float64,1}
  MVxzBTxz   :: SparseMatrixCSC{Float64,Int64}
  MVzxBVzx   :: Array{Float64,1}
  MVzxBTxz   :: SparseMatrixCSC{Float64,Int64}
  MVzzBVzz   :: Array{Float64,1}
  MVzzBTzz   :: SparseMatrixCSC{Float64,Int64}
  MTxxxBVx   :: SparseMatrixCSC{Float64,Int64}
  MTxxxBTxxx :: Array{Float64,1}
  MTxxzBVz   :: SparseMatrixCSC{Float64,Int64}
  MTxxzBTxxz :: Array{Float64,1}
  MTzzxBVx   :: SparseMatrixCSC{Float64,Int64}
  MTzzxBTzzx :: Array{Float64,1}
  MTzzzBVz   :: SparseMatrixCSC{Float64,Int64}
  MTzzzBTzzz :: Array{Float64,1}
  MTxzxBVz   :: SparseMatrixCSC{Float64,Int64}
  MTxzxBTxzx :: Array{Float64,1}
  MTxzzBVx   :: SparseMatrixCSC{Float64,Int64}
  MTxzzBTxzz :: Array{Float64,1}
end

function vel2lm(vp::Array{Float64,2}, vs::Array{Float64,2}, rho::Array{Float64,2})
    (m, n) = size(vp)
    lambda = zeros(typeof(vp[1]), m, n)
    mu     = zeros(typeof(vp[1]), m, n)
    for j = 1: n
        for i = 1: m
            mu[i,j] = rho[i,j] * vs[i,j]^2
            lambda[i,j] = rho[i,j] * vp[i,j]^2 - 2*mu[i,j]
        end
    end
    return lambda, mu
end

function lm2vel(lambda::Array{Float64,2}, mu::Array{Float64,2}, rho::Array{Float64,2})
    (m, n) = size(lambda)
    vp = Array(typeof(lambda[1]), m, n)
    vs = Array(typeof(lambda[1]), m, n)
    for j = 1: n
        for i = 1: m
            vp[i,j] = sqrt((lambda[i,j]+2*mu[i,j])/rho[i,j])
            vs[i,j] = sqrt(mu[i,j]/rho[i,j])
        end
    end
    return vp, vs
end

function DispStable!(vmax::Float64, vmin::Float64, f0::Float64, dz::Float64, dx::Float64, dt::Float64)
    h = vmin / 5 / f0
    d = minimum([dz dx])
    tt = 6*d / (7*sqrt(3)*vmax)
    if dx >= h || dz >= h || dt >= tt
       println("maximum acceptable grid size: $h")
       println("maximum acceptable time step: $tt")
       error("unstable or frequency dispersion")
    end
    return nothing
end

function MVxxBVxx(dVxx::SparseMatrixCSC{Float64,Int64}, rho::Array{Float64,2}, dt::Float64)
    (m,n) = size(dVxx)
    Mdiag = zeros(m*n)
    for ix = 1 : n-1
        for iz = 1 : m
            a3 = (rho[iz,ix+1]+rho[iz,ix])/2/dt + dVxx[iz,ix]/2
            tmp= (rho[iz,ix+1]+rho[iz,ix])/2/(dt*a3) - dVxx[iz,ix]/(2*a3)
            Mdiag[(ix-1)*m+iz] = tmp
        end
    end
    for iz = 1 : m
        a3 = rho[iz,n]/dt + dVxx[iz,n]/2
        tmp= rho[iz,n]/(dt*a3) - dVxx[iz,n]/(2*a3)
        Mdiag[(n-1)*m+iz] = tmp
    end
    M = Mdiag
    return M
end

function MVxxBTxx(dVxx::SparseMatrixCSC{Float64,Int64}, rho::Array{Float64,2}, dt::Float64, dx::Float64)
    a1 =  9/8
    a2 = -1/24
    (m, n) = size(dVxx)
    A3 = zeros(m*n)
    for ix = 1:n-1
        for iz = 1:m
            A3[(ix-1)*m+iz] = 1 / ((rho[iz,ix+1]+rho[iz,ix])/2/dt + dVxx[iz,ix]/2)
        end
    end
    for iz = 1 : m
        A3[(n-1)*m+iz] =  1 / (rho[iz,n]/dt + dVxx[iz,n]/2)
    end
    A = spzeros(n, n)
    A[1, 1] = -a1/dx; A[1, 2] = a1/dx; A[1, 3] = a2/dx;
    for ix = 2 : n-2
        A[ix, ix-1] = -a2/dx
        A[ix, ix  ] = -a1/dx
        A[ix, ix+1] =  a1/dx
        A[ix, ix+2] =  a2/dx
    end
    A[n-1, n-2] = -a2/dx;
    A[n-1, n-1] = -a1/dx; A[n-1, n  ] = a1/dx;
    A[n  , n-1] = -a2/dx; A[n  , n  ] = -a1/dx;
    Im = speye(m)
    M = kron(A, Im)
    M = spdiagm(A3) * M
    return M
end

function MVxzBVxz(dVxz::SparseMatrixCSC{Float64,Int64}, rho::Array{Float64,2}, dt::Float64)
    (m,n) = size(dVxz)
    Mdiag = zeros(m*n)
    for ix = 1 : n-1
        for iz = 1 : m
            a3 = (rho[iz,ix+1]+rho[iz,ix])/2/dt + dVxz[iz,ix]/2
            tmp= (rho[iz,ix+1]+rho[iz,ix])/2/(dt*a3) - dVxz[iz,ix]/(2*a3)
            Mdiag[(ix-1)*m+iz] = tmp
        end
    end
    for iz = 1 : m
        a3 = rho[iz,n]/dt + dVxz[iz,n]/2
        tmp= rho[iz,n]/(dt*a3) - dVxz[iz,n]/(2*a3)
        Mdiag[(n-1)*m+iz] = tmp
    end
    M = Mdiag
    return M
end

function MVxzBTxz(dVxz::SparseMatrixCSC{Float64, Int64}, rho::Array{Float64,2}, dt::Float64, dz::Float64)
    a1 =  9/8
    a2 = -1/24
    (m,n) = size(dVxz)
    A3 = zeros(m*n)
    for ix = 1:n-1
        for iz = 1:m
            A3[(ix-1)*m+iz] = 1 / ((rho[iz,ix+1]+rho[iz,ix])/2/dt + dVxz[iz,ix]/2)
        end
    end
    for iz = 1 : m
        A3[(n-1)*m+iz] =  1 / (rho[iz,n]/dt + dVxz[iz,n]/2)
    end
    A = spzeros(m, m)
    A[1, 1] = a1/dz; A[1, 2] = a2/dz;
    A[2, 1] =-a1/dz; A[2, 2] = a1/dz; A[2, 3] = a2/dz;
    for iz = 3 : m-1
        A[iz, iz-2] = -a2/dz
        A[iz, iz-1] = -a1/dz
        A[iz, iz  ] =  a1/dz
        A[iz, iz+1] =  a2/dz
    end
    A[m, m-2] = -a2/dz;
    A[m, m-1] = -a1/dz; A[m, m] = a1/dz;
    In = speye(n)
    M = kron(In, A)
    M = spdiagm(A3) * M
    return M
end

function MVzxBVzx(dVzx::SparseMatrixCSC{Float64,Int64}, rho::Array{Float64,2}, dt::Float64)
    (m,n) = size(dVzx)
    Mdiag = zeros(m*n)
    for ix = 1 : n
        for iz = 1 : m-1
            a3 = (rho[iz+1,ix]+rho[iz,ix])/2/dt + dVzx[iz,ix]/2
            tmp= (rho[iz+1,ix]+rho[iz,ix])/2/(dt*a3) - dVzx[iz,ix]/(2*a3)
            Mdiag[(ix-1)*m+iz] = tmp
        end
    end
    for ix = 1 : n
        a3 = rho[m,ix]/dt + dVzx[m,ix]/2
        tmp= rho[m,ix]/(dt*a3) - dVzx[m,ix]/(2*a3)
        Mdiag[(ix-1)*m+m] = tmp
    end
    M = Mdiag
    return M
end

function MVzxBTxz(dVzx::SparseMatrixCSC{Float64, Int64}, rho::Array{Float64,2}, dt::Float64, dx::Float64)
    a1 = 9/8
    a2 = -1/24
    (m,n) = size(dVzx)
    A3 = zeros(m*n)
    for ix = 1:n
        for iz = 1:m-1
            A3[(ix-1)*m+iz] = 1 / ((rho[iz+1,ix]+rho[iz,ix])/2/dt+dVzx[iz,ix]/2)
        end
    end
    for ix = 1 : n
        A3[(ix-1)*m+m] =  1 / (rho[m,ix]/dt + dVzx[m,ix]/2)
    end
    A = spzeros(n, n)
    A[1, 1] = a1/dx; A[1, 2] = a2/dx;
    A[2, 1] =-a1/dx; A[2, 2] = a1/dx;
    A[2, 3] = a2/dx;
    for ix = 3:n-1
        A[ix, ix-2] = -a2/dx
        A[ix, ix-1] = -a1/dx
        A[ix, ix  ] =  a1/dx
        A[ix, ix+1] =  a2/dx
    end
    A[n, n-2] =-a2/dx;
    A[n, n-1] =-a1/dx; A[n, n] = a1/dx;
    Im = speye(m)
    M = kron(A, Im)
    M = spdiagm(A3) * M
    return M
end

function MVzzBVzz(dVzz::SparseMatrixCSC{Float64,Int64}, rho::Array{Float64,2}, dt::Float64)
    (m,n) = size(dVzz)
    Mdiag = zeros(m*n)
    for ix = 1 : n
        for iz = 1 : m-1
            a3 = (rho[iz+1,ix]+rho[iz,ix])/2/dt + dVzz[iz,ix]/2
            tmp= (rho[iz+1,ix]+rho[iz,ix])/2/(dt*a3) - dVzz[iz,ix]/(2*a3)
            Mdiag[(ix-1)*m+iz] = tmp
        end
    end
    for ix = 1 : n
        a3 = rho[m,ix]/dt + dVzz[m,ix]/2
        tmp= rho[m,ix]/(dt*a3) - dVzz[m,ix]/(2*a3)
        Mdiag[(ix-1)*m+m] = tmp
    end
    M = Mdiag
    return M
end

function MVzzBTzz(dVzz::SparseMatrixCSC{Float64,Int64}, rho::Array{Float64,2}, dt::Float64, dz::Float64)
    a1 = 9/8
    a2 = -1/24
    (m,n) = size(dVzz)
    A3 = zeros(m*n)
    for ix = 1:n
        for iz = 1:m-1
            A3[(ix-1)*m+iz] = 1 / ((rho[iz+1,ix]+rho[iz,ix])/2/dt + dVzz[iz,ix]/2)
        end
    end
    for ix = 1 : n
        A3[(ix-1)*m+m] =  1 / (rho[m,ix]/dt + dVzz[m,ix]/2)
    end
    A = spzeros(m, m)
    A[1, 1] =-a1/dz; A[1, 2] = a1/dz; A[1, 3] = a2/dz;
    for iz = 2 : m-2
        A[iz, iz-1] = -a2/dz
        A[iz, iz  ] = -a1/dz
        A[iz, iz+1] =  a1/dz
        A[iz, iz+2] =  a2/dz
    end
    A[m-1, m-2] = -a2/dz;
    A[m-1, m-1] = -a1/dz; A[m-1, m] = a1/dz;
    A[m  , m-1] = -a2/dz; A[m  , m] = -a1/dz;
    In = speye(n)
    M = kron(In, A)
    M = spdiagm(A3) * M
    return M
end

function MTxxxBTxxx(dTxxx::SparseMatrixCSC{Float64,Int64}, dt::Float64)
    (m,n) = size(dTxxx)
    Mdiag = zeros(m*n)
    for ix = 1 : n
        for iz = 1 : m
            a3 = 1/dt + dTxxx[iz,ix]/2
            tmp= 1/(dt*a3) - dTxxx[iz,ix]/(2*a3)
            Mdiag[(ix-1)*m+iz] = tmp
        end
    end
    M = Mdiag
    return M
end

function MTxxxBVx(dTxxx::SparseMatrixCSC{Float64,Int64}, lambda::Array{Float64,2}, mu::Array{Float64,2}, dt::Float64, dx::Float64)
    a1 = 9/8
    a2 = -1/24
    (m,n) = size(dTxxx)
    A3 = zeros(m*n)
    for ix = 1 : n
        for iz = 1 : m
            A3[(ix-1)*m+iz] = (lambda[iz,ix]+2*mu[iz,ix]) / (1/dt+dTxxx[iz,ix]/2)
        end
    end
    A = spzeros(n, n)
    A[1, 1] = a1/dx; A[1, 2] = a2/dx;
    A[2, 1] =-a1/dx; A[2, 2] = a1/dx; A[2, 3] = a2/dx;
    for ix = 3:n-1
        A[ix, ix-2] = -a2/dx
        A[ix, ix-1] = -a1/dx
        A[ix, ix  ] =  a1/dx
        A[ix, ix+1] =  a2/dx
    end
    A[n, n-2] = -a2/dx;
    A[n, n-1] = -a1/dx; A[n, n] = a1/dx;
    Im = speye(m)
    M = kron(A, Im)
    M = spdiagm(A3) * M
    return M
end

function MTxxzBTxxz(dTxxz::SparseMatrixCSC{Float64,Int64}, dt::Float64)
    (m,n) = size(dTxxz)
    Mdiag = zeros(m*n)
    for ix = 1:n
        for iz = 1:m
            a3 = 1/dt + dTxxz[iz,ix]/2
            tmp= 1/(dt*a3) - dTxxz[iz,ix]/(2*a3)
            Mdiag[(ix-1)*m+iz] = tmp
        end
    end
    M = Mdiag
    return M
end

function MTxxzBVz(dTxxz::SparseMatrixCSC{Float64,Int64}, lambda::Array{Float64,2}, dt::Float64, dz::Float64)
    a1 = 9/8
    a2 = -1/24
    (m,n) = size(dTxxz)
    A3 = zeros(m*n)
    for ix = 1:n
        for iz = 1:m
            A3[(ix-1)*m+iz] = lambda[iz,ix] / (1/dt+dTxxz[iz,ix]/2)
        end
    end
    A = spzeros(m, m)
    A[1, 1] = a1/dz; A[1, 2] = a2/dz;
    A[2, 1] =-a1/dz; A[2, 2] = a1/dz; A[2, 3] = a2/dz;
    for iz = 3:m-1
        A[iz, iz-2] = -a2/dz
        A[iz, iz-1] = -a1/dz
        A[iz, iz  ] =  a1/dz
        A[iz, iz+1] =  a2/dz
    end
    A[m, m-2] = -a2/dz;
    A[m, m-1] = -a1/dz; A[m, m] = a1/dz;
    In = speye(n)
    M = kron(In, A)
    M = spdiagm(A3) * M
    return M
end

function MTzzxBTzzx(dTzzx::SparseMatrixCSC{Float64,Int64}, dt::Float64)
    (m,n) = size(dTzzx)
    Mdiag = zeros(m*n)
    for ix = 1 : n
        for iz = 1 : m
            a3 = 1/dt + dTzzx[iz,ix]/2
            tmp= 1/(dt*a3) - dTzzx[iz,ix]/(2*a3)
            Mdiag[(ix-1)*m+iz] = tmp
        end
    end
    M = Mdiag
    return M
end

function MTzzxBVx(dTzzx::SparseMatrixCSC{Float64,Int64}, lambda::Array{Float64,2}, dt::Float64, dx::Float64)
    a1 = 9/8
    a2 =-1/24
    (m, n) = size(dTzzx)
    A3 = zeros(m*n)
    for ix = 1:n
        for iz = 1:m
            A3[(ix-1)*m+iz] = lambda[iz,ix] / (1/dt+dTzzx[iz,ix]/2)
        end
    end
    A = spzeros(n, n)
    A[1, 1] = a1/dx; A[1, 2] = a2/dx;
    A[2, 1] =-a1/dx; A[2, 2] = a1/dx; A[2, 3] = a2/dx;
    for ix = 3:n-1
        A[ix, ix-2] = -a2/dx
        A[ix, ix-1] = -a1/dx
        A[ix, ix  ] =  a1/dx
        A[ix, ix+1] =  a2/dx
    end
    A[n, n-2] = -a2/dx;
    A[n, n-1] = -a1/dx; A[n, n] = a1/dx;
    Im = speye(m)
    M = kron(A, Im)
    M = spdiagm(A3) * M
    return M
end

function MTzzzBTzzz(dTzzz::SparseMatrixCSC{Float64,Int64}, dt::Float64)
    (m,n) = size(dTzzz)
    Mdiag = zeros(m*n)
    for ix = 1:n
        for iz = 1:m
            a3 = 1/dt + dTzzz[iz,ix]/2
            tmp= 1/(dt*a3) - dTzzz[iz,ix]/(2*a3)
            Mdiag[(ix-1)*m+iz] = tmp
        end
    end
    M = Mdiag
    return M
end

function MTzzzBVz(dTzzz::SparseMatrixCSC{Float64,Int64}, lambda::Array{Float64,2}, mu::Array{Float64,2}, dt::Float64, dz::Float64)
    a1 = 9/8
    a2 = -1/24
    (m,n) = size(dTzzz)
    A3 = zeros(m*n)
    for ix = 1:n
        for iz = 1:m
            A3[(ix-1)*m+iz] = (lambda[iz,ix]+2*mu[iz,ix]) / (1/dt+dTzzz[iz,ix]/2)
        end
    end
    A = spzeros(m, m)
    A[1, 1] = a1/dz; A[1, 2] = a2/dz;
    A[2, 1] =-a1/dz; A[2, 2] = a1/dz; A[2, 3] = a2/dz;
    for iz = 3 : m-1
        A[iz, iz-2] = -a2/dz
        A[iz, iz-1] = -a1/dz
        A[iz, iz  ] =  a1/dz
        A[iz, iz+1] =  a2/dz
    end
    A[m, m-2] = -a2/dz;
    A[m, m-1] = -a1/dz; A[m, m] = a1/dz;
    In = speye(n)
    M = kron(In, A)
    M = spdiagm(A3) * M
    return M
end

function MTxzxBTxzx(dTxzx::SparseMatrixCSC{Float64,Int64}, dt::Float64)
    (m,n) = size(dTxzx)
    Mdiag = zeros(m*n)
    for ix = 1 : n
        for iz = 1 : m
            a3 = 1/dt + dTxzx[iz,ix]/2
            tmp= 1/(dt*a3) - dTxzx[iz,ix]/(2*a3)
            Mdiag[(ix-1)*m+iz] = tmp
        end
    end
    M = Mdiag
    return M
end

function MTxzxBVz(dTxzx::SparseMatrixCSC{Float64,Int64}, mu::Array{Float64,2}, dt::Float64, dx::Float64)
    a1 = 9/8
    a2 =-1/24
    (m, n) = size(dTxzx)
    A3 = zeros(m*n)
    for ix = 1 : n-1
        for iz = 1 : m-1
            A3[(ix-1)*m+iz] = (mu[iz,ix]+mu[iz+1,ix]+mu[iz,ix+1]+mu[iz+1,ix+1])/4 / (1/dt+dTxzx[iz,ix]/2)
        end
    end
    for iz = 1 : m-1
        A3[(n-1)*m+iz] = (mu[iz,n]+mu[iz+1,n])/2 / (1/dt+dTxzx[iz,n]/2)
    end
    for ix = 1 : n-1
        A3[(ix-1)*m+m] = (mu[m,ix]+mu[m,ix+1])/2 / (1/dt+dTxzx[m,ix]/2)
    end
    A3[end] = mu[m,n] / (1/dt+dTxzx[m,n]/2)
    A = spzeros(n, n)
    A[1,1] =-a1/dx; A[1,2] = a1/dx; A[1,3] = a2/dx;
    for ix = 2 : n-2
        A[ix, ix-1] = -a2/dx
        A[ix, ix  ] = -a1/dx
        A[ix, ix+1] =  a1/dx
        A[ix, ix+2] =  a2/dx
    end
    A[n-1,n-2] =-a2/dx;
    A[n-1,n-1] =-a1/dx; A[n-1,n] = a1/dx;
    A[n  ,n-1] =-a2/dx; A[n  ,n] =-a1/dx;
    Im = speye(m)
    M = kron(A, Im)
    M = spdiagm(A3) * M
    return M
end

function MTxzzBTxzz(dTxzz::SparseMatrixCSC{Float64,Int64}, dt::Float64)
    (m, n) = size(dTxzz)
    Mdiag  = zeros(m*n)
    for ix = 1:n
        for iz = 1:m
            a3 = 1/dt + dTxzz[iz,ix]/2
            tmp= 1/(dt*a3) - dTxzz[iz,ix]/(2*a3)
            Mdiag[(ix-1)*m+iz] = tmp
        end
    end
    M = Mdiag
    return M
end

function MTxzzBVx(dTxzz::SparseMatrixCSC{Float64,Int64}, mu::Array{Float64,2}, dt::Float64, dz::Float64)
    a1 = 9/8
    a2 =-1/24
    (m,n) = size(dTxzz)
    A3 = zeros(m*n)
    for ix = 1 : n-1
        for iz = 1 : m-1
            A3[(ix-1)*m+iz] = (mu[iz,ix]+mu[iz+1,ix]+mu[iz,ix+1]+mu[iz+1,ix+1])/4 / (1/dt+dTxzz[iz,ix]/2)
        end
    end
    for iz = 1 : m-1
        A3[(n-1)*m+iz] = (mu[iz,n]+mu[iz+1,n])/2 / (1/dt+dTxzz[iz,n]/2)
    end
    for ix = 1 : n-1
        A3[(ix-1)*m+m] = (mu[m,ix]+mu[m,ix+1])/2 / (1/dt+dTxzz[m,ix]/2)
    end
    A3[end] = mu[m,n] / (1/dt+dTxzz[m,n]/2)
    A = spzeros(m, m)
    A[1,1] =-a1/dz; A[1,2] = a1/dz; A[1,3] = a2/dz;
    for iz = 2:m-2
        A[iz, iz-1] = -a2/dz
        A[iz, iz  ] = -a1/dz
        A[iz, iz+1] =  a1/dz
        A[iz, iz+2] =  a2/dz
    end
    A[m-1,m-2] =-a2/dz;
    A[m-1,m-1] =-a1/dz; A[m-1,m] = a1/dz;
    A[m  ,m-1] =-a2/dz; A[m  ,m] =-a1/dz;
    In = speye(n)
    M = kron(In, A)
    M = spdiagm(A3) * M
    return M
end

function CreateFidMtx(nz::Int64, nx::Int64, ext::Int64, iflag::Int64, vp::Array{Float64,2}, vs::Array{Float64,2}, rho::Array{Float64,2}, dz::Float64, dx::Float64, dt::Float64, f0::Float64)
    vmax = maximum(vp); vmin = minimum(vs);
    DispStable!(vmax, vmin, f0, dz, dx, dt)
    (lambda, mu) = vel2lm(vp, vs, rho)
    (lambda, mu, rho) = modelPadding(lambda, mu, rho, ext, iflag)
    dpCoef = DampBound(nz, nx, vmax, ext, rho, iflag, dz)
    MVxxBVxx1 = MVxxBVxx(dpCoef.dVxx, rho, dt)
    MVxxBTxx1 = MVxxBTxx(dpCoef.dVxx, rho, dt, dx)
    MVxzBVxz1 = MVxzBVxz(dpCoef.dVxz, rho, dt)
    MVxzBTxz1 = MVxzBTxz(dpCoef.dVxz, rho, dt, dz)
    MVzxBVzx1 = MVzxBVzx(dpCoef.dVzx, rho, dt)
    MVzxBTxz1 = MVzxBTxz(dpCoef.dVzx, rho, dt, dx)
    MVzzBVzz1 = MVzzBVzz(dpCoef.dVzz, rho, dt)
    MVzzBTzz1 = MVzzBTzz(dpCoef.dVzz, rho, dt, dz)
    MTxxxBVx1   =   MTxxxBVx(dpCoef.dTxxx, lambda, mu, dt, dx)
    MTxxxBTxxx1 = MTxxxBTxxx(dpCoef.dTxxx, dt)
    MTxxzBVz1   =   MTxxzBVz(dpCoef.dTxxz, lambda, dt, dz)
    MTxxzBTxxz1 = MTxxzBTxxz(dpCoef.dTxxz, dt)
    MTzzxBVx1   =   MTzzxBVx(dpCoef.dTzzx, lambda, dt, dx)
    MTzzxBTzzx1 = MTzzxBTzzx(dpCoef.dTzzx, dt)
    MTzzzBVz1   =   MTzzzBVz(dpCoef.dTzzz, lambda, mu, dt, dz)
    MTzzzBTzzz1 = MTzzzBTzzz(dpCoef.dTzzz, dt)
    MTxzxBVz1   =   MTxzxBVz(dpCoef.dTxzx, mu, dt, dx)
    MTxzxBTxzx1 = MTxzxBTxzx(dpCoef.dTxzx, dt)
    MTxzzBVx1   =   MTxzzBVx(dpCoef.dTxzz, mu, dt, dz)
    MTxzzBTxzz1 = MTxzzBTxzz(dpCoef.dTxzz, dt)
    fidMtx = FidMtx(MVxxBVxx1, MVxxBTxx1,   MVxzBVxz1, MVxzBTxz1,
                    MVzxBVzx1, MVzxBTxz1,   MVzzBVzz1, MVzzBTzz1,
                    MTxxxBVx1, MTxxxBTxxx1, MTxxzBVz1, MTxxzBTxxz1,
                    MTzzxBVx1, MTzzxBTzzx1, MTzzzBVz1, MTzzzBTzzz1,
                    MTxzxBVz1, MTxzxBTxzx1, MTxzzBVx1, MTxzzBTxzz1)
    return fidMtx
end

max_stars_repo_name: WneleiGao/ElasticWave


text:
# call order:
# TICON/test/runtests.jl
#  -> TICON/test/advection/runtests.jl

# This test suite is intended to test the advection physics module's
#   implementation of eqn_deepcopy

# 5 tests:
#   1. ensure eqn_copy matches eqn
#   2. ensure eqn_copy pointers do not equal eqn's
#   3. ensure that changes to eqn do not affect eqn_copy
#   4. ensure that changes to eqn_copy.q also change q_vec
#   5. ensure that changes to eqn_copy.res also change res_vec

global const test_eqn_copy_inputfile = "input_vals_channel_DG_foreqncopy.jl"

function test_eqn_copy(mesh, sbp, eqn, opts)

  # should we step through RK4 or CN a few steps before trying eqn_deepcopy?

  eqn_copy = eqn_deepcopy(mesh, sbp, eqn, opts)

  @testset "--- Testing eqn_deepcopy, DG q/q_vec & res/res_vec equality ---" begin

    @test ( mesh.isDG )== true
    @test ( (pointer(eqn.q) == pointer(eqn.q_vec)) )== true
    @test ( (pointer(eqn.res) == pointer(eqn.res_vec)) )== true
    @test ( (pointer(eqn_copy.q) == pointer(eqn_copy.q_vec)) )== true
    @test ( (pointer(eqn_copy.res) == pointer(eqn_copy.res_vec)) )== true

  end

  @testset "--- Testing eqn_deepcopy, copy phase, values ---" begin

    @test ( (eqn_copy.commsize == eqn.commsize) )== true
    @test ( (eqn_copy.myrank == eqn.myrank) )== true
    @test ( (eqn_copy.t == eqn.t) )== true
    @test ( (eqn_copy.res_type == eqn.res_type) )== true
    @test ( (eqn_copy.multiplyA0inv == eqn.multiplyA0inv) )== true
    @test ( (eqn_copy.params.LFalpha == eqn.params.LFalpha) )== true
    @test ( (eqn_copy.params.alpha_x == eqn.params.alpha_x) )== true
    @test ( (eqn_copy.params.alpha_y == eqn.params.alpha_y) )== true
    @test ( (eqn_copy.params.alpha_z == eqn.params.alpha_z) )== true
    # commented out since this is an artifact of some uadj adv BC code used to test.
    # @test ( (eqn_copy.params.sin_amplitude == eqn.params.sin_amplitude) )
    @test ( (eqn_copy.params.f == eqn.params.f) )== true
    @test ( (eqn_copy.params.time == eqn.params.time) )== true
    @test ( (eqn_copy.src_func == eqn.src_func) )== true
    @test ( (eqn_copy.flux_func == eqn.flux_func) )== true
    @test ( (eqn_copy.majorIterationCallback == eqn.majorIterationCallback) )== true

    # arrays
    @test ( (eqn_copy.q_face[1] == eqn.q_face[1]) )== true
    @test ( (eqn_copy.flux_parametric[1] == eqn.flux_parametric[1]) )== true
    @test ( (eqn_copy.flux_face[1] == eqn.flux_face[1]) )== true
    @test ( (eqn_copy.q_bndry[1] == eqn.q_bndry[1]) )== true
    @test ( (eqn_copy.bndryflux[1] == eqn.bndryflux[1]) )== true
    @test ( (eqn_copy.M[1] == eqn.M[1]) )== true
    @test ( (eqn_copy.Minv[1] == eqn.Minv[1]) )== true
    @test ( (eqn_copy.Minv3D[1] == eqn.Minv3D[1]) )== true
    @test ( (eqn_copy.q[1] == eqn.q[1]) )== true
    @test ( (eqn_copy.q_vec[1] == eqn.q_vec[1]) )== true
    @test ( (eqn_copy.res[1] == eqn.res[1]) )== true
    @test ( (eqn_copy.res_vec[1] == eqn.res_vec[1]) )== true

    # Zero length arrays?
    if length(eqn.aux_vars) > 0
      @test ( (eqn_copy.aux_vars[1] == eqn.aux_vars[1]) )== true
    end
    if length(eqn.res_edge) > 0
      @test ( (eqn_copy.res_edge[1] == eqn.res_edge[1]) )== true
    end
    # note: q_face_send and q_face_recv are now in ParallelData, see jcrean for details.
    #       Their purpose is now served by eqn.shared_data
    if length(eqn.shared_data) > 0
      @test ( (eqn_copy.shared_data[1] == eqn.shared_data[1]) )== true
    end
    if length(eqn.flux_sharedface) > 0
      @test ( (eqn_copy.flux_sharedface[1] == eqn.flux_sharedface[1]) )== true
    end

  end

  @testset "--- Testing eqn_deepcopy, copy phase, pointers of arrays ---" begin

    @test ( (pointer(eqn_copy.q_face) == pointer(eqn.q_face)) )== false
    @test ( (pointer(eqn_copy.aux_vars) == pointer(eqn.aux_vars)) )== false
    @test ( (pointer(eqn_copy.flux_parametric) == pointer(eqn.flux_parametric)) )== false
    @test ( (pointer(eqn_copy.flux_face) == pointer(eqn.flux_face)) )== false
    @test ( (pointer(eqn_copy.res_edge) == pointer(eqn.res_edge)) )== false
    @test ( (pointer(eqn_copy.q_bndry) == pointer(eqn.q_bndry)) )== false
    @test ( (pointer(eqn_copy.bndryflux) == pointer(eqn.bndryflux)) )== false
    @test ( (pointer(eqn_copy.M) == pointer(eqn.M)) )== false
    @test ( (pointer(eqn_copy.Minv) == pointer(eqn.Minv)) )== false
    @test ( (pointer(eqn_copy.Minv3D) == pointer(eqn.Minv3D)) )== false
    # note: q_face_send and q_face_recv are now in ParallelData, see jcrean for details.
    #       Their purpose is now served by eqn.shared_data
    @test ( (pointer(eqn_copy.shared_data) == pointer(eqn.shared_data)) )== false
    @test ( (pointer(eqn_copy.flux_sharedface) == pointer(eqn.flux_sharedface)) )== false
    @test ( (pointer(eqn_copy.q) == pointer(eqn.q)) )== false
    @test ( (pointer(eqn_copy.q_vec) == pointer(eqn.q_vec)) )== false
    @test ( (pointer(eqn_copy.res) == pointer(eqn.res)) )== false
    @test ( (pointer(eqn_copy.res_vec) == pointer(eqn.res_vec)) )== false

  end

  rand!(eqn.q_face)
  rand!(eqn.flux_parametric)
  rand!(eqn.flux_face)
  rand!(eqn.q_bndry)
  rand!(eqn.bndryflux)
  rand!(eqn.M)
  rand!(eqn.Minv)
  rand!(eqn.Minv3D)
  rand!(eqn.q)
  rand!(eqn.res)

  @testset "--- Testing eqn_deepcopy, eqn change phase ---" begin

    @test ( (eqn_copy.q_face[1] == eqn.q_face[1]) )== false
    @test ( (eqn_copy.flux_parametric[1] == eqn.flux_parametric[1]) )== false
    @test ( (eqn_copy.flux_face[1] == eqn.flux_face[1]) )== false
    @test ( (eqn_copy.q_bndry[1] == eqn.q_bndry[1]) )== false
    @test ( (eqn_copy.bndryflux[1] == eqn.bndryflux[1]) )== false
    @test ( (eqn_copy.M[1] == eqn.M[1]) )== false
    @test ( (eqn_copy.Minv[1] == eqn.Minv[1]) )== false
    @test ( (eqn_copy.Minv3D[1] == eqn.Minv3D[1]) )== false
    @test ( (eqn_copy.q[1] == eqn.q[1]) )== false
    @test ( (eqn_copy.q_vec[1] == eqn.q_vec[1]) )== false
    @test ( (eqn_copy.res[1] == eqn.res[1]) )== false
    @test ( (eqn_copy.res_vec[1] == eqn.res_vec[1]) )== false

  end

  rand!(eqn_copy.q)
  rand!(eqn_copy.res)

  @testset "--- Testing eqn_deepcopy, eqn_copy change phase ---" begin

    @test ( (eqn_copy.q[1] == eqn_copy.q_vec[1]) )== true
    @test ( (eqn_copy.res[1] == eqn_copy.res_vec[1]) )== true

  end

end   # end of function test_eqn_copy

add_func2!(AdvectionTests, test_eqn_copy, test_eqn_copy_inputfile, [TAG_SHORTTEST])

#=
Fields of advection eqn object:

  commsize
  myrank
  t
  res_type
  multiplyA0inv
  params.LFalpha
  params.alpha_x
  params.alpha_y
  params.alpha_z
  params.sin_amplitude
  params.f
  params.time
  src_func
  flux_func
  majorIterationCallback

  # arrays
  q_face
  aux_vars
  flux_parametric
  flux_face
  res_edge
  q_bndry
  bndryflux
  M
  Minv
  Minv3D
  q_face_send
  q_face_recv
  flux_sharedface
  q
  q_vec
  res
  res_vec

=#


text:
using DiffEqBase, DataStructures

import DiffEqBase: ODEProblem, solve

const warnkeywords =
    (:save_idxs, :d_discontinuities, :unstable_check, :save_everystep,
     :save_end, :initialize_save, :adaptive, :dt, :reltol, :dtmax,
     :dtmin, :force_dtmin, :internalnorm, :gamma, :beta1, :beta2,
     :qmax, :qmin, :qsteady_min, :qsteady_max, :qoldinit, :failfactor,
     :maxiters, :isoutofdomain, :unstable_check,
     :calck, :progress, :timeseries_steps, :tstops, :dense)

global warnlist = Set(warnkeywords)

abstract type TaylorAlgorithm <: DiffEqBase.DEAlgorithm end
struct TaylorMethod <: TaylorAlgorithm
    order::Int
end

TaylorMethod() = error("Maximum order must be specified for the Taylor method")

export TaylorMethod

function DiffEqBase.solve(
    prob::DiffEqBase.AbstractODEProblem{uType,tupType,isinplace},
    alg::AlgType,
    timeseries=[],ts=[],ks=[];
    verbose=true,
    saveat=eltype(tupType)[],
    abstol = 1e-6,
    save_start =  isempty(saveat) || saveat isa Number || prob.tspan[1] in saveat,
    save_end   =  isempty(saveat) || saveat isa Number || prob.tspan[2] in saveat,
    timeseries_errors=true, maxiters = 1_000_000,
    callback=nothing, kwargs...) where
        {uType, tupType, isinplace, AlgType <: TaylorAlgorithm}

    tType = eltype(tupType)

    if verbose
        warned = !isempty(kwargs) && check_keywords(alg, kwargs, warnlist)
        warned && warn_compat()
    end

    if haskey(prob.kwargs, :callback) || haskey(kwargs, :callback) || !isa(callback, Nothing)
        error("TaylorIntegration is not compatible with callbacks.")
    end

    sizeu = size(prob.u0)
    f = prob.f

    tspan = prob.tspan

    _, saveat_vec_, _ = tstop_saveat_disc_handling(1,saveat,1,tspan)
    saveat_vec = sort(saveat_vec_.valtree)

    if !isempty(saveat) && save_end && saveat_vec[end] != tspan[2]
        push!(saveat_vec,tspan[2])
    end
    if !isempty(saveat) && saveat_vec[1] != tspan[1]
        prepend!(saveat_vec,tspan[1])
    end

    if !isinplace && typeof(prob.u0) <: Vector{Float64}
        f! = (du, u, p, t) -> (du .= f(u, p, t); 0)
        if isempty(saveat_vec)
            t, vectimeseries = taylorinteg(f!, prob.u0, prob.tspan[1], prob.tspan[2],
                alg.order, abstol, prob.p, maxsteps=maxiters, parse_eqs=false)
        else
            t = saveat_vec
            vectimeseries = taylorinteg(f!, prob.u0, t,
                alg.order, abstol, prob.p, maxsteps=maxiters, parse_eqs=false)
        end

    elseif !isinplace && typeof(prob.u0) <: AbstractArray
        f! = (du, u, p, t) -> (du .= vec(f(reshape(u, sizeu), p, t)); 0)
        if isempty(saveat_vec)
            t, vectimeseries = taylorinteg(f!, prob.u0, prob.tspan[1], prob.tspan[2],
                alg.order, abstol, prob.p, maxsteps=maxiters, parse_eqs=false)
        else
            t = saveat_vec
            vectimeseries = taylorinteg(f!, prob.u0, t,
                alg.order, abstol, prob.p, maxsteps=maxiters, parse_eqs=false)
        end

    elseif isinplace && typeof(prob.u0) <: AbstractArray{eltype(uType),2}
        f! = (du, u, p, t) -> (
                dd = reshape(du, sizeu); uu = reshape(u, sizeu);
                f(dd, uu, p, t); u = vec(uu); du = vec(dd); 0)
        u0 = vec(prob.u0)
        if isempty(saveat_vec)
            t, vectimeseries = taylorinteg(f!, u0, prob.tspan[1], prob.tspan[2],
                alg.order, abstol, prob.p, maxsteps=maxiters, parse_eqs=false)
        else
            t = saveat_vec
            vectimeseries = taylorinteg(f!, prob.u0, t,
                alg.order, abstol, prob.p, maxsteps=maxiters, parse_eqs=false)
        end
    else
        if haskey(kwargs, :parse_eqs)
            parse_eqs = kwargs[:parse_eqs]
        else
            parse_eqs = true
        end
        if isempty(saveat_vec)
            t, vectimeseries = taylorinteg(f, prob.u0, prob.tspan[1],
                prob.tspan[2], alg.order, abstol, prob.p, maxsteps=maxiters,
                parse_eqs=parse_eqs)
        else
            t = saveat_vec
            vectimeseries = taylorinteg(f, prob.u0, t,
                alg.order, abstol, prob.p, maxsteps=maxiters,
                parse_eqs=parse_eqs)
        end
    end

    if save_start
      start_idx = 1
      _t = t
    else
      start_idx = 2
      _t = t[2:end]
    end

    if typeof(prob.u0) <: AbstractArray
      _timeseries = Vector{uType}(undef, 0)
      for i=start_idx:size(vectimeseries, 1)
          push!(_timeseries, reshape(view(vectimeseries, i, :, )', sizeu))
      end
    else
      _timeseries = vec(vectimeseries)[start_idx:end]
    end

    DiffEqBase.build_solution(prob,  alg, _t, _timeseries,
        timeseries_errors = timeseries_errors, retcode = :Success)
end

# This function is part of OrdinaryDiffEq.jl; MIT-licensed
# https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/8b5af68ee4126d3772fe690a7128468f1334b4d6/src/solve.jl#L411
function tstop_saveat_disc_handling(tstops, saveat, d_discontinuities, tspan)
    t0, tf = tspan
    tType = eltype(tspan)
    tdir = sign(tf - t0)

    tdir_t0 = tdir * t0
    tdir_tf = tdir * tf

    # time stops
    tstops_internal = BinaryMinHeap{tType}()
    if isempty(d_discontinuities) && isempty(tstops) # TODO: Specialize more
      push!(tstops_internal, tdir_tf)
    else
      for t in tstops
        tdir_t = tdir * t
        tdir_t0 < tdir_t ≤ tdir_tf && push!(tstops_internal, tdir_t)
      end

      for t in d_discontinuities
        tdir_t = tdir * t
        tdir_t0 < tdir_t ≤ tdir_tf && push!(tstops_internal, tdir_t)
      end

      push!(tstops_internal, tdir_tf)
    end

    # saving time points
    saveat_internal = BinaryMinHeap{tType}()
    if typeof(saveat) <: Number
      if (t0:saveat:tf)[end] == tf
        for t in (t0 + saveat):saveat:tf
          push!(saveat_internal, tdir * t)
        end
      else
        for t in (t0 + saveat):saveat:(tf - saveat)
          push!(saveat_internal, tdir * t)
        end
      end
    elseif !isempty(saveat)
      for t in saveat
        tdir_t = tdir * t
        tdir_t0 < tdir_t < tdir_tf && push!(saveat_internal, tdir_t)
      end
    end

    # discontinuities
    d_discontinuities_internal = BinaryMinHeap{tType}()
    sizehint!(d_discontinuities_internal.valtree, length(d_discontinuities))
    for t in d_discontinuities
      push!(d_discontinuities_internal, tdir * t)
    end

    tstops_internal, saveat_internal, d_discontinuities_internal
end


text:
using Documenter
using StaticArrays
using Test

using ScikitSpatial

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

include("types/base.jl")
include("computations/base.jl")


text:
function odeint(y0::Vector{Float64},t0::Float64,tf::Float64,sparams::WK3.SolverParams,
    mparams::WK3.ModelParams,k::Float64)
    # parameters for dependent variable scaling
    tiny = 1e-30;

    # integration options
    saveflag = true; # save at intervals ≈ dtsav

    # output variables
    tout = Float64[];
    yout = Vector{Float64}[];

    # initialization
    y = zeros(sparams.nvar);
    ys = zeros(sparams.nvar);
    dy = zeros(sparams.nvar);
    hdid = 0.;
    hnext = 0.;
    t = t0;
    y[:] = y0;
    nok = 0;
    nbad = 0;

    # make sure we're stepping in the right direction (forwards vs. backwards)
    h = (tf-t0) >= 0 ? (sparams.h0 >= 0 ? sparams.h0: -sparams.h0) : (sparams.h0 >= 0 ? -sparams.h0 : sparams.h0);

    # assure storage of first step
    if saveflag
        tsav = t-sparams.dtsav*2;
    end

    while true
        # evaluate derivatives
        # println("Derivative evaluation:")
        wk3odes(t,y,dy,mparams,k);
        # scaling to monitor accuracy
        # println("Scaling:")
        ys .= abs.(y) .+ abs.(dy.*h).+tiny;
        # store intermediate results
        if saveflag && abs.(t-tsav) > abs(sparams.dtsav)
            push!(tout,t)
            push!(yout,y[:])
            tsav = t;
        end
        # decrease stepsize if overshooting end of integration window
        if (t+h-tf)*(t+h-t0) > 0
            h = tf-t;
        end
        # attempt an integration step
        # println("RKQS:")
        t,hdid,hnext=rkqs(y,dy,t,h,hdid,hnext,ys,sparams,mparams,k);
        # check valve state, set ventricular flow accordingly
        E,~ = elastancefn(t*ts,mparams,k);
        Pv = E.*(y[2].*Vs.-mparams.V0)./Ps;
        if Pv < y[1] && y[3] <= 0.
            y[3] = 0.;
        end
        # track times step size was changed
        if hdid == h
            nok+=1;
        else
            nbad+=1;
        end
        # normal exit if finished (saving final step)
        if (t-tf)*(tf-t0) >= 0
            push!(yout,y)
            push!(tout,t)
            return tout,yout
        end
        if abs.(hnext) <= sparams.hmin
            push!(yout,y)
            push!(tout,t)
            println("Step size below minimum allowable in odeint. Aborting. hnext = $hnext. t = $t.")
            return tout,yout
        end
    end
end


text:
# TODO: Haven't made up my mind about whether or not to have a function to return
# a base header for both get_blob and put_blob

"""
	get_blob(blob, directory, container, storageaccount, storagekey)

Download file from blob storage using HTTP GET.
"""
function get_blob(blob, directory, container, storageaccount, storagekey)
	timestamp = now_http_date()
	signature = storage_signature(blob, directory, container, storageaccount, storagekey, timestamp)

    token = string("SharedKey ", storageaccount, ":", signature)

    header = [
        "Authorization" => token,
        "Content-Length" => "0",
        "x-ms-version" => X_MS_VERSION,
		"x-ms-date" => timestamp
    ]

	bloburl = blob_to_url(blob, directory, container, storageaccount)
    HTTP.get(bloburl, header)
end

"""
	put_blob(content, blob, directory, container, storageaccount,
			 storagekey, contenttype::String = "application/json")

Upload file to blob storage using HTTP PUT.
"""
function put_blob(content, blob, directory, container, storageaccount,
                  storagekey, contenttype::String = "application/json")
    # TODO: Check if content matches contenttype

	timestamp = now_http_date()
	signature = storage_signature(content, blob, directory, container, storageaccount, storagekey, timestamp, contenttype)

    token = string("SharedKey ", storageaccount, ":", signature)

    header = [
        "Authorization" => token,
		"Content-Length" => contentsize(content),
        "x-ms-version" => X_MS_VERSION,
        "x-ms-date" => timestamp,
        "x-ms-blob-type" => "Blockblob",
        "Content-type" => contenttype
    ]

	bloburl = blob_to_url(blob, directory, container, storageaccount)
    HTTP.put(bloburl, header, content)
end

"""
	blob_to_url(blob, directory, container, storageaccount)

Concatenate the components of a blob location to a URL.
"""
function blob_to_url(blob, directory, container, storageaccount)
    string(
		"https://", storageaccount,
		".blob.core.windows.net/", container, "/", directory, "/", blob
	)
end

max_stars_repo_name: robertdj/AzureBlob

max_issues_repo_name: robertdj/AzureBlob

max_forks_repo_name: robertdj/AzureBlob


text:
## Test pvalue adjustment methods ##
module Test_pval_adjustment

using MultipleTesting
using Test

@testset "p-Value adjustment" begin

    pi0 = 0.4

    # p-values without ties
    pval1 = [0.0, 0.0001, 0.001, 0.01, 0.05, 0.1, 0.2, 0.4, 0.75, 1.0]
    ref1 = Dict(
        Bonferroni         => [0.0, 0.001, 0.01, 0.1, 0.5, 1.0, 1.0, 1.0, 1.0, 1.0],
        Holm               => [0.0, 9e-4, 8e-3, 0.07, 0.3, 0.5, 0.8, 1.0, 1.0, 1.0],
        Hochberg           => [0.0, 9e-4, 8e-3, 0.07, 0.3, 0.5, 0.8, 1.0, 1.0, 1.0],
        Hommel             => [0.0, 9e-4, 8e-3, 0.07, 0.3, 0.5, 0.8, 1.0, 1.0, 1.0],
        BenjaminiHochberg  => [0.0, 0.0005, 0.003333333, 0.025, 0.1, 0.166666667, 0.285714286, 0.5, 0.833333333, 1.0],
        BenjaminiYekutieli => [0.0, 0.001464484, 0.009763228, 0.073224206, 0.292896825, 0.488161376, 0.836848073, 1.0, 1.0, 1.0],
        BenjaminiLiu       => [0.0, 0.0008096761, 0.0063776447, 0.0475542565, 0.1589448656, 0.2047550000, 0.2361600000, 0.2361600000, 0.2361600000, 0.2361600000],
        Sidak              => [0.0, 0.0009995501, 0.0099551198, 0.0956179250, 0.4012630608, 0.6513215599, 0.8926258176, 0.9939533824, 0.9999990463, 1.0000000000],
        ForwardStop        => [0.0, 0.0000500025, 0.0003668351, 0.0027877103, 0.0124888271, 0.0279674419, 0.0558497432, 0.1127217283, 0.2542297986, 1.0],
        BarberCandes       => [0.2857142857, 0.2857142857, 0.2857142857, 0.2857142857, 0.2857142857, 0.2857142857, 0.2857142857, 0.375, 1.0, 1.0]
    )

    # p-values with ties
    pval2 = [0.0001, 0.0001, 0.001, 0.01, 0.05, 0.1, 0.1, 0.4, 0.75, 1.0]
    ref2 = Dict(
        Bonferroni         => [0.001, 0.001, 0.01, 0.1, 0.5, 1.0, 1.0, 1.0, 1.0, 1.0],
        Holm               => [0.001, 0.001, 0.008, 0.07, 0.3, 0.5, 0.5, 1.0, 1.0, 1.0],
        Hochberg           => [9e-4, 9e-4, 8e-3, 0.07, 0.3, 0.4, 0.4, 1.0, 1.0, 1.0],
        Hommel             => [9e-4, 9e-4, 8e-3, 0.07, 0.25, 0.4, 0.4, 1.0, 1.0, 1.0],
        BenjaminiHochberg  => [0.0005, 0.0005, 0.003333333, 0.025, 0.1, 0.142857143, 0.142857143, 0.5, 0.833333333, 1.0],
        BenjaminiYekutieli => [0.001464484, 0.001464484, 0.009763228, 0.073224206, 0.292896825, 0.418424036, 0.418424036, 1.0, 1.0, 1.0],
        BenjaminiLiu       => [0.0009995501, 0.0009995501, 0.0063776447, 0.0475542565, 0.1589448656, 0.2047550000, 0.2047550000, 0.2352000000, 0.2352000000, 0.2352000000],
        Sidak              => [0.0009995501, 0.0009995501, 0.0099551198, 0.0956179250, 0.4012630608, 0.6513215599, 0.6513215599, 0.9939533824, 0.9999990463, 1.0000000000],
        ForwardStop        => [0.0001000050, 0.0001000050, 0.0004001701, 0.0028127115, 0.0125088281, 0.0279841094, 0.0390378817, 0.0980113495, 0.2411539063, 1.0],
        BarberCandes       => [0.2857142857, 0.2857142857, 0.2857142857, 0.2857142857, 0.2857142857, 0.2857142857, 0.2857142857, 0.375, 1.0, 1.0]
    )

    # smallest p-values from larger set
    k3 = 6
    n3 = length(pval1)
    pval3 = pval1[1:k3]
    pval3pad = [pval3; fill(1.0, n3-k3)]  # non-observed p-values equal to 1
    methods3 = [Bonferroni, Holm, Hochberg, Hommel, BenjaminiHochberg,
                BenjaminiYekutieli, BenjaminiLiu, Sidak, ForwardStop]

    @testset "pvalue adjustment $method" for method in keys(ref1)

        @test method <: PValueAdjustment
        @test typeof(method()) <: PValueAdjustment

        @test_throws MethodError method(0.1)

        ## no valid p-values as input
        @test_throws DomainError adjust([-1.0, 0.5], method())
        @test_throws DomainError adjust([0.5, 1.5], method())

        ## any single p-value is returned unchanged
        # this test is not valid for ForwardStop or Barber-Candès
        if !(method in [ForwardStop, BarberCandes])
            pval = rand(1)
            @test adjust(pval, method()) == pval
        end

        if (method == BarberCandes)
            pval = rand(1)
            @test adjust(pval, method()) == fill(1.0, size(pval))
        end

        ## compare with reference values
        @test isapprox( adjust(pval1, method()), ref1[method], atol = 1e-9 )
        @test isapprox( adjust(PValues(pval1), method()), ref1[method], atol = 1e-9 )

        # unsorted inputs
        for i in 1:10
            ord = MultipleTesting.unorder(pval1)
            @test isapprox( adjust(pval1[ord], method()), ref1[method][ord], atol = 1e-9 )
        end

        ## compare with reference values having ties
        @test isapprox( adjust(pval2, method()), ref2[method], atol = 1e-9 )
        @test isapprox( adjust(PValues(pval2), method()), ref2[method], atol = 1e-9 )

        # unsorted inputs
        # this test is not valid for ForwardStop
        if method != ForwardStop
            for i in 1:10
                ord = MultipleTesting.unorder(ref2[method])
                @test isapprox( adjust(pval2[ord], method()), ref2[method][ord], atol = 1e-9 ) # FIXME
            end
        end

        ## sorting order does not play a role
        for i in 1:10
            pval4 = sort(rand(10)) # all under H0
            ord = MultipleTesting.unorder(pval4)
            @test adjust(pval4[ord], method()) == adjust(pval4, method())[ord]
        end

        ## total number of tests explicitly specified
        if method in methods3
            # all p-values present: same reference values as for pval1
            @test adjust(pval3, length(pval3), method()) == adjust(pval3, method())
            @test adjust(PValues(pval3), length(pval3), method()) == adjust(pval3, method())

            # k smallest p-values present, total number n known
            @test adjust(pval3, n3, method()) == adjust(pval3pad, n3, method())[1:k3]

            # unsorted inputs
            for i in 1:10
                ord = MultipleTesting.unorder(pval3)
                @test adjust(pval3[ord], n3, method()) == (adjust(pval3pad, n3, method())[1:k3])[ord]
            end

            # k > n not allowed: test for any n in [1,k-1]
            @test_throws ArgumentError adjust(pval3, rand(1:k3-1), method())
        end

    end

    @testset "argument checking for number of tests" begin

        @test_throws ArgumentError MultipleTesting.check_number_tests(3, 2)
        @test MultipleTesting.check_number_tests(2, 2) == nothing
        @test MultipleTesting.check_number_tests(2, 3) == nothing

    end

    @testset "BarberCandès #2:" begin
        for k = 1:5
            pv = rand(BetaUniformMixtureModel(0.5, 0.5, 7.0), 40)
            @test isapprox(adjust(pv, BarberCandes()),
                           MultipleTesting.barber_candes_brute_force(pv), atol=1e-9)
        end
    end

    @testset "BarberCandès: All p-values < 0.5 (#87)" begin
        for pv in ([0.05, 0.1, 0.3], [0.01, 0.17, 0.25, 0.37, 0.47])
            n = length(pv)
            @test isapprox( adjust(PValues(pv), BarberCandes()), fill(1/n, n) )
        end
    end

    @testset "Step-up/down" begin

        x        = [1.0, 2.0, 2.0, 5.0, 4.0, 5.0, 7.0]
        ref_up   = [1.0, 2.0, 2.0, 4.0, 4.0, 5.0, 7.0]
        ref_down = [1.0, 2.0, 2.0, 5.0, 5.0, 5.0, 7.0]

        # step-up
        y = copy(x)
        MultipleTesting.stepup!(y)
        @test y == ref_up

        # step-down
        y = copy(x)
        MultipleTesting.stepdown!(y)
        @test y == ref_down

    end

end

end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module TestLQ

using Test, LinearAlgebra, Random
using LinearAlgebra: BlasComplex, BlasFloat, BlasReal, rmul!, lmul!

n = 10

# Split n into 2 parts for tests needing two matrices
n1 = div(n, 2)
n2 = 2*n1

srand(1234321)

areal = randn(n,n)/2
aimg  = randn(n,n)/2
a2real = randn(n,n)/2
a2img  = randn(n,n)/2
breal = randn(n,2)/2
bimg  = randn(n,2)/2

# helper functions to unambiguously recover explicit forms of an LQPackedQ
squareQ(Q::LinearAlgebra.LQPackedQ) = (n = size(Q.factors, 2); lmul!(Q, Matrix{eltype(Q)}(I, n, n)))
rectangularQ(Q::LinearAlgebra.LQPackedQ) = convert(Array, Q)

@testset for eltya in (Float32, Float64, ComplexF32, ComplexF64)
    a = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(areal, aimg) : areal)
    a2 = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(a2real, a2img) : a2real)
    asym = a' + a                 # symmetric indefinite
    apd  = a' * a                 # symmetric positive-definite
    ε = εa = eps(abs(float(one(eltya))))

    @testset for eltyb in (Float32, Float64, ComplexF32, ComplexF64, Int)
        b = eltyb == Int ? rand(1:5, n, 2) : convert(Matrix{eltyb}, eltyb <: Complex ? complex.(breal, bimg) : breal)
        εb = eps(abs(float(one(eltyb))))
        ε = max(εa,εb)

        α = rand(eltya)
        aα = fill(α,1,1)
        @test lq(α).L*lq(α).Q ≈ lq(aα).L*lq(aα).Q
        @test abs(lq(α).Q[1,1]) ≈ one(eltya)
        tab = promote_type(eltya,eltyb)

        for i = 1:2
            let a = i == 1 ? a : view(a, 1:n - 1, 1:n - 1), b = i == 1 ? b : view(b, 1:n - 1), n = i == 1 ? n : n - 1
                lqa   = lq(a)
                l,q   = lqa.L, lqa.Q
                qra   = qr(a)
                @testset "Basic ops" begin
                    @test size(lqa,1) == size(a,1)
                    @test size(lqa,3) == 1
                    @test size(lqa.Q,3) == 1
                    @test_throws ErrorException lqa.Z
                    @test Array(copy(adjoint(lqa))) ≈ a'
                    @test q*squareQ(q)' ≈ Matrix(I, n, n)
                    @test l*q ≈ a
                    @test Array(lqa) ≈ a
                    @test Array(copy(lqa)) ≈ a
                    lstring = sprint(show,l)
                    qstring = sprint(show,q)
                    @test sprint(show,lqa) == "$(typeof(lqa)) with factors L and Q:\n$lstring\n$qstring"
                end
                @testset "Binary ops" begin
                    @test a*(lqa\b) ≈ b atol=3000ε
                    @test lqa*b ≈ qra.Q*qra.R*b atol=3000ε
                    @test (sq = size(q.factors, 2); *(Matrix{eltyb}(I, sq, sq), adjoint(q))*squareQ(q)) ≈ Matrix(I, n, n) atol=5000ε
                    if eltya != Int
                        @test Matrix{eltyb}(I, n, n)*q ≈ convert(AbstractMatrix{tab},q)
                    end
                    @test q*b ≈ squareQ(q)*b atol=100ε
                    @test transpose(q)*b ≈ transpose(squareQ(q))*b atol=100ε
                    @test q'*b ≈ squareQ(q)'*b atol=100ε
                    @test a*q ≈ a*squareQ(q) atol=100ε
                    @test a*transpose(q) ≈ a*transpose(squareQ(q)) atol=100ε
                    @test a*q' ≈ a*squareQ(q)' atol=100ε
                    @test a'*q ≈ a'*squareQ(q) atol=100ε
                    @test a'*q' ≈ a'*squareQ(q)' atol=100ε
                    @test_throws DimensionMismatch q*b[1:n1 + 1]
                    @test_throws DimensionMismatch adjoint(q) * Matrix{eltya}(undef,n+2,n+2)
                    @test_throws DimensionMismatch Matrix{eltyb}(undef,n+2,n+2)*q
                    if isa(a, DenseArray) && isa(b, DenseArray)
                        # use this to test 2nd branch in mult code
                        pad_a = vcat(I, a)
                        pad_b = hcat(I, b)
                        @test pad_a*q ≈ pad_a*squareQ(q) atol=100ε
                        @test transpose(q)*pad_b ≈ transpose(squareQ(q))*pad_b atol=100ε
                        @test q'*pad_b ≈ squareQ(q)'*pad_b atol=100ε
                    end
                end
            end
        end

        @testset "Matmul with LQ factorizations" begin
            lqa = lq(a[:,1:n1])
            l,q = lqa.L, lqa.Q
            @test rectangularQ(q)*rectangularQ(q)' ≈ Matrix(I, n1, n1)
            @test squareQ(q)'*squareQ(q) ≈ Matrix(I, n1, n1)
            @test_throws DimensionMismatch rmul!(Matrix{eltya}(I, n+1, n+1),q)
            @test lmul!(adjoint(q), rectangularQ(q)) ≈ Matrix(I, n1, n1)
            @test_throws DimensionMismatch rmul!(Matrix{eltya}(I, n+1, n+1), adjoint(q))
            @test_throws BoundsError size(q,-1)
        end
    end
end

@testset "getindex on LQPackedQ (#23733)" begin
    local m, n
    function getqs(F::LinearAlgebra.LQ)
        implicitQ = F.Q
        sq = size(implicitQ.factors, 2)
        explicitQ = lmul!(implicitQ, Matrix{eltype(implicitQ)}(I, sq, sq))
        return implicitQ, explicitQ
    end

    m, n = 3, 3 # reduced Q 3-by-3, full Q 3-by-3
    implicitQ, explicitQ = getqs(lq(randn(m, n)))
    @test implicitQ[1, 1] == explicitQ[1, 1]
    @test implicitQ[m, 1] == explicitQ[m, 1]
    @test implicitQ[1, n] == explicitQ[1, n]
    @test implicitQ[m, n] == explicitQ[m, n]

    m, n = 3, 4 # reduced Q 3-by-4, full Q 4-by-4
    implicitQ, explicitQ = getqs(lq(randn(m, n)))
    @test implicitQ[1, 1] == explicitQ[1, 1]
    @test implicitQ[m, 1] == explicitQ[m, 1]
    @test implicitQ[1, n] == explicitQ[1, n]
    @test implicitQ[m, n] == explicitQ[m, n]
    @test implicitQ[m+1, 1] == explicitQ[m+1, 1]
    @test implicitQ[m+1, n] == explicitQ[m+1, n]

    m, n = 4, 3 # reduced Q 3-by-3, full Q 3-by-3
    implicitQ, explicitQ = getqs(lq(randn(m, n)))
    @test implicitQ[1, 1] == explicitQ[1, 1]
    @test implicitQ[n, 1] == explicitQ[n, 1]
    @test implicitQ[1, n] == explicitQ[1, n]
    @test implicitQ[n, n] == explicitQ[n, n]
end

@testset "size on LQPackedQ (#23780)" begin
    # size(Q::LQPackedQ) yields the shape of Q's full/square form
    for ((mA, nA), nQ) in (
        ((3, 3), 3), # A 3-by-3 => full/square Q 3-by-3
        ((3, 4), 4), # A 3-by-4 => full/square Q 4-by-4
        ((4, 3), 3) )# A 4-by-3 => full/square Q 3-by-3
        @test size(lq(randn(mA, nA)).Q) == (nQ, nQ)
    end
end

@testset "postmultiplication with / right-application of LQPackedQ (#23779)" begin
    function getqs(F::LinearAlgebra.LQ)
        implicitQ = F.Q
        explicitQ = lmul!(implicitQ, Matrix{eltype(implicitQ)}(I, size(implicitQ)...))
        return implicitQ, explicitQ
    end
    # for any shape m-by-n of LQ-factored matrix, where Q is an LQPackedQ
    # A_mul_B*(C, Q) (Ac_mul_B*(C, Q)) operations should work for
    # *-by-n (n-by-*) C, which we test below via n-by-n C
    for (mA, nA) in ((3, 3), (3, 4), (4, 3))
        implicitQ, explicitQ = getqs(lq(randn(mA, nA)))
        C = randn(nA, nA)
        @test *(C, implicitQ) ≈ *(C, explicitQ)
        @test *(C, adjoint(implicitQ)) ≈ *(C, adjoint(explicitQ))
        @test *(adjoint(C), implicitQ) ≈ *(adjoint(C), explicitQ)
        @test *(adjoint(C), adjoint(implicitQ)) ≈ *(adjoint(C), adjoint(explicitQ))
    end
    # where the LQ-factored matrix has at least as many rows m as columns n,
    # Q's full/square and reduced/rectangular forms have the same shape (n-by-n). hence we expect
    # _only_ *-by-n (n-by-*) C to work in A_mul_B*(C, Q) (Ac_mul_B*(C, Q)) ops.
    # and hence the n-by-n C tests above suffice.
    #
    # where the LQ-factored matrix has more columns n than rows m,
    # Q's full/square form is n-by-n whereas its reduced/rectangular form is m-by-n.
    # hence we need also test *-by-m C with
    # A*_mul_B(C, Q) ops, as below via m-by-m C.
    mA, nA = 3, 4
    implicitQ, explicitQ = getqs(lq(randn(mA, nA)))
    C = randn(mA, mA)
    zeroextCright = hcat(C, zeros(eltype(C), mA))
    zeroextCdown = vcat(C, zeros(eltype(C), (1, mA)))
    @test *(C, implicitQ) ≈ *(zeroextCright, explicitQ)
    @test *(adjoint(C), implicitQ) ≈ *(adjoint(zeroextCdown), explicitQ)
    @test_throws DimensionMismatch C * adjoint(implicitQ)
    @test_throws DimensionMismatch adjoint(C) * adjoint(implicitQ)
end

end # module TestLQ

max_stars_repo_name: ninjin/julia

max_issues_repo_name: ninjin/julia

max_forks_repo_name: ninjin/julia


text:
"""
Wrapper for single-objective function with domain type VT
"""
mutable struct SingleObjective
    f # JuMP-expression TODO: use JuMPTypes or something?
    sense::MOI.OptimizationSense
    initialvalue::Dict{String,Any} # Variable => Initial value
    # TODO: implement bound in algorithm
    bound::Float64 # Hard lower/upper depending on sense
end

SingleObjective(; sense=MOI.MIN_SENSE) = SingleObjective(Any, sense, Dict{String,Any}(), NaN)
function SingleObjective(f::Union{AbstractJuMPScalar,MOI.AbstractScalarFunction,NonlinearExpression}; sense::MOI.OptimizationSense=MOI.MIN_SENSE,
                         iv::Dict{String,Any}=Dict{String,Any}(),
                         bound::Float64=NaN)
    SingleObjective(f, sense, iv, bound)
end
value(obj::SingleObjective) = value(obj.f)

"""
Orients the objective in the minimization sense
"""
function sensevalue(obj::SingleObjective)
    # To get Φ and Fstar with right signs
    if obj.sense == MOI.MAX_SENSE
        return -value(obj.f)
    else
        return value(obj.f)
    end
end

"""
Stores extension data for multi-objective
optimization inside JuMP Model
"""
mutable struct MultiData{Tx,To}
    objectives::Vector{SingleObjective}

    linear::Bool
    # Number of points in each direction
    # of Pareto submanifold
    pointsperdim::Int
    #
    # stored values
    utopiavarvalues::Vector{Vector{Tx}}
    paretovarvalues::Vector{Vector{To}}
    paretofront

    Phi::Array{To,2}
    Fstar::Vector{To}
    inequality::Bool
end

MultiData(; pointsperdim=10, linear=false, inequality=false) =
    MultiData(Array{SingleObjective}(undef, 0),
              linear, pointsperdim,
              Vector{Float64}[], Vector{Float64}[], Any[],
              Array{Float64}(undef, 2, 2),
              Array{Float64}(undef, 2), inequality)

function multi_model(optimizer; linear=false)
    m = Model(optimizer)
    return multi_model(m; linear=linear)
end

function multi_model(m::JuMP.AbstractModel; linear=false)
    m.optimize_hook = solvehook # defined in methods.jl
    m.ext[:Multi] = MultiData(;linear=linear)
    return m
end

"""
Gets the MultiData struct from a model if it exists
"""
function get_multidata(m::Model)
    if haskey(m.ext, :Multi)
        return m.ext[:Multi]::MultiData
    else
        error("This functionality is only available for MultiJuMP models")
    end
end

"""
Sets up the tree of all possible
convex combinations of the objectives
"""
function betas(levels, parts)
    map(Combinatorics.combinations(1:(levels + parts - 1), (levels - 1))) do x
        ([x;levels + parts] .- [0;x] .- 1) / parts
    end
end


text:
const DATA_SESSION = typedb.protocol.Session_Type.DATA
const SCHEMA_SESSION = typedb.protocol.Session_Type.SCHEMA

function dbconnect(f::Base.Callable, host::AbstractString, port::Int = 1729)
    client = CoreClient(host)
    try
        f(client)
    finally
        # close(client)
    end
end

function Base.open(
    f::Function,
    client::AbstractCoreClient,
    schema::AbstractString,
    type = DATA_SESSION
)
    session = Session(client, schema, type)
    try
        f(session)
    finally
        close(session)
    end
end

function readwrite(
    f::Function,
    session::AbstractCoreSession,
    direction::Integer,
    options::TypeDBOptions
)
    transaction = Transaction(session, session.sessionID, direction, options)
    try
        f(transaction)
    finally
        close(transaction)
    end
end

function Base.read(
    f::Function,
    session::AbstractCoreSession,
    options::TypeDBOptions = TypeDBOptions()
)
    readwrite(f, session, typedb.protocol.Transaction_Type.READ, options)
end

function Base.write(
    f::Function,
    session::AbstractCoreSession,
    options::TypeDBOptions = TypeDBOptions()
)
    readwrite(f, session, typedb.protocol.Transaction_Type.WRITE, options)
end


text:
"""
$(TYPEDEF)
Terms are part of the multi-sorted algebra that is part of High-Level Petri Net
markings, inscriptions and conditions.

See also [`AbstractDeclaration`](@ref).
"""
abstract type AbstractTerm end
abstract type AbstractOperator <: AbstractTerm end
abstract type AbstractSort end

"""
Declarations are the core of high-level Petri Net.
They define objects/names that are used for conditions, inscriptions, markings.
They are attached to PNML nets and pages.

$(TYPEDEF)
"""
abstract type AbstractDeclaration <: HLAnnotation end

pid(decl::AbstractDeclaration) = decl.id
name(decl::AbstractDeclaration) = decl.name

abstract type SortDeclaration <: AbstractDeclaration end
abstract type OperatorDeclaration <: AbstractDeclaration end

"""
$(TYPEDEF)
$(TYPEDFIELDS)
"""
struct UnknownDeclaration  <: AbstractDeclaration
    id::Symbol
    name::String
    nodename::String
    content::Vector{AnyElement}
    #sort::S
    #com::ObjectCommon
end

"""
$(TYPEDEF)
$(TYPEDFIELDS)
"""
struct VariableDeclaration{S}  <: AbstractDeclaration
    id::Symbol
    name::String
    sort::S
    #com::ObjectCommon
    #xml::XMLNode
end

VariableDeclaration(pdict::PnmlDict, xml::XMLNode) =
    VariableDeclaration(PnmlLabel(pdict, xml), ObjectCommon(pdict), xml)

"""
$(TYPEDEF)
$(TYPEDFIELDS)
"""
struct NamedSort{S} <: SortDeclaration
    id::Symbol
    name::String
    def::S # BuiltInSort, MultisetSort, ProductSort, UserSort
end

"""
$(TYPEDEF)
$(TYPEDFIELDS)

User-declared sort.
"""
struct Partition{S,PE} <: SortDeclaration
    id::Symbol
    name::String
    def::S # Refers to a NamedSort
    element::PE # 0 or more PartitionElements.
end

"""
$(TYPEDEF)
$(TYPEDFIELDS)
"""
struct NamedOperator{V,T} <: OperatorDeclaration
    id::Symbol
    name::String
    parameter::Vector{V}
    def::T # sort of term
end

"""
$(TYPEDEF)
$(TYPEDFIELDS)

Example input: <variable refvariable="varx"/>
"""
struct Variable <: AbstractTerm
    variableDecl::Symbol
end

struct BuiltInOperator <: AbstractOperator end
struct BuiltInConst <: AbstractOperator end
struct MultiSetOperator <: AbstractOperator end
struct PnmlTuple <: AbstractOperator end

"""
$(TYPEDEF)
$(TYPEDFIELDS)
"""
struct UserOperator <: AbstractOperator
    declaration::Symbol # operator declaration id
end

"""
Label of a net or page that holds zero or more [`AbstractDeclaration`](@ref).

$(TYPEDEF)
$(TYPEDFIELDS)
"""
struct Declaration <: HLAnnotation
    declarations::Vector{AbstractDeclaration} #TODO type stability?
    com::ObjectCommon
    #TODO attach XML node?
end

Declaration() = Declaration(Vector{AbstractDeclaration}[], ObjectCommon())

declarations(d::Declaration) = d.declarations
Base.length(d::Declaration) = length(declarations(d))

# Flattening pages combines declarations into the first page.
function Base.append!(l::Declaration, r::Declaration)
    append!(declarations(l), declarations(r))
    append!(l.com, r.com) # Only merges collections.
end

function Base.empty!(d::Declaration)
    empty!(declarations(d))
    empty!(d.com) #! JET will bark?
end


text:
mutable struct IVertex{T} <: Vertex{T}
  value::T
  inputs::Vector{IVertex{T}}

  IVertex{T}(x) where T = new(x, [])
end

IVertex(x) = IVertex{typeof(x)}(x)

value(v::IVertex) = v.value
inputs(v::IVertex) = v.inputs
outputs(v::IVertex) = ()
nout(v::IVertex) = 0

function thread!(to::IVertex, from::IVertex)
  push!(inputs(to), from)
  return to
end

function prethread!(to::IVertex, from::IVertex)
  pushfirst!(inputs(to), from)
  return to
end

il(v::Vertex) = convert(IVertex, v)

vertex(a...) = IVertex{Any}(a...)

vertex(x::Vertex) = convert(IVertex{Any}, x)

struct Stop
  x
end

stop(x) = Stop(x)

unwrap_stop(x) = false, x
unwrap_stop(x::Stop) = true, x.x

function walk!(v::IVertex, pre, post, cache = ODict())
  haskey(cache, v) && return cache[v]::typeof(v)
  stop, v′ = unwrap_stop(pre(v))
  cache[v] = v′
  stop || map!(v -> walk!(v, pre, post, cache), v′.inputs, v′.inputs)
  cache[v] = post(v′)
end

prewalk!(f, v::IVertex) = walk!(v, f, identity)
postwalk!(f, v::IVertex) = walk!(v, identity, f)

copy1(v::IVertex) = typeof(v)(v.value, v.inputs...)

precopy(v) = copy1(v)
precopy(v::Stop) = Stop(v.x)

walk(v::IVertex, pre, post) = walk!(v, v -> precopy(pre(v)), post)

prewalk(f, v::IVertex) = walk(v, f, identity)
postwalk(f, v::IVertex) = walk(v, identity, f)

copy(v::IVertex) = walk(v, identity, identity)

Base.map(f, v::IVertex) = prewalk(v -> typeof(v)(f(value(v)), inputs(v)...), v)

Base.replace(v::IVertex, pat, r) = prewalk(v -> v == pat ? r : v, v)

prefor(f, v) = prewalk!(v -> (f(v); v), v)

# TODO: check we don't get equivalent hashes for different graphs

function hash(v::IVertex, h::UInt = UInt(0), seen = OSet())
  h = hash(value(v), h)
  v in seen ? (return h) : push!(seen, v)
  for n in inputs(v)
    h ⊻= hash(n, h, seen)
  end
  return h
end

function iscyclic(v::IVertex)
  is = false
  prefor(v -> is |= (v < v), v)
  return is
end

function contains_(haystack::IVertex, needle)
  result = false
  map(v -> result |= v == needle, haystack)
  return result
end


text:
module TableWidgets

using InteractBase, Widgets, CSSUtil, Observables, WebIO
using Tables
using IterTools
import Observables: AbstractObservable, @map, @map!, @on
import Widgets: AbstractWidget, components, widgettype, layout!

import InteractBulma
using DataStructures
import DataStructures: reset!

using MacroTools

export categoricalselector, rangeselector, selector, selectors
export dataeditor, addfilter

const examplefolder = joinpath(@__DIR__, "..", "examples")

include("utils.jl")
include("selector.jl")
include("selectors.jl")
include("table.jl")
include("edit.jl")
include("filter.jl")

end # module


text:
module NFFTTools

using LinearAlgebra
using AbstractFFTs, FFTW  # FFTW just because of FFTW.ESTIMATE
using AbstractNFFTs

export sdc
export calculateToeplitzKernel, calculateToeplitzKernel!, convolveToeplitzKernel!

include("samplingDensity.jl")
include("Toeplitz.jl")

end


text:
"""
    Mosaic(T::Type, S::Type)

Data structure for generating a mosaic plot, a comparison between two categorical variables.

# Example

    using OnlineStats, Plots
    x = [rand() > .8 for i in 1:10^5]
    y = rand([1,2,2,3,3,3], 10^5)
    o = fit!(Mosaic(Bool, Int), zip(x, y))
    plot(o)
"""
struct Mosaic{T <: CountMap} <: OnlineStat{TwoThings}
    cm::T
end
Mosaic(T::Type, S::Type) = Mosaic(CountMap{Tuple{T,S}}())
value(o::Mosaic) = value(o.cm)
nobs(o::Mosaic) = nobs(o.cm)
_fit!(o::Mosaic, xy) = fit!(o.cm, Tuple(xy))
_merge!(a::Mosaic, b::Mosaic) = merge!(a.cm, b.cm)

function add_zero_counts!(o::Mosaic)
    v = value(o)
    akeys = [k[1] for k in keys(v)]
    bkeys = [k[2] for k in keys(v)]
    for ai in akeys, bi in bkeys
        get!(v, (ai, bi), 0)
    end
    sort!(v)
    o
end

function split_countmaps(o::Mosaic)
    v = value(o)
    akeys, bkeys = unique(first.(keys(v))), unique(last.(keys(v)))
    a = OrderedDict(k => sum(last, filter(x -> x[1][1] == k, v)) for k in akeys)
    b = OrderedDict(k => sum(last, filter(x -> x[1][2] == k, v)) for k in bkeys)
    sort!(a), sort!(b)
end

@recipe function f(o::Mosaic)
    n = nobs(o)
    a, b = split_countmaps(o)
    d = value(add_zero_counts!(o))

    x = vcat(0.0, cumsum([av / n for av in values(a)]))

    seriestype := :bar
    bar_width := diff(x)
    ylims --> (0,1)
    xlims --> (0, 1)
    label --> permutedims(string.(collect(keys(b))))
    xticks --> (midpoints(x), string.(collect(keys(a))))
    linewidth --> 0.5

    y = reverse(cumsum([d[(av,bv)] / a[av] for av in keys(a), bv in keys(b)], dims=2), dims=2)
    midpoints(x), y
end


text:
module ResidueClass

import Base: print, println,  *, /, -, +, ==, divrem, div, rem, copy, one, zero, sign, convert

using Polynomial

export Residue

### TO DO

##
# residue class ring with value a mod p^i
immutable Residue{T}

	val::T
	## generator of ideal
	p::T
	# power of ideal
	i::Int64
	# Ideal m = p^i, saved so that its computational more efficient
	m::T

	function Residue(val::T, p::T, i::Int, m::T)
		sign(val) >= 0 ? new(val,p,i,m) : new(val+m,p,i,m)
	end
end

## Factory ##
Residue{T}(a::T, p::T) = Residue{T}(rem(a,p), p, 1, p)
Residue{T}(a::T, p::T, i::Int) = (m = p^i; Residue{T}(a % m, p, i, m))
Residue{T}(a::T, p::T, i::Int, m::T) = Residue{T}(a % m, p, i, m)
Residue{T}(newVal::T, a::Residue{T}) = Residue{T}(newVal % a.m, a.p, a.i, a.m)
ResidueInRange{T}(newVal::T, a::Residue{T}) = Residue{T}(newVal, a.p, a.i, a.m)
# erzeugt aus einer Liste von Polynomen die entsprechenden Polynome über dem Restklassenring
function Residue{T}(a::Vector{Poly{T}}, p::T, i::Int)
	m = p^i
	array = Array{Poly{Residue{T}}}(size(a))
	for k=1:length(a)
		array[k] = Poly(Residue{T}(a[k].a, p, i, m), a[k].var)
	end
	return array
end

function Residue{T}(a::Vector{T}, p::T, i::Int, m::T)
	array = Array{Residue{T}}(size(a))
	for k=1:length(a)
		array[k] = Residue{T}(a[k] % m, p, i, m)
	end
	return array
end

Residue{T}(a::Vector{T}, p::T, i::Int) = Residue(a, p, i, p^i)
Residue{T}(a::Vector{T}, p::T) = Residue(a, p, 1, p)

convert{T}(t::Type{Residue{T}}, a::Residue{T}) = a

copy(a::Residue) = Residue(a.val, a)
one{T}(a::Residue{T}) = Residue(one(T), a.p, a.i, a.m)
zero{T}(a::Residue{T}) = Residue(zero(T), a)

# REVIEW
one{T}(t::Type{Residue{T}}) = one(T)
zero{T}(t::Type{Residue{T}}) = zero(T)

# ## Printing ##

print(io::IO, a::Residue) = (a.i == 1) ? print(io, "$(a.val) mod $(a.p)") :  print(io, "$(a.val) mod $(a.p)^$(a.i)=$(a.m)")

# ## Equates ##
==(a::Residue, b::Residue) = ((a.p != b.p) || (a.i != b.i)) ? false : a.val == b.val
==(a::WeakRef,b::Residue{WeakRef}) = (Residue(a % b.m, b) == b)
==(a::Residue{WeakRef},b::WeakRef) = ==(b,a)

## Arithmetics ##

sign{T}(a::Residue{T}) = sign(a.val)

# plus #
+{T}(a::Residue{T}, b::Residue{T}) = (a.m != b.m) ? error("Can't add residues!") : Residue(a.val+b.val, a)
+{T}(a::Residue{T}, b::T) = Residue(a.val+b, a)
+{T}(a::T, b::Residue{T}) = +(b,a)

# minus #
-{T}(a::Residue{T}, b::Residue{T}) = (a.m != b.m) ? error("Can't add residues!") : Residue(a.val - b.val, a)
-{T}(a::Residue{T}, b::T) = Residue(rem(a.val - b,a.m), a)
-{T}(a::T, b::Residue{T}) = -(b,a)

# unary minus #
-{T}(a::Residue{T}) = Residue(-a.val, a)

# # multiply #
*{T}(a::T, b::Residue{T}) = Residue(rem(a * b.val, b.m), b)
*{T}(a::Residue{T}, b::T) = *(b,a)
*{T}(a::Residue{T}, b::Residue{T}) = (a.m != b.m) ? error("Can't multiply! Quotients ($(a.m) and $(b.m))") : Residue(a.val*b.val, b)
function *{T}(a::Vector{Residue{T}}, b::Residue{T})
	c = Array{Residue{T}}(size(a))
	for i=1:length(a)
		c[i] = a[i] * b
	end
	return c
end
 *{T}(a::Residue{T}, b::Vector{Residue{T}}) = *(b,a)

function ^{T}(a::Residue{T}, n::Int)

    if n == 0
    	return ResidueInRange(one(T), a)
    end
    # a^-n = (a^-1)^n
    a_new = n > 0 ? a : inverse(a)
    b = copy(a_new)
    for k in reverse(digits(abs(n),2))[2:end]
        b = (k == 1) ? (b * b) * a_new : b * b
    end
    return b
end

# a is invertibale iff gcd(a, a.m) == 1, bezout coefficient s is then the inverse
function inverse{T}(a::Residue{T})
	(gcd, s, t) = gcdx(a.val, a.m)
    if gcd != one(T)
    	error("Value has no Inverse")
    end
    return Residue(s, a)
end

# division #

function divrem{T}(a::Residue{T}, b::Residue{T})
	if (a.m != b.m)
		error("Quotients ($(a.m) and $(b.m))")
	end
	(div, rem) = divrem(a.val, b.val)
	return (Residue(div, a), Residue(rem, a))
end

div{T}(a::Residue{T}, b::Residue{T}) = divrem(a,b)[1]
rem{T}(a::Residue{T}, b::Residue{T}) = divrem(a,b)[2]
/{T}(a::Residue{T}, b::Residue{T}) = a * inverse(b)

# function gcdx{T}(a::Residue{T}, b::Residue{T})
# 	(gcd, s, t) = gcdx(a.val, b.val)
# 	return (Residue(gcd, a), Residue(s, a), Residue(t, a))
# end

# gcd{T}(a::Residue{T}, b::Residue{T}) = Residue(gcd(a.val, b.val), a)

end

max_stars_repo_name: saschatimme/Polynomial

max_issues_repo_name: saschatimme/Polynomial

max_forks_repo_name: saschatimme/Polynomial


text:
abstract type MajFlowNorm end

struct OneNorm <: MajFlowNorm end

struct InfNorm <: MajFlowNorm end

export SortedProbVecMult

struct SortedProbVecMult{T, VT <: AbstractVector{T}} <: AbstractVector{T}
    distinct_entries::VT
    multiplicities::Vector{Int}
end

Base.size(r::SortedProbVecMult) = (sum(r.multiplicities),)
Base.eltype(r::SortedProbVecMult{T}) where {T} = T

Base.@propagate_inbounds function Base.getindex(r::SortedProbVecMult, i::Int)
    @boundscheck (1 <= i <= length(r)) || throw(BoundsError(r,i))
    ks = r.multiplicities
    seen_so_far = 0
    for j = eachindex(ks)
        seen_so_far += ks[j]
        if seen_so_far >= i
            return r.distinct_entries[j]
        end
    end
end

function SortedProbVecMult(v::AbstractVector)
    v_sorted = sort(v; rev = true)
    distinct_entries = eltype(v)[]
    multiplicities = Int[]
    for elt in v_sorted
        if isempty(distinct_entries) || distinct_entries[end] != elt
            push!(distinct_entries, elt)
            push!(multiplicities, 1)
        else
            multiplicities[end] += 1
        end

    end
    SortedProbVecMult{eltype(distinct_entries), typeof(distinct_entries)}(distinct_entries, multiplicities)
end

function Base.collect(spvm::SortedProbVecMult)
    @unpack distinct_entries, multiplicities = spvm
    v = eltype(distinct_entries)[]
    sizehint!(v, length(spvm))
    for (i, elt) in enumerate(distinct_entries)
        for j = 1:multiplicities[i]
            push!(v, elt)
        end
    end
    return v
end

# Default
majmin(v::SortedProbVecMult, ϵ) = majmin(OneNorm(), v, ϵ)

function majmin(nrm::MajFlowNorm, spvm::SortedProbVecMult, ϵ)
    majmin!(nrm, deepcopy(spvm), ϵ)
end

function majmin!(nrm::OneNorm, spvm::SortedProbVecMult, ϵ)
    @unpack distinct_entries, multiplicities = spvm
    @assert ϵ >= 0
    if length(distinct_entries) == 1
        return spvm
    elseif length(distinct_entries) == 2
        _majmin_onenorm_2entries!(spvm::SortedProbVecMult, ϵ)
    else
        kp = first(multiplicities)
        km = last(multiplicities)

        δp = kp * (distinct_entries[1] - distinct_entries[2])
        δm = km * (distinct_entries[end-1] - distinct_entries[end])
        gap = min(δp, δm)
        if ϵ < gap
            distinct_entries[1] -= ϵ / kp
            distinct_entries[end] += ϵ / km
            return spvm
        else
            if δp < δm
                # evolve for time δp and recurse
                deleteat!(distinct_entries, 1)
                deleteat!(multiplicities, 1)
                multiplicities[1] += kp
                distinct_entries[end] += δp / km
                return majmin!(nrm, spvm, ϵ - δp)
            elseif δm < δp
                # evolve for time δm and recurse
                pop!(distinct_entries)
                pop!(multiplicities)
                multiplicities[end] += km
                distinct_entries[1] -= δm / kp
                return majmin!(nrm, spvm, ϵ - δm)
            else
                # evolve for time gap=δm=δp and recurse
                deleteat!(distinct_entries, 1)
                deleteat!(multiplicities, 1)
                multiplicities[1] += kp
                pop!(distinct_entries)
                pop!(multiplicities)
                multiplicities[end] += km
                return majmin!(nrm, spvm, ϵ - gap)
            end
        end
    end
end

# With only 2 entries, the next crossing isn't just determined by
# gap = min(δp, δm).
function _majmin_onenorm_2entries!(spvm::SortedProbVecMult, ϵ)
    @unpack distinct_entries, multiplicities = spvm
    @assert length(distinct_entries) == length(multiplicities) == 2
    kp, km = multiplicities
    μp, μm = distinct_entries
    # the larger entry is decreasing at a rate t/kp
    # the smaller entry is increasing at a rate t/km
    # Need to solve   `μp - t/kp == μm + t/km` for `t`
    # t = (μp - μm) / ( 1/km + 1/kp )
    t = (μp - μm) * (km * kp) / (km + kp)
    if ϵ < t
        distinct_entries[1] -= ϵ / kp
        distinct_entries[2] += ϵ / km
        return spvm
    else
        pop!(distinct_entries)
        pop!(multiplicities)
        distinct_entries[] -= t / kp
        multiplicities[] += km
        return spvm
    end
end


text:
int_or_bigint(x) = (y=mod(x,Int); x==y ? convert(Int,x) : convert(BigInt,y))

mutable struct UnpickleState
    stack :: Stack
    proto :: Int
    pos :: Int
    memo :: Memo
end
UnpickleState() = UnpickleState(Stack(), HIGHEST_PROTOCOL, 1, Memo())

function read_opcode(io::IO, state::UnpickleState)
    ans = read(io, OpCode)
    state.pos += sizeof(OpCode)
    ans
end

function read_prim_l(io::IO, state::UnpickleState, ::Type{T}) where {T}
    ans = ltoh(read(io, T))
    state.pos += sizeof(T)
    ans
end

function read_prim_b(io::IO, state::UnpickleState, ::Type{T}) where {T}
    ans = ntoh(read(io, T))
    state.pos += sizeof(T)
    ans
end

read_u1(io::IO, state::UnpickleState) = read_prim_l(io, state, UInt8)
read_u2(io::IO, state::UnpickleState) = read_prim_l(io, state, UInt16)
read_u4(io::IO, state::UnpickleState) = read_prim_l(io, state, UInt32)
read_u8(io::IO, state::UnpickleState) = read_prim_l(io, state, UInt64)
read_s4(io::IO, state::UnpickleState) = read_prim_l(io, state, Int32)
read_f8(io::IO, state::UnpickleState) = read_prim_b(io, state, Float64)

function read_bytes(io::IO, state::UnpickleState, sz::Integer)
    ans = read(io, sz)
    state.pos += length(ans)
    length(ans) == sz || error("unexpected end of file")
    ans
end

function read_line(io::IO, state::UnpickleState)
    ans = readuntil(io, '\n', keep=false)
    state.pos += sizeof(ans) + 1
    ans
end

function read_stringnl_noescape(io::IO, state::UnpickleState)
    read_line(io, state)
end

function read_decimalnl_short(io::IO, state::UnpickleState)
    line = read_line(io, state)
    if line == "00"
        return false
    elseif line == "01"
        return true
    end
    return int_or_bigint(parse(BigInt, line))
end

function read_floatnl(io::IO, state::UnpickleState)
    line = read_line(io, state)
    return parse(Float64, line)
end

function unpickle(io::IO)
    state = UnpickleState()
    # optional first PROTO opcode
    op = read_opcode(io, state)
    if op == OP_PROTO
        state.proto = Int(read_u1(io, state))
        stop = false
    else
        stop = doop(op, state, io)
    end
    # remaining opcodes
    while !stop
        op = read_opcode(io, state)
        stop = doop(op, state, io)
    end
    # done
    len = length(state.stack)
    if len == 1
        return pop!(state.stack)
    else
        error("unexpected STOP")
    end
end

function doop(op::OpCode, state::UnpickleState, io::IO)
    if op == OP_MARK
        mark!(state.stack)

    elseif op == OP_STOP
        return true

    elseif op == OP_POP
        pop!(state.stack)

    elseif op == OP_POP_MARK
        poptomark!(state.stack)

    elseif op == OP_DUP
        push!(state.stack, top(state.stack))

    elseif op == OP_FLOAT
        val = read_floatnl(io, state)
        push!(state.stack, val)

    elseif op == OP_INT
        val = read_decimalnl_short(io, state)
        push!(state.stack, val)

    elseif op == OP_BININT
        val = int_or_bigint(read_s4(io, state))
        push!(state.stack, val)

    elseif op == OP_BININT1
        val = int_or_bigint(read_u1(io, state))
        push!(state.stack, val)

    elseif op == OP_LONG
        error("opcode not implemented: $op")

    elseif op == OP_BININT2
        val = int_or_bigint(read_u2(io, state))
        push!(state.stack, val)

    elseif op == OP_NONE
        val = nothing
        push!(state.stack, val)

    elseif op == OP_PERSID
        error("opcode not implemented: $op")

    elseif op == OP_BINPERSID
        error("opcode not implemented: $op")

    elseif op == OP_REDUCE
        args = pop!(state.stack)::PyTuple
        func = pop!(state.stack)
        val = PyFuncCall(func, args.values)
        push!(state.stack, val)

    elseif op == OP_STRING
        error("opcode not implemented: $op")

    elseif op == OP_BINSTRING
        error("opcode not implemented: $op")

    elseif op == OP_SHORT_BINSTRING
        error("opcode not implemented: $op")

    elseif op == OP_UNICODE
        error("opcode not implemented: $op")

    elseif op == OP_BINUNICODE
        sz = read_u4(io, state)
        val = String(read_bytes(io, state, sz))
        push!(state.stack, val)

    elseif op == OP_APPEND
        x = pop!(state.stack)
        list = top(state.stack)::PyList
        push!(list.values, x)

    elseif op == OP_BUILD
        arg = pop!(state.stack)
        obj = top(state.stack)
        pysetstate!(obj, arg)

    elseif op == OP_GLOBAL
        mod = read_stringnl_noescape(io, state)
        attr = read_stringnl_noescape(io, state)
        val = PyGlobal(mod, attr)
        push!(state.stack, val)

    elseif op == OP_DICT
        kvs = poptomark!(state.stack)
        iseven(length(kvs)) || error("odd number of keys and values")
        val = PyDict(Pair{Any,Any}[Pair{Any,Any}(kvs[i], kvs[i+1]) for i in 1:2:length(kvs)])
        push!(state.stack, val)

    elseif op == OP_EMPTY_DICT
        val = PyDict()
        push!(state.stack, val)

    elseif op == OP_APPENDS
        xs = poptomark!(state.stack)
        list = top(state.stack)::PyList
        append!(list.values, xs)

    elseif op == OP_GET
        error("opcode not implemented: $op")

    elseif op == OP_BINGET
        idx = read_u1(io, state)
        val = state.memo[idx]
        push!(state.stack, val)

    elseif op == OP_INST
        error("opcode not implemented: $op")

    elseif op == OP_LONG_BINGET
        idx = read_u4(io, state)
        val = state.memo[idx]
        push!(state.stack, val)

    elseif op == OP_LIST
        val = PyList(poptomark!(state.stack))
        push!(state.stack, val)

    elseif op == OP_EMPTY_LIST
        val = PyList()
        push!(state.stack, val)

    elseif op == OP_OBJ
        error("opcode not implemented: $op")

    elseif op == OP_PUT
        idx = read_decimalnl_short(io, state)
        state.memo[idx] = top(state.stack)

    elseif op == OP_BINPUT
        idx = read_u1(io, state)
        state.memo[idx] = top(state.stack)

    elseif op == OP_LONG_BINPUT
        idx = read_u4(io, state)
        state.memo[idx] = top(state.stack)

    elseif op == OP_SETITEM
        v = pop!(state.stack)
        k = pop!(state.stack)
        dict = top(state.stack)::PyDict
        push!(dict.items, Pair{Any,Any}(k, v))

    elseif op == OP_TUPLE
        val = PyTuple(poptomark!(state.stack))
        push!(state.stack, val)

    elseif op == OP_EMPTY_TUPLE
        val = PyTuple()
        push!(state.stack, val)

    elseif op == OP_SETITEMS
        kvs = poptomark!(state.stack)
        iseven(length(kvs)) || error("odd number of keys and values")
        dict = top(state.stack)::PyDict
        for i in 1:2:length(kvs)
            push!(dict.items, Pair{Any,Any}(kvs[i], kvs[i+1]))
        end

    elseif op == OP_BINFLOAT
        val = read_f8(io, state)
        push!(state.stack, val)

    elseif op == OP_PROTO
        error("unexpected op: $op")

    elseif op == OP_NEWOBJ
        args = pop!(state.stack)::PyTuple
        cls = pop!(state.stack)
        val = PyNewObj(cls, args.values)
        push!(state.stack, val)

    elseif op == OP_EXT1
        error("opcode not implemented: $op")

    elseif op == OP_EXT2
        error("opcode not implemented: $op")

    elseif op == OP_EXT4
        error("opcode not implemented: $op")

    elseif op == OP_TUPLE1
        x1 = pop!(state.stack)
        val = PyTuple(Any[x1])
        push!(state.stack, val)

    elseif op == OP_TUPLE2
        x2 = pop!(state.stack)
        x1 = pop!(state.stack)
        val = PyTuple(Any[x1, x2])
        push!(state.stack, val)

    elseif op == OP_TUPLE3
        x3 = pop!(state.stack)
        x2 = pop!(state.stack)
        x1 = pop!(state.stack)
        val = PyTuple(Any[x1, x2, x3])
        push!(state.stack, val)

    elseif op == OP_NEWTRUE
        val = true
        push!(state.stack, val)

    elseif op == OP_NEWFALSE
        val = false
        push!(state.stack, val)

    elseif op == OP_LONG1
        error("opcode not implemented: $op")

    elseif op == OP_LONG4
        error("opcode not implemented: $op")

    elseif op == OP_BINBYTES
        sz = read_u4(io, state)
        val = PyBytes(read_bytes(io, state, sz))
        push!(state.stack, val)

    elseif op == OP_SHORT_BINBYTES
        sz = read_u1(io, state)
        val = PyBytes(read_bytes(io, state, sz))
        push!(state.stack, val)

    elseif op == OP_SHORT_BINUNICODE
        sz = read_u1(io, state)
        val = String(read_bytes(io, state, sz))
        push!(state.stack, val)

    elseif op == OP_BINUNICODE8
        sz = read_u8(io, state)
        val = String(read_bytes(io, state, sz))
        push!(state.stack, val)

    elseif op == OP_BINBYTES8
        sz = read_u8(io, state)
        val = PyBytes(read_bytes(io, state, sz))
        push!(state.stack, val)

    elseif op == OP_EMPTY_SET
        val = PySet()
        push!(state.stack, val)

    elseif op == OP_ADDITEMS
        xs = poptomark!(state.stack)
        set = top(state.stack)::PySet
        append!(set.values, xs)

    elseif op == OP_FROZENSET
        val = PyFrozenSet(poptomark!(state.stack))
        push!(state.stack, val)

    elseif op == OP_NEWOBJ_EX
        kwargs = pop!(state.stack)::PyDict
        args = pop!(state.stack)::PyTuple
        cls = pop!(state.stack)
        val = PyNewObj(cls, args.values, kwargs.items)
        push!(state.stack, val)

    elseif op == OP_STACK_GLOBAL
        attr = pop!(state.stack)::String
        mod = pop!(state.stack)::String
        val = PyGlobal(mod, attr)
        push!(state.stack, val)

    elseif op == OP_MEMOIZE
        push!(state.memo, top(state.stack))

    elseif op == OP_FRAME
        # TODO: we could reuse `frameio` to avoid some allocations
        sz = read_u8(io, state)
        pos = state.pos
        frameio = IOBuffer(read_bytes(io, state, sz))
        while !eof(frameio)
            op = read_opcode(frameio, state)
            stop = doop(op, state, frameio)
            stop && return true
        end
        state.pos = pos + sz

    elseif op == OP_BYTEARRAY8
        sz = read_u8(io, state)
        val = PyByteArray(read_bytes(io, state, sz))
        push!(state.stack, val)

    elseif op == OP_NEXT_BUFFER
        error("opcode not implemented: $op")

    elseif op == OP_READONLY_BUFFER
        error("opcode not implemented: $op")

    else
        error("opcode not implemented: $op")
    end

    return false
end