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

# The order in which files are included is important,
# since certain things need to be defined before others use them

## Interval type

if haskey(ENV, "IA_VALID") == true
    const validity_check = true
else
    const validity_check = false
end

abstract type AbstractInterval{T} <: Real end

struct Interval{T<:Real} <: AbstractInterval{T}
    lo :: T
    hi :: T

    function Interval{T}(a::Real, b::Real) where T<:Real

        a = _normalisezero(a)
        b = _normalisezero(b)

        if validity_check

            if is_valid_interval(a, b)
                return new(a, b)

            else
                @warn "Invalid input, empty interval is returned"
                return new(T(Inf), T(-Inf))
            end

        end

        new(a, b)

    end
end

@inline _normalisezero(a::Real) = ifelse(iszero(a) && signbit(a), copysign(a, 1), a)

## Outer constructors

Interval(a::T, b::T) where T<:Real = Interval{T}(a, b)
Interval(a::T) where T<:Real = Interval(a, a)
Interval(a::Tuple) = Interval(a...)
Interval(a::T, b::S) where {T<:Real, S<:Real} = Interval(promote(a,b)...)

## Concrete constructors for Interval, to effectively deal only with Float64,
# BigFloat or Rational{Integer} intervals.
Interval(a::T, b::T) where T<:Integer = Interval(float(a), float(b))

# Constructors for Irrational
# Single argument Irrational constructor are in IntervalArithmetic.jl
# as generated functions need to be define last.
Interval{T}(a::Irrational, b::Irrational) where {T<:Real} = Interval{T}(T(a, RoundDown), T(b, RoundUp))
Interval{T}(a::Irrational, b::Real) where {T<:Real} = Interval{T}(T(a, RoundDown), b)
Interval{T}(a::Real, b::Irrational) where {T<:Real} = Interval{T}(a, T(b, RoundUp))

Interval(a::Irrational, b::Irrational) = Interval{Float64}(a, b)
Interval(a::Irrational, b::Real) = Interval{Float64}(a, b)
Interval(a::Real, b::Irrational) = Interval{Float64}(a, b)

Interval(x::Interval) = x
Interval(x::Complex) = Interval(real(x)) + im*Interval(imag(x))

Interval{T}(x) where T = Interval(convert(T, x))

Interval{T}(x::Interval) where T = atomic(Interval{T}, x)

size(x::Interval) = (1,)

"""
    is_valid_interval(a::Real, b::Real)

Check if `(a, b)` constitute a valid interval
"""
function is_valid_interval(a::Real, b::Real)

    # println("isvalid()")

    if isnan(a) || isnan(b)
        return false
    end

    a > b && return false

    if a == Inf || b == -Inf
        return false
    end

    return true
end

"""
    interval(a, b)

`interval(a, b)` checks whether [a, b] is a valid `Interval`, using the (non-exported) `is_valid_interval` function. If so, then an `Interval(a, b)` object is returned; if not, a warning is printed and the empty interval is returned.
"""
function interval(a::T, b::S) where {T<:Real, S<:Real}
    if !is_valid_interval(a, b)
        @warn "Invalid input, empty interval is returned"
        return emptyinterval(promote_type(T, S))
    end

    return Interval(a, b)
end

interval(a::Real) = interval(a, a)
interval(a::Interval) = a

"Make an interval even if a > b"
function force_interval(a, b)
    a > b && return interval(b, a)
    return interval(a, b)
end

## Include files
include("special.jl")
include("macros.jl")
include("rounding_macros.jl")
include("rounding.jl")
include("conversion.jl")
include("precision.jl")
include("set_operations.jl")
include("arithmetic.jl")
include("functions.jl")
include("trigonometric.jl")
include("hyperbolic.jl")
include("complex.jl")

# Syntax for intervals

function ..(a::T, b::S) where {T, S}
    if !is_valid_interval(a, b)
        @warn "Invalid input, empty interval is returned"
        return emptyinterval(promote_type(T, S))
    end
    Interval(atomic(Interval{T}, a).lo, atomic(Interval{S}, b).hi)
end

function ..(a::T, b::Irrational{S}) where {T, S}
    if !is_valid_interval(a, b)
        @warn "Invalid input, empty interval is returned"
        return emptyinterval(promote_type(T, Irrational{S}))
    end
    R = promote_type(T, Irrational{S})
    Interval(atomic(Interval{R}, a).lo, R(b, RoundUp))
end

function ..(a::Irrational{T}, b::S) where {T, S}
    if !is_valid_interval(a, b)
        @warn "Invalid input, empty interval is returned"
        return emptyinterval(promote_type(Irrational{T}, S))
    end
    R = promote_type(Irrational{T}, S)
    return Interval(R(a, RoundDown), atomic(Interval{R}, b).hi)
end

function ..(a::Irrational{T}, b::Irrational{S}) where {T, S}
    return interval(a, b)
end

# ..(a::Integer, b::Integer) = interval(a, b)
# ..(a::Integer, b::Real) = interval(a, nextfloat(float(b)))
# ..(a::Real, b::Integer) = interval(prevfloat(float(a)), b)
#
# ..(a::Real, b::Real) = interval(prevfloat(float(a)), nextfloat(float(b)))

macro I_str(ex)  # I"[3,4]"
    @interval(ex)
end

a ± b = (a-b)..(a+b)
±(a::Interval, b) = (a.lo - b)..(a.hi + b)

"""
Computes the integer hash code for an `Interval` using the method for composite types used in `AutoHashEquals.jl`
"""
hash(x::Interval, h::UInt) = hash(x.hi, hash(x.lo, hash(:Interval, h)))


text:
module EwaldSum

include("Madelung.jl")

end # module

ext: jl


text:
##############################################################################################
module TestIEEE39bus
##############################################################################################
##############################################################################################
# imports and exports
using DataAssimilationBenchmarks.DeSolvers, DataAssimilationBenchmarks.L96
using JLD, Statistics

##############################################################################################
##############################################################################################
# Test to see if the model without noise reaches the synchronous steady state

function test_synchrony()
    try
        # load the observations
        path = joinpath(@__DIR__, "../src/data/time_series/")
        obs = load(path * "IEEE39bus_time_series_seed_0000_diff_0.000_tanl_0.01" *
                   "_nanl_05000_spin_1500_h_0.010.jld2")
        obs = obs["obs"][:, 3001:end]

        # take the standard deviation of the model state after warm up
        sd = std(obs, dims=2)
        if sum(sd .< 0.01) == 20
            true
        else
            false
        end
    catch
        false
    end
end

##############################################################################################
# end module

end

max_stars_repo_name: sukhreens/DataAssimilationBenchmarks

max_issues_repo_name: sukhreens/DataAssimilationBenchmarks

max_forks_repo_name: sukhreens/DataAssimilationBenchmarks


text:
abstract type Builder end

"""
    build(project, builder[, docs])

Build every document named in `docs` using `builder`. `docs` defaults
to `keys(project.sources)`, meaning every document will be built. If you only
want to rebuild previously built files, use [`rebuild`](#).
"""
function build(project::Project, builder::Builder)
    return project

    # Build all documents
    dirtypaths = addfiles!(project, project.sources)

    # Save to disk
    build(builder, project, dirtypaths)

    return builder
end

"""
    build(project)

Build project to a temporary directory with [`HTML`](#) format.
"""
build(project) = build(project, FileBuilder(HTML(), Path(mktempdir())))

function fullbuild(project, builder)
    paths = rewritesources!(project)
    build(builder, project, paths)
end

"""
    rebuild(project, builder)

Build previously built documents in projects using `builder`. Equivalent to
`build(project, builder, keys(project.outputs))`.
"""
rebuild(project, builder) = build(project, builder, keys(project.outputs))

"""
    struct FileBuilder <: Builder

Build every document to a file in `dir` using output `format`.
"""
struct FileBuilder <: Builder
    format::Format
    dir::AbstractPath
end

function build(builder::FileBuilder, project::Project, dirtypaths = keys(project.outputs))
    # TODO: make threadable for performance
    for p in collect(dirtypaths)
        buildtofile(project.outputs[p], p, builder.dir, builder.format)
    end

    # TODO: make threadable for performance
    for rewriter in project.rewriters
        postbuild(rewriter, project, builder)
    end
end

function buildtofile(xtree, p, dir, format)
    fullpath = withext(joinpath(dir, p), formatextension(format))
    try
        mkpath(parent(fullpath))
    catch end
    render!(fullpath, xtree, format)
end


text:
"""
Subroutine dicf

Given a sparse symmetric matrix A in compressed row storage,
this subroutine computes an incomplete Cholesky factorization.

Implementation of dicf is based on the Jones-Plassmann code.
Arrays indf and list define the data structure.
At the beginning of the computation of the j-th column,

    For k < j, indf[k] is the index of A for the first
    nonzero l[i,k] in the k-th column with i >= j.

    For k < j, list[i] is a pointer to a linked list of column
    indices k with i = L.rowval[indf[k]].

For the computation of the j-th column, the array indr records
the row indices. Hence, if nlj is the number of nonzeros in the
j-th column, then indr[1],...,indr[nlj] are the row indices.
Also, for i > j, indf[i] marks the row indices in the j-th
column so that indf[i] = 1 if l[i,j] is not zero.

MINPACK-2 Project. May 1998.
Argonne National Laboratory.
Chih-Jen Lin and Jorge J. More'.
"""
function dicf(n, nnz, L::TronSparseMatrixCSC, p, indr, indf, list, w)
    zero = 0.0
    insortf = 20

    info = 0
    @inbounds for j=1:n
        indf[j] = 0
        list[j] = 0
    end

    # Make room for L by moving A to the last n*p positions in a.

    np = n*p
    @inbounds for j=1:n+1
        L.colptr[j] = L.colptr[j] + np
    end
    @inbounds for j=nnz:-1:1
        L.rowval[np+j] = L.rowval[j]
        L.tril_vals[np+j] = L.tril_vals[j]
    end

    # Compute the incomplete Cholesky factorization.

    isj = L.colptr[1]
    L.colptr[1] = 1
    @inbounds for j=1:n

        # Load column j into the array w. The first and last elements
        # of the j-th column of A are a(isj) and a(iej).

        nlj = 0
        iej = L.colptr[j+1] - 1

        @inbounds for ip in isj:iej
            i = L.rowval[ip]
            w[i] = L.tril_vals[ip]
            nlj = nlj + 1
            indr[nlj] = i
            indf[i] = 1
        end

        # Exit if the current pivot is not positive.
        if L.diag_vals[j] <= zero
            info = -j
            return info
        end
        L.diag_vals[j] = sqrt(L.diag_vals[j])

        # Update column j using the previous columns.
        k = list[j]
        while k != 0
            isk = indf[k]
            iek = L.colptr[k+1]-1

            # Set lval to l[j,k].
            lval = L.tril_vals[isk]

            # Update indf and list.

            newk = list[k]
            isk = isk + 1
            if isk < iek
                indf[k] = isk
                list[k] = list[L.rowval[isk]]
                list[L.rowval[isk]] = k
            end
            k = newk

            # Compute the update a[i,i] <- a[i,j] - l[i,k]*l[j,k].
            # In this loop we pick up l[i,k] for k < j and i > j.
            @inbounds for ip=isk:iek
                i = L.rowval[ip]
                if indf[i] != 0
                    w[i] = w[i] - lval*L.tril_vals[ip]
                else
                    indf[i] = 1
                    nlj = nlj + 1
                    indr[nlj] = i
                    w[i] = -lval*L.tril_vals[ip]
                end
            end
        end

        # Compute the j-th column of L.
        @inbounds for k=1:nlj
            w[indr[k]] = w[indr[k]]/L.diag_vals[j]
        end

        # Set mlj to the number of nonzeros to be retained.
        mlj = min(iej-isj+1+p,nlj)
        kth = nlj - mlj + 1

        if nlj >= 1
            # Determine the kth smallest elements in the current column
            dsel2(nlj,w,indr,kth)

            # Sort the row indices of the selected elements. Insertion
            # sort is used for small arrays, and heap sort for larger
            # arrays. The sorting of the row indices is required so that
            # we can retrieve l[i,k] with i > k from indf[k].
            if mlj <= insortf
                insort(mlj,view(indr, kth:n))
            else
                ihsort(mlj,view(indr, kth:n))
            end
        end

        # Store the largest elements in L. The first and last elements
        # of the j-th column of L are L.tril_vals[newisj] and L.tril_vals[newiej].
        newisj = L.colptr[j]
        newiej = newisj + mlj - 1
        @inbounds for k in newisj:newiej
            L.tril_vals[k] = w[indr[k-newisj+kth]]
            L.rowval[k] = indr[k-newisj+kth]
        end

        # Update the diagonal elements.
        @inbounds for k=kth:nlj
            L.diag_vals[indr[k]] = L.diag_vals[indr[k]] - w[indr[k]]^2
        end

        # Update indr and list for the j-th column.
        if newisj < newiej
            indf[j] = newisj
            list[j] = list[L.rowval[newisj]]
            list[L.rowval[newisj]] = j
        end

        # Clear out elements j+1,...,n of the array indf.
        @inbounds for k=1:nlj
            indf[indr[k]] = 0
        end

        # Update isj and L.colptr.
        isj = L.colptr[j+1]
        L.colptr[j+1] = newiej + 1
    end

    return info
end

function dicf(n, nnz, L::TronDenseMatrix, p, indr, indf, list, w)
    zero = 0.0
    info = 0

    # We perform left-looking Cholesky factorization.
    @inbounds for j=1:n
        if L.vals[j,j] <= zero
            info = -j
            return info
        end
        L.vals[j,j] = sqrt(L.vals[j,j])

        # Update column j using the previous columns.
        @inbounds for k=1:j-1,i=j+1:n
            L.vals[i,j] = L.vals[i,j] - L.vals[i,k]*L.vals[j,k]
        end

        # Compute the j-th column of L.
        @inbounds for i=j+1:n
            L.vals[i,j] = L.vals[i,j] / L.vals[j,j]
            L.vals[i,i] = L.vals[i,i] - L.vals[i,j]^2
        end
    end

    return info
end

#=
# Right-looking Cholesky
function dicf(n::Int,L::CuDeviceArray{Float64})
    tx = threadIdx().x
    ty = threadIdx().y

    zero = 0.0
    Lmn = L[tx,ty]
    @inbounds for j=1:n
        # Update the diagonal.
        if (tx == j) && (ty == j)
            if Lmn > zero
                L[j,j] = sqrt(Lmn)
            else
                L[j,j] = -1.0
            end
        end
        CUDA.sync_threads()

        Ljj = L[j,j]
        if Ljj <= zero
            CUDA.sync_threads()
            return -1
        end

        # Update the jth column.
        if ty == j
            L[tx,j] = Lmn / Ljj
        end
        CUDA.sync_threads()

        # Update the trailing submatrix. To avoid if-conditional,
        # we update the whole matrix, but only the trailing part
        # will be saved in next iteration.
        Lmn = Lmn - L[tx,j] * L[ty,j]
    end

    if tx > ty
        L[ty,tx] = L[tx,ty]
    end
    CUDA.sync_threads()

    return 0
end
=#

# Left-looking Cholesky
@inline function dicf(n::Int,L::CuDeviceArray{Float64,2})
    tx = threadIdx().x
    ty = threadIdx().y

    @inbounds for j=1:n
        # Apply the pending updates.
        if j > 1
            if tx >= j && tx <= n && ty == 1
                for k=1:j-1
                    L[tx,j] -= L[tx,k] * L[j,k]
                end
            end
        end
        CUDA.sync_threads()

        if (L[j,j] <= 0)
            CUDA.sync_threads()
            return -1
        end

        Ljj = sqrt(L[j,j])
        if tx >= j && tx <= n && ty == 1
            L[tx,j] /= Ljj
        end
        CUDA.sync_threads()
    end

    if tx <= n && ty == 1
        @inbounds for j=1:n
            if tx > j
                L[j,tx] = L[tx,j]
            end
        end
    end
    CUDA.sync_threads()

    return 0
end


text:
module DFViews

using DataFrames, Gtk, Dates, Stella, Labels

export dfview, dfdesc

function eltype3(df::DataFrame,v::Symbol)
    if typeof(df[v]) <: CategoricalArray
        return eltype(df[v].pool.index)
    end
    return Missings.T(eltype(df[v]))
end

function colmissing(df)
    return [size(df,1) - sum(completecases(df[[s]])) for s in names(df)]
end

"""
    dfview(df::DataFrame; start=1,maxrows=1000)

Displays the content of a DataFrame in a separate window. This is a very basic
data viewer that does not have any other functionality. The default maximum
number of rows is 1000.

## Options
- start: starting row (default: 1)
- maxrows: number of rows to display (default: 1000)

"""
function dfview(df::DataFrame; start=1, maxrows=1000)

    # df - dataframe name
    (nrow,ncol) = size(df)

    # variable names
    varnames = names(df)

    # data types
    etypes = Vector{DataType}()
    for v in names(df)
        push!(etypes,String)
    end

    # define a listStore
    ls = GtkListStore(Int64,etypes...)

    # add data to the listStore
    #row_array = Array{Any,1}(ncol)
    aa = Array{Any}(undef,ncol)

    # number of rows to output
    finish = start + (nrow > maxrows ? maxrows : nrow) - 1
    for i in start:finish
        for j in 1:ncol
            aa[j] = ismissing(df[i,j]) ? "" : string(df[i,j])
        end
        push!(ls,tuple(i, aa...))
    end

    # TreeView
    tv = GtkTreeView(GtkTreeModel(ls))

    # let's put borders around the cells
    set_gtk_property!(tv,:enable_grid_lines,3)

    # variable names on the column headings
    r = GtkCellRendererText()

    # setproperty!(r,:xalign, 1.) # left justification
    c = GtkTreeViewColumn("Obs", r, Dict([("text",0)]))
    push!(tv,c)
    for (i,v) in enumerate(varnames)
        if eltype3(df,v) <: Number
            set_gtk_property!(r,:xalign,1)
        else
            set_gtk_property!(r,:xalign,0)
        end
        push!(tv,GtkTreeViewColumn(string(v), r, Dict([("text",i)])))
    end

	# add a Frame with scollbars
	sw = GtkScrolledWindow(tv)

    # create a frame and add it to the window
    nrows = start + 999
    win = GtkWindow(sw,"DataFrame - $nrow rows, $ncol columns - Rows $start - $nrows" , 600, 800)

    Gtk.showall(win)
end

"""
    dfdesc(df::DataFrame; labels::Union{Nothing,Label} = nothing)

Displays a list of all variables of a DataFrame in a separate windows. If `labels`
is specified, label names and variable labels are displayed.
"""
function dfdesc(df::DataFrame,labels::Union{Nothing,Label} = nothing)

    # df - dataframe name
    (nrow,ncol) = size(df)

    # variable names
    varnames = names(df)

    # define a listStore
    if labels == nothing
        # Row, Variable Name, Atype, Eltype, N Missing, P Missing
        ls = GtkListStore(Int,String,String,String,Int,String)
    else
        # Row, Variable Name, Atype, Eltype, N Missing, P Missing, Lblname, Description
        ls = GtkListStore(Int,String,String,String,Int,String,String,String)
    end

    # add data to the listStore
    aa = Array{Any}(undef,labels == nothing ? 6 : 8)

    # number of missing values in each column - DataFrames.colmissing is not working
    nmiss = colmissing(df)

    for i in 1:ncol
        # row
        aa[1] = i

        # variable name
        aa[2] = string(varnames[i])

        # array type
        aa[3] = Stella.atype(df,varnames[i])

        # eltype
        aa[4] = Stella.etype(df,varnames[i])

        # number of rows with missing values
        aa[5] = nmiss[i]

        # percent missing
        aa[6] = string( round(100 * nmiss[i] / size(df,1), digits = 2), "%")

        if labels != nothing

            # label name
            lname = Labels.lblname(labels,varnames[i])
            aa[7] = lname == nothing ? "" : string(lname)

            # variable description
            aa[8] = Labels.varlab(labels,varnames[i])
        end

        push!(ls,tuple(aa...))
    end

    # TreeView
    tv = GtkTreeView(GtkTreeModel(ls))

    # let's put borders around the cells
    set_gtk_property!(tv,:enable_grid_lines,3)

    # variable names on the column headings
    r = GtkCellRendererText()

    # setproperty!(r,:xalign, 1.) # left justification
    set_gtk_property!(r,:xalign,1)
    push!(tv,GtkTreeViewColumn("Num", r, Dict([("text",0)])))
    for (i,v) in enumerate(["Variable Name","Array Type","Eltype","N Miss","% Miss","Lblname","Description"])
        if labels == nothing && i > 5
            continue
        end
        r = GtkCellRendererText()
        if i in [4,5]
            set_gtk_property!(r,:xalign,1)
        else
            set_gtk_property!(r,:xalign,0)
        end
        push!(tv,GtkTreeViewColumn(v, r, Dict([("text",i)])))
    end

	# add a Frame with scollbars
	sw = GtkScrolledWindow(tv)

    # create a frame and add it to the window
    w = labels == nothing ? 450 : 750
    h = ncol > 25 ? 800 : 30 * ncol
    win = GtkWindow(sw,"DataFrame Table of Contents" , w, h)
    # push!(win,sw)

    Gtk.showall(win)
end

end # module


text:
# REST API adapters
# the adapters returns the raw JSON in Dict{String, Any}
export AuthAPI, ServiceAPI, ProjectAPI, JobAPI

abstract type AbstractAPI end

endpoint(api::AbstractAPI) = api.endpoint
api_layers(::AbstractAPI; kw...) = HTTP.stack(;kw...)
function api_uri(api::AbstractAPI, path::String, query=nothing)
    uri = endpoint(api)::URI
    # NOTE: this directly set the query to query
    if query === nothing
        return joinpath(uri, path)
    else
        isempty(uri.query) || error("non empty query in API endpoint is not supported")
        return URI(joinpath(uri, path); query=query)
    end
end

function to_json(response)
    body = HTTP.payload(response, String)
    return JSON.parse(body)
end

"abstract REST API type for IBM Q"
abstract type IBMQAPI <: AbstractAPI end

"""
    create_headers(api::IBMQAPI[, headers=HTTP.Header[]; kw...])

Create headers for a IBMQ REST API object.
"""
function create_headers(::IBMQAPI, headers=HTTP.Header[]; kw...)
    headers = HTTP.mkheaders(headers)
    # insert headers
    # TODO: use version
    HTTP.setkv(headers, "X-Qx-Client-Application", "IBMQClient.jl")
    HTTP.setkv(headers, "Connection", "keep-alive")
    HTTP.setkv(headers, "Accept", "*/*")
    access_token = Base.get(kw, :access_token, nothing)
    if access_token === nothing
        HTTP.rmkv(headers, "X-Access-Token")
    else
        HTTP.setkv(headers, "X-Access-Token", access_token)
    end
    return headers
end

function request(api::IBMQAPI, method::String, path::String, body=HTTP.nobody; headers=HTTP.Header[], query=nothing, kw...)
    headers = create_headers(api, headers; kw...)

    if body === nothing
        body = HTTP.nobody
    elseif body isa Dict
        # we convert Julia Dict to JSON in IBMQ REST API
        HTTP.setkv(headers, "Content-Type", "application/json")
        if body isa Dict
            body = JSON.json(body)
        end
    elseif body isa Vector{UInt8}
    elseif body isa AbstractString
    else
        error("invalid content type: $(typeof(body))")
    end

    uri = api_uri(api, path, query)
    layers = api_layers(api; kw...)
    response = HTTP.request(layers, method, uri, HTTP.mkheaders(headers), body; kw...)
    return to_json(response)
end

function get(api::AbstractAPI, path::String, body=HTTP.nobody; kw...)
    request(api, "GET", path, body; kw...)
end

post(api::AbstractAPI, path::String, body=HTTP.nobody; kw...) = request(api, "POST", path, body; kw...)
put(api::AbstractAPI, path::String, body=HTTP.nobody; kw...) = request(api, "PUT", path, body; kw...)
patch(api::AbstractAPI, path::String, body=HTTP.nobody; kw...) = request(api, "PATCH", path, body; kw...)
head(api::AbstractAPI, path::String, body=HTTP.nobody; kw...) = request(api, "HEAD", path, body; kw...)

"""
    AuthAPI <: IBMQAPI

IBM Q Authentication REST API.
"""
Base.@kwdef struct AuthAPI <: IBMQAPI
    endpoint::URI = URI("https://auth.quantum-computing.ibm.com/api")
end

"""
    login(::AuthAPI, token::String)

Login with user token.
"""
function login(api::AuthAPI, token::String)
    return post(api, "users/loginWithToken", Dict("apiToken" => token))
end

"""
    user_info(::AuthAPI, access_token::String)

Get user info of given IBM Q account access token.
"""
function user_info(api::AuthAPI, access_token::String)
    return get(api, "users/me"; access_token=access_token)
end

"""
    user_urls(::AuthAPI, access_token::String)

Get user urls to create services etc.
"""
function user_urls(api::AuthAPI, access_token::String)
    return user_info(api, access_token)["urls"]
end

"""
    ServiceAPI <: IBMQAPI

IBM Q service REST API.
"""
struct ServiceAPI <: IBMQAPI
    endpoint::URI
end

"""
    ServiceAPI([uri="https://api.quantum-computing.ibm.com/api"])

Create IBM Q service REST API object.
"""
ServiceAPI() = ServiceAPI("https://api.quantum-computing.ibm.com/api")
ServiceAPI(uri::String) = ServiceAPI(URI(uri))

"""
    ServiceAPI(::AuthAPI, access_token::String)

Create IBM Q service REST API by querying authentication server.
"""
function ServiceAPI(auth::AuthAPI, access_token::String)
    return ServiceAPI(URI(user_urls(auth, access_token)["http"]))
end

"""
    hubs(api::ServiceAPI, access_token::String)

Get alll IBM hubs.
"""
function hubs(api::ServiceAPI, access_token::String)
    return get(api, "Network"; access_token=access_token)
end

"""
    user_hubs(api::ServiceAPI, access_token::String)

Get given users' hubs.
"""
function user_hubs(api::ServiceAPI, access_token::String)
    ret_hubs = Dict{String, String}[]
    for hub in hubs(api, access_token)
        hub_name = hub["name"]
        for (group_name, group) in hub["groups"], (project_name, project) in group["projects"]
            entry = Dict(
                "hub" => hub_name,
                "group" => group_name,
                "project" => project_name
            )

            if Base.get(project, "isDefault", false)
                insert!(ret_hubs, 1, entry)
            else
                push!(ret_hubs, entry)
            end
        end
    end
    return ret_hubs
end

# this is similar to Account class in qiskit
"""
    ProjectAPI <: IBMQAPI

IBM Q Project REST API.
"""
struct ProjectAPI <: IBMQAPI
    endpoint::URI
end

"""
    ProjectAPI(base, hub::String, group::String, project::String)

Create IBM Q Project REST API from given base uri, `hub`, `group` and project name `project`.
"""
function ProjectAPI(base::URI, hub::String, group::String, project::String)
    return ProjectAPI(base.uri, hub, group, project)
end

function ProjectAPI(base::String, hub::String, group::String, project::String)
    return ProjectAPI(URI(base * ibmq_template_hubs(hub, group, project)))
end

function ibmq_template_hubs(hub::String, group::String, project::String)
    return "/Network/$hub/Groups/$group/Projects/$project"
end

"""
    devices(api::ProjectAPI, access_token::String; timeout = 0)

Query available devices.
"""
function devices(api::ProjectAPI, access_token::String; timeout = 0)
    return get(api, "devices/v/1"; readtimeout = timeout, access_token = access_token)
end

"""
    jobs(api::ProjectAPI, access_token::String; descending::Bool=true, limit::Int=10, skip::Int=0, extra_filter=nothing)

Query available jobs.

## Args
- `limit`: Maximum number of items to return.
- `skip`: Offset for the items to return.
- `descending`: Whether the jobs should be in descending order.
- `extra_filter`: Additional filtering passed to the query.
"""
function jobs(api::ProjectAPI, access_token::String; descending::Bool=true, limit::Int=10, skip::Int=0, extra_filter=nothing)
    order = descending ? "DESC" : "ASC"
    query = Dict{String, Any}(
        "order" => "creationDate " * order,
        "limit" => limit,
        "skip" => skip,
    )

    if !isnothing(extra_filter)
        query["where"] = extra_filter
    end

    get(api, "Jobs/status/v/1"; access_token=access_token, query=Dict("filter" => JSON.json(query)))
end

"""
    create_remote_job(api::ProjectAPI, device::IBMQDevice, access_token::String; job_name=nothing, job_share_level=nothing, job_tags=nothing)

Create a job instance on the remote server.

## Args

- `device_name`: A IBM Q device name.
- `job_name`: Custom name to be assigned to the job.
- `job_share_level`: Level the job should be shared at.
- `job_tags`: Tags to be assigned to the job.
"""
function create_remote_job(api::ProjectAPI, device_name::String, access_token::String; job_name=nothing, job_share_level=nothing, job_tags=nothing)
    payload = Dict{String, Any}(
        "backend" => Dict{String, Any}(
            "name" => device_name,
        ),
        "allowObjectStorage" => true,
    )

    if !isnothing(job_name)
        payload["name"] = job_name
    end

    if !isnothing(job_share_level)
        payload["shareLevel"] = job_share_level
    end

    if !isnothing(job_tags)
        payload["tags"] = job_tags
    end

    return post(api, "Jobs", payload; access_token=access_token)
end

struct JobAPI <: IBMQAPI
    endpoint::URI
end

function JobAPI(api::ProjectAPI, job_id::String)
    JobAPI(joinpath(api.endpoint, "Jobs/$job_id"))
end

function upload_url(api::JobAPI, access_token::String)
    get(api, "jobUploadUrl"; access_token=access_token)
end

function retreive_job_info(api::JobAPI, access_token::String)
    get(api, "v/1"; access_token=access_token)
end

function cancel(api::JobAPI, access_token::String)
    post(api, "cancel"; access_token=access_token)
end

function properties(api::JobAPI, access_token::String)
    get(api, "properties"; access_token=access_token)
end

function result_url(api::JobAPI, access_token::String)
    get(api, "resultDownloadUrl"; access_token=access_token)
end

function download_url(api::JobAPI, access_token::String)
    get(api, "jobDownloadUrl"; access_token=access_token)
end

function status(api::JobAPI, access_token::String)
    get(api, "status/v/1"; access_token=access_token)
end

function put_object_storage(api::JobAPI, url::String, qobj::Schema.Qobj, access_token::String)
    put_object_storage(api, url, to_dict(qobj; include_defaults=true, exclude_nothing=true), access_token)
end

function put_object_storage(api::JobAPI, url::String, qobj::AbstractDict{String, Any}, access_token::String)
    put_object_storage(api, url, JSON.json(qobj), access_token)
end

function put_object_storage(api::JobAPI, url::String, qobj::String, access_token::String)
    HTTP.put(url, create_headers(api; access_token=access_token), qobj; readtimeout=600)
end

function get_object_storage(api::JobAPI, url::String, access_token::String)
    response = HTTP.get(url, create_headers(api; access_token=access_token); readtimeout=600)
    return to_json(response)
end

function callback_upload(api::JobAPI, access_token::String)
    post(api, "jobDataUploaded"; access_token=access_token)
end

function callback_download(api::JobAPI, access_token::String)
    post(api, "resultDownloaded"; access_token=access_token)
end


text:
using OrdinaryDiffEq, ODE, ODEInterfaceDiffEq, DiffEqDevTools
tspan = (0.0,1.0)
# Linear ODE
f = (u,p,t) -> (1.01*u)
(::typeof(f))(::Type{Val{:analytic}},u₀,p,t) = u₀*exp(1.01*t)
probnum = ODEProblem(f,1/2,tspan)

# 2D Linear ODE
f = (du,u,p,t) -> begin
  @inbounds for i in eachindex(u)
    du[i] = 1.01*u[i]
  end
end
(::typeof(f))(::Type{Val{:analytic}},u₀,p,t) = u₀*exp(1.01*t)
prob = ODEProblem(f,rand(100,100),tspan)
using Plots; gr()

setups = [Dict(:alg=>RK4());Dict(:alg=>RK4(),:dense=>false);Dict(:alg=>RK4(),:dense=>false,:save_everystep=>false);Dict(:alg=>ode4())]
names = ["RK4 Continuous";"RK4 Non-Dense";"RK4 Save End";"ODE.jl rk4"]
shoot = Shootout(probnum,setups;dt=1/2^(6),names=names);

@show shoot.times # Times
@show shoot.errors # Errors
@show shoot.effs # Efficiencies
@show shoot.effratios[shoot.bestidx,:] # Efficiencies

plot(shoot)

setups = [Dict(:alg=>RK4());Dict(:alg=>RK4(),:dense=>false);Dict(:alg=>RK4(),:dense=>false,:save_everystep=>false);Dict(:alg=>ode4())]
names = ["RK4 Continuous";"RK4 Non-Dense";"RK4 Save End";"ode4"]
shoot = Shootout(prob,setups;dt=1/2^(6),names=names)
println(shoot.effratios[shoot.bestidx,:])
plot(shoot)

setups = [Dict(:alg=>DP5())
          Dict(:abstol=>1e-3,:reltol=>1e-6,:alg=>ode45()); # Fix ODE to be normal
          Dict(:alg=>dopri5())]
shoot = Shootout(probnum,setups)
println(shoot.times)
println(shoot.errors)
println(shoot.effratios[shoot.bestidx,:])
plot(shoot)

using DiffEqBenchmarks
DiffEqBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])


text:
using AutoHashEquals
using TermInterface

"""
Abstract type representing a pattern used in all the various pattern matching backends.
You can use the `Pattern` constructor to recursively convert an `Expr` (or any type satisfying [`Metatheory.TermInterface`](@ref)) to a [`Pattern`](@ref).
"""
abstract type Pattern end

import Base.==
==(a::Pattern, b::Pattern) = false
TermInterface.arity(p::Pattern) = 0
function isground(p::Pattern)
    false
end

"""
Pattern variables will first match on any subterm
and instantiate the substitution to that subterm.
"""
mutable struct PatVar <: Pattern
    name::Symbol
    idx::Int
end
function ==(a::PatVar, b::PatVar)
    # (a.name == b.name)
    a.idx == b.idx
end
PatVar(var) = PatVar(var, -1)

"""
A pattern literal will match only against an instance of itself.
Example:
```julia
PatLiteral(2)
```
Will match only against values that are equal (using `Base.(==)`) to 2.

```julia
PatLiteral(:a)
```
Will match only against instances of the literal symbol `:a`.
"""
@auto_hash_equals struct PatLiteral{T} <: Pattern
    val::T
end
function isground(p::PatLiteral)
    true
end

"""
Type assertions on a [`PatVar`](@ref), will match if and only if
the type of the matched term for the pattern variable `var` is a subtype
of `type`.
Type assertions are supported in the left hand of rules
to match and access literal values both when using classic
rewriting and EGraph based rewriting.
To use a type assertion pattern, add `::T` after
a pattern variable in the `left_hand` of a rule.
"""
struct PatTypeAssertion <: Pattern
    var::PatVar
    type::Type
    hash::Ref{UInt}
    PatTypeAssertion(v,t) = new(v, t, Ref{UInt}(0))
end
function Base.hash(t::PatTypeAssertion, salt::UInt)
    !iszero(salt) && return hash(hash(t, zero(UInt)), salt)
    h = t.hash[]
    !iszero(h) && return h
    h′ = hash(t.var,  hash(t.type, salt))
    t.hash[] = h′
    return h′
end

struct PatSplatVar <: Pattern
    var::PatVar
    hash::Ref{UInt}
    PatSplatVar(v) = new(v, Ref{UInt}(0))
end
function Base.hash(t::PatSplatVar, salt::UInt)
    !iszero(salt) && return hash(hash(t, zero(UInt)), salt)
    h = t.hash[]
    !iszero(h) && return h
    h′ = hash(t.var, salt)
    t.hash[] = h′
    return h′
end

"""
This type of pattern will match if and only if
the two subpatterns exist in the same equivalence class,
in the e-graph on which the matching is performed.
**Can be used only in the e-graphs backend**
"""
struct PatEquiv <: Pattern
    left::Pattern
    right::Pattern
    hash::Ref{UInt}
    PatEquiv(l,r) = new(l,r, Ref{UInt}(0))
end

function ==(a::PatEquiv, b::PatEquiv)
    a.left == b.left && a.right == b.right
end

function Base.hash(t::PatEquiv, salt::UInt)
    !iszero(salt) && return hash(hash(t, zero(UInt)), salt)
    h = t.hash[]
    !iszero(h) && return h
    h′ = hash(t.left,  hash(t.right, salt))
    t.hash[] = h′
    return h′
end

function isground(p::PatEquiv)
   isground(p.left) && isground(p.right)
end

"""
Term patterns will match
on terms of the same `arity` and with the same
function symbol (`head`).
"""
struct PatTerm <: Pattern
    head::Any
    args::Vector{Pattern}
    hash::Ref{UInt}
    PatTerm(h,args) = new(h,args, Ref{UInt}(0))
end
TermInterface.gethead(p::PatTerm) = p.head
TermInterface.arity(p::PatTerm) = length(p.args)

function Base.hash(t::PatTerm, salt::UInt)
    !iszero(salt) && return hash(hash(t, zero(UInt)), salt)
    h = t.hash[]
    !iszero(h) && return h
    h′ = hash(t.head,  hash(t.args, salt))
    t.hash[] = h′
    return h′
end

function isground(p::PatTerm)
    mapreduce(isground, (&), p.args)
end

"""
This pattern type matches on a function application
but instead of strictly matching on a head symbol,
it has a pattern variable as head. It can be used for
example to match arbitrary function calls.
"""
@auto_hash_equals struct PatAllTerm <: Pattern
    head::PatVar
    args::Vector{Pattern}
end
TermInterface.arity(p::PatAllTerm) = length(p.args)

# ==============================================
# ================== PATTERN VARIABLES =========
# ==============================================

"""
Collects pattern variables appearing in a pattern into a vector of symbols
"""
patvars(p::PatLiteral, s) = s
patvars(p::PatVar, s) = push!(s, p.name)
patvars(p::PatTypeAssertion, s) = patvars(p.var, s)
patvars(p::PatSplatVar, s) = patvars(p.var, s)
function patvars(p::PatEquiv, s)
    patvars(p.left, s)
    patvars(p.right, s)
end

function patvars(p::PatTerm, s)
    for x ∈ p.args
        patvars(x, s)
    end
    return s
end

function patvars(p::PatAllTerm, s)
    push!(s, p.head.name)
    for x ∈ p.args
        patvars(x, s)
    end
    return s
end

function patvars(p::Pattern)
    unique(patvars(p, Symbol[]))
end

# ==============================================
# ================== DEBRUJIN INDEXING =========
# ==============================================

Base.setindex!(p::PatLiteral, pvars) = nothing
function Base.setindex!(p::PatVar, pvars)
    p.idx = findfirst((==)(p.name), pvars)
end
Base.setindex!(p::PatTypeAssertion, pvars) = setindex!(p.var, pvars)
Base.setindex!(p::PatSplatVar, pvars) = setindex!(p.var, pvars)

function Base.setindex!(p::PatEquiv, pvars)
    setindex!(p.left, pvars)
    setindex!(p.right, pvars)
end

function Base.setindex!(p::PatTerm, pvars)
    for x ∈ p.args
        setindex!(x, pvars)
    end
end

function Base.setindex!(p::PatAllTerm, pvars)
    setindex!(p.head, pvars)
    for x ∈ p.args
        setindex!(x, pvars)
    end
end


text:
## Methods for root finding which use a bracket

## Bisection_method for floats and int

# type to throw on succesful convergence
type StateConverged
    x0::Real
end

# type to throw on failure
type ConvergenceFailed
    reason::AbstractString
end

## From Jason Merrill https://gist.github.com/jwmerrill/9012954
## cf. http://squishythinking.com/2014/02/22/bisecting-floats/
# Alternative "mean" definition that operates on the binary representation
# of a float. Using this definition, bisection will never take more than
# 64 steps.
function _middle(x::Float64, y::Float64)
  # Use the usual float rules for combining non-finite numbers
  if !isfinite(x) || !isfinite(y)
    return x + y
  end

  # Always return 0.0 when inputs have opposite sign
  if sign(x) != sign(y) && x != 0.0 && y != 0.0
    return 0.0
  end

  negate = x < 0.0 || y < 0.0

  xint = reinterpret(UInt64, abs(x))
  yint = reinterpret(UInt64, abs(y))
  unsigned = reinterpret(Float64, (xint + yint) >> 1)

  return negate ? -unsigned : unsigned
end

"""

Find zero using modified bisection method for Float64 arguments.

This is guaranteed to take no more than 64 steps. The `a42` alternative usually has
fewer iterations, but this seems to find the value with fewer function evaluations.

Terminates with `x1` when the bracket length of `[x0,x2]` is `<= max(xtol, xtolrel*abs(x1))` where `x1` is the midpoint . The tolerances can be set to 0, in which case, the termination occurs when `nextfloat(x0) = x2`.

The bracket `[a,b]` must be bounded.

"""
function find_zero(f, a::Float64, b::Float64; xtol::Real=0.0, xtolrel::Real=0.0, verbose::Bool=false)

    if (xtol < 0.0) | (xtolrel < 0.0)
        error("Tolerances must be non-negative")
    end

    x0, y0 = a, f(a)
    x2, y2 = b, f(b)

    y0 == 0 && return a
    y2 == 0.0 && return b
    if sign(y0) * sign(y2) > 0
        error("[a,b] is not a bracket")
    end

    x1 = _middle(x0, x2)
    y1 = f(x1)

    while x0 < x1 && x1 < x2
        if sign(y0) == sign(y1)
            x0, x2 = x1, x2
            y0, y2 = y1, y2
        else
            x0, x2 = x0, x1
            y0, y2 = y0, y1
        end

        x1 = _middle(x0, x2)
        y1 = f(x1)

        sign(y1) == 0 && return x1
        abs(x2 - x0) <= max(xtol, xtolrel*abs(x1)) && return(x1)

        verbose && println(@sprintf("xi =%18.15f,\t f(xi) = %18.15f", x1, f(x1)))
    end

    return abs(y0) < abs(y2) ? x0 : x2
end

"""

find_zero using Algorithm 4.2  of Alefeld, Potra and Shi for Big numbers

"""
typealias BigSomething @compat  Union{BigFloat, BigInt}
function find_zero(f, a::BigSomething, b::BigSomething;
                   xtol=zero(a),
                   maxeval::Int=100,
                   verbose::Bool=false
                   )
    a42(f, a, b; xtol=xtol, maxeval=maxeval, verbose=verbose)
end

"""

Finds the root of a continuous function within a provided
interval [a, b], without requiring derivatives. It is based on algorithm 4.2
described in: 1. G. E. Alefeld, F. A. Potra, and Y. Shi, "Algorithm 748:
enclosing zeros of continuous functions," ACM Trans. Math. Softw. 21,
327–344 (1995).

input:
    `f`: function to find the root of
    `a`, `b`: the initial bracket, with: a < b, f(a)*f(b) < 0
    `xtol`: acceptable error (it's safe to set zero for machine precision)
    `maxeval`:  maximum number of iterations

output:
    an estimate of the zero of f

By John Travers

"""
function a42(f, a, b;
                   xtol=zero(float(a)),
                   maxeval::Int=15,
                   verbose::Bool=false
                   )
    if a >= b || sign(f(a))*sign(f(b)) >= 0
        error("on input a < b and f(a)f(b) < 0 must both hold")
    end
    if xtol < 0.0
        error("tolerance must be >= 0.0")
    end

    c = secant(f, a, b)
    a42a(f, a, b, c,
         xtol=xtol, maxeval=maxeval, verbose=verbose)
end

"""

Split Alefeld, F. A. Potra, and Y. Shi algorithm 4.2 into a function
where `c` is passed in.

Solve f(x) = 0 over bracketing interval [a,b] starting at c, with a < c < b

"""
function a42a(f, a, b, c=(0.5)*(a+b);
                   xtol=zero(float(a)),
                   maxeval::Int=15,
                   verbose::Bool=false
              )

    try
        # re-bracket and check termination
        a, b, d = bracket(f, a, b, c, xtol)
        ee = d
        for n = 2:maxeval
            # use either a cubic (if possible) or quadratic interpolation
            if n > 2 && distinct(f, a, b, d, ee)
                c = ipzero(f, a, b, d, ee)
            else
                c = newton_quadratic(f, a, b, d, 2)
            end
            # re-bracket and check termination
            ab, bb, db = bracket(f, a, b, c, xtol)
            eb = d
            # use another cubic (if possible) or quadratic interpolation
            if distinct(f, ab, bb, db, eb)
                cb = ipzero(f, ab, bb, db, eb)
            else
                cb = newton_quadratic(f, ab, bb, db, 3)
            end
            # re-bracket and check termination
            ab, bb, db = bracket(f, ab, bb, cb, xtol)
            # double length secant step; if we fail, use bisection
            u = abs(f(ab)) < abs(f(bb)) ? ab : bb
            cb = u - 2*f(u)/(f(bb) - f(ab))*(bb - ab)
            ch = abs(cb - u) > (bb - ab)/2 ? ab + (bb - ab)/2 : cb
            # re-bracket and check termination
            ah, bh, dh = bracket(f, ab, bb, ch, xtol)
            # if not converging fast enough bracket again on a bisection
            if bh - ah < 0.5*(b - a)
                a = ah
                b = bh
                d = dh
                ee = db
            else
                ee = dh
                a, b, d = bracket(f, ah, bh, ah + (bh - ah)/2,
                                  xtol)
            end

            verbose && println(@sprintf("a=%18.15f, n=%s", a, n))

            if nextfloat(ch) * prevfloat(ch) <= 0
                throw(StateConverged(ch))
            end
            if nextfloat(a) >= b
                throw(StateConverged(a))
            end
        end
        throw(ConvergenceFailed("More than $maxeval iterations before convergence"))
    catch ex
        if isa(ex, StateConverged)
            return ex.x0
        else
            rethrow(ex)
        end
    end
    throw(ConvergenceFailed("More than $maxeval iterations before convergence"))
end

# calculate a scaled tolerance
# based on algorithm on page 340 of [1]
function tole(a, b, fa, fb, tol)
    u = abs(fa) < abs(fb) ? abs(a) : abs(b)
    2u*eps(1.0) + tol
end

# bracket the root
# inputs:
#     - f: the function
#     - a, b: the current bracket with a < b, f(a)f(b) < 0
#     - c within (a,b): current best guess of the root
#     - tol: desired accuracy
#
# if root is not yet found, return
#     ab, bb, d
# with:
#     - [ab, bb] a new interval within [a, b] with f(ab)f(bb) < 0
#     - d a point not inside [ab, bb]; if d < ab, then f(ab)f(d) > 0,
#       and f(d)f(bb) > 0 otherwise
#
# if the root is found, throws a StateConverged instance with x0 set to the
# root.
#
# based on algorithm on page 341 of [1]
function bracket(f, a, b, c, tol)

    fa = f(a)
    fb = f(b)

    if !(a <= c <= b)
        error("c must be in (a,b)")
    end
    delta = 0.7*tole(a, b, fa, fb, tol)
    if b - a <= 4delta
        c = (a + b)/2
    elseif c <= a + 2delta
        c = a + 2delta
    elseif c >= b - 2delta
        c = b - 2delta
    end
    fc = f(c)
    if fc == 0
        throw(StateConverged(c))
    elseif sign(fa)*sign(fc) < 0
        aa = a
        bb = c
        db = b
    else
        aa = c
        bb = b
        db = a
    end
    faa = f(aa)
    fbb = f(bb)
    if bb - aa < 2*tole(aa, bb, faa, fbb, tol)
        x0 = abs(faa) < abs(fbb) ? aa : bb
        throw(StateConverged(x0))
    end
    aa, bb, db
end

# take a secant step, if the resulting guess is very close to a or b, then
# use bisection instead
function secant(f, a, b)
    c = a - f(a)/(f(b) - f(a))*(b - a)
    tol = eps(1.0)*5
    if c <= a + abs(a)*tol || c >= b - abs(b)*tol
        return a + (b - a)/2
    end
    return c
end

# approximate zero of f using quadratic interpolation
# if the new guess is outside [a, b] we use a secant step instead
# based on algorithm on page 330 of [1]
function newton_quadratic(f, a, b, d, k::Int)
    fa = f(a)
    fb = f(b)
    fd = f(d)
    B = (fb - fa)/(b - a)
    A = ((fd - fb)/(d - b) - B)/(d - a)
    if A == 0
        return secant(f, a, b)
    end
    r = A*fa > 0 ? a : b
    for i = 1:k
        r -= (fa + (B + A*(r - b))*(r - a))/(B + A*(2*r - a - b))
    end
    if isnan(r) || (r <= a || r >= b)
        r = secant(f, a, b)
    end
    return r
end

# approximate zero of f using inverse cubic interpolation
# if the new guess is outside [a, b] we use a quadratic step instead
# based on algorithm on page 333 of [1]
function ipzero(f, a, b, c, d)
    fa = f(a)
    fb = f(b)
    fc = f(c)
    fd = f(d)
    Q11 = (c - d)*fc/(fd - fc)
    Q21 = (b - c)*fb/(fc - fb)
    Q31 = (a - b)*fa/(fb - fa)
    D21 = (b - c)*fc/(fc - fb)
    D31 = (a - b)*fb/(fb - fa)
    Q22 = (D21 - Q11)*fb/(fd - fb)
    Q32 = (D31 - Q21)*fa/(fc - fa)
    D32 = (D31 - Q21)*fc/(fc - fa)
    Q33 = (D32 - Q22)*fa/(fd - fa)
    c = a + (Q31 + Q32 + Q33)
    if (c <= a) || (c >= b)
        c = newton_quadratic(f, a, b, d, 3)
    end
    return c
end

# floating point comparison function
function almost_equal(x, y)
    const min_diff = realmin(Float64)*32
    abs(x - y) < min_diff
end

# check that all interpolation values are distinct
function distinct(f, a, b, d, e)
    f1 = f(a)
    f2 = f(b)
    f3 = f(d)
    f4 = f(e)
    !(almost_equal(f1, f2) || almost_equal(f1, f3) || almost_equal(f1, f4) ||
      almost_equal(f2, f3) || almost_equal(f2, f4) || almost_equal(f3, f4))
end

"""
Searches for zeros  of `f` in an interval [a, b].

Basic algorithm used:

* split interval [a,b] into `no_pts` subintervals.
* For each bracketing interval finds a bracketed zero.
* For other subintervals does a quick search with a derivative-free method.

If there are many zeros relative to the number of points, the process
is repeated with more points, in hopes of finding more zeros for
oscillating functions.

"""
function find_zeros(f, a::Real, b::Real, args...;
                    no_pts::Int=100,
                    ftol::Real=10*eps(), reltol::Real=10*eps(),
                    kwargs...)

    a, b = a < b ? (a,b) : (b,a)
    rts = eltype(promote(float(a),b))[]
    xs = vcat(a, a + (b-a)*sort(rand(no_pts)), b)

    ## Look in [ai, bi)
    for i in 1:(no_pts+1)
        ai,bi=xs[i:i+1]
        if isapprox(f(ai), 0.0, rtol=reltol, atol=ftol)
            push!(rts, ai)
        elseif sign(f(ai)) * sign(f(bi)) < 0
            push!(rts, fzero(f, ai, bi))
        else
            try
                x = fzero(f, ai + (0.5)* (bi-ai), order=8, maxeval=10, ftol=ftol, reltol=reltol)
                if ai < x < bi
                    push!(rts, x)
                end
            catch e
            end
        end
    end
    ## finally, b?
    isapprox(f(b), 0.0, rtol=reltol, atol=ftol) && push!(rts, b)

    ## redo if it appears function oscillates alot in this interval...
    if length(rts) > (1/4) * no_pts
        return(find_zeros(f, a, b, args...; no_pts = 10*no_pts, ftol=ftol, reltol=reltol, kwargs...))
    else
        return(sort(rts))
    end
end


text:
module Images

using PyCall: PyObject

using DockerPy: Collection, docker, @pyinterface

export Image, history, reload, save, tag, build, push, pull, search

mutable struct Image
    o::PyObject
end
Image() = Image(docker.models.Image())

history(x::Image) = PyObject(x).history()
reload(x::Image) = PyObject(x).reload()
save(x::Image; chunk_size = 2097152, named = false) =
    PyObject(x).save(chunk_size = chunk_size, named = named)
tag(x::Image, repository; tag = nothing, kwargs...) =
    PyObject(x).tag(repository, tag = tag, kwargs...)
Base.show(io::IO, x::Image) = print(
    io,
    "Image(tag = \"$(PyObject(x).tags[1])\", short_id = \"$(PyObject(x).short_id)\")",
)

const ImageCollection = Collection{Image}

function build(
    x::ImageCollection;
    path::AbstractString,
    fileobj,
    tag::AbstractString,
    quiet::Bool,
    nocache::Bool,
    rm::Bool = false,
    timeout::Integer,
    custom_context::Bool,
    encoding::AbstractString,
    pull::Bool,
    forcerm::Bool,
    dockerfile::AbstractString,
    buildargs,
    container_limits,
    shmsize::Integer,
    labels,
    cache_from,
    target::AbstractString,
    network_mode::AbstractString,
    squash::Bool,
    extra_hosts,
    platform::AbstractString,
    isolation::AbstractString,
    use_config_proxy::Bool,
)
    result = PyObject(x).build(
        path = path,
        fileobj = fileobj,
        tag = tag,
        quiet = quiet,
        nocache = nocache,
        rm = rm,
        timeout = timeout,
        custom_context = custom_context,
        encoding = encoding,
        pull = pull,
        forcerm = forcerm,
        dockerfile = dockerfile,
        buildargs = buildargs,
        container_limits = container_limits,
        shmsize = shmsize,
        labels = labels,
        cache_from = cache_from,
        target = target,
        network_mode = network_mode,
        squash = squash,
        extra_hosts = extra_hosts,
        platform = platform,
        isolation = isolation,
        use_config_proxy = use_config_proxy,
    )
    return Image(result[1]), string(collect(result[2]))
end # function build
pull(x::ImageCollection, repository; tag = nothing, kwargs...) =
    Image.(PyObject(x).pull(repository, tag = tag, kwargs...))
push(x::ImageCollection, repository; tag = nothing, kwargs...) =
    PyObject(x).pull(repository, tag = tag, kwargs...)
search(x::ImageCollection, term::AbstractString) = PyObject(x).search(term)
Base.rm(x::ImageCollection, image; force::Bool = false, noprune::Bool = false) =
    PyObject(x).remove(image, force = force, noprune = noprune)
Base.collect(x::ImageCollection) = map(Image, PyObject(x).list(all = true))
Base.get(x::ImageCollection, name::AbstractString) = Image(PyObject(x).get(name))
Base.empty!(x::ImageCollection, filters = nothing) = PyObject(x).prune(filters = filters)

@pyinterface Image

end # module Images


text:
#Declaring the abstract type
abstract type PokeType end
#Declaring the sub-types
struct Normal <: PokeType end
struct Fire <: PokeType end
struct Water <: PokeType end
struct Electric <: PokeType end
struct Grass <: PokeType end
struct Ice <: PokeType end
struct Fighting <: PokeType end
struct Poison <: PokeType end
struct Ground <: PokeType end
#Declaing constants that will not change
const NO_EFFECT = 0.0
const NOT_VERY_EFFECTIVE = 0.5
const NORMAL_EFFECTIVE = 1.0
const SUPER_EFFECTIVE = 2.0

#Defining some catch-all interactions
eff(atk::T1, def::T2) where {T1 <: PokeType, T2 <: PokeType} = NORMAL_EFFECTIVE
eff(atk::T, def::T) where {T <: PokeType} = NOT_VERY_EFFECTIVE

#Defining specific interactions between types
eff(atk:: Fighting, def:: Fighting) = NORMAL_EFFECTIVE
eff(atk:: Ground, def:: Ground) = NORMAL_EFFECTIVE
eff(atk::Fire,     def::Grass)    = SUPER_EFFECTIVE
eff(atk::Fire,     def::Ice)      = SUPER_EFFECTIVE
eff(atk::Water,    def::Fire)     = SUPER_EFFECTIVE
eff(atk::Water,    def::Ground)   = SUPER_EFFECTIVE
eff(atk::Electric, def::Water)    = SUPER_EFFECTIVE
eff(atk::Grass,    def::Water)    = SUPER_EFFECTIVE
eff(atk::Grass,    def::Ground)   = SUPER_EFFECTIVE
eff(atk::Ice,      def::Grass)    = SUPER_EFFECTIVE
eff(atk::Ice,      def::Ground)   = SUPER_EFFECTIVE
eff(atk::Fighting, def::Normal)   = SUPER_EFFECTIVE
eff(atk::Fighting, def::Ice)      = SUPER_EFFECTIVE
eff(atk::Poison,   def::Grass)    = SUPER_EFFECTIVE
eff(atk::Ground,   def::Fire)     = SUPER_EFFECTIVE
eff(atk::Ground,   def::Electric) = SUPER_EFFECTIVE
eff(atk::Ground,   def::Poison)   = SUPER_EFFECTIVE
eff(atk::Fire,     def::Water)  = NOT_VERY_EFFECTIVE
eff(atk::Water,    def::Grass)  = NOT_VERY_EFFECTIVE
eff(atk::Electric, def::Grass)  = NOT_VERY_EFFECTIVE
eff(atk::Grass,    def::Fire)   = NOT_VERY_EFFECTIVE
eff(atk::Grass,    def::Poison) = NOT_VERY_EFFECTIVE
eff(atk::Ice,      def::Fire)  = NOT_VERY_EFFECTIVE
eff(atk::Ice,      def::Water)  = NOT_VERY_EFFECTIVE
eff(atk::Fighting, def::Poison) = NOT_VERY_EFFECTIVE
eff(atk::Poison,   def::Ground) = NOT_VERY_EFFECTIVE
eff(atk::Ground,   def::Grass)  = NOT_VERY_EFFECTIVE
eff(atk::Electric, def::Ground) = NO_EFFECT

#Reads value of effectiveness and returns string
function eff_string(effectiveness)
    if effectiveness == 0.0
        return "not effective."
    elseif effectiveness == 0.5
        return "not very effective."
    elseif effectiveness == 1.0
        return "normally effective."
    elseif effectiveness == 2.0
        return "super effective!"
    end
end

#Function for when one pokemon attacks another
#Reads two variables of Type <: PokeType
#Checks the effectiveness of attacker and defender combo
#Prints the result of interaction
function attack(atk, def)
    effectiveness = eff(atk, def)
    println("A Pokemon used a $(typeof(atk)) attack")
    println("against a $(typeof(def)) Pokemon")
    println("it was $(eff_string(effectiveness))")
end

#Defining the Type of the attacker and defender
attacker = Fire()
defender = Fire()

#Checks outcome using attack()
attack(attacker, defender)

max_stars_repo_path: learning exercises/Poketypes.jl

max_stars_repo_name: ian-double-u/NNfromScratch

max_issues_repo_path: learning exercises/Poketypes.jl

max_issues_repo_name: ian-double-u/NNfromScratch

max_forks_repo_path: learning exercises/Poketypes.jl

max_forks_repo_name: ian-double-u/NNfromScratch


text:
using AWSSDK.Batch: describe_job_queues

struct JobQueue
    arn::String

    function JobQueue(queue::AbstractString)
        arn = job_queue_arn(queue)
        arn === nothing && error("No job queue ARN found for: $queue")
        new(arn)
    end
end

Base.:(==)(a::JobQueue, b::JobQueue) = a.arn == b.arn
describe(queue::JobQueue) = describe_job_queue(queue)
describe_job_queue(queue::JobQueue) = describe_job_queue(queue.arn)
max_vcpus(queue::JobQueue) = sum(max_vcpus(ce) for ce in compute_environments(queue))

function compute_environments(queue::JobQueue)
    ce_order = describe(queue)["computeEnvironmentOrder"]

    compute_envs = Vector{ComputeEnvironment}(undef, length(ce_order))
    for ce in ce_order
        i, arn = ce["order"], ce["computeEnvironment"]
        compute_envs[i] = ComputeEnvironment(arn)
    end

    return compute_envs
end

function job_queue_arn(queue::AbstractString)
    startswith(queue, "arn:") && return queue
    json = describe_job_queue(queue)
    isempty(json) ? nothing : json["jobQueueArn"]
end

function describe_job_queue(queue::AbstractString)::OrderedDict
    json = @mock describe_job_queues(Dict("jobQueues" => [queue]))
    queues = json["jobQueues"]
    len = length(queues)::Int
    @assert len <= 1
    return len == 1 ? first(queues) : OrderedDict()
end


text:
module InteractBase

using WebIO, OrderedCollections, Observables, CSSUtil, Colors, JSExpr
import Observables: ObservablePair, AbstractObservable
import JSExpr: JSString
using Random
using Dates
using Base64: stringmime
using JSON
using Knockout
using Knockout: js_lambda
using Widgets
import Widgets:
    observe,
    AbstractWidget,
    div,
    Widget,
    widget,
    widgettype,
    @layout!,
    components,
    input,
    spinbox,
    textbox,
    textarea,
    autocomplete,
    datepicker,
    timepicker,
    colorpicker,
    checkbox,
    toggle,
    filepicker,
    opendialog,
    savedialog,
    slider,
    rangeslider,
    rangepicker,
    button,
    dropdown,
    radiobuttons,
    checkboxes,
    toggles,
    togglebuttons,
    tabs,
    entry,
    latex,
    alert,
    highlight,
    notifications,
    confirm,
    togglecontent,
    tabulator,
    accordion,
    mask,
    tooltip!,
    wdglabel,
    slap_design!,
    @manipulate,
    manipulatelayout,
    triggeredby,
    onchange

import Observables: throttle, _val

export observe, Widget, widget

export @manipulate

export filepicker, opendialog, savedialog, datepicker, timepicker, colorpicker, spinbox

export autocomplete, input, dropdown, checkbox, textbox, textarea, button, toggle, togglecontent

export slider, rangeslider, rangepicker

export radiobuttons, togglebuttons, tabs, checkboxes, toggles

export latex, alert, confirm, highlight, notifications, accordion, tabulator, mask

export onchange

export settheme!, resettheme!, gettheme, availablethemes, NativeHTML

export slap_design!

abstract type WidgetTheme<:Widgets.AbstractBackend; end
struct NativeHTML<:WidgetTheme; end

const font_awesome = joinpath(@__DIR__, "..", "assets", "all.js")
const prism_js = joinpath(@__DIR__, "..", "assets", "prism.js")
const prism_css = joinpath(@__DIR__, "..", "assets", "prism.css")
const highlight_css = joinpath(@__DIR__, "..", "assets", "highlight.css")
const nouislider_min_js = joinpath(@__DIR__, "..", "assets", "nouislider.min.js")
const nouislider_min_css = joinpath(@__DIR__, "..", "assets", "nouislider.min.css")
const style_css = joinpath(@__DIR__, "..", "assets", "style.css")

include("classes.jl")
include("backends.jl")
include("utils.jl")
include("input.jl")
include("slider.jl")
include("optioninput.jl")
include("layout.jl")
include("output.jl")
include("modifiers.jl")

function __init__()
    if Widgets.get_backend() === Widgets.DummyBackend()
        Widgets.set_backend!(NativeHTML())
    end
end

end # module


text:
module MixtureCollapsedGibbs
using Distributions
using Plots
using StatsBase

include("gibbs.jl")

abstract type ModelType end
abstract type Parametric <: ModelType end
abstract type NonParametric <: ModelType end

model_name(::Type{Parametric}) = "parametric"
model_name(::Type{NonParametric}) = "nonparametric"

function log_predictive_likelihood(::Type{Parametric}, data_id, cluster_id, state)
    """Predictive likelihood of the data at data_id is generated
    by cluster_id given the currenbt state.

    From Section 2.4 of
    http://www.cs.ubc.ca/~murphyk/Papers/bayesGauss.pdf
    """
    ss = state.suffstats[cluster_id]
    hp_mean = state.hp.λ1
    hp_var = state.hp.λ2
    param_var = state.v
    x = state.data[data_id]
    return _helper_log_predictive_likelihood(ss, hp_mean, hp_var, param_var, x)
end

function _helper_log_predictive_likelihood(ss, hp_mean, hp_var, param_var, x)
    posterior_sigma2 = 1 / (ss.count * 1.0 / param_var + 1.0 / hp_var)
    predictive_mu =
        posterior_sigma2 * (hp_mean * 1.0 / hp_var + ss.count * ss.μ * 1.0 / param_var)
    predictive_sigma2 = param_var + posterior_sigma2
    predictive_sd = sqrt(predictive_sigma2)
    return logpdf(Normal(predictive_mu, predictive_sd), x)
end

function log_cluster_assign_score(::Type{Parametric}, cluster_id, state)
    """Log-likelihood that a new point generated will
    be assigned to cluster_id given the current state.
    """
    current_cluster_size = state.suffstats[cluster_id].count
    num_clusters = length(state.cluster_ids)
    α = state.hp.α
    return log(current_cluster_size + α * 1.0 / num_clusters)
end

function cluster_assignment_distribution(::Type{Parametric}, data_id, state)
    """Compute the marginal distribution of cluster assignment
    for each cluster.
    """
    labels = []
    scores = []
    for (cid, _) in state.suffstats
        score = log_predictive_likelihood(Parametric, data_id, cid, state)
        score += log_cluster_assign_score(Parametric, cid, state)
        score = exp(score)
        push!(labels, cid)
        push!(scores, score)
    end
    scores /= sum(scores)
    return labels, scores
end

function add_datapoint_to_suffstats(x, ss)
    """Add datapoint to sufficient stats for normal component
    """
    nμ = (ss.μ * (ss.count) + x) / (ss.count + 1)
    ncount = ss.count + 1
    return SS(nμ, ncount)
end

function remove_datapoint_from_suffstats(x, ss)
    """Remove datapoint from sufficient stats for normal component
    """
    nμ = (ss.μ * (ss.count) - x * 1.0) / (ss.count - 1)
    ncount = ss.count - 1
    return SS(nμ, ncount)
end

function gibb_step!(::Type{Parametric}, state::GibbsState)
    both = zip(state.data, state.assignment)
    for (data_id, (datapoint, cid)) in enumerate(both)
        state.suffstats[cid] = remove_datapoint_from_suffstats(
            datapoint, state.suffstats[cid]
        )
        labels, scores = cluster_assignment_distribution(Parametric, data_id, state)
        cid = labels[rand(Categorical(scores))]
        state.assignment[data_id] = cid
        state.suffstats[cid] = add_datapoint_to_suffstats(
            state.data[data_id], state.suffstats[cid]
        )
    end

    return nothing
end

function log_cluster_assign_score(::Type{NonParametric}, cluster_id, state)
    """Log-likelihood that a new point generated will
    be assigned to cluster_id given the current state.
    """
    if cluster_id == typemax(Int)
        return log(state.hp.α)
    else
        return log(state.suffstats[cluster_id].count)
    end
end

function log_predictive_likelihood(::Type{NonParametric}, data_id, cluster_id, state)
    """Predictive likelihood of the data at data_id is generated
    by cluster_id given the currenbt state.

    From Section 2.4 of
    http://www.cs.ubc.ca/~murphyk/Papers/bayesGauss.pdf
    """
    if cluster_id == typemax(Int)
        ss = SS(0, 0)
    else
        ss = state.suffstats[cluster_id]
    end

    hp_mean = state.hp.λ1
    hp_var = state.hp.λ2
    param_var = state.v
    x = state.data[data_id]
    return _helper_log_predictive_likelihood(ss, hp_mean, hp_var, param_var, x)
end

function cluster_assignment_distribution(::Type{NonParametric}, data_id, state)
    """Compute the marginal distribution of cluster assignment
    for each cluster.
    """
    labels = []
    scores = []
    for cid in [keys(state.suffstats)..., typemax(Int)]
        score = log_predictive_likelihood(NonParametric, data_id, cid, state)
        score += log_cluster_assign_score(NonParametric, cid, state)
        score = exp(score)
        push!(labels, cid)
        push!(scores, score)
    end
    scores /= sum(scores)
    return labels, scores
end

function create_cluster!(state)
    cluster_id = maximum(keys(state.suffstats)) + 1
    state.suffstats[cluster_id] = SS(0, 0)
    push!(state.cluster_ids, cluster_id)
    return cluster_id
end

function destroy_cluster!(state, cluster_id)
    delete!(state.suffstats, cluster_id)
    filter!(x -> x != cluster_id, state.cluster_ids)
    return nothing
end

function prune_clusters!(state)
    for cid in state.cluster_ids
        if state.suffstats[cid].count == 0
            destroy_cluster!(state, cid)
        end
    end
end

function sample_assignment(data_id, state)
    """Sample new assignment from marginal distribution.
    If cluster is "`new`", create a new cluster.
    """
    labels, scores = cluster_assignment_distribution(NonParametric, data_id, state)
    cid = labels[rand(Categorical(scores))]
    if cid == typemax(Int)
        return create_cluster!(state)
    else
        return cid
    end
end

function gibb_step!(::Type{NonParametric}, state)
    """Collapsed Gibbs sampler for Dirichlet Process Mixture Model
    """
    both = zip(state.data, state.assignment)
    for (data_id, (datapoint, cid)) in enumerate(both)
        state.suffstats[cid] = remove_datapoint_from_suffstats(
            datapoint, state.suffstats[cid]
        )
        prune_clusters!(state)
        cid = sample_assignment(data_id, state)
        state.assignment[data_id] = cid
        state.suffstats[cid] = add_datapoint_to_suffstats(
            state.data[data_id], state.suffstats[cid]
        )
    end
end

function plot_cluster!(state::GibbsState; opts...)
    data = state.data
    assignment = state.assignment
    histogram!(data; group=assignment, alpha=0.5, opts...)
    ylims!(0, 150)
    xlims!(-8, 8)
    return nothing
end

function run(T, n_clusters)
    Random.seed!(0)
    data = generate(1000)
    state = GibbsState(data, n_clusters; α=0.01)
    n_steps = 5
    fig = plot(; legend=false, layout=(n_steps, 1), size=(400, 600))
    for i in 1:n_steps
        plot_cluster!(state; subplot=i)
        annotate!(-5, 120, "step = $((i-1))"; subplot=i)
        gibb_step!(T, state)
    end
    display(fig)

    return nothing
end

function animation(T, n_clusters)
    Random.seed!(0)
    data = generate(1000)
    state = GibbsState(data, n_clusters; α=0.1)

    anim = @animate for i in 1:50
        plot(; legend=false)
        plot_cluster!(state)
        annotate!(-5, 120, "step = $((i-1))")
        gibb_step!(T, state)
    end

    return gif(anim, "assets/$(model_name(T))_mixture_collapsed_gibbs.gif"; fps=1)
end

end


text:
using Flux,Functors
const tol = 1e-18
invr = (r; p = 1, kwargs...) -> r == 0 ? 0.0 : 1 / r^p

DIFF_COEFFS = ((0.0, 0.5), (-2.0, 1.0))
function δn(x, dx, n, dim)
    i = round(Int, x / dx)
    if abs(x - i) < tol
        return i < length(DIFF_COEFFS[n]) ? DIFF_COEFFS[n][i+1] / dx^dim : 0.0
    end
    return 0.0
end

"""
radial function
"""
struct Radfunc
    # c::AbstractVector
    # μ::AbstractVector
    # σ::AbstractVector
    p::AbstractVector
    f
    rmin::AbstractFloat
    rmax::AbstractFloat
end
@functor Radfunc

# Flux.trainable(m::Radfunc) = [m.p]
Flux.trainable(m::Radfunc) = [m.f,m.p]
# Flux.trainable(m::Radfunc) = [m.c, m.μ, m.σ]

function Radfunc(; rmin = 0.0, rmax = nothing, n = 8)
    # r = Array(LinRange(rmin, rmax, n))
    # c = exp.(-r / rmax)
    # μ = 1.0r
    # σ = 1.0r.+rmax/n
    # Radfunc(c, μ, σ, rmin, rmax)
    p=ones(4)
    f=Chain(Dense(1,n,leakyrelu),Dense(n,1))
    Radfunc(p,f, rmin, rmax)
end

const T=Float64
function (m::Radfunc)(r)::T
    @unpack f, rmin, rmax ,p= m
    # @unpack f, rmin, rmax = m
    # @unpack c, μ, σ, rmin, rmax = m
    if r < m.rmin - tol || r > m.rmax + tol
        return convert(T,0.0)
    end
    # dot(c, exp.(-((r .- μ) ./ σ) .^ 2))
    # (f([r])-f([rmax]))[1]
    c,k=p[1:2]
    # a=p[3:6]
    # c*exp(-abs(k)*r)*dot(a,r.^(0:3))
    convert(T,c*exp(-abs(k)*r)*f([r])[1])
end

max_stars_repo_name: aced-differentiate/equiv-docs

max_issues_repo_name: aced-differentiate/equiv-docs

max_forks_repo_name: aced-differentiate/equiv-docs


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule GAP_pkg_normalizinterface_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("GAP_pkg_normalizinterface")
JLLWrappers.@generate_main_file("GAP_pkg_normalizinterface", UUID("74f4ad05-5f81-5a78-8197-4e83b84c9c17"))
end  # module GAP_pkg_normalizinterface_jll


text:
# define test kernel
function test_kernel!(
    b  :: Matrix{Float64},
    x  :: Float64,
    dx :: Float64
    )  :: Nothing

    for i in eachindex(b)
        b[i] += dx * x^i * exp(-x^2) * sin(2.0 * pi * x)
    end

    return nothing
end

# define benchmark kernel
function bench_kernel!(
    b :: Matrix{Float64},
    x :: Float64
    ) :: Nothing

    for i in eachindex(b)
        b[i] = x^i * exp(-x^2) * sin(2.0 * pi * x)
    end

    return nothing
end

"""
    test_flow() :: Nothing
Run consistency checks for flow equations by testing integrators.
"""
function test_flow() :: Nothing

    # benchmark handcrafted integrators against QuadGK
    test_integrators()

    return nothing
end

"""
    test_integrators() :: Nothing
Run consistency checks for integrators by computing test integrals and comparing to QuadGK.
"""
function test_integrators() :: Nothing

    # init test dummys
    b1 = zeros(Float64, 10, 10)
    b2 = (copy(b1), copy(b1), copy(b1))
    b3 = (copy(b1), copy(b1), copy(b1))

    @testset "quadrature" begin
        # compute integral with QuadGK
        quadgk!((b, x) -> bench_kernel!(b, x), b1, 1.0, 5.0, atol = 1e-8, rtol = 1e-8)

        # compute integral with simps!
        integrate_lin!((b, x, dx) -> test_kernel!(b, x, dx), b2, 1.0, 5.0, 100, 1e-8, 1e-8)
        integrate_log!((b, x, dx) -> test_kernel!(b, x, dx), b3, 1.0, 5.0, 100, 1e-8, 1e-8)

        @test b1 ≈ b2[1]
        @test b1 ≈ b3[1]
    end

    return nothing
end