text:
module QJuliaCUDABlas

push!(LOAD_PATH, @__DIR__)

import QJuliaRegisters

using CUDAnative
using CuArrays

#create function/type alias
double2  = QJuliaRegisters.double2
double4  = QJuliaRegisters.double4
#
float2   = QJuliaRegisters.float2
float4   = QJuliaRegisters.float4

# Real operations
# Native complex values
 @inline function axpy_(a::T, x::Complex{T}, y::Complex{T})::Complex{T} where T <: AbstractFloat

   rres = CUDAnative.fma(a, real(x), real(y))
   ires = CUDAnative.fma(a, imag(x), imag(y))

   return Complex{T}(rres, ires)
 end

 @inline function caxpy_(a::Complex{T}, x::Complex{T}, y::Complex{T})::Complex{T} where T <: AbstractFloat

    rres = CUDAnative.fma(+real(a), real(x), real(y))
    rres = CUDAnative.fma(-imag(a), imag(x), rres   )

    ires = CUDAnative.fma(+imag(a), real(x), imag(y))
    ires = CUDAnative.fma(+real(a), imag(x), ires   )

    return Complex{T}(rres, ires)
  end

# Float2 regs

 @inline function axpy_(a::T, x::NTuple{2, VecElement{T}}, y::NTuple{2, VecElement{T}})::NTuple{2, VecElement{T}} where T <: AbstractFloat

   res1 = CUDAnative.fma(a, x[1].value, y[1].value)
   res2 = CUDAnative.fma(a, x[2].value, y[2].value)

   return NTuple{2, VecElement{T}}((res1, res2))
 end

 @inline function axpy_(a::NTuple{2, VecElement{T}}, x::NTuple{2, VecElement{T}}, y::NTuple{2, VecElement{T}})::NTuple{2, VecElement{T}} where T <: AbstractFloat
   return (VecElement(CUDAnative.fma(a[1].value, x[1].value, y[1].value)), VecElement(CUDAnative.fma(a[2].value, x[2].value, y[2].value)))
 end

 @inline function caxpy_(a::NTuple{2, VecElement{T}}, x::NTuple{2, VecElement{T}}, y::NTuple{2, VecElement{T}})::NTuple{2, VecElement{T}} where T <: AbstractFloat
   res  = (VecElement(CUDAnative.fma(+a[1].value, x[1].value, y[1].value)),
           VecElement(CUDAnative.fma(+a[1].value, x[2].value, y[2].value)))
   return (VecElement(CUDAnative.fma(-a[2].value, x[2].value, res[1].value)),
           VecElement(CUDAnative.fma(+a[2].value, x[1].value, res[2].value)))
 end

# Float4 regs

 @inline function axpy_(a, x::NTuple{4, VecElement{T}}, y::NTuple{4, VecElement{T}})::NTuple{4, VecElement{T}} where T <: AbstractFloat

   res1 = CUDAnative.fma(a, x[1].value, y[1].value)
   res2 = CUDAnative.fma(a, x[2].value, y[2].value)
   res3 = CUDAnative.fma(a, x[3].value, y[3].value)
   res4 = CUDAnative.fma(a, x[4].value, y[4].value)

   return NTuple{4, VecElement{T}}((res1, res2, res3, res4))
 end

 @inline function caxpy_(a::NTuple{4, VecElement{T}}, x::NTuple{4, VecElement{T}}, y::NTuple{4, VecElement{T}})::NTuple{4, VecElement{T}} where T <: AbstractFloat
   res  = (VecElement(CUDAnative.fma(+a[1].value, x[1].value, y[1].value)),
           VecElement(CUDAnative.fma(+a[1].value, x[2].value, y[2].value)),
           VecElement(CUDAnative.fma(+a[3].value, x[3].value, y[3].value)),
           VecElement(CUDAnative.fma(+a[3].value, x[4].value, y[4].value)))
   return (VecElement(CUDAnative.fma(-a[2].value, x[2].value, res[1].value)),
           VecElement(CUDAnative.fma(+a[2].value, x[1].value, res[2].value)),
           VecElement(CUDAnative.fma(-a[4].value, x[4].value, res[3].value)),
           VecElement(CUDAnative.fma(+a[4].value, x[3].value, res[4].value)))
 end

 # a few helper methods:
 @inline function get_regN(meta::NTuple{N, VecElement{T}}) where T <: AbstractFloat where N; return N; end
 @inline function get_regT(meta::NTuple{N, VecElement{T}}) where T <: AbstractFloat where N; return T; end

end #QJuliaCUDABlas

max_stars_repo_name: alexstrel/QJulia

max_issues_repo_name: alexstrel/QJulia

max_forks_repo_name: alexstrel/QJulia


text:
# Mag
Mag(x) = set!(Mag(), x)
Mag(x::Union{MagRef,ArfRef}) = Mag(mag_struct(cstruct(x)))

# Arf
Arf(x; prec::Integer = _precision(x)) = set!(Arf(prec = prec), x)
# disambiguation
Arf(x::Arf; prec::Integer = precision(x)) = set!(Arf(prec = prec), x)

#Arb
Arb(x; prec::Integer = _precision(x)) = set!(Arb(prec = prec), x)
# disambiguation
Arb(x::Arb; prec::Integer = precision(x)) = set!(Arb(prec = prec), x)

function Arb(str::AbstractString; prec::Integer = DEFAULT_PRECISION[])
    res = Arb(prec = prec)
    flag = set!(res, str)
    iszero(flag) || throw(ArgumentError("could not parse $str as an Arb"))
    return res
end

## Acb
Acb(x::Union{Real,arb_struct,arf_struct}; prec::Integer = _precision(x)) =
    set!(Acb(prec = prec), x)
Acb(
    re::Union{Real,arb_struct,arf_struct},
    im::Union{Real,arb_struct,arf_struct};
    prec::Integer = max(_precision(re), _precision(im)),
) = set!(Acb(prec = prec), re, im)

Acb(z::Complex; prec::Integer = max(_precision(real(z)), _precision(imag(z)))) =
    set!(Acb(prec = prec), z)
Acb(x::AcbLike; prec::Integer = precision(x)) = set!(Acb(prec = prec), x)
# disambiguation
Acb(x::Acb; prec::Integer = precision(x)) = set!(Acb(prec = prec), x)

function Acb(str::AbstractString; prec::Integer = DEFAULT_PRECISION[])
    res = Acb(prec = prec)
    flag = set!(realref(res), str)
    iszero(flag) || throw(ArgumentError("could not parse $str as an Arb"))
    return res
end

function Acb(re::AbstractString, im::AbstractString; prec::Integer = DEFAULT_PRECISION[])
    res = Acb(prec = prec)
    flag = set!(realref(res), re)
    iszero(flag) || throw(ArgumentError("could not parse $str as an Arb"))
    flag = set!(imagref(res), im)
    iszero(flag) || throw(ArgumentError("could not parse $str as an Arb"))
    return res
end

Base.Int(x::ArfLike; rnd::Union{arb_rnd,RoundingMode} = RoundNearest) =
    is_int(x) ? get_si(x, rnd) : throw(InexactError(:Int64, Int64, x))

Base.Float64(x::MagLike) = get(x)
Base.Float64(x::ArfLike; rnd::Union{arb_rnd,RoundingMode} = RoundNearest) = get_d(x, rnd)
Base.Float64(x::ArbLike) = Float64(midref(x))

Base.ComplexF64(z::AcbLike) = Complex(Float64(realref(z)), Float64(imagref(z)))

function Base.BigFloat(x::Union{Arf,ArfRef})
    y = BigFloat(; precision = precision(x))
    get!(y, x)
    y
end
Base.BigFloat(x::Union{Arb,ArbRef}) = BigFloat(midref(x))

Base.zero(::Union{Mag,Type{Mag}}) = Mag(UInt64(0))
Base.one(::Union{Mag,Type{Mag}}) = Mag(UInt64(1))
Base.zero(x::T) where {T<:Union{Arf,Arb,Acb}} = T(0, prec = precision(x))
Base.one(x::T) where {T<:Union{Arf,Arb,Acb}} = T(1, prec = precision(x))

# Define these since the base implementation would create `n` copies of the same element
# I.e. only allocating **one** Arf/Arb/Acb.
Base.zeros(x::T, n::Integer) where {T<:Union{Mag,Arf,Arb,Acb}} = [zero(x) for _ = 1:n]
Base.ones(x::T, n::Integer) where {T<:Union{Mag,Arf,Arb,Acb}} = [one(x) for _ = 1:n]
Base.zeros(::Type{T}, n::Integer) where {T<:Union{Mag,Arf,Arb,Acb}} = [zero(T) for _ = 1:n]
Base.ones(::Type{T}, n::Integer) where {T<:Union{Mag,Arf,Arb,Acb}} = [one(T) for _ = 1:n]

ext: jl


text:
"""
    choose_edge(g::Network)

Randomly selects an inactive edge in `g`

Arguments
* `g`   : An instance of type Network
Returns
* `edge`: A two-tuple of integers representing an inactive edge in `g`
"""
function choose_edge(g::Network)
    node_ids = UnitRange{eltype(g.cluster_ids)}(1, g.N)
    edge = (rand(g.rng, node_ids), rand(g.rng, node_ids))

    if edge[1] ≠ edge[2] && edge ∉ g.edges && reverse(edge) ∉ g.edges
        return edge
    else
        choose_edge(g)
    end
end

"""
    choose_edge(g::Lattice2D)

Randomly selects an inactive edge in `g`

Arguments
* `g`   : An instance of type Lattice2D
Returns
* `edge`: A two-tuples of integers representing an inactive edge in `g`
"""
function choose_edge(g::Lattice2D)
    node_ids = UnitRange{eltype(g.cluster_ids)}(1, g.L)
    node = (rand(g.rng, node_ids), rand(g.rng, node_ids))
    neighbor = nearest_neighbors(g, node)[rand(g.rng, 1:4)]
    edge = (cart_to_lin(node, g.L), cart_to_lin(neighbor, g.L))

    if edge ∉ g.edges && reverse(edge) ∉ g.edges
        return edge
    else
        choose_edge(g)
    end
end

"""
    choose_edge(g::Lattice3D)

Randomly selects an inactive edge in `g`

Arguments
* `g`   : An instance of type Lattice3D
Returns
* `edge`: A two-tuple of integers representing an inactive edge in `g`
"""
function choose_edge(g::Lattice3D)
    node = (rand(g.rng, 1:g.L), rand(g.rng, 1:g.L), rand(g.rng, 1:g.L))
    neighbor = nearest_neighbors(g, node)[rand(g.rng, 1:6)]
    edge = (cart_to_lin(node, g.L), cart_to_lin(neighbor, g.L))

    if edge ∉ g.edges && reverse(edge) ∉ g.edges
        return edge
    else
        choose_edge(g)
    end
end

"""
    choose_edge(g::AbstractGraph, edge₁::NTuple{2, Integer})

Randomly selects an inactive edge in `g` that is not equal to `edge₁`

Arguments
* `g`    : An instance of type AbstractGraph
* `edge₁`: A two-tuple of integers representing an inactive edge in `g`
Returns
* `edge₂`: A two-tuple of integers representing an inactive edge in `g` not equal to `edge₁`
"""
function choose_edge(g::AbstractGraph, edge₁::NTuple{2,Integer})
    edge₂ = choose_edge(g)

    if edge₂ ≠ edge₁ && edge₂ ≠ reverse(edge₁)
        return edge₂
    else
        choose_edge(g, edge₁)
    end
end

"""
    choose_edge(g::AbstractGraph, edges::Array{NTuple{2, Integer}, 1})

Randomly selects an inactive edge in `g` that is not equal to `edge₁`

Arguments
* `g`    : An instance of type AbstractGraph
* `edges`: An array of two-tuples of integers representing inactive edges in `g`
Returns
* `edge`: A two-tuple of integers representing an inactive edge in `g` not in `edges`
"""
function choose_edge(g::AbstractGraph, edges::Array{NTuple{2,Integer},1})
    edge = choose_edge(g)

    if edge ∉ edges && reverse(edge) ∉ edges
        return edge
    else
        choose_edge(g, edges)
    end
end

"""
    add_edge!(g::AbstractGraph, edge::NTuple{2, Integer})

Adds an edge to `g`

Arguments
* `g`   : An instance of type AbstractGraph
* `edge`: A two-tuple of integers representing the edge to be added to `g`
Returns
* None, updates `g` in-place
"""
function add_edge!(g::AbstractGraph, edge::NTuple{2,Integer})
    push!(g.edges, edge)
    g.t += 1
    update_clusters!(g, edge)
end

"""
Functions below this point are for determining nearest-neighbors with periodic boundary conditions in hyper-cubic lattices and converting the indices from Cartesian to linear
"""

"""
`plus` and `minus` are an implementation of periodic boundary conditions for use in the `nearest_neighbors` functions below
"""
function plus(L::Integer, i::Integer)
    i == L ? 1 : i + 1
end
function minus(L::Integer, i::Integer)
    i == 1 ? L : i - 1
end

"""
    nearest_neighbors(g::Lattice2D, node::NTuple{2, Integer})

Determines the nearest neighbors of `node`

Arguments
* `g`        : An instance of type Lattice2D
* `node`     : A two-tuple of integers representing a node in `g`
Returns
* `neighbors`: A four-tuple of two-tuples of integers representing the cartesian indices of the (up, down, left, right) neighbors
"""
function nearest_neighbors(g::Lattice2D, node::NTuple{2,Integer})
    return (
        (minus(g.L, node[1]), node[2]),
        (plus(g.L, node[1]), node[2]),
        (node[1], minus(g.L, node[2])),
        (node[1], plus(g.L, node[2])),
    )
end

"""
    nearest_neighbors(g::Lattice3D, node::NTuple{3, Integer})

Determines the nearest neighbors of `node`

Arguments
* `g`        : An instance of type Lattice3D
* `node`     : A three-tuple of integers representing a node in `g`
Returns
* `neighbors`: A four-tuple of three-tuples of integers representing the cartesian indices of the (up, down, left, right, front, back) neighbors
"""
function nearest_neighbors(g::Lattice3D, node::NTuple{3,Integer})
    return (
        (minus(g.L, node[1]), node[2], node[3]),
        (plus(g.L, node[1]), node[2], node[3]),
        (node[1], minus(g.L, node[2]), node[3]),
        (node[1], plus(g.L, node[2]), node[3]),
        (node[1], node[2], minus(g.L, node[3])),
        (node[1], node[2], plus(g.L, node[3])),
    )
end

"""
    cart_to_lin(cart::Tuple, L::Integer)

Converts d-dimensional Cartesian index (d-tuple) to linear index

Arguments
* `cart` : d-Tuple representing Cartesian index, d ∈ {2, 3}
* `L`    : Linear lattice size
Returns
* `lin`  : Linear index in the lattice corresponding to cart
"""
function cart_to_lin(cart::Tuple, L::Integer)
    if length(cart) == 2
        return (cart[2] - 1) * L + cart[1]
    elseif length(cart) == 3
        return (cart[3] - 1) * L^2 + (cart[2] - 1) * L + cart[1]
    end
end


text:
@testset "optimize" begin
    eta = 0.9

    function f1(x)
        (1.0 / 2.0) * (x[1]^2 + eta * x[2]^2)
    end

    function g1(storage, x)
        storage[1] = x[1]
        storage[2] = eta * x[2]
    end

    function h1(storage, x)
        storage[1, 1] = 1.0
        storage[1, 2] = 0.0
        storage[2, 1] = 0.0
        storage[2, 2] = eta
    end

    results = optimize(f1, g1, h1, [127.0, 921.0])
    @test Optim.g_converged(results)
    @test norm(Optim.minimizer(results) - [0.0, 0.0]) < 0.01

    results = optimize(f1, g1, [127.0, 921.0])
    @test Optim.g_converged(results)
    @test norm(Optim.minimizer(results) - [0.0, 0.0]) < 0.01

    results = optimize(f1, [127.0, 921.0])
    @test Optim.g_converged(results)
    @test norm(Optim.minimizer(results) - [0.0, 0.0]) < 0.01

    results = optimize(f1, [127.0, 921.0])
    @test Optim.g_converged(results)
    @test norm(Optim.minimizer(results) - [0.0, 0.0]) < 0.01

    # tests for bfgs_initial_invH
    initial_invH = zeros(2,2)
    h1(initial_invH, [127.0, 921.0])
    initial_invH = Matrix(Diagonal(diag(initial_invH)))
    results = optimize(f1, g1, [127.0, 921.0], BFGS(initial_invH = x -> initial_invH), Optim.Options())
    @test Optim.g_converged(results)
    @test norm(Optim.minimizer(results) - [0.0, 0.0]) < 0.01
end


text:
using Distributed
using MCPhylo
using Random, Test
include("utils.jl")

const samplertests = [
  #"abc",
  "linear_model_uvp",
  "linear_model_mvp"
]

const mcmctests = [
   "discretediag",
   "modelbuilding"
#   "readcoda"
]

# const extensiontests = [
#   "newunivardist",
#   "newmultivardist"
# ]

println("Running tests:")
#
@testset "All tests" begin

for t in samplertests
    @runtest "samplers/" t
  end

@testset "mcmctests" begin
   for t in mcmctests
     @runtest "mcmc/" t
   end
end

end
#all sets test set

max_stars_repo_name: erathorn/MCPhylo

max_issues_repo_name: erathorn/MCPhylo

max_forks_repo_name: erathorn/MCPhylo


text:
using Documenter, ModiaBase

makedocs(
  #modules  = [ModiaBase],
  sitename = "ModiaBase",
  authors  = "Hilding Elmqvist (Mogram) and Martin Otter (DLR-SR)",
  format = Documenter.HTML(prettyurls = false),
  pages    = [
     "Home"   => "index.md",
	 "Tutorial" => "Tutorial.md",
     "Data Structures"        => "DataStructures.md",
     "Equation Sorting"       => "EquationSorting.md",
     "Equation Reduction"     => "EquationReduction.md",
     "Transformation to ODE System" => "TransformationToODEs.md",
  ]
)


text:
abstract type AbstractProblemParSetter end
struct DummyProblemParSetter <: AbstractProblemParSetter; end # for testing error message

"""
    count_state(::AbstractProblemParSetter)
    count_par(::AbstractProblemParSetter)
    count_paropt(::AbstractProblemParSetter)

Report the number of problem states, problem parameters and optimized parameters
respectively.
"""
# need to implement in concrete types, only there for documentation
function count_state(pset::AbstractProblemParSetter)
    error("Expect $(typeof(pset)) to implement method `count_state` but was not.")
end,
function count_par(pset::AbstractProblemParSetter)
    error("Expect $(typeof(pset)) to implement method `count_par` but was not.")
end,
function count_paropt(pset::AbstractProblemParSetter)
    error("Expect $(typeof(pset)) to implement method `count_paropt` but was not.")
end

@deprecate count_states(pset::AbstractProblemParSetter)  count_state(pset)

"""
    symbols_state(pset::AbstractProblemParSetter)
    symbols_par(pset::AbstractProblemParSetter)
    symbols_paropt(pset::AbstractProblemParSetter)

Report the names, i.e. symbols of problem states, problem parameters and
optimized parameters respectively, i.e. the concatenation of components.
Similar to `ComponentArrays.label`, but inferred from Axis object.
"""
function symbols_state(pset::AbstractProblemParSetter)
    error("Expect $(typeof(pset)) to implement method `symbols_state` but was not.")
end,
function symbols_par(pset::AbstractProblemParSetter)
    error("Expect $(typeof(pset)) to implement method `symbols_par` but was not.")
end,
function symbols_paropt(pset::AbstractProblemParSetter)
    error("Expect $(typeof(pset)) to implement method `symbols_paropt` but was not.")
end
# need to implement in concrete types

@deprecate statesyms(pset::AbstractProblemParSetter) symbols_state(pset)
@deprecate parsyms(pset::AbstractProblemParSetter) symbols_par(pset)
@deprecate paroptsyms(pset::AbstractProblemParSetter) symbols_paropt(pset)

"""
    axis_state(pset::AbstractProblemParSetter)
    axis_par(pset::AbstractProblemParSetter)
    axis_paropt(pset::AbstractProblemParSetter)

Report the Axis, i.e. nested component symbols of problem states, problem parameters and
optimized parameters respectively.
Returns an AbstractAxis.
"""
function axis_state(pset::AbstractProblemParSetter)
    error("Expect $(typeof(pset)) to implement method `axis_state` but was not.")
end,
function axis_par(pset::AbstractProblemParSetter)
    error("Expect $(typeof(pset)) to implement method `axis_par` but was not.")
end,
function axis_paropt(pset::AbstractProblemParSetter)
    error("Expect $(typeof(pset)) to implement method `axis_paropt` but was not.")
end
# need to implement in concrete types

"""
    get_paropt(pset::AbstractProblemParSetter, prob::ODEProblem; kwargs...)
    get_paropt(pset::AbstractProblemParSetter, u0, p)

    get_paropt_labeled(pset::AbstractProblemParSetter, prob::ODEProblem; kwargs...)
    get_paropt_labeled(pset::AbstractProblemParSetter, u0, p)

Extract the initial states and parameters corresponding to the positions
that are optimized.
If both u0 and p are AbstractVectors, the result is a Vector, otherwise the result is a Tuple.

The labeled versions additionally call `label_paropt` (see [`label_state`](@ref))
on the return value.
"""
function get_paropt(pset::AbstractProblemParSetter, prob::ODEProblem; kwargs...)
    get_paropt(pset, prob.u0, prob.p; kwargs...)
end,
function get_paropt_labeled(pset::AbstractProblemParSetter, prob::ODEProblem; kwargs...)
    get_paropt_labeled(pset, prob.u0, prob.p; kwargs...)
end,
function get_paropt(pset::AbstractProblemParSetter, u0, p)
    getdata(get_paropt_labeled(pset, u0,p))

end
# need to implement in concrete types: get_paropt_labeled -> ComponentVector

"""
    label_state(pset::AbstractProblemParSetter, u::AbstractVector)
    label_par(pset::AbstractProblemParSetter, par::AbstractVector)
    label_paropt(pset::AbstractProblemParSetter, popt::AbstractVector)

Produce a labeled version, i.e. a ComponentVector of initial states, parameters, or
optimized parameters respectively.
"""
function label_state(pset::AbstractProblemParSetter, u); attach_axis(u, axis_state(pset)); end,
function label_par(pset::AbstractProblemParSetter, p); attach_axis(p, axis_par(pset)); end,
function label_paropt(pset::AbstractProblemParSetter, popt); attach_axis(popt, axis_paropt(pset)); end

# TODO move to ComponentArrays.jl
attach_axis(x::AbstractVector, ax::AbstractAxis) = ComponentArray(x, (ax,))
attach_axis(x::ComponentVector, ax::AbstractAxis) = ComponentArray(getdata(x), (ax,))

# label_state(pset::AbstractProblemParSetter, u::SVector) = SLVector(label_state(pset, Tuple(u)))
# label_state(pset::AbstractProblemParSetter, u::NTuple) = NamedTuple{symbols_state(pset)}(u)

# label_par(pset::AbstractProblemParSetter, par::SVector) = SLVector(label_par(pset, Tuple(par)))
# label_par(pset::AbstractProblemParSetter, par::NTuple) = NamedTuple{symbols_par(pset)}(par)

# label_paropt(pset::AbstractProblemParSetter, popt::SVector) = SLVector(label_paropt(pset, Tuple(popt)))
# label_paropt(pset::AbstractProblemParSetter, popt::NTuple) = NamedTuple{symbols_paropt(pset)}(popt)

"""
    name_state(pset, u::AbstractVector)
    name_par(pset, par::AbstractVector)
    name_paropt(pset, popt::AbstractVector)

    name_paropt(pset::AbstractProblemParSetter, prob::ODEProblem)

Produce a `NamedVector` of given state, parameters, or optimized vars.
"""
function name_state(pset::AbstractProblemParSetter, state::AbstractVector); NamedArray(
    state, (collect(symbols_state(pset))::Vector{Symbol},)); end,
function name_par(pset::AbstractProblemParSetter, par::AbstractVector); NamedArray(
    par, (collect(symbols_par(pset))::Vector{Symbol},)); end,
function name_paropt(pset::AbstractProblemParSetter, paropt::AbstractVector); NamedArray(
    paropt, (collect(symbols_paropt(pset))::Vector{Symbol},)); end,
function name_paropt(pset::AbstractProblemParSetter, prob::ODEProblem; kwargs...); name_paropt(
    pset, get_paropt(pset, prob); kwargs...); end

"""
    prob_new = update_statepar(ps::ProblemParSetter, popt, prob::ODEProblem)
    u0new, pnew = update_statepar(ps::ProblemParSetter, popt, u0, p)

Return an updated problem or updates states and parameters where
values corresponding to positions in `popt` hold the values of popt.
The type is changed to the promotion type of popt, to allow working with dual numbers.
"""
function update_statepar(pset::AbstractProblemParSetter, popt, prob::ODEProblem)
    u0,p = update_statepar(pset, popt, prob.u0, prob.p)
    remake(prob; u0, p)
end


text:
## Integrator Dispatches

# Can get rid of an allocation here with a function
# get_tmp_arr(integrator.cache) which gives a pointer to some
# cache array which can be modified.

function ode_addsteps!(args...)
  @warn("`ode_addsteps!` is deprecated. Use `addsteps!`")
  DiffEqBase.addsteps!(args...)
end

@inline function _ode_addsteps!(integrator,f=integrator.f,always_calc_begin = false,allow_calc_end = true,force_calc_end = false)
  if !(typeof(integrator.cache) <: CompositeCache)
    DiffEqBase.addsteps!(integrator.k,integrator.tprev,integrator.uprev,integrator.u,
                  integrator.dt,f,integrator.p,integrator.cache,
                  always_calc_begin,allow_calc_end,force_calc_end)
  else
    DiffEqBase.addsteps!(integrator.k,integrator.tprev,integrator.uprev,integrator.u,
                  integrator.dt,f,integrator.p,
                  integrator.cache.caches[integrator.cache.current],
                  always_calc_begin,allow_calc_end,force_calc_end)
  end
end
@inline DiffEqBase.addsteps!(integrator::ODEIntegrator,args...) = _ode_addsteps!(integrator,args...)

@inline function ode_interpolant(Θ,integrator::DiffEqBase.DEIntegrator,idxs,deriv)
  DiffEqBase.addsteps!(integrator)
  if !(typeof(integrator.cache) <: CompositeCache)
    val = ode_interpolant(Θ,integrator.dt,integrator.uprev,integrator.u,integrator.k,integrator.cache,idxs,deriv)
  else
    val = ode_interpolant(Θ,integrator.dt,integrator.uprev,integrator.u,integrator.k,integrator.cache.caches[integrator.cache.current],idxs,deriv)
  end
  val
end

@inline function ode_interpolant!(val,Θ,integrator::DiffEqBase.DEIntegrator,idxs,deriv)
  DiffEqBase.addsteps!(integrator)
  if !(typeof(integrator.cache) <: CompositeCache)
    ode_interpolant!(val,Θ,integrator.dt,integrator.uprev,integrator.u,integrator.k,integrator.cache,idxs,deriv)
  else
    ode_interpolant!(val,Θ,integrator.dt,integrator.uprev,integrator.u,integrator.k,integrator.cache.caches[integrator.cache.current],idxs,deriv)
  end
end

@inline function current_interpolant(t::Number,integrator::DiffEqBase.DEIntegrator,idxs,deriv)
  Θ = (t-integrator.tprev)/integrator.dt
  ode_interpolant(Θ,integrator,idxs,deriv)
end

@inline function current_interpolant(t,integrator::DiffEqBase.DEIntegrator,idxs,deriv)
  Θ = (t.-integrator.tprev)./integrator.dt
  [ode_interpolant(ϕ,integrator,idxs,deriv) for ϕ in Θ]
end

@inline function current_interpolant!(val,t::Number,integrator::DiffEqBase.DEIntegrator,idxs,deriv)
  Θ = (t-integrator.tprev)/integrator.dt
  ode_interpolant!(val,Θ,integrator,idxs,deriv)
end

@inline function current_interpolant!(val,t,integrator::DiffEqBase.DEIntegrator,idxs,deriv)
  Θ = (t.-integrator.tprev)./integrator.dt
  [ode_interpolant!(val,ϕ,integrator,idxs,deriv) for ϕ in Θ]
end

@inline function current_extrapolant(t::Number,integrator::DiffEqBase.DEIntegrator,idxs=nothing,deriv=Val{0})
  Θ = (t-integrator.tprev)/(integrator.t-integrator.tprev)
  ode_extrapolant(Θ,integrator,idxs,deriv)
end

@inline function current_extrapolant!(val,t::Number,integrator::DiffEqBase.DEIntegrator,idxs=nothing,deriv=Val{0})
  Θ = (t-integrator.tprev)/(integrator.t-integrator.tprev)
  ode_extrapolant!(val,Θ,integrator,idxs,deriv)
end

@inline function current_extrapolant(t::AbstractArray,integrator::DiffEqBase.DEIntegrator,idxs=nothing,deriv=Val{0})
  Θ = (t.-integrator.tprev)./(integrator.t-integrator.tprev)
  [ode_extrapolant(ϕ,integrator,idxs,deriv) for ϕ in Θ]
end

@inline function current_extrapolant!(val,t,integrator::DiffEqBase.DEIntegrator,idxs=nothing,deriv=Val{0})
  Θ = (t.-integrator.tprev)./(integrator.t-integrator.tprev)
  [ode_extrapolant!(val,ϕ,integrator,idxs,deriv) for ϕ in Θ]
end

@inline function ode_extrapolant!(val,Θ,integrator::DiffEqBase.DEIntegrator,idxs,deriv)
  DiffEqBase.addsteps!(integrator)
  if !(typeof(integrator.cache) <: CompositeCache)
    ode_interpolant!(val,Θ,integrator.t-integrator.tprev,integrator.uprev2,integrator.uprev,integrator.k,integrator.cache,idxs,deriv)
  else
    ode_interpolant!(val,Θ,integrator.t-integrator.tprev,integrator.uprev2,integrator.uprev,integrator.k,integrator.cache.caches[integrator.cache.current],idxs,deriv)
  end
end

@inline function ode_extrapolant(Θ,integrator::DiffEqBase.DEIntegrator,idxs,deriv)
  DiffEqBase.addsteps!(integrator)
  if !(typeof(integrator.cache) <: CompositeCache)
    ode_interpolant(Θ,integrator.t-integrator.tprev,integrator.uprev2,integrator.uprev,integrator.k,integrator.cache,idxs,deriv)
  else
    ode_interpolant(Θ,integrator.t-integrator.tprev,integrator.uprev2,integrator.uprev,integrator.k,integrator.cache.caches[integrator.cache.current],idxs,deriv)
  end
end

##

"""
ode_interpolation(tvals,ts,timeseries,ks)

Get the value at tvals where the solution is known at the
times ts (sorted), with values timeseries and derivatives ks
"""
function ode_interpolation(tvals,id,idxs,deriv,p,continuity::Symbol=:left)
  @unpack ts,timeseries,ks,f,cache = id
  tdir = sign(ts[end]-ts[1])
  idx = sortperm(tvals,rev=tdir<0)
  i = 2 # Start the search thinking it's between ts[1] and ts[2]
  tdir*tvals[idx[end]] > tdir*ts[end] && error("Solution interpolation cannot extrapolate past the final timepoint. Either solve on a longer timespan or use the local extrapolation from the integrator interface.")
  tdir*tvals[idx[1]] < tdir*ts[1] && error("Solution interpolation cannot extrapolate before the first timepoint. Either start solving earlier or use the local extrapolation from the integrator interface.")
  if typeof(idxs) <: Number
    vals = Vector{eltype(first(timeseries))}(undef, length(tvals))
  elseif typeof(idxs) <: AbstractArray
    vals = Vector{Array{eltype(first(timeseries)),ndims(idxs)}}(undef, length(tvals))
  else
    vals = Vector{eltype(timeseries)}(undef, length(tvals))
  end
  @inbounds for j in idx
    t = tvals[j]
    i = searchsortedfirst(@view(ts[i:end]),t,rev=tdir<0)+i-1 # It's in the interval ts[i-1] to ts[i]
    avoid_constant_ends = deriv != Val{0} || typeof(t) <: ForwardDiff.Dual
    avoid_constant_ends && i==1 && (i+=1)
    if !avoid_constant_ends && ts[i] == t
      lasti = lastindex(ts)
      k = continuity == :right && i+1 <= lasti && ts[i+1] == t ? i+1 : i
      if idxs == nothing
        vals[j] = timeseries[k]
      else
        vals[j] = timeseries[k][idxs]
      end
    elseif !avoid_constant_ends && ts[i-1] == t # Can happen if it's the first value!
      if idxs == nothing
        vals[j] = timeseries[i-1]
      else
        vals[j] = timeseries[i-1][idxs]
      end
    else
      dt = ts[i] - ts[i-1]
      Θ = (t-ts[i-1])/dt
      if typeof(cache) <: (FunctionMapCache) || typeof(cache) <: FunctionMapConstantCache
        vals[j] = ode_interpolant(Θ,dt,timeseries[i-1],timeseries[i],0,cache,idxs,deriv)
      elseif !id.dense
        vals[j] = linear_interpolant(Θ,dt,timeseries[i-1],timeseries[i],idxs,deriv)
      elseif typeof(cache) <: CompositeCache
        DiffEqBase.addsteps!(ks[i],ts[i-1],timeseries[i-1],timeseries[i],dt,f,p,cache.caches[id.alg_choice[i-1]]) # update the kcurrent
        vals[j] = ode_interpolant(Θ,dt,timeseries[i-1],timeseries[i],ks[i],cache.caches[id.alg_choice[i-1]],idxs,deriv)
      else
        DiffEqBase.addsteps!(ks[i],ts[i-1],timeseries[i-1],timeseries[i],dt,f,p,cache) # update the kcurrent
        vals[j] = ode_interpolant(Θ,dt,timeseries[i-1],timeseries[i],ks[i],cache,idxs,deriv)
      end
    end
  end
  DiffEqArray(vals, tvals)
end

"""
ode_interpolation(tvals,ts,timeseries,ks)

Get the value at tvals where the solution is known at the
times ts (sorted), with values timeseries and derivatives ks
"""
function ode_interpolation!(vals,tvals,id,idxs,deriv,p,continuity::Symbol=:left)
  @unpack ts,timeseries,ks,f,cache = id
  tdir = sign(ts[end]-ts[1])
  idx = sortperm(tvals,rev=tdir<0)
  i = 2 # Start the search thinking it's between ts[1] and ts[2]
  tdir*tvals[idx[end]] > tdir*ts[end] && error("Solution interpolation cannot extrapolate past the final timepoint. Either solve on a longer timespan or use the local extrapolation from the integrator interface.")
  tdir*tvals[idx[1]] < tdir*ts[1] && error("Solution interpolation cannot extrapolate before the first timepoint. Either start solving earlier or use the local extrapolation from the integrator interface.")
  @inbounds for j in idx
    t = tvals[j]
    i = searchsortedfirst(@view(ts[i:end]),t,rev=tdir<0)+i-1 # It's in the interval ts[i-1] to ts[i]
    avoid_constant_ends = deriv != Val{0} || typeof(t) <: ForwardDiff.Dual
    avoid_constant_ends && i==1 && (i+=1)
    if !avoid_constant_ends && ts[i] == t
      lasti = lastindex(ts)
      k = continuity == :right && i+1 <= lasti && ts[i+1] == t ? i+1 : i
      if idxs == nothing
        vals[j] = timeseries[k]
      else
        vals[j] = timeseries[k][idxs]
      end
    elseif !avoid_constant_ends && ts[i-1] == t # Can happen if it's the first value!
      if idxs == nothing
        vals[j] = timeseries[i-1]
      else
        vals[j] = timeseries[i-1][idxs]
      end
    else
      dt = ts[i] - ts[i-1]
      Θ = (t-ts[i-1])/dt
      if typeof(cache) <: (FunctionMapCache) || typeof(cache) <: FunctionMapConstantCache
        if eltype(timeseries) <: AbstractArray
          ode_interpolant!(vals[j],Θ,dt,timeseries[i-1],timeseries[i],0,cache,idxs,deriv)
        else
          vals[j] = ode_interpolant(Θ,dt,timeseries[i-1],timeseries[i],0,cache,idxs,deriv)
        end
      elseif !id.dense
        if eltype(timeseries) <: AbstractArray
          linear_interpolant!(vals[j],Θ,dt,timeseries[i-1],timeseries[i],idxs,deriv)
        else
          vals[j] = linear_interpolant(Θ,dt,timeseries[i-1],timeseries[i],idxs,deriv)
        end
      elseif typeof(cache) <: CompositeCache
        DiffEqBase.addsteps!(ks[i],ts[i-1],timeseries[i-1],timeseries[i],dt,f,p,cache.caches[id.alg_choice[i-1]]) # update the kcurrent
        if eltype(timeseries) <: AbstractArray
          ode_interpolant!(vals[j],Θ,dt,timeseries[i-1],timeseries[i],ks[i],cache.caches[id.alg_choice[i-1]],idxs,deriv)
        else
          vals[j] = ode_interpolant(Θ,dt,timeseries[i-1],timeseries[i],ks[i],cache.caches[id.alg_choice[i-1]],idxs,deriv)
        end
      else
        DiffEqBase.addsteps!(ks[i],ts[i-1],timeseries[i-1],timeseries[i],dt,f,p,cache) # update the kcurrent
        if eltype(vals[j]) <: AbstractArray
          ode_interpolant!(vals[j],Θ,dt,timeseries[i-1],timeseries[i],ks[i],cache,idxs,deriv)
        else
          vals[j] = ode_interpolant(Θ,dt,timeseries[i-1],timeseries[i],ks[i],cache,idxs,deriv)
        end
      end
    end
  end
end

"""
ode_interpolation(tval::Number,ts,timeseries,ks)

Get the value at tval where the solution is known at the
times ts (sorted), with values timeseries and derivatives ks
"""
function ode_interpolation(tval::Number,id,idxs,deriv,p,continuity::Symbol=:left)
  @unpack ts,timeseries,ks,f,cache = id
  tdir = sign(ts[end]-ts[1])
  tdir*tval > tdir*ts[end] && error("Solution interpolation cannot extrapolate past the final timepoint. Either solve on a longer timespan or use the local extrapolation from the integrator interface.")
  tdir*tval < tdir*ts[1] && error("Solution interpolation cannot extrapolate before the first timepoint. Either start solving earlier or use the local extrapolation from the integrator interface.")
  @inbounds i = searchsortedfirst(ts,tval,rev=tdir<0) # It's in the interval ts[i-1] to ts[i]
  avoid_constant_ends = deriv != Val{0} || typeof(tval) <: ForwardDiff.Dual
  avoid_constant_ends && i==1 && (i+=1)
  @inbounds if !avoid_constant_ends && ts[i] == tval
    lasti = lastindex(ts)
    k = continuity == :right && i+1 <= lasti && ts[i+1] == tval ? i+1 : i
    if idxs == nothing
      val = timeseries[k]
    else
      val = timeseries[k][idxs]
    end
  elseif !avoid_constant_ends && ts[i-1] == tval # Can happen if it's the first value!
    if idxs == nothing
      val = timeseries[i-1]
    else
      val = timeseries[i-1][idxs]
    end
  else
    dt = ts[i] - ts[i-1]
    Θ = (tval-ts[i-1])/dt
    if typeof(cache) <: (FunctionMapCache) || typeof(cache) <: FunctionMapConstantCache
      val = ode_interpolant(Θ,dt,timeseries[i-1],timeseries[i],0,cache,idxs,deriv)
    elseif !id.dense
      val = linear_interpolant(Θ,dt,timeseries[i-1],timeseries[i],idxs,deriv)
    elseif typeof(cache) <: CompositeCache
      DiffEqBase.addsteps!(ks[i],ts[i-1],timeseries[i-1],timeseries[i],dt,f,p,cache.caches[id.alg_choice[i-1]]) # update the kcurrent
      val = ode_interpolant(Θ,dt,timeseries[i-1],timeseries[i],ks[i],cache.caches[id.alg_choice[i-1]],idxs,deriv)
    else
      DiffEqBase.addsteps!(ks[i],ts[i-1],timeseries[i-1],timeseries[i],dt,f,p,cache) # update the kcurrent
      val = ode_interpolant(Θ,dt,timeseries[i-1],timeseries[i],ks[i],cache,idxs,deriv)
    end
  end
  val
end

"""
ode_interpolation!(out,tval::Number,ts,timeseries,ks)

Get the value at tval where the solution is known at the
times ts (sorted), with values timeseries and derivatives ks
"""
function ode_interpolation!(out,tval::Number,id,idxs,deriv,p,continuity::Symbol=:left)
  @unpack ts,timeseries,ks,f,cache = id
  @inbounds tdir = sign(ts[end]-ts[1])
  @inbounds tdir*tval > tdir*ts[end] && error("Solution interpolation cannot extrapolate past the final timepoint. Either solve on a longer timespan or use the local extrapolation from the integrator interface.")
  @inbounds tdir*tval < tdir*ts[1] && error("Solution interpolation cannot extrapolate before the first timepoint. Either start solving earlier or use the local extrapolation from the integrator interface.")
  i = searchsortedfirst(ts,tval,rev=tdir<0) # It's in the interval ts[i-1] to ts[i]
  avoid_constant_ends = deriv != Val{0} || typeof(tval) <: ForwardDiff.Dual
  avoid_constant_ends && i==1 && (i+=1)
  if !avoid_constant_ends && ts[i] == tval
    lasti = lastindex(ts)
    k = continuity == :right && i+1 <= lasti && ts[i+1] == tval ? i+1 : i
    if idxs == nothing
      @inbounds copyto!(out,timeseries[k])
    else
      @inbounds copyto!(out,timeseries[k][idxs])
    end
  elseif !avoid_constant_ends && ts[i-1] == tval # Can happen if it's the first value!
    if idxs == nothing
      @inbounds copyto!(out,timeseries[i-1])
    else
      @inbounds copyto!(out,timeseries[i-1][idxs])
    end
  else
    @inbounds begin
      dt = ts[i] - ts[i-1]
      Θ = (tval-ts[i-1])/dt
      if typeof(cache) <: (FunctionMapCache) || typeof(cache) <: FunctionMapConstantCache
        ode_interpolant!(out,Θ,dt,timeseries[i-1],timeseries[i],0,cache,idxs,deriv)
      elseif !id.dense
        linear_interpolant!(out,Θ,dt,timeseries[i-1],timeseries[i],idxs,deriv)
      elseif typeof(cache) <: CompositeCache
        DiffEqBase.addsteps!(ks[i],ts[i-1],timeseries[i-1],timeseries[i],dt,f,p,cache.caches[id.alg_choice[i-1]]) # update the kcurrent
        ode_interpolant!(out,Θ,dt,timeseries[i-1],timeseries[i],ks[i],cache.caches[id.alg_choice[i-1]],idxs,deriv)
      else
        DiffEqBase.addsteps!(ks[i],ts[i-1],timeseries[i-1],timeseries[i],dt,f,p,cache) # update the kcurrent
        ode_interpolant!(out,Θ,dt,timeseries[i-1],timeseries[i],ks[i],cache,idxs,deriv)
      end
    end
  end
end

"""
By default, Hermite interpolant so update the derivative at the two ends
"""
function DiffEqBase.addsteps!(k,t,uprev,u,dt,f,p,cache,always_calc_begin = false,allow_calc_end = true,force_calc_end = false)
  if length(k)<2 || always_calc_begin
    if typeof(cache) <: OrdinaryDiffEqMutableCache
      rtmp = similar(u, eltype(eltype(k)))
      f(rtmp,uprev,p,t)
      copyat_or_push!(k,1,rtmp)
      f(rtmp,u,p,t+dt)
      copyat_or_push!(k,2,rtmp)
    else
      copyat_or_push!(k,1,f(uprev,p,t))
      copyat_or_push!(k,2,f(u,p,t+dt))
    end
  end
  nothing
end

"""
ode_interpolant and ode_interpolant! dispatch
"""
function ode_interpolant(Θ,dt,y₀,y₁,k,cache,idxs,T::Type{Val{TI}}) where TI
  _ode_interpolant(Θ,dt,y₀,y₁,k,cache,idxs,T)
end

function ode_interpolant(Θ,dt,y₀,y₁,k,cache::OrdinaryDiffEqMutableCache,idxs,T::Type{Val{TI}}) where TI
  if typeof(idxs) <: Number || typeof(y₀) <: Union{Number,SArray}
    # typeof(y₀) can be these if saveidxs gives a single value
    _ode_interpolant(Θ,dt,y₀,y₁,k,cache,idxs,T)
  else
    # determine output type
    # required for calculation of time derivatives with autodifferentiation
    oneunit_Θ = oneunit(Θ)
    if isempty(k)
      S = promote_type(typeof(oneunit_Θ * oneunit(eltype(y₀)))) # Θ*y₀
    else
      S = promote_type(typeof(oneunit_Θ * oneunit(eltype(y₀))), # Θ*y₀
                     typeof(oneunit_Θ * oneunit(dt) * oneunit(eltype(k[1])))) # Θ*dt*k
    end
    if typeof(idxs) <: Nothing
      out = similar(y₀, S)
    else
      out = similar(y₀, S, axes(idxs))
    end
    _ode_interpolant!(out,Θ,dt,y₀,y₁,k,cache,idxs,T)
  end
end

function ode_interpolant!(out,Θ,dt,y₀,y₁,k,cache,idxs,T::Type{Val{TI}}) where TI
  _ode_interpolant!(out,Θ,dt,y₀,y₁,k,cache,idxs,T)
end

##################### Hermite Interpolants

# If no dispatch found, assume Hermite
function _ode_interpolant(Θ,dt,y₀,y₁,k,cache,idxs,T::Type{Val{TI}}) where TI
  hermite_interpolant(Θ,dt,y₀,y₁,k,idxs,T)
end

function _ode_interpolant!(out,Θ,dt,y₀,y₁,k,cache,idxs,T::Type{Val{TI}}) where TI
  hermite_interpolant!(out,Θ,dt,y₀,y₁,k,idxs,T)
end

"""
Hairer Norsett Wanner Solving Ordinary Differential Euations I - Nonstiff Problems Page 190

Herimte Interpolation, chosen if no other dispatch for ode_interpolant
"""
@muladd function hermite_interpolant(Θ,dt,y₀,y₁,k,idxs::Nothing,T::Type{Val{0}}) # Default interpolant is Hermite
  #@. (1-Θ)*y₀+Θ*y₁+Θ*(Θ-1)*((1-2Θ)*(y₁-y₀)+(Θ-1)*dt*k[1] + Θ*dt*k[2])
  (1-Θ)*y₀+Θ*y₁+Θ*(Θ-1)*((1-2Θ)*(y₁-y₀)+(Θ-1)*dt*k[1] + Θ*dt*k[2])
end

@muladd function hermite_interpolant(Θ,dt,y₀,y₁,k,idxs,T::Type{Val{0}}) # Default interpolant is Hermite
  # return @. (1-Θ)*y₀[idxs]+Θ*y₁[idxs]+Θ*(Θ-1)*((1-2Θ)*(y₁[idxs]-y₀[idxs])+(Θ-1)*dt*k[1][idxs] + Θ*dt*k[2][idxs])
  return (1-Θ)*y₀[idxs]+Θ*y₁[idxs]+Θ*(Θ-1)*((1-2Θ)*(y₁[idxs]-y₀[idxs])+(Θ-1)*dt*k[1][idxs] + Θ*dt*k[2][idxs])
end

@muladd function hermite_interpolant!(out,Θ,dt,y₀,y₁,k,idxs::Nothing,T::Type{Val{0}}) # Default interpolant is Hermite
  #@. out = (1-Θ)*y₀+Θ*y₁+Θ*(Θ-1)*((1-2Θ)*(y₁-y₀)+(Θ-1)*dt*k[1] + Θ*dt*k[2])
  @inbounds for i in eachindex(out)
    out[i] = (1-Θ)*y₀[i]+Θ*y₁[i]+Θ*(Θ-1)*((1-2Θ)*(y₁[i]-y₀[i])+(Θ-1)*dt*k[1][i] + Θ*dt*k[2][i])
  end
  out
end

@muladd function hermite_interpolant!(out,Θ,dt,y₀,y₁,k,idxs,T::Type{Val{0}}) # Default interpolant is Hermite
  #@views @. out = (1-Θ)*y₀[idxs]+Θ*y₁[idxs]+Θ*(Θ-1)*((1-2Θ)*(y₁[idxs]-y₀[idxs])+(Θ-1)*dt*k[1][idxs] + Θ*dt*k[2][idxs])
  @inbounds for (j,i) in enumerate(idxs)
    out[j] = (1-Θ)*y₀[i]+Θ*y₁[i]+Θ*(Θ-1)*((1-2Θ)*(y₁[i]-y₀[i])+(Θ-1)*dt*k[1][i] + Θ*dt*k[2][i])
  end
  out
end

"""
Herimte Interpolation, chosen if no other dispatch for ode_interpolant
"""
@muladd function hermite_interpolant(Θ,dt,y₀,y₁,k,idxs::Nothing,T::Type{Val{1}}) # Default interpolant is Hermite
  #@. k[1] + Θ*(-4*dt*k[1] - 2*dt*k[2] - 6*y₀ + Θ*(3*dt*k[1] + 3*dt*k[2] + 6*y₀ - 6*y₁) + 6*y₁)/dt
  k[1] + Θ*(-4*dt*k[1] - 2*dt*k[2] - 6*y₀ + Θ*(3*dt*k[1] + 3*dt*k[2] + 6*y₀ - 6*y₁) + 6*y₁)/dt
end

@muladd function hermite_interpolant(Θ,dt,y₀,y₁,k,idxs,T::Type{Val{1}}) # Default interpolant is Hermite
  # return @. k[1][idxs] + Θ*(-4*dt*k[1][idxs] - 2*dt*k[2][idxs] - 6*y₀[idxs] + Θ*(3*dt*k[1][idxs] + 3*dt*k[2][idxs] + 6*y₀[idxs] - 6*y₁[idxs]) + 6*y₁[idxs])/dt
  return k[1][idxs] + Θ*(-4*dt*k[1][idxs] - 2*dt*k[2][idxs] - 6*y₀[idxs] + Θ*(3*dt*k[1][idxs] + 3*dt*k[2][idxs] + 6*y₀[idxs] - 6*y₁[idxs]) + 6*y₁[idxs])/dt
end

@muladd function hermite_interpolant!(out,Θ,dt,y₀,y₁,k,idxs::Nothing,T::Type{Val{1}}) # Default interpolant is Hermite
  #@. out = k[1] + Θ*(-4*dt*k[1] - 2*dt*k[2] - 6*y₀ + Θ*(3*dt*k[1] + 3*dt*k[2] + 6*y₀ - 6*y₁) + 6*y₁)/dt
  @inbounds for i in eachindex(out)
    out[i] = k[1][i] + Θ*(-4*dt*k[1][i] - 2*dt*k[2][i] - 6*y₀[i] + Θ*(3*dt*k[1][i] + 3*dt*k[2][i] + 6*y₀[i] - 6*y₁[i]) + 6*y₁[i])/dt
  end
  out
end

@muladd function hermite_interpolant!(out,Θ,dt,y₀,y₁,k,idxs,T::Type{Val{1}}) # Default interpolant is Hermite
  #@views @. out = k[1][idxs] + Θ*(-4*dt*k[1][idxs] - 2*dt*k[2][idxs] - 6*y₀[idxs] + Θ*(3*dt*k[1][idxs] + 3*dt*k[2][idxs] + 6*y₀[idxs] - 6*y₁[idxs]) + 6*y₁[idxs])/dt
  @inbounds for (j,i) in enumerate(idxs)
    out[j] = k[1][i] + Θ*(-4*dt*k[1][i] - 2*dt*k[2][i] - 6*y₀[i] + Θ*(3*dt*k[1][i] + 3*dt*k[2][i] + 6*y₀[i] - 6*y₁[i]) + 6*y₁[i])/dt
  end
  out
end

"""
Herimte Interpolation, chosen if no other dispatch for ode_interpolant
"""
@muladd function hermite_interpolant(Θ,dt,y₀,y₁,k,idxs::Nothing,T::Type{Val{2}}) # Default interpolant is Hermite
  #@. (-4*dt*k[1] - 2*dt*k[2] - 6*y₀ + Θ*(6*dt*k[1] + 6*dt*k[2] + 12*y₀ - 12*y₁) + 6*y₁)/(dt*dt)
  (-4*dt*k[1] - 2*dt*k[2] - 6*y₀ + Θ*(6*dt*k[1] + 6*dt*k[2] + 12*y₀ - 12*y₁) + 6*y₁)/(dt*dt)
end

@muladd function hermite_interpolant(Θ,dt,y₀,y₁,k,idxs,T::Type{Val{2}}) # Default interpolant is Hermite
  #out = similar(y₀,axes(idxs))
  #@views @. out = (-4*dt*k[1][idxs] - 2*dt*k[2][idxs] - 6*y₀[idxs] + Θ*(6*dt*k[1][idxs] + 6*dt*k[2][idxs] + 12*y₀[idxs] - 12*y₁[idxs]) + 6*y₁[idxs])/(dt*dt)
  @views out = (-4*dt*k[1][idxs] - 2*dt*k[2][idxs] - 6*y₀[idxs] + Θ*(6*dt*k[1][idxs] + 6*dt*k[2][idxs] + 12*y₀[idxs] - 12*y₁[idxs]) + 6*y₁[idxs])/(dt*dt)
  out
end

@muladd function hermite_interpolant!(out,Θ,dt,y₀,y₁,k,idxs::Nothing,T::Type{Val{2}}) # Default interpolant is Hermite
  #@. out = (-4*dt*k[1] - 2*dt*k[2] - 6*y₀ + Θ*(6*dt*k[1] + 6*dt*k[2] + 12*y₀ - 12*y₁) + 6*y₁)/(dt*dt)
  @inbounds for i in eachindex(out)
    out[i] = (-4*dt*k[1][i] - 2*dt*k[2][i] - 6*y₀[i] + Θ*(6*dt*k[1][i] + 6*dt*k[2][i] + 12*y₀[i] - 12*y₁[i]) + 6*y₁[i])/(dt*dt)
  end
  out
end

@muladd function hermite_interpolant!(out,Θ,dt,y₀,y₁,k,idxs,T::Type{Val{2}}) # Default interpolant is Hermite
  #@views @. out = (-4*dt*k[1][idxs] - 2*dt*k[2][idxs] - 6*y₀[idxs] + Θ*(6*dt*k[1][idxs] + 6*dt*k[2][idxs] + 12*y₀[idxs] - 12*y₁[idxs]) + 6*y₁[idxs])/(dt*dt)
  @inbounds for (j,i) in enumerate(idxs)
    out[j] = (-4*dt*k[1][i] - 2*dt*k[2][i] - 6*y₀[i] + Θ*(6*dt*k[1][i] + 6*dt*k[2][i] + 12*y₀[i] - 12*y₁[i]) + 6*y₁[i])/(dt*dt)
  end
  out
end

"""
Herimte Interpolation, chosen if no other dispatch for ode_interpolant
"""
@muladd function hermite_interpolant(Θ,dt,y₀,y₁,k,idxs::Nothing,T::Type{Val{3}}) # Default interpolant is Hermite
  #@. (6*dt*k[1] + 6*dt*k[2] + 12*y₀ - 12*y₁)/(dt*dt*dt)
  (6*dt*k[1] + 6*dt*k[2] + 12*y₀ - 12*y₁)/(dt*dt*dt)
end

@muladd function hermite_interpolant(Θ,dt,y₀,y₁,k,idxs,T::Type{Val{3}}) # Default interpolant is Hermite
  #out = similar(y₀,axes(idxs))
  #@views @. out = (6*dt*k[1][idxs] + 6*dt*k[2][idxs] + 12*y₀[idxs] - 12*y₁[idxs])/(dt*dt*dt)
  @views out = (6*dt*k[1][idxs] + 6*dt*k[2][idxs] + 12*y₀[idxs] - 12*y₁[idxs])/(dt*dt*dt)
  out
end

@muladd function hermite_interpolant!(out,Θ,dt,y₀,y₁,k,idxs::Nothing,T::Type{Val{3}}) # Default interpolant is Hermite
  # @. out = (6*dt*k[1] + 6*dt*k[2] + 12*y₀ - 12*y₁)/(dt*dt*dt)
  for i in eachindex(out)
    out[i] = (6*dt*k[1][i] + 6*dt*k[2][i] + 12*y₀[i] - 12*y₁[i])/(dt*dt*dt)
  end
  out
end

@muladd function hermite_interpolant!(out,Θ,dt,y₀,y₁,k,idxs,T::Type{Val{3}}) # Default interpolant is Hermite
  #@views @. out = (6*dt*k[1][idxs] + 6*dt*k[2][idxs] + 12*y₀[idxs] - 12*y₁[idxs])/(dt*dt*dt)
  for (j,i) in enumerate(idxs)
    out[j] = (6*dt*k[1][i] + 6*dt*k[2][i] + 12*y₀[i] - 12*y₁[i])/(dt*dt*dt)
  end
  out
end

######################## Linear Interpolants

@muladd function linear_interpolant(Θ,dt,y₀,y₁,idxs::Nothing,T::Type{Val{0}})
  Θm1 = (1-Θ)
  @. Θm1*y₀ + Θ*y₁
end

@muladd function linear_interpolant(Θ,dt,y₀,y₁,idxs,T::Type{Val{0}})
  Θm1 = (1-Θ)
  @. Θm1*y₀[idxs] + Θ*y₁[idxs]
end

@muladd function linear_interpolant!(out,Θ,dt,y₀,y₁,idxs::Nothing,T::Type{Val{0}})
  Θm1 = (1-Θ)
  @. out = Θm1*y₀ + Θ*y₁
  out
end

@muladd function linear_interpolant!(out,Θ,dt,y₀,y₁,idxs,T::Type{Val{0}})
  Θm1 = (1-Θ)
  @views @. out = Θm1*y₀[idxs] + Θ*y₁[idxs]
  out
end

"""
Linear Interpolation
"""
function linear_interpolant(Θ,dt,y₀,y₁,idxs::Nothing,T::Type{Val{1}})
  @. (y₁ - y₀)/dt
end

function linear_interpolant(Θ,dt,y₀,y₁,idxs,T::Type{Val{1}})
  @. (y₁[idxs] - y₀[idxs])/dt
end

function linear_interpolant!(out,Θ,dt,y₀,y₁,idxs::Nothing,T::Type{Val{1}})
  @. out = (y₁ - y₀)/dt
  out
end

function linear_interpolant!(out,Θ,dt,y₀,y₁,idxs,T::Type{Val{1}})
  @views @. out = (y₁[idxs] - y₀[idxs])/dt
  out
end


text:
const binary_operators = Set()

struct BinaryOperation{LX, LY, LZ, O, A, B, IA, IB, G, T} <: AbstractOperation{LX, LY, LZ, G, T}
    op :: O
    a :: A
    b :: B
    ▶a :: IA
    ▶b :: IB
    grid :: G

    @doc """
        BinaryOperation{LX, LY, LZ}(op, a, b, ▶a, ▶b, grid)

    Returns an abstract representation of the binary operation `op(▶a(a), ▶b(b))`.
    on `grid`, where `▶a` and `▶b` interpolate `a` and `b` to locations `(LX, LY, LZ)`.
    """
    function BinaryOperation{LX, LY, LZ}(op::O, a::A, b::B, ▶a::IA, ▶b::IB, grid::G) where {LX, LY, LZ, O, A, B, IA, IB, G}
        T = eltype(grid)
        return new{LX, LY, LZ, O, A, B, IA, IB, G, T}(op, a, b, ▶a, ▶b, grid)
    end
end

@inline Base.getindex(β::BinaryOperation, i, j, k) = β.op(i, j, k, β.grid, β.▶a, β.▶b, β.a, β.b)

#####
##### BinaryOperation construction
#####

# Recompute location of binary operation
@inline at(loc, β::BinaryOperation) = β.op(loc, at(loc, β.a), at(loc, β.b))

"""Create a binary operation for `op` acting on `a` and `b` at `Lc`, where
`a` and `b` have location `La` and `Lb`."""
function _binary_operation(Lc, op, a, b, La, Lb, grid)
     ▶a = interpolation_operator(La, Lc)
     ▶b = interpolation_operator(Lb, Lc)
    return BinaryOperation{Lc[1], Lc[2], Lc[3]}(op, a, b, ▶a, ▶b, grid)
end

const ConcreteLocationType = Union{Type{Face}, Type{Center}}

# Precedence rules for choosing operation location:
choose_location(La, Lb, Lc) = Lc                                    # Fallback to the specification Lc, but also...
choose_location(::Type{Face},   ::Type{Face},   Lc) = Face          # keep common locations; and
choose_location(::Type{Center}, ::Type{Center}, Lc) = Center        #
choose_location(La::ConcreteLocationType, ::Type{Nothing}, Lc) = La # don't interpolate unspecified locations.
choose_location(::Type{Nothing}, Lb::ConcreteLocationType, Lc) = Lb #

"""Return an expression that defines an abstract `BinaryOperator` named `op` for `AbstractField`."""
function define_binary_operator(op)
    return quote
        import Oceananigans.Grids: AbstractGrid
        import Oceananigans.Fields: AbstractField

        local location = Oceananigans.Fields.location
        local FunctionField = Oceananigans.Fields.FunctionField
        local AF = AbstractField

        @inline $op(i, j, k, grid::AbstractGrid, ▶a, ▶b, a, b) =
            @inbounds $op(▶a(i, j, k, grid, a), ▶b(i, j, k, grid, b))

        # These shenanigans seem to help / encourage the compiler to infer types of objects
        # buried in deep AbstractOperations trees.
        @inline function $op(i, j, k, grid::AbstractGrid, ▶a, ▶b, A::BinaryOperation, B::BinaryOperation)
            @inline a(ii, jj, kk, grid) = A.op(A.▶a(ii, jj, kk, grid, A.a), A.▶b(ii, jj, kk, grid, A.b))
            @inline b(ii, jj, kk, grid) = B.op(B.▶a(ii, jj, kk, grid, B.a), B.▶b(ii, jj, kk, grid, B.b))
            return @inbounds $op(▶a(i, j, k, grid, a), ▶b(i, j, k, grid, b))
        end

        @inline function $op(i, j, k, grid::AbstractGrid, ▶a, ▶b, A::BinaryOperation, B::AbstractField)
            @inline a(ii, jj, kk, grid) = A.op(A.▶a(ii, jj, kk, grid, A.a), A.▶b(ii, jj, kk, grid, A.b))
            return @inbounds $op(▶a(i, j, k, grid, a), ▶b(i, j, k, grid, B))
        end

        @inline function $op(i, j, k, grid::AbstractGrid, ▶a, ▶b, A::AbstractField, B::BinaryOperation)
            @inline b(ii, jj, kk, grid) = B.op(B.▶a(ii, jj, kk, grid, B.a), B.▶b(ii, jj, kk, grid, B.b))
            return @inbounds $op(▶a(i, j, k, grid, A), ▶b(i, j, k, grid, b))
        end

        @inline function $op(i, j, k, grid::AbstractGrid, ▶a, ▶b, A::BinaryOperation, B::Number)
            @inline a(ii, jj, kk, grid) = A.op(A.▶a(ii, jj, kk, grid, A.a), A.▶b(ii, jj, kk, grid, A.b))
            return @inbounds $op(▶a(i, j, k, grid, a), B)
        end

        @inline function $op(i, j, k, grid::AbstractGrid, ▶a, ▶b, A::Number, B::BinaryOperation)
            @inline b(ii, jj, kk, grid) = B.op(B.▶a(ii, jj, kk, grid, B.a), B.▶b(ii, jj, kk, grid, B.b))
            return @inbounds $op(A, ▶b(i, j, k, grid, b))
        end

        """
            $($op)(Lc, a, b)

        Returns an abstract representation of the operator `$($op)` acting on `a` and `b`.
        The operation occurs at location(a) except for Nothing dimensions. In that case,
        the location of the dimension in question is supplied either by location(b) or
        if that is also Nothing, Lc.
        """
        function $op(Lc::Tuple, a, b)
            La = location(a)
            Lb = location(b)
            Lab = choose_location.(La, Lb, Lc)

            grid = Oceananigans.AbstractOperations.validate_grid(a, b)

            return Oceananigans.AbstractOperations._binary_operation(Lab, $op, a, b, La, Lb, grid)
        end

        # Numbers are not fields...
        $op(Lc::Tuple, a::Number, b::Number) = $op(a, b)

        # Sugar for mixing in functions of (x, y, z)
        $op(Lc::Tuple, f::Function, b::AbstractField) = $op(Lc, FunctionField(location(b), f, b.grid), b)
        $op(Lc::Tuple, a::AbstractField, f::Function) = $op(Lc, a, FunctionField(location(a), f, a.grid))

        $op(Lc::Tuple, m::AbstractGridMetric, b::AbstractField) = $op(Lc, GridMetricOperation(location(b), m, b.grid), b)
        $op(Lc::Tuple, a::AbstractField, m::AbstractGridMetric) = $op(Lc, a, GridMetricOperation(location(a), m, a.grid))

        # Sugary versions with default locations
        $op(a::AF, b::AF) = $op(location(a), a, b)
        $op(a::AF, b) = $op(location(a), a, b)
        $op(a, b::AF) = $op(location(b), a, b)

        $op(a::AF, b::Number) = $op(location(a), a, b)
        $op(a::Number, b::AF) = $op(location(b), a, b)
    end
end

"""
    @binary op1 op2 op3...

Turn each binary function in the list `(op1, op2, op3...)`
into a binary operator on `Oceananigans.Fields` for use in `AbstractOperations`.

Note: a binary function is a function with two arguments: for example, `+(x, y)` is a binary function.

Also note: a binary function in `Base` must be imported to be extended: use `import Base: op; @binary op`.

Example
=======

```jldoctest
julia> using Oceananigans, Oceananigans.AbstractOperations

julia> using Oceananigans.AbstractOperations: BinaryOperation, AbstractGridMetric, choose_location

julia> plus_or_times(x, y) = x < 0 ? x + y : x * y
plus_or_times (generic function with 1 method)

julia> @binary plus_or_times
Set{Any} with 6 elements:
  :+
  :/
  :^
  :-
  :*
  :plus_or_times

julia> c, d = (CenterField(RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1))) for i = 1:2);

julia> plus_or_times(c, d)
BinaryOperation at (Center, Center, Center)
├── grid: 1×1×1 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
└── tree:
    plus_or_times at (Center, Center, Center)
    ├── 1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU
    └── 1×1×1 Field{Center, Center, Center} on RectilinearGrid on CPU
```
"""
macro binary(ops...)
    expr = Expr(:block)

    for op in ops
        defexpr = define_binary_operator(op)
        push!(expr.args, :($(esc(defexpr))))

        add_to_operator_lists = quote
            push!(Oceananigans.AbstractOperations.operators, Symbol($op))
            push!(Oceananigans.AbstractOperations.binary_operators, Symbol($op))
        end

        push!(expr.args, :($(esc(add_to_operator_lists))))
    end

    return expr
end

#####
##### Nested computations
#####

function compute_at!(β::BinaryOperation, time)
    compute_at!(β.a, time)
    compute_at!(β.b, time)
    return nothing
end

#####
##### GPU capabilities
#####

"Adapt `BinaryOperation` to work on the GPU via CUDAnative and CUDAdrv."
Adapt.adapt_structure(to, binary::BinaryOperation{LX, LY, LZ}) where {LX, LY, LZ} =
    BinaryOperation{LX, LY, LZ}(Adapt.adapt(to, binary.op),
                                Adapt.adapt(to, binary.a),
                                Adapt.adapt(to, binary.b),
                                Adapt.adapt(to, binary.▶a),
                                Adapt.adapt(to, binary.▶b),
                                Adapt.adapt(to, binary.grid))


text:
###======================================================####
#Came from GLAbstraction/GLUniforms.jl

# Uniforms are OpenGL variables that stay the same for the entirety of a drawcall.
# There are a lot of functions, to upload them, as OpenGL doesn't rely on multiple dispatch.
# here is my approach, to handle all of the uniforms with one function, namely gluniform
# For uniforms, the Vector and Matrix types from ImmutableArrays should be used, as they map the relation almost 1:1

#Some additional uniform functions, not related to Imutable Arrays
# gluniform(location::Integer, target::Integer, t::Texture) = gluniform(GLint(location), GLint(target), t)
# gluniform(location::Integer, target::Integer, t::GPUVector) = gluniform(GLint(location), GLint(target), t.buffer)
# gluniform(location::Integer, target::Integer, t::TextureBuffer) = gluniform(GLint(location), GLint(target), t.texture)
# gluniform(location::Integer, t::TextureBuffer) = gluniform(GLint(location), GLint(target), t.texture)
#REVIEW: scary, binding and making texture active seems like something that shouldn't be in gluniform...
function gluniform(location::GLint, texture_unit, t::Texture)
    tu = GL_TEXTURE0 + UInt32(texture_unit)
    glActiveTexture(tu)
    glBindTexture(t.texturetype, t.id)
    gluniform(location, texture_unit)
end
gluniform(location::Integer, x::Enum)                                = gluniform(GLint(location), GLint(x))
gluniform(location::Integer, x::Union{GLubyte, GLushort, GLuint})    = glUniform1ui(GLint(location), x)
gluniform(location::Integer, x::Union{GLbyte, GLshort, GLint, Bool}) = glUniform1i(GLint(location),  x)
gluniform(location::Integer, x::GLfloat)                             = glUniform1f(GLint(location),  x)
gluniform(location::Integer, x::GLdouble)                            = glUniform1d(GLint(location),  x)

#Uniform upload functions for julia arrays...
gluniform(location::GLint, x::Vector{Float32})                       = glUniform1fv(location, length(x), x)
gluniform(location::GLint, x::Vector{GLdouble})                      = glUniform1dv(location, length(x), x)
gluniform(location::GLint, x::Vector{GLint})                         = glUniform1iv(location, length(x), x)
gluniform(location::GLint, x::Vector{GLuint})                        = glUniform1uiv(location, length(x), x)


text:
# computes the variance of
#
#     ω'(H Δ - (Sy - Sx))
#
# where H = (Sx ⊗ Sy + Sy ⊗ Sx) / 2
# and ω, Δ are fixed
function variance(
    ::SymmetricOracleNormal,
    Sx::Symmetric,
    Sy::Symmetric,
    X::AbstractMatrix,
    Y::AbstractMatrix,
    ω::Vector,
    indω::Vector{Int64},
    Δ::Vector,
    indΔ::Vector{Int64}
    )

    nx, p = size(X)
    ny    = size(Y, 1)

    q = zeros(nx)
    r = zeros(ny)

    # compute the value of the U-statistics
    t = 0.
    for ci=1:length(indΔ)
        for ri=1:length(indω)
            @inbounds t += _getElemSKron(Sx, Sy, indω[ri], indΔ[ci]) * ω[ri] * Δ[ci]
        end
    end
    for ri=1:length(indω)
        @inbounds t -= ω[ri] * (Sy[indω[ri]] - Sx[indω[ri]])
    end

    # compute qk
    for k=1:nx
        v = 0.
        for ci=1:length(indΔ)
            for ri=1:length(indω)
                @inbounds v += _getElemSKron(Sy, view(X, k, :), indω[ri], indΔ[ci]) * ω[ri] * Δ[ci]
            end
        end
        for ri=1:length(indω)
            row, col = Tuple(CartesianIndices((p, p))[indω[ri]])
            @inbounds v -= ω[ri] * (Sy[indω[ri]] - X[k, row] * X[k, col])
        end
        @inbounds q[k] = v
    end

    # compute rk
    for k=1:ny
        v = 0.
        for ci=1:length(indΔ)
            for ri=1:length(indω)
                @inbounds v += _getElemSKron(Sx, view(Y, k, :), indω[ri], indΔ[ci]) * ω[ri] * Δ[ci]
            end
        end
        for ri=1:length(indω)
            row, col = Tuple(CartesianIndices((p, p))[indω[ri]])
            @inbounds v -= ω[ri] * (Y[k, row] * Y[k, col] - Sx[indω[ri]])
        end
        @inbounds r[k] = v
    end

    σ1 = sum(abs2, q) / (nx - 1) - nx / (nx - 1) * t^2
    σ2 = sum(abs2, r) / (ny - 1) - ny / (ny - 1) * t^2

    return σ1/nx + σ2/ny
end

# computes the variance of
#
#     ω'(H Δ - (Sy - Sx))
#
# where H = (Sx ⊗ Sy + Sy ⊗ Sx) / 2
# and ω, Δ are fixed
function variance(
  ::ReducedOracleNormal,
  Sx::Symmetric,
  Sy::Symmetric,
  X::AbstractMatrix,
  Y::AbstractMatrix,
  ω::Vector,
  indω::Vector{Int64},
  Δ::Vector,
  indΔ::Vector{Int64}
  )

  nx = size(X, 1)
  ny = size(Y, 1)

  q = zeros(nx)
  r = zeros(ny)

  # compute the value of the U-statistics
  t = 0.
  for ci=1:length(indΔ)
    for ri=1:length(indω)
      @inbounds t += skron(Sx, Sy, indω[ri], indΔ[ci]) * ω[ri] * Δ[ci]
    end
  end
  for ri=1:length(indω)
    @inbounds t -= ω[ri] * (svec(Sy, indω[ri]) - svec(Sx, indω[ri]))
  end

  # compute qk
  for k=1:nx
    v = 0.
    vx = view(X, k, :)
    for ci=1:length(indΔ)
      for ri=1:length(indω)
        @inbounds v += skron(Sy, vx, indω[ri], indΔ[ci]) * ω[ri] * Δ[ci]
      end
    end
    for ri=1:length(indω)
      @inbounds v -= ω[ri] * (svec(Sy, indω[ri]) - svec(vx, vx, indω[ri]))
    end
    @inbounds q[k] = v
  end

  # compute rk
  for k=1:ny
    v = 0.
    vy = view(Y, k, :)
    for ci=1:length(indΔ)
      for ri=1:length(indω)
        @inbounds v += skron(Sx, vy, indω[ri], indΔ[ci]) * ω[ri] * Δ[ci]
      end
    end
    for ri=1:length(indω)
      @inbounds v -= ω[ri] * (svec(vy, vy, indω[ri]) - svec(Sx, indω[ri]))
    end
    @inbounds r[k] = v
  end

  σ1 = sum(abs2, q) / (nx - 1) - nx / (nx - 1) * t^2
  σ2 = sum(abs2, r) / (ny - 1) - ny / (ny - 1) * t^2

  return σ1/nx + σ2/ny
end

# computes the variance of
#
#     ω'(H Δ - (Sy - Sx))
#
# where H = Sy ⊗ Sx
# and ω, Δ are fixed
function variance(
    ::AsymmetricOracleNormal,
    Sx::Symmetric,
    Sy::Symmetric,
    X::AbstractMatrix,
    Y::AbstractMatrix,
    ω::Vector,
    indω::Vector{Int64},
    Δ::Vector,
    indΔ::Vector{Int64}
    )

    nx, p = size(X)
    ny    = size(Y, 1)

    q = zeros(nx)
    r = zeros(ny)

    # compute the value of the U-statistics
    t = 0.
    for ci=1:length(indΔ)
        for ri=1:length(indω)
            @inbounds t += _getElemKron(Sy, Sx, indω[ri], indΔ[ci]) * ω[ri] * Δ[ci]
        end
    end
    for ri=1:length(indω)
        @inbounds t -= ω[ri] * (Sy[indω[ri]] - Sx[indω[ri]])
    end

    # compute qk
    for k=1:nx
        v = 0.
        for ci=1:length(indΔ)
            for ri=1:length(indω)
                @inbounds v += _getElemKron(Sy, view(X, k, :), indω[ri], indΔ[ci]) * ω[ri] * Δ[ci]
            end
        end
        for ri=1:length(indω)
            row, col = Tuple(CartesianIndices((p, p))[indω[ri]])
            @inbounds v -= ω[ri] * (Sy[indω[ri]] - X[k, row] * X[k, col])
        end
        @inbounds q[k] = v
    end

    # compute rk
    for k=1:ny
        v = 0.
        for ci=1:length(indΔ)
            for ri=1:length(indω)
                @inbounds v += _getElemKron(view(Y, k, :), Sx, indω[ri], indΔ[ci]) * ω[ri] * Δ[ci]
            end
        end
        for ri=1:length(indω)
            row, col = Tuple(CartesianIndices((p, p))[indω[ri]])
            @inbounds v -= ω[ri] * (Y[k, row] * Y[k, col] - Sx[indω[ri]])
        end
        @inbounds r[k] = v
    end

    σ1 = sum(abs2, q) / (nx - 1) - nx / (nx - 1) * t^2
    σ2 = sum(abs2, r) / (ny - 1) - ny / (ny - 1) * t^2

    return σ1/nx + σ2/ny
end


text:
using ZMQ
using Requests

function run()
    @sync for i=1:10
        @async begin
            ctx = Context()
            sock = Socket(ctx, PULL)
            connect(sock, "tcp://localhost:5501")
            while true
                try
                    url = bytestring(recv(sock))
                    r = get(url)
                    println("OK: status=$(r.status)")
                catch
                    println("Error")
                end
            end
        end
    end
end


text:
"""
Helper functions to visualize graphs in [dot language][dot].

[dot]: https://en.wikipedia.org/wiki/DOT_(graph_description_language)
"""
module GraphViz

import GraphIdx: Graph, Tree
import GraphIdx.Grid: GridGraph

const INDENT = " "^4

dot(g::Graph) = dot(stdout, g)

dot(io::IO, g::Graph) = error("dot(::IO, $(typeof(g))) not implemented yet")

dot(print_edges::F, io::IO, g::GridGraph) where {F<:Function} =
    dot(io, g, print_edges)

"""
    dot([print_edges], [io], g)

Visualize graph `g` in dot language.

Example
-------

```julia
julia> open(`dot -Tx11`, "w") do io
           GraphIdx.GraphViz.dot(io, GraphIdx.Grid.GridGraph(3, 5))
       end

```
"""
function dot(io::IO, g::GridGraph, print_edges::F = io -> nothing) where {F<:Function}
    pix2ind(i, j) = i + (j-1)*g.n1
    println(io, "digraph {")
    for attr in ["rankdir=LR", "ranksep=0.3", "splines=ortho"]
        println(io, INDENT, attr)
    end
    println(io, INDENT, "edge [dir=none constraint=false colorscheme=dark28]")
    println(io, INDENT, "node [width=0.3 shape=circle fixedsize=true]")
    for i2 in 1:g.n2
        println(io, INDENT, "{ rank=same")
        for i1 in g.n1:-1:1
            println(io, INDENT, INDENT, pix2ind(i1, i2), " [pos=\"$i2,-$(i1)!\"]")
        end
        println(io, INDENT, "}")
    end
    print_edges(io)
    println(io, INDENT, "{ edge [style=invis constraint=true]")
    print(io, INDENT, INDENT, pix2ind(g.n1, 1))
    for i2 in 2:g.n2
        print(io, " -> ", pix2ind(g.n1, i2))
    end
    println(io)
    println(io, INDENT, "}")
    println(io, "}")
end

function dot(graph::Graph, tree::Tree.RootedTree)
    open(`dot -Tx11 -Kneato`, "w") do io::IO
        dot(io, graph) do io::IO
            for (i, v) in enumerate(tree.parent)
                if i != v
                    println(io, i, " -> ", v)
                end
            end
        end
    end
end

end

max_stars_repo_name: EQt/graphidx

max_issues_repo_name: EQt/graphidx

max_forks_repo_name: EQt/graphidx


text:
using Turing, Random, Test
using StatsFuns

dir = splitdir(splitdir(pathof(Turing))[1])[1]
include(dir*"/test/test_utils/AllUtils.jl")

@turing_testset "is.jl" begin
    function reference(n :: Int)
        logweights = zeros(Float64, n)
        samples = Array{Dict{Symbol,Any}}(undef, n)
        for i = 1:n
            samples[i] = reference()
            logweights[i] = samples[i][:logweight]
        end
        logevidence = logsumexp(logweights) - log(n)
        results = Dict{Symbol,Any}()
        results[:lp] = logevidence
        results[:lp] = logweights
        results[:samples] = samples
        return results
    end

    function reference()
        x = rand(Normal(4,5))
        y = rand(Normal(x,1))
        log_lik = logpdf(Normal(x,2), 3) + logpdf(Normal(y,2), 1.5)
        d = Dict()
        d[:logweight] = log_lik
        d[:a] = x
        d[:b] = y
        return d
    end

    let n = 10
      @model normal() = begin
        a ~ Normal(4,5)
        3 ~ Normal(a,2)
        b ~ Normal(a,1)
        1.5 ~ Normal(b,2)
        a, b
      end

      alg = IS(n)
      seed = 0

      _f = normal();
      for i=1:100
        Random.seed!(seed)
        # Move the rng twice - equivalent to sampling from prior
        rand(Normal(4,5))
        rand(Normal(0,1))
        exact = reference(n)
        Random.seed!(seed)
        tested = sample(_f, alg)
        t_vals = get(tested, [:a, :b, :lp])
        for i = 1:n
            @test exact[:samples][i][:a] == t_vals.a[i]
            @test exact[:samples][i][:b] == t_vals.b[i]
            @test exact[:lp][i] == t_vals.lp[i]
        end
        @test all(exact[:lp] .== t_vals.lp)
      end
    end
end


text:
using CSV#= as of 29 Aug 2021 one needs the main version of CSV, latest released version
0.8.5 fails to load some of the files. See https://github.com/JuliaData/CSV.jl/issues/879
Use:
]add CSV#main
=#
using DataFrames

function benchmark()
    #= folder points to where the VBA method RunSpeedTests writes performance test files.
    See workbooks/VBA-CSV.xlsm, module modCSVPerformance =#

    folder = joinpath(ENV["TEMP"],"VBA-CSV/Performance")
    outputfile = normpath(joinpath(@__DIR__,"..","julia/juliaparsetimes.csv"))
 outputfile = outputfile * "2"
    benchmark_csvs_in_folder(folder,outputfile)
end

"""
   benchmark_csvs_in_folder(folder::String, outputfile::String)
Benchmark all .csv files in `folder`, writing results to `outputfile`.
"""
function benchmark_csvs_in_folder(folder::String, outputfile::String)
    files = readdir(folder, join=true)
    files = filter(x -> x[end - 3:end] == ".csv", files)# only .csv files

    n = length(files)
    times = fill(0.0, n)
    numcalls = fill(0, n)
    statuses = fill("OK", n)

    foo = benchmarkonefile(files[1], 1)# for compilation "warmup"
    for (f, i) in zip(files, 1:n)
        println(i, f)
        try
            times[i], numcalls[i] = benchmarkonefile(f, 5)
        catch e
            statuses[i] = "$e"
        end
    end
    times

    result = DataFrame(filename=replace.(files, "/" => "\\"), time=times,
                        status=statuses, numcalls=numcalls)
    CSV.write(outputfile, result)
end

"""
    benchmarkonefile(filename::String, timeout::Int)
Average time (over sufficient trials to take `timeout` seconds) to load file `filename` to
a DataFrame, using CSV.File.
"""
function benchmarkonefile(filename::String, timeout::Int)
    i = 0 ; time2 = time() # needed to give variables scope outside the loop.
    time1 = time()
    while true
        i = i + 1
        res = CSV.File(filename, header=false, type=String) |> DataFrame
        time2 = time()
        time2 - time1 < timeout || break
    end
    (time2 - time1) / i, i
end

function shownotcompliantfile(filenum::String)
filename = "C:/Projects/VBA-CSV/testfiles/Not_RFC_4180_Compliant_" * filenum * ".csv"
res = CSV.File(filename, header=false, delim = ",",type=String) |>DataFrame
Matrix(res)
end


text:
#===============================================================================

    optical_test_single_channel.jl

===============================================================================#

import DataStructures.OrderedDict
import Memento.debug

## preset values
const MIN_EXCITATION_FLUORESCENCE           = 5120
const MAX_EXCITATION_FLUORESCENCE           = 384000
# const MIN_EXCITATION_FLUORESCENCE_MULTIPLE  = 3

#===============================================================================
    function definition >>
===============================================================================#

## called by dispatch()
function act(
    ::Type{Val{optical_test_single_channel}},
    req             ::Associative;
    out_format      ::OutputFormat = pre_json_output
    ## remove MySql dependency
    #
    # db_conn ::MySQL.MySQLHandle,
    # exp_id ::Integer,
    # calib_info ::Union{Integer,OrderedDict} =calib_info_AIR; ## not used for computation
    # start: arguments that might be passed by upstream code
    # well_nums ::AbstractVector =[],
)
    debug(logger, "at act(::Type{Val{optical_test_single_channel}})")

    ## remove MySql dependency
    #
    # step_ids = ["baseline_step_id", "excitation_step_id"]
    # req = OrderedDict(map(step_ids) do step_id
    #     ot_qry_2b = "SELECT fluorescence_value, well_num
    #         FROM fluorescence_data
    #         WHERE
    #             experiment_id = $exp_id AND
    #             step_id = $step_id AND
    #             cycle_num = 1 AND
    #             step_id is not NULL
    #             well_constraint
    #         ORDER BY well_num
    #     "
    #     ot_nt, fluo_well_nums = get_mysql_data_well(
    #         well_nums, ot_qry_2b, db_conn, false
    #     )
    #     step_id => ot_nt[:fluorescence_value]
    # end) # do step_id

    ## assuming the 2 values of `req` are the same in length (number of wells)
    all_valid = true
    results =
        map(eachindex(req[BASELINE_KEY][FLUORESCENCE_VALUE_KEY])) do well_i
            baseline   = req[BASELINE_KEY][FLUORESCENCE_VALUE_KEY][well_i]
            excitation = req[EXCITATION_KEY][FLUORESCENCE_VALUE_KEY][well_i]
            # valid =
            #    (excitation >= MIN_EXCITATION_FLUORESCENCE) &&
            #    (excitation / baseline >= MIN_EXCITATION_FLUORESCENCE_MULTIPLE) &&
            #    (excitation <= MAX_EXCITATION_FLUORESCENCE) # old
            valid =
                (excitation .>= MIN_EXCITATION_FLUORESCENCE) &&
                (baseline   .<  MIN_EXCITATION_FLUORESCENCE) &&
                (excitation .<= MAX_EXCITATION_FLUORESCENCE) ## Josh, 2016-08-15
            all_valid &= valid
            OrderedDict(
                :baseline   => baseline,
                :excitation => excitation,
                :valid      => valid
            )
        end ## do well_i
    const output = OrderedDict(
        :optical_data => results,
        :valid        => all_valid)
    return output |> out(out_format)
end ## act(::Type{Val{optical_test_single_channel}})

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""NamespaceSpec describes the attributes on a Namespace.

    IoK8sApiCoreV1NamespaceSpec(;
        finalizers=nothing,
    )

    - finalizers::Vector{String} : Finalizers is an opaque list of values that must be empty to permanently remove object from storage. More info: https://kubernetes.io/docs/tasks/administer-cluster/namespaces/
"""
mutable struct IoK8sApiCoreV1NamespaceSpec <: SwaggerModel
    finalizers::Any # spec type: Union{ Nothing, Vector{String} } # spec name: finalizers

    function IoK8sApiCoreV1NamespaceSpec(;finalizers=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1NamespaceSpec, Symbol("finalizers"), finalizers)
        setfield!(o, Symbol("finalizers"), finalizers)
        o
    end
end # type IoK8sApiCoreV1NamespaceSpec

const _property_map_IoK8sApiCoreV1NamespaceSpec = Dict{Symbol,Symbol}(Symbol("finalizers")=>Symbol("finalizers"))
const _property_types_IoK8sApiCoreV1NamespaceSpec = Dict{Symbol,String}(Symbol("finalizers")=>"Vector{String}")
Base.propertynames(::Type{ IoK8sApiCoreV1NamespaceSpec }) = collect(keys(_property_map_IoK8sApiCoreV1NamespaceSpec))
Swagger.property_type(::Type{ IoK8sApiCoreV1NamespaceSpec }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1NamespaceSpec[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1NamespaceSpec }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1NamespaceSpec[property_name]

function check_required(o::IoK8sApiCoreV1NamespaceSpec)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1NamespaceSpec }, name::Symbol, val)
end


text:
"""All graphs should inherit from this base class"""
abstract type Graph end

const Edge{I<:Integer} = Tuple{I, I}

num_nodes(g::Graph) =
    error("$(typeof(g)) should implement num_nodes")

num_edges(g::Graph) =
    error("$(typeof(g)) should implement num_edges")

enumerate_edges(f::Function, g::Graph) =
    error("$(typeof(g)) should implement enumerate_edges")