text:
"""
    reduce(step, xf, reducible; [init, simd, basesize, stoppable]) :: T

Thread-based parallelization of [`foldl`](@ref).  The "bottom"
reduction function `step(::T, ::T) :: T` must be associative and
`init` must be its identity element.

Transducers composing `xf` must be stateless (e.g., [`Map`](@ref),
[`Filter`](@ref), [`Cat`](@ref), etc.) except for [`ScanEmit`](@ref).
Note that [`Scan`](@ref) is not supported (although possible in
theory).  Early termination requires Julia ≥ 1.3.

Use [`tcollect`](@ref) or [`tcopy`](@ref) to collect results into a
container.

See also: [Parallel processing tutorial](@ref tutorial-parallel),
[`foldl`](@ref), [`dreduce`](@ref).

# Keyword Arguments
- `basesize::Integer = amount(reducible) ÷ nthreads()`: A size of
  chunk in `reducible` that is processed by each worker.  A smaller
  size may be required when:
    * computation time for processing each item fluctuates a lot
    * computation can be terminated by [`reduced`](@ref) or
      transducers using it, such as [`ReduceIf`](@ref)
- `stoppable::Bool`: [This option usually does not have to be set
  manually.]  Transducers.jl's `reduce` executed in the "stoppable"
  mode used for optimizing reduction with [`reduced`](@ref) has a
  slight overhead if `reduced` is not used.  This mode can be disabled
  by passing `stoppable = false`.  It is usually automatically
  detected and set appropriately.  Note that this option is purely for
  optimization and does not affect the result value.
- For other keyword arguments, see [`foldl`](@ref).

!!! compat "Transducers.jl 0.4.23"

    Keyword option `stoppable` requires at least Transducers.jl 0.4.23.

# Examples
```jldoctest
julia> using Transducers

julia> reduce(+, Map(exp) |> Map(log), 1:3)
6.0

julia> using BangBang: append!!

julia> reduce(append!!, Map(x -> 1:x), 1:2; basesize=1, init=Union{}[])
3-element Array{Int64,1}:
 1
 1
 2
```
"""
Base.reduce

"""
    mapreduce(xf, step, reducible; init, simd)

!!! warning

    `mapreduce` exists primary for backward compatibility.  It is
    recommended to use `reduce`.

Like [`reduce`](@ref) but `step` is _not_ automatically wrapped by
[`Completing`](@ref).
"""
Base.mapreduce

struct SizedReducible{T} <: Reducible
    reducible::T
    basesize::Int
end

foldable(reducible::SizedReducible) = reducible.reducible

"""
    Transducers.issmall(reducible, basesize) :: Bool

Check if `reducible` collection is considered small compared to
`basesize` (an integer).  Fold functions such as [`reduce`](@ref)
switches to sequential `__foldl__` when `issmall` returns `true`.

Default implementation is `amount(reducible) <= basesize`.
"""
issmall

issmall(reducible, basesize) = amount(reducible) <= basesize

issmall(reducible::SizedReducible) =
    issmall(reducible.reducible, max(reducible.basesize, 1))

function _halve(reducible::SizedReducible)
    left, right = halve(reducible.reducible)
    return (
        SizedReducible(left, reducible.basesize),
        SizedReducible(right, reducible.basesize),
    )
end

struct TaskContext
    listening::Vector{Threads.Atomic{Bool}}
    cancellables::Vector{Threads.Atomic{Bool}}
end

TaskContext() = TaskContext([], [])

function splitcontext(ctx::TaskContext)
    c = Threads.Atomic{Bool}(false)
    return (
        fg = TaskContext(ctx.listening, vcat(ctx.cancellables, c)),
        bg = TaskContext(vcat(ctx.listening, c), ctx.cancellables),
    )
end

function should_abort(ctx::TaskContext)
    for c in ctx.listening
        c[] && return true
    end
    return false
end

function cancel!(ctx::TaskContext)
    for c in ctx.cancellables
        c[] = true
    end
end

struct DummyTask end
Base.schedule(::DummyTask) = nothing

function transduce_assoc(
    xform::Transducer,
    step::F,
    init,
    coll0;
    simd::SIMDFlag = Val(false),
    basesize::Integer = amount(coll0) ÷ Threads.nthreads(),
    stoppable::Union{Bool,Nothing} = nothing,
) where {F}
    rf0 = _reducingfunction(xform, step; init = init)
    rf, coll = retransform(rf0, coll0)
    if stoppable === nothing
        stoppable = _might_return_reduced(rf, init, coll)
    end
    acc = @return_if_reduced _transduce_assoc_nocomplete(
        maybe_usesimd(rf, simd),
        init,
        coll,
        basesize,
        stoppable,
    )
    result = complete(rf, acc)
    if unreduced(result) isa DefaultInitOf
        throw(EmptyResultError(rf))
        # See how `transduce(rf, init, coll)` is implemented in ./processes.jl
    end
    return result
end

if VERSION >= v"1.3-alpha"
    maybe_collect(coll) = coll
else
    maybe_collect(coll::AbstractArray) = coll
    maybe_collect(coll) = collect(coll)
end

function _transduce_assoc_nocomplete(
    rf::F,
    init,
    coll,
    basesize,
    stoppable = true,
) where {F}
    reducible = SizedReducible(maybe_collect(coll), basesize)
    @static if VERSION >= v"1.3-alpha"
        return _reduce(TaskContext(), stoppable, DummyTask(), rf, init, reducible)
    else
        return _reduce_threads_for(rf, init, reducible)
    end
end

@noinline _reduce_basecase(rf::F, init::I, reducible) where {F,I} =
    restack(foldl_nocomplete(rf, start(rf, init), foldable(reducible)))
# `restack` here is crucial when using heap-allocated accumulator.
# See `ThreadsX.unique` and the MWE extracted from it:
# https://github.com/tkf/Restacker.jl/blob/master/benchmark/bench_unique.jl

function _reduce(
    ctx,
    stoppable,
    next_task,
    rf::R,
    init::I,
    reducible::Reducible,
) where {R,I}
    if should_abort(ctx)
        # As other tasks may be calling `fetch` on `next_task`, it
        # _must_ be scheduled at some point to avoid dead lock:
        stoppable && schedule(next_task)
        # Maybe use `error=false`?  Or pass something and get it via `yieldto`?
        return init
    end
    if issmall(reducible)
        stoppable && schedule(next_task)
        acc = _reduce_basecase(rf, init, reducible)
        if acc isa Reduced
            cancel!(ctx)
        end
        return acc
    else
        left, right = _halve(reducible)
        fg, bg = splitcontext(ctx)
        task = nonsticky!(@task _reduce(bg, stoppable, next_task, rf, init, right))
        stoppable || schedule(task)
        a0 = _reduce(fg, stoppable, task, rf, init, left)
        b0 = fetch(task)
        a = @return_if_reduced a0
        should_abort(ctx) && return a  # slight optimization
        b0 isa Reduced && return combine_right_reduced(rf, a, b0)
        return combine(rf, a, b0)
    end
end

combine_right_reduced(rf, a, b0::Reduced) =
    reduced(combine(_realbottomrf(rf), a, unreduced(b0)))

function _reduce_threads_for(rf, init, reducible::SizedReducible{<:AbstractArray})
    arr = reducible.reducible
    basesize = reducible.basesize
    nthreads = max(
        1,
        basesize <= 1 ? length(arr) : length(arr) ÷ basesize
    )
    if nthreads == 1
        return foldl_nocomplete(rf, start(rf, init), arr)
    else
        w = length(arr) ÷ nthreads
        results = Vector{Any}(undef, nthreads)
        Threads.@threads for i in 1:nthreads
            if i == nthreads
                chunk = @view arr[(i - 1) * w + 1:end]
            else
                chunk = @view arr[(i - 1) * w + 1:i * w]
            end
            results[i] = foldl_nocomplete(rf, start(rf, init), chunk)
        end
        # It can be done in `log2(n)` for loops but it's not clear if
        # `combine` is compute-intensive enough so that launching
        # threads is worth enough.  Let's merge the `results`
        # sequentially for now.
        return combine_all(rf, results)
    end
end

function combine_all(rf, results)
    step = combine_step(rf)
    return transduce(ensurerf(Completing(step)), Init(step), results)
end

combine_step(rf) =
    asmonoid() do a0, b0
        a = @return_if_reduced a0
        b0 isa Reduced && return combine_right_reduced(rf, a, b0)
        return combine(rf, a, b0)
    end

# The output of `reduce` is correct regardless of the value of
# `stoppable`.  Thus, we can use `return_type` here purely for
# optimization.
_might_return_reduced(rf, init, coll) =
    Base.typeintersect(
        Core.Compiler.return_type(
            _reduce_dummy,  # simulate the output type of `_reduce`
            typeof((rf, init, coll)),
        ),
        Reduced,
    ) !== Union{}

_reduce_dummy(rf, init, coll) =
    __reduce_dummy(rf, init, SizedReducible(maybe_collect(coll), 1))

function __reduce_dummy(rf, init, reducible)
    if issmall(reducible)
        return _reduce_basecase(rf, init, reducible)
    else
        left, right = halve(reducible)
        a = _reduce_dummy(rf, init, left)
        b = _reduce_dummy(rf, init, right)
        a isa Reduced && return a
        b isa Reduced && return combine_right_reduced(rf, a, b)
        return combine(rf, a, b)
    end
end

# Disambiguation:
Base.mapreduce(
    xform::Transducer,
    step::F,
    itr::AbstractArray;
    init = DefaultInit,
    kwargs...,
) where {F} = unreduced(transduce_assoc(xform, step, init, itr; kwargs...))
Base.mapreduce(
    xform::Transducer,
    step::F,
    itr::AbstractArrayOrBroadcasted;
    init = DefaultInit,
    kwargs...,
) where {F} = unreduced(transduce_assoc(xform, step, init, itr; kwargs...))
Base.mapreduce(
    xform::Transducer,
    step::F,
    itr::Base.SkipMissing{<:AbstractArray};
    init = DefaultInit,
    kwargs...,
) where {F} = unreduced(transduce_assoc(xform, step, init, itr; kwargs...))
Base.mapreduce(
    xform::Transducer,
    step::F,
    itr::Number;
    init = DefaultInit,
    kwargs...,
) where {F} = unreduced(transduce_assoc(xform, step, init, itr; kwargs...))

Base.mapreduce(xform::Transducer, step::F, itr;
               init = DefaultInit, kwargs...) where {F} =
    unreduced(transduce_assoc(xform, step, init, itr; kwargs...))

Base.reduce(step::F, xform::Transducer, itr; kwargs...) where {F} =
    mapreduce(xform, Completing(step), itr; kwargs...)

Base.reduce(step::F, foldable::Foldable; kwargs...) where {F} =
    reduce(step, extract_transducer(foldable)...; kwargs...)

"""
    tcopy(xf::Transducer, T, reducible; basesize) :: Union{T, Empty{T}}
    tcopy(xf::Transducer, reducible::T; basesize) :: Union{T, Empty{T}}
    tcopy([T,] itr; basesize) :: Union{T, Empty{T}}

Thread-based parallel version of [`copy`](@ref).
Keyword arguments are passed to [`reduce`](@ref).

See also: [Parallel processing tutorial](@ref tutorial-parallel)
(especially [Example: parallel `collect`](@ref tutorial-parallel-collect)).

!!! compat "Transducers.jl 0.4.5"

    New in version 0.4.5.

!!! compat "Transducers.jl 0.4.8"

    `tcopy` now accepts iterator comprehensions and eductions.

# Examples
```jldoctest
julia> using Transducers

julia> tcopy(Map(x -> x => x^2), Dict, 2:2)
Dict{Int64,Int64} with 1 entry:
  2 => 4

julia> using TypedTables

julia> @assert tcopy(Map(x -> (a=x,)), Table, 1:1) == Table(a=[1])

julia> using StructArrays

julia> @assert tcopy(Map(x -> (a=x,)), StructVector, 1:1) == StructVector(a=[1])
```

`tcopy` works with iterator comprehensions and eductions (unlike
[`copy`](@ref), there is no need for manual conversion with
[`eduction`](@ref)):

```jldoctest; setup = :(using Transducers, StructArrays, DataFrames)
julia> table = StructVector(a = [1, 2, 3], b = [5, 6, 7]);

julia> @assert tcopy(
           (A = row.a + 1, B = row.b - 1) for row in table if isodd(row.a)
       ) == StructVector(A = [2, 4], B = [4, 6])

julia> @assert tcopy(
           DataFrame,
           (A = row.a + 1, B = row.b - 1) for row in table if isodd(row.a)
       ) == DataFrame(A = [2, 4], B = [4, 6])

julia> @assert tcopy(eduction(
           Filter(row -> isodd(row.a)) |> Map(row -> (A = row.a + 1, B = row.b - 1)),
           table,
       )) == StructVector(A = [2, 4], B = [4, 6])
```

If you have [`Cat`](@ref) or [`MapCat`](@ref) at the end of the
transducer, consider using [`reduce`](@ref) directly:

```jldoctest
julia> using Transducers
       using DataFrames

julia> @assert tcopy(
           Map(x -> DataFrame(a = [x])) |> MapCat(eachrow),
           DataFrame,
           1:2;
           basesize = 1,
       ) == DataFrame(a = [1, 2])

julia> using BangBang: Empty, append!!

julia> @assert reduce(
           append!!,
           Map(x -> DataFrame(a = [x])),
           1:2;
           basesize = 1,
           # init = Empty(DataFrame),
       ) == DataFrame(a = [1, 2])
```

Note that above snippet assumes that it is OK to mutate the dataframe
returned by the transducer.  Use `init = Empty(DataFrame)` if this is
not the case.

This approach of using `reduce` works with other containers; e.g.,
with `TypedTables.Table`:

```jldoctest; setup = :(using Transducers)
julia> using TypedTables

julia> @assert reduce(
           append!!,
           Map(x -> Table(a = [x])),
           1:2;
           basesize = 1,
           # init = Empty(Table),
       ) == Table(a = [1, 2])
```
"""
tcopy(xf, T, reducible; kwargs...) =
    reduce(append!!, xf |> Map(SingletonVector), reducible; init = Empty(T), kwargs...)
tcopy(xf, reducible; kwargs...) = tcopy(xf, _materializer(reducible), reducible; kwargs...)

function tcopy(::Type{T}, itr; kwargs...) where {T}
    xf, foldable = extract_transducer(itr)
    return tcopy(xf, T, foldable; kwargs...)
end

function tcopy(itr; kwargs...)
    xf, foldable = extract_transducer(itr)
    return tcopy(xf, foldable; kwargs...)
end

tcopy(xf, T::Type{<:AbstractSet}, reducible; kwargs...) =
    reduce(union!!, xf |> Map(SingletonVector), reducible; init = Empty(T), kwargs...)

function tcopy(
    ::typeof(Map(identity)),
    T::Type{<:AbstractSet},
    array::PartitionableArray;
    basesize::Integer = max(1, length(array) ÷ Threads.nthreads()),
    kwargs...,
)
    @argcheck basesize >= 1
    return reduce(
        union!!,
        Map(identity),
        Iterators.partition(array, basesize);
        init = Empty(T),
        basesize = 1,
        kwargs...,
    )
end

"""
    tcollect(xf::Transducer, reducible; basesize) :: Union{Vector, Empty{Vector}}
    tcollect(itr; basesize) :: Union{Vector, Empty{Vector}}

Thread-based parallel version of [`collect`](@ref).
This is just a short-hand notation of `tcopy(xf, Vector, reducible)`.
Use [`tcopy`](@ref) to get a container other than a `Vector`.

See also: [Parallel processing tutorial](@ref tutorial-parallel)
(especially [Example: parallel `collect`](@ref tutorial-parallel-collect)).

!!! compat "Transducers.jl 0.4.5"

    New in version 0.4.5.

!!! compat "Transducers.jl 0.4.8"

    `tcollect` now accepts iterator comprehensions and eductions.

# Examples
```jldoctest
julia> using Transducers

julia> tcollect(Map(x -> x^2), 1:2)
2-element Array{Int64,1}:
 1
 4

julia> tcollect(x^2 for x in 1:2)
2-element Array{Int64,1}:
 1
 4
```
"""
tcollect(xf, reducible; kwargs...) = tcopy(xf, Vector, reducible; kwargs...)
tcollect(itr; kwargs...) = tcollect(extract_transducer(itr)...; kwargs...)


text:
"""
  Model(expr, parts::Vector{FormulaTerm})

Representation of a multipart model composed of a vector of FormulaTerm's.
It is inspired by the [https://github.com/matthieugomez/FixedEffectModels.jl](FixedEffectModels.jl) package
"""
struct Model
  expr
  parts::Vector{FormulaTerm}
end

"""
    @model(formula1, formula2, ...)

Parses a sequence of comma-separated FormulaTerm's representing a multipart model.
"""
macro model(args...)
  Expr(:call, :model_helper, (esc(Base.Meta.quot(a)) for a in args)...)
end

function model_helper(args...)
  expr = args
  parts = FormulaTerm[]
  for part in args
    push!(parts, @eval(@formula($(part.args[2]) ~ $(part.args[3]))))
  end
  Model(expr,parts)
end

function Base.show(io::IO, m::Model)
  print(io, "$(length(m.parts))-part model: [")
  notfirst = false
  for p in m.parts
    if notfirst
      print(io, ", ")
    end
    print(io, p)
    notfirst = true
  end
  print(io, "]")
  nothing
end

function mergerhsterms(a, b)
  terms = union(a, b)

  ina = findall((in)(a), terms)
  inb = findall((in)(b), terms)

  (terms...,), ina, inb
end

"maps inzero and inpos to ModelMatrix columns (important with factor variables)"
function mapins(inzero, inpos, appliedschema)
  mappedinzero = Int[]
  mappedinpos = Int[]
  from = 1
  to = 0
  for t=1:length(appliedschema)
    to = from + width(appliedschema[t]) - 1
    if t in inzero
      [push!(mappedinzero, i) for i in from:to]
    end
    if t in inpos
      [push!(mappedinpos, i) for i in from:to]
    end
    from = to + 1
  end
  mappedinzero, mappedinpos
end

"""
  getrhsterms(m, lhs)

Returns the right-hand-side Terms of the part of model `m` with left-hand-side
symbol `lhs`.
"""
function getrhsterms(m::Model, lhs::Symbol)
  ix = findfirst(p->p.lhs==term(lhs),m.parts)
  @assert ix > 0 "The model is missing a formula with $lhs on its left-hand-side."
  f = m.parts[ix]
  f.rhs
end

StatsModels.missing_omit(data::T, formula::TupleTerm) where T<:ColumnTable =
    missing_omit(NamedTuple{tuple(termvars(formula)...)}(data))

function StatsModels.modelcols(trms::TupleTerm, data, counts::AbstractMatrix;
    model::Type{M}=DCR, contrasts=Dict{Symbol,Any}()) where M

  cols = columntable(data)
  cols, nonmissing = missing_omit(cols, trms)
  s = schema(trms, cols, contrasts)
  as = apply_schema(trms, s, M)
  modelcols(as, nonmissing, cols, counts)
end

function StatsModels.modelcols(as::TupleTerm, nonmissing, cols, counts::AbstractMatrix)

  covars = modelcols(MatrixTerm(as), cols)

  if !all(nonmissing)
    counts = counts[nonmissing,:]
  end

  covars, counts, as
end

"Similar to StatsModels.TableRegressionModel, but also holds a counts matrix and Model"
struct TableCountsRegressionModel{M,D,C} <: RegressionModel
  model::M    # actual fitted model (result of fit())
  data::D     # data table used
  counts::C   # counts matrix used
  f::Model    # @model specification
  schema      # applied schema
  sprojdir::Union{Nothing,Symbol}   # lhs variable in forward regression
end
TableCountsRegressionModel(model, data, counts, f, schema) = TableCountsRegressionModel(model, data, counts, f, schema, nothing)

function Base.show(io::IO, model::TableCountsRegressionModel)
    println(io, typeof(model))
    println(io)
    println(io, model.f)
end

ext: jl


text:
using BNPRegressionGGA2021
using Documenter

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

makedocs(;
    modules=[BNPRegressionGGA2021],
    authors="Iván Gutiérrez <ivangutierrez1988@gmail.com> and contributors",
    repo="https://github.com/igutierrezm/BNPRegressionGGA2021.jl/blob/{commit}{path}#{line}",
    sitename="BNPRegressionGGA2021.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://igutierrezm.github.io/BNPRegressionGGA2021.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/igutierrezm/BNPRegressionGGA2021.jl",
)


text:
module ServiceSolicitation

using ThreadingUtilities, VectorizationBase
using VectorizationBase: num_threads, cache_linesize
using StrideArraysCore: object_and_preserve
using Requires

export batch, num_threads

include("request.jl")
include("batch.jl")
include("unsignediterator.jl")

# reset_workers!() = WORKERS[] = UInt128((1 << (num_threads() - 1)) - 1)
reset_workers!() = WORKERS[] = UInt128((1 << (Threads.nthreads() - 1)) - 1)
function __init__()
    reset_workers!()
    resize!(STATES, num_threads() * cache_linesize())
    STATES .= 0x00
    @require ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" include("forwarddiff.jl")
end

end


text:
getMinimumVolume(r::Reservoir, c::MinimumReservoirLevelCondition) = computeVolume(bathymetry(r), getMinimumLevel(c))
getMinimumVolume(r::Reservoir, c::DamOutputCondition) = error("Minimum reservoir volume has no meaning for the condition.")
minVolume(r::Reservoir, c::MinimumReservoirLevelCondition) = getMinimumVolume(r, c)


text:
using GaussianMarkovRandomFields

using Random
using BenchmarkTools
using SparseArrays, LinearAlgebra
using LDLFactorizations
using Distributions, PDMats

using ForwardDiff
using ReverseDiff
# using Zygote # broken
using FiniteDiff
using Turing

using Test

function make_precision(ρ, k)
    return spdiagm(-1 => -ρ*ones(k-1), 0 => ones(k), 1 => -ρ*ones(k-1))
end

function test_variables()
    k = 50
    ρ = 0.4
    μ = zeros(k)
    Q = make_precision(ρ, k)
    return k, ρ, μ, Q
end

@testset "Constructors" begin
    k = 10
    T = Float32
    Q = sprandn(k, k, 0.1)
    Q = Q'Q + I
    μ = zeros(T, k)
    L, P = cholesky_ldl(Q)
    F = ldl(Q)

    g1 = GMRF(Q)
    g2 = GMRF(μ, Q)
    g3 = GMRF(μ, Q, L)
    g4 = GMRF(μ, F)
    g5 = GMRF(F)
    @test g1 == g2
    @test g1 == g3
    @test mean(g1) == μ
    @test prec(g1) == Q
    @test eltype(g2) == T
    for gf in [g4, g5]
        @test all(mean(gf) .≈ mean(g1))
        @test all(prec(gf) .≈ prec(g1))
    end
    #  @test all(var(g1) .≈ 1 ./ diag(Q))
    #  Tests for variance below
    @test_throws ErrorException cov(g1)
end

@testset "Linear Algebra" begin
    k, ρ, μ, Q = test_variables()

    L, P = cholesky_ldl(Q)
    L1 = sparse(cholesky(Q).L)
    @test all(Matrix(L * L') .≈ Matrix(Q)[P, P])
    @test prod(diag(L).^2) .≈ det(Q)
    @test all(Matrix(L) .≈ Matrix(sparse(cholesky(Q).L)))
    @test logdet_ldl(Q) .≈ log(det(Matrix(Q)))
    @test logdet_ldl(Q) .≈ logdet(Matrix(Q))
end

@testset "Likelihood" begin
    k, ρ, μ, Q = test_variables()
    Qd = Matrix(Q)

    d = GMRF(μ, Q)
    dd = MvNormalCanon(PDMat(Qd))
    dd1 = MvNormal(zeros(k), inv(Symmetric(Matrix(Q))))
    Random.seed!(1)
    x = rand(d)
    @test logpdf(d, x) ≈ logpdf(dd, x)
    @test logpdf(d, x) ≈ logpdf(dd1, x)
end

@testset "Performance" begin
    Random.seed!(1)
    k, ρ, μ, Q = test_variables()
    Qd = Matrix(Q)

    t_logdet_ldl = @benchmark logdet_ldl($Q)
    t_logdet_ldl_d = @benchmark logdet_ldl($Qd)
    t_logdet_d = @benchmark logdet($Qd)
    println(" ")
    println("logdet_ldl (sparse): $(minimum(t_logdet_ldl))")
    println("logdet_ldl (dense): $(minimum(t_logdet_ldl_d))")
    println("logdet (dense): $(minimum(t_logdet_d))")

    @test minimum(t_logdet_ldl.times) .< minimum(t_logdet_ldl_d.times)
    @test minimum(t_logdet_ldl.times) .< minimum(t_logdet_d.times)

    d = GMRF(μ, Q)
    d1 = MvNormalCanon(PDSparseMat(Q))
    dd = MvNormalCanon(PDMat(Qd))
    dd1 = MvNormal(zeros(k), inv(Symmetric(Matrix(Q))))
    x = rand(d)

    t_logpdf_s = @benchmark logpdf($d, $x)
    t_logpdf_s1 = @benchmark logpdf($d1, $x)
    t_logpdf_d = @benchmark logpdf($dd, $x)
    t_logpdf_d1 = @benchmark logpdf($dd1, $x)
    println(" ")
    println("logpdf (sparse, LDL): $(minimum(t_logpdf_s))")
    println("logpdf (sparse, MvNormalCanon): $(minimum(t_logpdf_s1))")
    println("logpdf (dense, MvNormalCanon): $(minimum(t_logpdf_d))")
    println("logdet (dense, MvNormal): $(minimum(t_logpdf_d1))")

    @test_broken minimum(t_logpdf_s.times) .< minimum(t_logpdf_s1.times)
    @test minimum(t_logpdf_s.times) .< minimum(t_logpdf_d.times)
    # @test minimum(t_logpdf_s.times) .< minimum(t_logpdf_d1.times)
end

@testset "AD" begin
    Random.seed!(1)
    k, ρ, μ, Q = test_variables()
    d = GMRF(μ, Q)
    x = rand(d)

    logistic(x) = 1 / (1 + exp(-x))

    function test_loglik(θ)
        σ = exp(θ[1])
        ρ = σ/2 * logistic(θ[2])
        μ = θ[3:end]
        Q = make_precision(ρ, length(μ))
        d = GMRF(μ, Q)
        return logpdf(d, x)
    end

    θ = randn(k + 2)
    Δfinite = FiniteDiff.finite_difference_gradient(test_loglik, θ)
    @test all(ForwardDiff.gradient(test_loglik, θ) .≈ Δfinite)
    @test all(ReverseDiff.gradient(test_loglik, θ) .≈ Δfinite)
    # Zygote broken for some reason
    # @test all(Zygote.gradient(test_loglik, θ) .≈ Δfinite)
end

@testset "Turing" begin
    Random.seed!(1)
    k, ρ, μ, Q = test_variables()
    d = GMRF(μ, Q)
    x = rand(d)

    @model function GMRFTest(x)
        n = length(x)
        ρ ~ Uniform(0, 0.5)
        Q = make_precision(ρ, n)
        x ~ GMRF(Q)
    end

    m = GMRFTest(x)
    Turing.setadbackend(:forwarddiff)
    @test typeof(sample(m, NUTS(), 50)) <: Turing.Chains
    Turing.setadbackend(:reversediff)
    @test typeof(sample(m, NUTS(), 50)) <: Turing.Chains
end

@testset "Marginal variance" begin
    ### Setup AR1 example
    function ar1_precision(n::I, ρ::T) where {I<:Integer,T}
        offdiagvec = -ρ * ones(T, n - 1)
        diagvec = T[one(ρ); ones(T, n - 2) .+ ρ^2; one(ρ)]
        return  spdiagm(-1 => offdiagvec, 0 => diagvec, 1 => offdiagvec)
    end
    function ar1_precision(n, ρ, σ)
        ρ, σ = promote(ρ, σ)
        return 1 / σ^2 * ar1_precision(n, ρ)
    end
    #  function ar1_precision(n, ρ::T, σ::AbstractVector{T}) where T
        #  scalemat = Diagonal(σ)
        #  return scalemat * ar1_precision(n, ρ) * scalemat
    #  end

    k, ρ, μ, Q = test_variables()

    Q1 = ar1_precision(k, ρ)
    Q2 = ar1_precision(k, ρ, 2.0)
    # Qns = ar1_precision(k, ρ, range(0.1, 1.9, length = k))

    n = 10_000   # number of measurements in time series
    σ = 2.0      # innovation variance
    ρ = 0.4      # autocorrelation coefficient

    Qbig = ar1_precision(n, ρ, σ)
    dbig = GMRF(Qbig)

    # Compare to the analytic marginal variance for an AR(1) process and ensure
    # that we don't overflow the stack on large precision matrices.
    @test all(var(dbig) .≈ σ^2 / (1 - ρ^2))

    # Test that all elements are computed accurately
    Q_small = ar1_precision(k, ρ)
    # Convert to generic sparse matrix
    L_small = dropzeros(sparse(cholesky(Matrix(Q_small)).L))
    # Invert directly to get full covariance matrix for comparison
    Σ_small = inv(Matrix(Q_small))

    Σ_elem = similar(Σ_small)
    for j in k:-1:1, i in k:-1:1
        Σ_elem[i, j] = GaussianMarkovRandomFields.cov_element(L_small, i, j)
    end

    @test all(Σ_elem .≈ Σ_small)
end


text:
@testset "Check field membership" begin
    R, (x, y, z) = Nemo.PolynomialRing(Nemo.QQ, ["x", "y", "z"])

    @test check_field_membership(
        [[R(1), x + y], [R(1), x * y], [z, (x + y)^2]],
        [(x^2 + y^2) // R(1), (x^3 + y^3) // (z - x * y), R(1) // (z + x + y), z // x],
        0.99
    ) == [true, true, true, false]

    @test check_field_membership(
        [[x + y + z, x^2 + y^2 + z^2, (x + y + z)^2, x^3 + y^3 + z^3]],
        [x + y + z, x * y * z, x + y + 2 * z, x // (y + z)],
        0.99
    ) == [true, true, false, false]
end


text:
using Test, Random
import LuxurySparse: PermMatrix, pmrand

Random.seed!(2)
p1 = PermMatrix([1,4,2,3],[0.1, 0.2, 0.4im, 0.5])
p2 = PermMatrix([2,1,4,3],[0.1, 0.2, 0.4, 0.5])
#p3 = PermMatrix([4,1,2,3],[0.5, 0.4im, 0.3, 0.2])
p3 = pmrand(4)
sp = sprand(4, 4, 0.3)
v = [0.5, 0.3im, 0.2, 1.0]

@testset "basic" begin
    @test p1==copy(p1)
    @test eltype(p1) == ComplexF64
    @test eltype(p2) == Float64
    @test eltype(p3) == Float64
    @test size(p1) == (4, 4)
    @test size(p3) == (4, 4)
    @test size(p1, 1) == size(p1, 2) == 4
    @test Matrix(p1) == [0.1 0 0 0; 0 0 0 0.2; 0 0.4im 0 0; 0 0 0.5 0]
end

@testset "sparse" begin
    @test nnz(p1) == 4
    @test nonzeros(p1) == p1.vals
    @test dropzeros!(p1) == p1
end

@testset "linalg" begin
    @test inv(p1) * p1 ≈ Matrix(I, 4, 4)
    @test p1*transpose(p1) == diagm(0=>p1.vals.^2)
    #@test p1*adjoint(p1) == diagm(0=>abs.(p1.vals).^2)
    #@test all(isapprox.(adjoint(p3), transpose(conj(Matrix(p3)))))
    @test p1*p1' == diagm(0=>abs.(p1.vals).^2)
    @test all(isapprox.(p3', transpose(conj(Matrix(p3)))))
end

@testset "mul" begin
    @test p3*p2 == SparseMatrixCSC(p3)*p2 == Matrix(p3)*p2

    # Multiply vector
    @test p3*v == Matrix(p3)*v
    @test v'*p3 == v'*Matrix(p3)
    @test p3*collect(1:4) == p3.perm.*p3.vals

    # Diagonal matrices
    Dv = Diagonal(v)
    @test p3*Dv == Matrix(p3)*Dv
    @test Dv*p3 == Dv*Matrix(p3)
end

@testset "elementary" begin
    @test all(isapprox.(conj(p1), conj(Matrix(p1))))
    @test all(isapprox.(real(p1), real(Matrix(p1))))
    @test all(isapprox.(imag(p1), imag(Matrix(p1))))
end

@testset "basicmath" begin
    @test p1*2 == Matrix(p1)*2
    @test p1/2 == Matrix(p1)/2
end

@testset "memorysafe" begin
    @test p1 == PermMatrix([1,4,2,3],[0.1, 0.2, 0.4im, 0.5])
    @test p2 == PermMatrix([2,1,4,3],[0.1, 0.2, 0.4, 0.5])
    @test v == [0.5, 0.3im, 0.2, 1.0]
end


text:
module TestLB

using LibCloud.LB
using Base.Test
using PyCall

function common(lb)
    println("Regions:\n", list_regions(lb))
    println("Algorithms:\n", list_supported_algorithms(lb))
    println("Protocols:\n", list_protocols(lb))

    for balancer in list_balancers(lb)
        @test isa(balancer, LoadBalancer)
        println(balancer)
        try
            for member in list_members(balancer)
                @test isa(member, Member)
                println("    - ", member)
            end
        end
    end
end

function elb(key, secret, region="us-east-1")
    lb = LBDriver(LBProvider.ELB, key, secret, region)
    common(lb)
end

function gce(service_account_id, key, project)
	lb = LBDriver(LBProvider.GCE, service_account_id, key, project=project)
    common(lb)
end

end # module TestLB


text:
using DesignAutomation

@testset "Cyclic Coordinate Search" begin
    fn(x) = (x[1]-5)^2+2*(x[2]-5)^2
    x0 = [0., 0.]
    res = ccs(fn, x0)

    @test res.minimizer ≈ [5.0, 5.0]
    @test res.minimum + 1 ≈ 1 # answer is approx 0
    @test res.iterations == 3
    @test res.converged == true

    fn(x) = (2*x[2]-3*x[1]^2)^4 + (5-2*x[1])^2
    x0 = [0., 0.]
    res = ccs(fn, x0)

    @test res.minimizer ≈ [2.5, 9.375]
    @test res.minimum + 1 ≈ 1 # answer is approx 0
    @test res.iterations == 4
    @test res.converged == true
end

max_stars_repo_name: camirmas/DesignAutomation

max_issues_repo_name: camirmas/DesignAutomation

max_forks_repo_name: camirmas/DesignAutomation


text:
for T in subtypes(AbstractPointwise)
    @eval begin
        function (ρ::$T)(node::AbstractNode)
            return Node(ρ(node.data), node.ranges)
        end
    end
end


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

JLLWrappers.@generate_main_file_header("librealsense")
JLLWrappers.@generate_main_file("librealsense", UUID("65edd895-dc1e-518e-896c-dc343a3c7f2e"))
end  # module librealsense_jll


text:
# Command line args
!isempty(ARGS) || error("No config supplied.")
isfile(ARGS[1]) || error("Cannot read '$(ARGS[1])'")
isabspath(ARGS[1]) || error("Please use an absolute path for the config.")
println("Config supplied: '$(ARGS[1])'")
using SVDD, OneClassActiveLearning
include(ARGS[1])

# Load packages
using Serialization, Distributed, Logging

# setup workers
localhost = filter(x -> x[1] == "localhost", worker_list)
remote_servers = filter(x -> x[1] != "localhost", worker_list)

length(localhost) > 0 && addprocs(localhost[1][2], exeflags=exeflags)
length(remote_servers) > 0 && addprocs(remote_servers, sshflags=sshflags, exeflags=exeflags)

# validate package versions
@everywhere function get_git_hash(path)
    cmd = `git -C $path rev-parse HEAD`
    (gethostname(), strip(read(cmd, String)))
end

function setup_julia_environment()
    local_githash = get_git_hash(JULIA_ENV)[2]
    for id in workers()
        remote_name, remote_githash = remotecall_fetch(get_git_hash, id, JULIA_ENV)
        @assert remote_githash == local_githash "Host: $remote_name has version mismatch." *
                                                    "Hash is '$remote_githash' instead of '$local_githash'."
    end
end

# Load remote packages and functions
@info "Loading packages on all workers."
@everywhere using Pkg
setup_julia_environment()
@everywhere using SVDD, OneClassActiveLearning, Memento, Gurobi, Random
@everywhere import SVDD: RandomOCClassifier

@everywhere fmt_string = "[{name} | {date} | {level}]: {msg}"
@everywhere loglevel = "debug"

@everywhere function setup_logging(experiment)
    setlevel!(getlogger("root"), "error")
    setlevel!(getlogger(OneClassActiveLearning), loglevel)
    setlevel!(getlogger(SVDD), loglevel)

    exp_logfile = joinpath(experiment[:log_dir], "experiment", "$(experiment[:hash]).log")
    worker_logfile = joinpath(experiment[:log_dir], "worker", "$(gethostname())_$(getpid()).log")

    WORKER_LOGGER = Memento.config!("runner", "debug"; fmt=fmt_string)

    exp_handler = DefaultHandler(exp_logfile, DefaultFormatter(fmt_string))
    push!(getlogger(OneClassActiveLearning), exp_handler, experiment[:hash])
    push!(getlogger(SVDD), exp_handler, experiment[:hash])
    push!(WORKER_LOGGER, exp_handler, experiment[:hash])

    worker_handler = DefaultHandler(worker_logfile, DefaultFormatter(fmt_string))
    setlevel!(gethandlers(WORKER_LOGGER)["console"], "error")
    push!(WORKER_LOGGER, worker_handler)

    return WORKER_LOGGER
end

@everywhere function cleanup_logging(worker_logger::Logger, experiment_hash)
    delete!(gethandlers(getlogger(OneClassActiveLearning)), experiment_hash)
    delete!(gethandlers(getlogger(SVDD)), experiment_hash)
    delete!(gethandlers(worker_logger), experiment_hash)
    return nothing
end

@everywhere function Memento.warn(logger::Logger, error::Gurobi.GurobiError)
    Memento.warn(logger, "GurobiError(exit_code=$(error.code), msg='$(error.msg)')")
end

@everywhere function Memento.debug(logger::Logger, error::Gurobi.GurobiError)
    Memento.warn(logger, "GurobiError(exit_code=$(error.code), msg='$(error.msg)')")
end

@everywhere function Memento.debug(logger::Logger, error::ErrorException)
    Memento.warn(logger, "Caught ErrorException, msg='$(error.msg)')")
end

@everywhere function run_experiment(experiment::Dict)
    # Make experiments deterministic
    Random.seed!(0)

    WORKER_LOGGER = setup_logging(experiment)
    info(WORKER_LOGGER, "host: $(gethostname()) worker: $(getpid()) exp_hash: $(experiment[:hash])")
    if isfile(experiment[:output_file])
        warn(WORKER_LOGGER, "Aborting experiment because the output file already exists. Filename: $(experiment[:output_file])")
        cleanup_logging(WORKER_LOGGER, experiment[:hash])
        return nothing
    end

    res = Result(experiment)
    errorfile = joinpath(experiment[:log_dir], "worker", "$(gethostname())_$(getpid())")
    try
        time_exp = @elapsed res = OneClassActiveLearning.active_learn(experiment)
        res.al_summary[:runtime] = Dict(:time_exp => time_exp)
    catch e
        res.status[:exit_code] = Symbol(typeof(e))
        @warn "Experiment $(experiment[:hash]) finished with unkown error."
        @warn e
        @warn stacktrace(catch_backtrace())
        warn(WORKER_LOGGER, "Experiment $(experiment[:hash]) finished with unkown error.")
        warn(WORKER_LOGGER, e)
    finally
        if res.status[:exit_code] != :success
            info(WORKER_LOGGER, "Writing error hash to $errorfile.error.")
            open("$errorfile.error", "a") do f
                print(f, "$(experiment[:hash])\n")
            end
        end
        info(WORKER_LOGGER, "Writing result to $(experiment[:output_file]).")
        OneClassActiveLearning.write_result_to_file(experiment[:output_file], res)
        cleanup_logging(WORKER_LOGGER, experiment[:hash])
    end

    return nothing
end

# load and run experiments
all_experiments = []
for s in readdir(data_output_root)
    if occursin(".csv", s)
        continue
    end
    @info "Running experiments in directory $s"
    exp_dir = joinpath(data_output_root, s)
    @info "Loading experiments.jld"
    # load experiments
    experiments = deserialize(open(joinpath(exp_dir, "experiments.jser")))
    append!(all_experiments, experiments)
end
@info "Running $(length(all_experiments)) experiments."
@info "Running experiments..."
pmap(run_experiment, all_experiments)
@info "Done."

# cleanup
rmprocs(workers())

max_stars_repo_name: kit-dbis/ocal-evaluation

max_issues_repo_name: kit-dbis/ocal-evaluation

max_forks_repo_name: kit-dbis/ocal-evaluation


text:
"""
function beam_test()
Test the time-dependent and steady state beam problem.

"""
function beam_test()
    #
    dt = 0.02
    n = 20
    stepnum = 5
    (t, se, xe, fhist, fhistt) = ivpBeam(n, dt, stepnum)
    beamtdout = (length(fhist) == 6) && (norm(fhistt, Inf) < 5.e-5)
    beamtdout || println("error in beam_test beamtdout")
    (pout, nout) = ptcBeam(10, 100)
    pout2 = ptcBeam(10,100; jknowsdt=true)
    kdtdiff=norm(pout.solution-pout2.solution,Inf) +
     norm(pout2.history-pout.history,Inf)
    kdtok=(kdtdiff < 1.e-14)
    kdtok || println("error is knowsdt test")
    nsolp = norm(pout.solution)
    nsoln = norm(nout.solution)
    itp = length(pout.history)
    pnormok = (nsolp > 5.0) && (nsoln < 1.e-15)
    pnormok || println("error in beam_test pnromok")
    presok = (itp < 100) && (pout.history[itp] < 1.e-10)
    presok || println("error in beam_test presok")
    return beamtdout && pnormok && presok && kdtok
end


text:
"""
	init_eye(S)

Initialize `O[i]` by an `di x r` diagonal matrix with ones on diagonal.
"""
function init_eye(
    S :: Matrix{Matrix{T}},
    r :: Integer
    ) where T <: BlasReal
    [Matrix{T}(I, size(S[i, i], 1), r) for i in 1:size(S, 1)]
end

"""
	init_tb(S)

Compute initial point following Ten Berge's second upper bound
(<https://doi.org/10.1007/BF02294053>, p. 273). This is same as the
Liu-Wang-Wang starting point strategy 2 (<https://doi.org/10.1137/15M100883X>,
Algorithm 3.2, p. 1495). Take the eigenvectors corresponding to the `r` largest
eigenvalues of `S`. This is a `D x r` orthogonal matrix, ``D=\\sum_{i=1}^m d_i``.
For each `di x r` block, project to the Stiefel manifold.
These blocks constitute an initial point.
"""
function init_tb(
    S :: Matrix{Matrix{T}},
    r :: Integer
    ) where T <: BlasReal
    m = size(S, 1)
    # only need evecs corresponding to largest r evals
	V = (eigen!(Symmetric(cat(S))).vectors)[:, end:-1:end-r+1]
    rowoffset = 0
    O = Vector{Matrix{T}}(undef, m)
	for i in 1:m
		di     = size(S[i, i], 1)
		rowidx = rowoffset+1:rowoffset+di
		Visvd  = svd!(V[rowidx, :])
		O[i]   = Visvd.U * Visvd.Vt
		rowoffset += di
	end
	return O
end

"""
	init_sb(S, r)

Compute initial point by Shapiro-Botha (<https://doi.org/10.1137/0609032>, p. 380).
Also see the extension by Won-Zhou-Lange paper (Lemma B.1). Replace diagonal
blocks `S[i, i]` by ``S_{i,i} - \\sum_{j} P_{ij} * D_{ij} * P_{ij}``, where
``P_{ij} * D_{ij} * Q_{ij}`` is the SVD of ``S_{ij}``. The resulting matrix is
negative semidefinite. Take the eigenvectors corresponding to the `r` largest
eigenvalues. This is a `D x r` orthogonal matrix, ``D=\\sum_{i=1}^m d_i``.
For each `di x r` block, project to the Stiefel manifold.
These blocks constitute an initial point.
"""
function init_sb(
    S :: Matrix{Matrix{T}},
    r :: Integer
    ) where T <: BlasReal
	m    = size(S, 1)
    SS   = similar(S)
    Ssvd = [svd(S[i, j]) for i in 1:m, j in 1:m] # wasteful, not using symmetry
    @inbounds for i in 1:m
        SS[i, i] = copy(S[i, i])
        for j in 1:m
            (j ≠ i) && (SS[i, j] = copy(S[i, j]))
            svdij      = Ssvd[i, j]
            SS[i, i] .-= svdij.U * Diagonal(svdij.S) * transpose(svdij.U)
		end
	end
	# only need evecs corresponding to largest r evals
	V = (eigen!(Symmetric(cat(SS))).vectors)[:, end:-1:end-r+1]
    O = Vector{Matrix{T}}(undef, m)
    rowoffset = 0
	@inbounds for i in 1:m
		di     = size(SS[i, i], 1)
		rowidx = rowoffset+1:rowoffset+di
        Visvd  = svd!(V[rowidx, :])
        O[i]   = Visvd.U * Visvd.Vt
		rowoffset += di
	end
	return O
end

"""
	init_lww1(S, r)

Compute initial point following Liu-Wang-Wang starting point strategy 1
(<https://doi.org/10.1137/15M100883X>, Algorithm 3.1, p. 1494). Set `O[1]` to
the top `r` eigenvectors of `S[1, 1]`. Then `Ok = Uk * Qk`, where `Uk` is the
top `r` eigenvectors of `S[k, k]` and `Qk` is the Q factor in the QR decomposition
of `Uk * sum_{j<k} S[k, j] * O[j]`.
"""
function init_lww1(
    S :: Matrix{Matrix{T}},
    r :: Integer
    ) where T <: BlasReal
    m = size(S, 1)
    O = Vector{Matrix{T}}(undef, m)
    for k in 1:m
        if k == 1
            O[k] = eigen(Symmetric(S[k, k])).vectors[:, end:-1:end-r+1]
        else
            Uk = eigen(Symmetric(S[k, k])).vectors[:, end:-1:end-r+1]
            @views O[k] = Uk * qr(transpose(Uk) * (S[k:k, 1:k-1] * O[1:k-1])[1]).Q
        end
    end
    return O
end


text:
using Soss

gps = @model speed,n,t begin
    x ~ For(1:n) do i
        Normal(μ = speed * t[i])
    end
    return x
end

radar = @model speed,n begin
    v ~ For(1:n) do i
        Normal(speed, 5.0)
    end
end

fullmodel = @model n,t begin
    speed ~ Normal(0.0, 100.0)
    x ~ gps(;speed, n, t)
    v ~ radar(;speed, n)
end

n = 10
t = 1:10

truth = rand(fullmodel(;n, t))
v = truth.v

post = fullmodel(;n, t) |(v=truth.v,)

dynamicHMC(post)

logdensity(post, truth)

speed = truth.speed
x = truth.x
v = truth.v

_ℓ = 0.0
_ℓ += logdensity(Normal(0.0, 100.0), speed)
speed = predict(Normal(0.0, 100.0), speed)
_ℓ += logdensity(gps(; speed, n, t), x)
x = predict(gps(; speed, n, t), x)
_ℓ += logdensity(radar(; speed, n), v)
v = predict(radar(; speed, n), v)
_ℓ


text:
################################################################################
#                   CORE DISPATCHVARIABLEREF METHOD EXTENSIONS
################################################################################
# Extend dispatch_variable_ref
function dispatch_variable_ref(
    model::InfiniteModel,
    index::DerivativeIndex
    )::DerivativeRef
    return DerivativeRef(model, index)
end

# Extend _add_data_object
function _add_data_object(
    model::InfiniteModel,
    object::VariableData{<:Derivative}
    )::DerivativeIndex
    return MOIUC.add_item(model.derivatives, object)
end

# Extend _data_dictionary (type based)
function _data_dictionary(model::InfiniteModel, ::Type{Derivative})
    return model.derivatives
end

# Extend _data_dictionary (reference based)
function _data_dictionary(dref::DerivativeRef)
    return JuMP.owner_model(dref).derivatives
end

# Extend _data_object
function _data_object(dref::DerivativeRef)
    object = get(_data_dictionary(dref), JuMP.index(dref), nothing)
    object === nothing && error("Invalid derivative reference, cannot find " *
    "corresponding derivative in the model. This is likely " *
    "caused by using the reference of a deleted derivative.")
    return object
end

# Extend _core_variable_object
function _core_variable_object(dref::DerivativeRef)
    return _data_object(dref).variable
end

# Define getter function for deriv.is_vector_start
function _is_vector_start(dref::DerivativeRef)::Bool
    return _core_variable_object(dref).is_vector_start
end

"""
    derivative_argument(dref::DerivativeRef)::GeneralVariableRef

Returns the infinite variable/derivative reference that is the input the
differential operator (i.e., the dependent variable of the derivative).

**Example**
```julia-repl
julia> derivative_argument(dref)
x(t)
```
"""
function derivative_argument(dref::DerivativeRef)::GeneralVariableRef
    return _core_variable_object(dref).variable_ref
end

"""
    operator_parameter(dref::DerivativeRef)::GeneralVariableRef

Returns the infinite parameter reference that is what the differential operator
is operating with respect to (i.e., the independent  variable of the derivative).

**Example**
```julia-repl
julia> operator_parameter(dref)
t
```
"""
function operator_parameter(dref::DerivativeRef)::GeneralVariableRef
    return _core_variable_object(dref).parameter_ref
end

"""
    derivative_method(dref::DerivativeRef)::AbstractDerivativeMethod

Returns the evaluation method employed by `dref` that determines the numerical
computation scheme that will be used to evaluate the derivative. Note that this
is set on by the infinite parameter with respect to which the derivative is
defined.

**Example**
```julia-repl
julia> derivative_method(dref)
FiniteDifference(Backward, true)
```
"""
function derivative_method(dref::DerivativeRef)
    return derivative_method(operator_parameter(dref))
end

# Extend _object_numbers
function _object_numbers(dref::DerivativeRef)::Vector{Int}
    return _object_numbers(derivative_argument(dref))
end

# Extend _parameter_numbers
function _parameter_numbers(dref::DerivativeRef)::Vector{Int}
    return _parameter_numbers(derivative_argument(dref))
end

# Get derivative index for variable-parameter pair if one exists
function _existing_derivative_index(
    vref::GeneralVariableRef,
    pref::GeneralVariableRef
    )::Union{DerivativeIndex, Nothing}
    model = JuMP.owner_model(vref)
    return get(model.deriv_lookup, (vref, pref), nothing)
end

# Get access the derivative constraint indices
function _derivative_constraint_dependencies(
    dref::DerivativeRef
    )::Vector{InfOptConstraintIndex}
    return _data_object(dref).deriv_constr_indices
end

"""
    has_derivative_constraints(dref::DerivativeRef)::Bool

Return a `Bool` whether `dref` has been evaluated within the `InfiniteModel` and
has derivative constraints that have been added to the `InfiniteModel`. Note this
does not indicate if such constraints have been added to the optimizer model. Thus,
with normal usage (i.e., not using `evaluate`) this should always return `false`.
"""
function has_derivative_constraints(dref::DerivativeRef)::Bool
    return !isempty(_derivative_constraint_dependencies(dref))
end

## Set helper methods for adapting data_objects with parametric changes
# No change needed
function _adaptive_data_update(
    dref::DerivativeRef,
    d::D,
    data::VariableData{D}
    )::Nothing where {D <: Derivative}
    data.variable = d
    return
end

# Reconstruction is necessary
function _adaptive_data_update(
    dref::DerivativeRef,
    d::D1,
    data::VariableData{D2}
    )::Nothing  where {D1, D2}
    new_data = VariableData(d, data.name, data.lower_bound_index,
                            data.upper_bound_index, data.fix_index,
                            data.zero_one_index, data.integrality_index,
                            data.measure_indices, data.constraint_indices,
                            data.in_objective, data.point_var_indices,
                            data.semi_infinite_var_indices,
                            data.derivative_indices, data.deriv_constr_indices)
    _data_dictionary(dref)[JuMP.index(dref)] = new_data
    return
end

# Extend _set_core_variable_object for DerivativeRefs
function _set_core_variable_object(
    dref::DerivativeRef,
    d::Derivative
    )::Nothing
    _adaptive_data_update(dref, d, _data_object(dref))
    return
end

################################################################################
#                          DEFINTION METHODS
################################################################################
"""
    build_derivative(_error::Function, info::JuMP.VariableInfo,
                     argument_ref::GeneralVariableRef,
                     parameter_ref::GeneralVariableRef
                     )::Derivative

Constructs and returns a [`Derivative`](@ref) with a differential operator that
depends on `parameter_ref` and operates on `argument_ref`. Variable `info` can also
be provided to associate this derivative with bounds and a starting value function
like that of infinite variables. Errors when `argument_ref` is not an
infinite/semi-infinite variable or derivative that depends on `parameter_ref`.

**Example**
```julia-repl
julia> @infinite_parameter(m, t in [0, 1]); @infinite_variable(m, x(t));

julia> info = VariableInfo(false, 0, false, 0, false, 0, false, 0, false, false);

julia> build_derivative(error, info, x, t)
Derivative{GeneralVariableRef}(VariableInfo{Float64,Float64,Float64,Function}(false, 0.0, false, 0.0, false, 0.0, false, start_func, false, false), true, x(t), t)
````
"""
function build_derivative(
    _error::Function,
    info::JuMP.VariableInfo,
    argument_ref::GeneralVariableRef,
    parameter_ref::GeneralVariableRef
    )
    # check the derivative numerator and denominator
    if !(_index_type(parameter_ref) <: InfiniteParameterIndex)
        _error("Derivatives must be with respect to an infinite parameter, but " *
               "$parameter_ref was given.")
    elseif !(_parameter_number(parameter_ref) in _parameter_numbers(argument_ref))
        _error("Derivative arguments must dependent on the infinite parameter " *
               "it is with respect to. Here $argument_ref does not depend on " *
               "$parameter_ref.")
    elseif isequal(argument_ref, parameter_ref)
        _error("Cannot specify the operator parameter as the derivative argument.")
    end
    # check and format the info correctly
    prefs = raw_parameter_refs(argument_ref)
    new_info, is_vect_func = _check_and_format_infinite_info(_error, info, prefs)
    # make the derivative and return it
    return Derivative(new_info, is_vect_func, argument_ref, parameter_ref)
end

"""
    Deriv{V, P} <: InfOptVariableType

A `DataType` to assist in making derivative variables. This can be passed as an
extra argument to `@variable` to make such a variable:
```julia
@variable(model, var_expr, Deriv(inf_var, inf_par), kwargs...)
```
Here `inf_var` is the infinite variable that is being operated on and `inf_par`
is the infinite parameter that the derivative is defined with respect to.

**Fields**
- `argument::V`: The infinite variable being operated on.
- `operator_parameter::P:` The infinite parameter that determines the derivative.
"""
struct Deriv{V, P} <: InfOptVariableType
    argument::V
    operator_parameter::P
end

"""
    JuMP.build_variable(_error::Function, info::JuMP.VariableInfo,
                        var_type::Deriv)::InfiniteVariable{GeneralVariableRef}

Build and return a first order derivative based on `info` and `var_type`. Errors
if the information in `var_type` is invalid. See [`Deriv`](@ref) for more
information.

**Example**
```julia-repl
julia> info = VariableInfo(false, 0, false, 0, false, 0, true, 0, false, false);

julia> deriv_var = build_variable(error, info, Deriv(y, t));
```
"""
function JuMP.build_variable(
    _error::Function,
    info::JuMP.VariableInfo,
    var_type::Deriv;
    extra_kw_args...
    )
    # Error extra keyword arguments
    for (kwarg, _) in extra_kw_args
        _error("Keyword argument $kwarg is not for use with derivatives.")
    end
    # check for valid inputs
    ivref = var_type.argument
    if !(ivref isa GeneralVariableRef)
        _error("Expected an infinite variable reference dependency ",
               "of type `GeneralVariableRef`, but got an argument of type ",
               "$(typeof(ivref)).")
    end
    pref = var_type.operator_parameter
    if !(pref isa GeneralVariableRef)
        _error("Expected an infinite parameter reference dependency ",
               "of type `GeneralVariableRef`, but got an argument of type ",
               "$(typeof(pref)).")
    end
    return build_derivative(_error, info, ivref, pref)
end

"""
    add_derivative(model::InfiniteModel, d::Derivative,
                   [name::String = ""])::GeneralVariableRef

Adds a derivative `d` to `model` and returns a `GeneralVariableRef` that points
to it. Errors if the derivative dependencies do not belong to `model`. Note that
`d` should be built using [`build_derivative`](@ref) to avoid nuance internal
errors.

**Example**
```julia-repl
julia> @infinite_parameter(m, t in [0, 1]); @variable(m, x, Infinite(t));

julia> info = VariableInfo(false, 0, false, 0, false, 0, false, 0, false, false);

julia> d = build_derivative(error, info, x, t);

julia> dref = add_derivative(m, d)
∂/∂t[x(t)]
```
"""
function add_derivative(
    model::InfiniteModel,
    d::Derivative,
    name::String = ""
    )::GeneralVariableRef
    # check the validity
    vref = dispatch_variable_ref(d.variable_ref)
    pref = dispatch_variable_ref(d.parameter_ref)
    JuMP.check_belongs_to_model(vref, model)
    JuMP.check_belongs_to_model(pref, model)
    # check if we already have a derivative
    existing_index = _existing_derivative_index(d.variable_ref, d.parameter_ref)
    if existing_index === nothing
        # add it to the model and make the reference
        data_object = VariableData(d, name)
        dindex = _add_data_object(model, data_object)
        dref = DerivativeRef(model, dindex)
        # update the mappings
        push!(_derivative_dependencies(vref), dindex)
        push!(_derivative_dependencies(pref), dindex)
        # update the derivative lookup dict
        model.deriv_lookup[(d.variable_ref, d.parameter_ref)] = dindex
        # add the info constraints
        gvref = _make_variable_ref(model, dindex)
        _set_info_constraints(d.info, gvref, dref)
    else
        dref = DerivativeRef(model, existing_index)
        gvref = _make_variable_ref(model, existing_index)
        _update_info_constraints(d.info, gvref, dref)
        if !isempty(name)
            set_name(dref, name)
        end
    end
    return gvref
end

# Extend JuMP.add_variable to enable macro definition
function JuMP.add_variable(
    model::InfiniteModel,
    d::Derivative,
    name::String = ""
    )::GeneralVariableRef
    return add_derivative(model, d, name)
end

## Define function that build derivative expressions following rules of calculus
# GeneralVariableRef
function _build_deriv_expr(
    vref::GeneralVariableRef,
    pref
    )::Union{GeneralVariableRef, Float64}
    if isequal(vref, pref)
        return 1.0
    elseif _parameter_number(pref) in _parameter_numbers(vref)
        dindex = _existing_derivative_index(vref, pref)
        model = JuMP.owner_model(vref)
        if dindex === nothing
            info = VariableInfo(false, NaN, false, NaN, false, NaN, false,
                                s -> NaN, false, false)
            d = Derivative(info, true, vref, pref)
            return add_derivative(model, d)
        else
            return _make_variable_ref(model, dindex)
        end
    else
        return 0.0
    end
end

# AffExpr
function _build_deriv_expr(aff::JuMP.GenericAffExpr, pref)
    return _MA.@rewrite(sum(c * _build_deriv_expr(v, pref)
                        for (c, v) in JuMP.linear_terms(aff)))
end

# Quad Expr (implements product rule)
function _build_deriv_expr(quad::JuMP.GenericQuadExpr, pref)
    return _MA.@rewrite(sum(c * (_build_deriv_expr(v1, pref) * v2 +
                        v1 * _build_deriv_expr(v2, pref))
                        for (c, v1, v2) in JuMP.quad_terms(quad)) +
                        _build_deriv_expr(quad.aff, pref))
end

# Real number
function _build_deriv_expr(expr::Real, pref)::Float64
    return 0.0
end

# Fallback
function _build_deriv_expr(expr, pref)
    error("Unsupported expression type $(typeof(expr)) for derivative definition.")
end

# Define recursive function to call _build_deriv_expr
function _recursive_deriv_build(expr, prefs
    )::Union{JuMP.AbstractJuMPScalar, Float64}
    if isempty(prefs)
        return expr
    else
        new_expr = _build_deriv_expr(expr, first(prefs))
        return _recursive_deriv_build(new_expr, prefs[2:end])
    end
end

"""
    deriv(expr::JuMP.AbstractJuMPScalar, pref1::GeneralVariableRef[, ....]
          )::Union{JuMP.AbstractJuMPScalar, Float64}

Apply appropriate calculus methods to define and return the derivative expression of `expr`
with respect to the infinite parameter(s) `pref1`, pref2`, etc. in that respective
order. This will implicilty build and add individual [`Derivative`](@ref)s as
appropriate. Errors if no infinite parameter is given or if the parameters are
not infinite.

**Example**
```julia-repl
julia> @infinite_parameter(m, t in [0, 1])
t

julia> @variable(m, x, Infinite(t))
x(t)

julia> @variable(m, z)
z

julia> deriv_expr = deriv(x^2 + z, t, t)
2 ∂/∂t[∂/∂t[x(t)]]*x(t) + 2 ∂/∂t[x(t)]²
```
"""
function deriv(expr, prefs::GeneralVariableRef...)
    # Check inputs
    if !all(_index_type(pref) <: InfiniteParameterIndex for pref in prefs)
        error("Can only take derivative with respect to infinite parameters.")
    elseif length(prefs) == 0
        error("Must specify at least one infinite parameter with which to define " *
              "the derivative operator.")
    end
    # Build the derivative expression
    return _recursive_deriv_build(expr, prefs)
end

"""
    @deriv(expr, pref_expr1[, ...]
           )::Union{JuMP.AbstractJuMPScalar, Float64}

The macro variant of [`deriv`](@ref) that is more efficient for expression building
and enables symbolic differential operator parameter defintions via `pref_expr`s.
Like `deriv` expr can be any InfiniteOpt expression and the appropriate calculus
rules will applied to `expr` to take its derivative with respect to the indicated
infinite parameters detailed by the `pref_expr`s. The resulting derivative
expression will contain individual derivatives that were created and added to the
InfiniteModel as needed. Here each `pref_expr` arugment can be of the form:
- `pref::GeneralVariableRef`: An indiviudal infinite parameter reference
- `(pref::GeneralVariableRef)^(p::Int)`: An infinite parameter applied `p` times.
Thus, the syntax `@deriv(expr, pref^2)` is equivalent to `@deriv(expr, pref, pref)`.

This will error if `pref_expr` is an unrecongnized syntax, no infinite parameter
is given, or if any of the specified parameters are not infinite.

**Example**
```julia-repl
julia> @infinite_parameter(m, t in [0, 1])
t

julia> @variable(m, x, Infinite(t))
x(t)

julia> @variable(m, z)
z

julia> deriv_expr = @deriv(x^2 + z, t^2)
2 ∂/∂t[∂/∂t[x(t)]]*x(t) + 2 ∂/∂t[x(t)]²
```
"""
macro deriv(expr, args...)
    # process the arugments
    extra, kwargs, _, _ = _extract_kwargs(args)
    # error if kwargs are given
    _error(str...) = _macro_error(:deriv, (expr, args...), __source__, str...)
    isempty(kwargs) || _error("Invalid keyword argument given.")
    # expand the parameter references as needed with powers
    pref_exprs = []
    for p in extra
        if isexpr(p, :call) && p.args[1] == :^
            append!(pref_exprs, repeat([esc(p.args[2])], p.args[3]))
        else
            push!(pref_exprs, esc(p))
        end
    end
    # prepare the code to call deriv
    expression = _MA.rewrite_and_return(expr)
    return :( deriv($expression, $(pref_exprs...)) )
end

"""
    ∂(expr::JuMP.AbstractJuMPScalar, pref1::GeneralVariableRef[, ....]
      )::Union{JuMP.AbstractJuMPScalar, Float64}

This serves as a convenient unicode wrapper for [`deriv`](@ref). The `∂` is
produced via `\\partial`.
"""
function ∂(
    expr,
    prefs::GeneralVariableRef...
    )::Union{JuMP.AbstractJuMPScalar, Float64}
    return deriv(expr, prefs...)
end

"""
    @∂(expr, pref_expr1[, ...])::Union{JuMP.AbstractJuMPScalar, Float64}

This serves as a convenient unicode wrapper for [`@deriv`](@ref). The `∂` is
produced via `\\partial`.
"""
macro ∂(expr, args...)
    return esc(:( @deriv($expr, $(args...)) ))
end

################################################################################
#                           PARAMETER REFERENCES
################################################################################
"""
    raw_parameter_refs(dref::DerivativeRef)::VectorTuple

Return the raw [`VectorTuple`](@ref InfiniteOpt.Collections.VectorTuple) of the
parameter references that `dref` depends on. This is primarily an internal method
where [`parameter_refs`](@ref parameter_refs(::DerivativeRef))
is intended as the preferred user function.
"""
function raw_parameter_refs(dref::DerivativeRef)
    return raw_parameter_refs(derivative_argument(dref))
end

"""
    parameter_refs(dref::DerivativeRef)::Tuple

Return the parameter references associated with the infinite derivative `dref`. This
is formatted as a Tuple of containing the parameter references as they inputted
to define `dref`.

**Example**
```julia-repl
julia> parameter_refs(deriv)
(t,)
```
"""
function parameter_refs(dref::DerivativeRef)
    return parameter_refs(derivative_argument(dref))
end

"""
    parameter_list(dref::DerivativeRef)::Vector{GeneralVariableRef}

Return a vector of the parameter references that `dref` depends on. This is
primarily an internal method where [`parameter_refs`](@ref parameter_refs(::DerivativeRef))
is intended as the preferred user function.
"""
function parameter_list(dref::DerivativeRef)::Vector{GeneralVariableRef}
    return raw_parameter_refs(dref).values
end

################################################################################
#                           VARIABLE INFO METHODS
################################################################################
# Set info for infinite derivatives
function _update_variable_info(
    dref::DerivativeRef,
    info::JuMP.VariableInfo
    )::Nothing
    vref = derivative_argument(dref)
    pref = operator_parameter(dref)
    is_vect_func = _is_vector_start(dref)
    new_deriv = Derivative(_format_infinite_info(info), is_vect_func, vref, pref)
    _set_core_variable_object(dref, new_deriv)
    return
end

"""
    set_start_value_function(dref::DerivativeRef,
                             start::Union{Real, Function})::Nothing

Set the start value function of `dref`. If `start::Real` then a function is
generated to such that the start value will be `start` for the entire infinite
domain. If `start::Function` then this function should map to a scalar start value
given a support value arguments matching the format of the parameter elements in
`parameter_refs(dref)`.

**Example**
```julia-repl
julia> set_start_value_function(dref, 1) # all start values will be 1

julia> set_start_value_function(dref, my_func) # each value will be made via my_func
```
"""
function set_start_value_function(
    dref::DerivativeRef,
    start::Union{Real, Function}
    )::Nothing
    info = _variable_info(dref)
    set_optimizer_model_ready(JuMP.owner_model(dref), false)
    prefs = raw_parameter_refs(dref)
    temp_info = JuMP.VariableInfo(info.has_lb, info.lower_bound, info.has_ub,
                                 info.upper_bound, info.has_fix, info.fixed_value,
                                 true, start, info.binary, info.integer)
    new_info, is_vect_func = _check_and_format_infinite_info(error, temp_info, prefs)
    vref = derivative_argument(dref)
    pref = operator_parameter(dref)
    new_deriv = Derivative(new_info, is_vect_func, vref, pref)
    _set_core_variable_object(dref, new_deriv)
    return
end

"""
    reset_start_value_function(dref::DerivativeRef)::Nothing

Remove the existing start value function and return to the default. Generally,
this is triggered by deleting an infinite parameter that `dref` depends on.

**Example**
```julia-repl
julia> reset_start_value_function(dref)
```
"""
function reset_start_value_function(dref::DerivativeRef)::Nothing
    info = _variable_info(dref)
    set_optimizer_model_ready(JuMP.owner_model(dref), false)
    new_info = JuMP.VariableInfo(info.has_lb, info.lower_bound, info.has_ub,
                                 info.upper_bound, info.has_fix, info.fixed_value,
                                 false, s -> NaN, info.binary, info.integer)
    vref = derivative_argument(dref)
    pref = operator_parameter(dref)
    new_deriv = Derivative(new_info, true, vref, pref)
    _set_core_variable_object(dref, new_deriv)
    return
end

# TODO maybe throw errors for binary and integer?

################################################################################
#                              MODEL QUERIES
################################################################################
"""
    num_derivatives(model::InfiniteModel)::Int

Returns the number of derivatives that have been defined in `model`. Note that
nested derivatives will be counted in accordance with their components (e.g.,
``\\frac{d^2 x(t)}{dt^2} = ``\\frac{d}{dt}\\left(\\frac{d x(t)}{dt} \\right)``
will count as 2 derivatives.)

**Example**
```julia-repl
julia> num_derivatives(model)
12
```
"""
function num_derivatives(model::InfiniteModel)::Int
    return length(_data_dictionary(model, Derivative))
end

"""
    all_derivatives(model::InfiniteModel)::Vector{GeneralVariableRef}

Returns a list of all the individual derivatives stored in `model`.

**Example**
```julia-repl
julia> all_derivatives(model)
3-element Array{GeneralVariableRef,1}:
 ∂/∂t[T(x, t)]
 ∂/∂x[T(x, t)]
 ∂/∂x[∂/∂x[T(x, t)]]
```
"""
function all_derivatives(model::InfiniteModel)::Vector{GeneralVariableRef}
    vrefs_list = Vector{GeneralVariableRef}(undef, num_derivatives(model))
    for (i, (index, _)) in enumerate(_data_dictionary(model, Derivative))
        vrefs_list[i] = _make_variable_ref(model, index)
    end
    return vrefs_list
end

################################################################################
#                      EVALUATION METHOD MODIFICATIONS/QUERIES
################################################################################
# Fallback (needs to be done via the infinite parameter)
function set_derivative_method(::DerivativeRef, ::AbstractDerivativeMethod)
    error("Cannot specify the derivative method in terms of the derivative, must " *
          " do so in terms of the operator parameter by calling `set_derivative_method` " *
          " using the operator parameter.")
end

"""
    derivative_constraints(dref::DerivativeRef)::Vector{InfOptConstraintRef}

Return a list of the derivative evaluation constraints for `dref` that have been
added directly to the `InfiniteModel` associated with `dref`. An empty vector is
returned is there are no such constraints.
"""
function derivative_constraints(dref::DerivativeRef)::Vector{InfOptConstraintRef}
    return [_make_constraint_ref(JuMP.owner_model(dref), idx)
            for idx in _derivative_constraint_dependencies(dref)]
end

"""
    delete_derivative_constraints(dref::DerivativeRef)::Nothing

Delete any derivative constraints of `dref` that have been directly added to the
`InfiniteModel`.
"""
function delete_derivative_constraints(dref::DerivativeRef)::Nothing
    model = JuMP.owner_model(dref)
    for cref in derivative_constraints(dref)
        JuMP.delete(model, cref)
    end
    empty!(_derivative_constraint_dependencies(dref))
    pref = dispatch_variable_ref(operator_parameter(dref))
    if all(!has_derivative_constraints(DerivativeRef(model, idx))
            for idx in _derivative_dependencies(pref))
        _set_has_derivative_constraints(pref, false)
    end
    return
end

################################################################################
#                                 DELETION
################################################################################
# Extend _delete_variable_dependencies (for use with JuMP.delete)
function _delete_variable_dependencies(dref::DerivativeRef)::Nothing
    # remove variable info constraints associated with dref
    _delete_info_constraints(dref)
    # update variable and parameter mapping
    vref = derivative_argument(dref)
    pref = operator_parameter(dref)
    filter!(e -> e != JuMP.index(dref), _derivative_dependencies(vref))
    filter!(e -> e != JuMP.index(dref), _derivative_dependencies(pref))
    model = JuMP.owner_model(dref)
    # delete the variable-parameter pair from the anti-duplication dictionary
    delete!(model.deriv_lookup, (vref, pref))
    # delete any derivative constraints associated with this derivative
    delete_derivative_constraints(dref)
    # delete associated point variables and mapping
    for index in _point_variable_dependencies(dref)
        JuMP.delete(model, dispatch_variable_ref(model, index))
    end
    # delete associated semi-infinite variables and mapping
    for index in _semi_infinite_variable_dependencies(dref)
        JuMP.delete(model, dispatch_variable_ref(model, index))
    end
    # delete associated derivative variables and mapping
    for index in _derivative_dependencies(dref)
        JuMP.delete(model, dispatch_variable_ref(model, index))
    end
    return
end


text:
using StatsFuns
using Base.Test
import StatsFuns.RFunctions
using Compat

function check_rmath(fname, statsfun, rmathfun, params, aname, a, isprob, rtol)
    v = statsfun(params..., a)
    rv = rmathfun(params..., a)
    if isprob
        rd = abs(v / rv - 1.0)
        if rd > rtol
            error("$fname deviates too much from Rmath at " *
                "params = $params, $aname = $a:\n" *
                "  v = $v (rv = $rv)\n  |v/rv - 1| = $rd > $rtol.")
        end
    else
        τ = (1.0 + abs(rv)) * rtol
        ad = abs(v - rv)
        if ad > τ
            error("$fname deviates too much from Rmath at " *
                "params = $params, $aname = $a:\n" *
                "  v = $v (rv = $rv)\n  |v - rv| = $ad > $τ.")
        end
    end
end

get_statsfun(fname) = eval(Symbol(fname))
get_rmathfun(fname) = eval(parse(string("RFunctions.", fname)))

function rmathcomp(basename, params, X::AbstractArray, rtol=100eps(float(one(eltype(X)))))
    pdf     = string(basename, "pdf")
    logpdf  = string(basename, "logpdf")
    cdf     = string(basename, "cdf")
    ccdf    = string(basename, "ccdf")
    logcdf  = string(basename, "logcdf")
    logccdf = string(basename, "logccdf")
    invcdf     = string(basename, "invcdf")
    invccdf    = string(basename, "invccdf")
    invlogcdf  = string(basename, "invlogcdf")
    invlogccdf = string(basename, "invlogccdf")
    rand       = string(basename, "rand")

    stats_pdf     = get_statsfun(pdf)
    stats_logpdf  = get_statsfun(logpdf)
    stats_cdf     = get_statsfun(cdf)
    stats_ccdf    = get_statsfun(ccdf)
    stats_logcdf  = get_statsfun(logcdf)
    stats_logccdf = get_statsfun(logccdf)
    stats_invcdf     = get_statsfun(invcdf)
    stats_invccdf    = get_statsfun(invccdf)
    stats_invlogcdf  = get_statsfun(invlogcdf)
    stats_invlogccdf = get_statsfun(invlogccdf)

    rmath_pdf     = get_rmathfun(pdf)
    rmath_logpdf  = get_rmathfun(logpdf)
    rmath_cdf     = get_rmathfun(cdf)
    rmath_ccdf    = get_rmathfun(ccdf)
    rmath_logcdf  = get_rmathfun(logcdf)
    rmath_logccdf = get_rmathfun(logccdf)
    rmath_invcdf     = get_rmathfun(invcdf)
    rmath_invccdf    = get_rmathfun(invccdf)
    rmath_invlogcdf  = get_rmathfun(invlogcdf)
    rmath_invlogccdf = get_rmathfun(invlogccdf)

    # tackle rand specially
    has_rand = true
    if basename == "nbeta" || basename == "nfdist" || basename == "ntdist"
        has_rand = false
    end
    rmath_rand = has_rand ? get_rmathfun(rand) : nothing

    for i = 1:length(X)
        x = X[i]
        check_rmath(pdf, stats_pdf, rmath_pdf,
            params, "x", x, true, rtol)
        check_rmath(logpdf, stats_logpdf, rmath_logpdf,
            params, "x", x, false, rtol)
        check_rmath(cdf, stats_cdf, rmath_cdf,
            params, "x", x, true, rtol)
        check_rmath(ccdf, stats_ccdf, rmath_ccdf,
            params, "x", x, true, rtol)
        check_rmath(logcdf, stats_logcdf, rmath_logcdf,
            params, "x", x, false, rtol)
        check_rmath(logccdf, stats_logccdf, rmath_logccdf,
            params, "x", x, false, rtol)

        p = rmath_cdf(params..., x)
        cp = rmath_ccdf(params..., x)
        lp = rmath_logcdf(params..., x)
        lcp = rmath_logccdf(params..., x)

        check_rmath(invcdf, stats_invcdf, rmath_invcdf,
            params, "q", p, false, rtol)
        check_rmath(invccdf, stats_invccdf, rmath_invccdf,
            params, "q", cp, false, rtol)
        check_rmath(invlogcdf, stats_invlogcdf, rmath_invlogcdf,
            params, "lq", lp, false, rtol)
        check_rmath(invlogccdf, stats_invlogccdf, rmath_invlogccdf,
            params, "lq", lcp, false, rtol)

        # make sure that rand works
        if has_rand
            rmath_rand(params...)
        end
    end
end

function rmathcomp_tests(basename, configs)
    println("\ttesting $basename ...")
    for (params, data) in configs
        rmathcomp(basename, params, data)
    end
end

### Test cases

rmathcomp_tests("beta", [
    ((1.0, 1.0), 0.01:0.01:0.99),
    ((2.0, 3.0), 0.01:0.01:0.99),
    ((10.0, 2.0), 0.01:0.01:0.99),
])

rmathcomp_tests("binom", [
    ((1, 0.5), 0:1),
    ((1, 0.7), 0:1),
    ((8, 0.6), 0:8),
    ((20, 0.1), 0:20),
    ((20, 0.9), 0:20),
])

rmathcomp_tests("chisq", [
    ((1,), 0.0:0.1:8.0),
    ((4,), 0.0:0.1:8.0),
    ((9,), 0.0:0.1:8.0),
])

rmathcomp_tests("fdist", [
    ((1, 1), (0.0:0.1:5.0)),
    ((2, 1), (0.0:0.1:5.0)),
    ((5, 2), (0.0:0.1:5.0)),
    ((10, 1), (0.0:0.1:5.0)),
    ((10, 3), (0.0:0.1:5.0)),
])

rmathcomp_tests("gamma", [
    ((1.0, 1.0), (0.05:0.05:12.0)),
    ((0.5, 1.0), (0.05:0.05:12.0)),
    ((3.0, 1.0), (0.05:0.05:12.0)),
    ((9.0, 1.0), (0.05:0.05:12.0)),
    ((2.0, 3.0), (0.05:0.05:12.0)),
])

rmathcomp_tests("hyper", [
    ((2, 3, 4), 0:4)
])

rmathcomp_tests("nbeta", [
    ((1.0, 1.0, 0.0), 0.01:0.01:0.99),
    ((2.0, 3.0, 0.0), 0.01:0.01:0.99),
    ((1.0, 1.0, 2.0), 0.01:0.01:0.99),
    ((3.0, 4.0, 2.0), 0.01:0.01:0.99),
])

rmathcomp_tests("nbinom", [
    ((1, 0.5), 0:20),
    ((3, 0.5), 0:20),
    ((3, 0.2), 0:20),
    ((3, 0.8), 0:20),
])

rmathcomp_tests("nchisq", [
    ((2, 1), 0.0:0.2:8.0),
    ((2, 3), 0.0:0.2:8.0),
    ((4, 1), 0.0:0.2:8.0),
    ((4, 3), 0.0:0.2:8.0),
])

rmathcomp_tests("nfdist", [
    ((1.0, 1.0, 0.0), 0.1:0.1:10.0),
    ((1.0, 1.0, 2.0), 0.1:0.1:10.0),
    ((2.0, 3.0, 1.0), 0.1:0.1:10.0),
])

rmathcomp_tests("norm", [
    ((0.0, 1.0), -6.0:0.01:6.0),
    ((2.0, 1.0), -3.0:0.01:7.0),
    ((0.0, 0.5), -3.0:0.01:3.0),
    ((0, 1), -3:3),
    ((0, 1), -Float16(3):Float16(0.01):Float16(3))
])

rmathcomp_tests("ntdist", [
    ((0, 1), -4.0:0.1:10.0),
    ((0, 4), -4.0:0.1:10.0),
    ((2, 1), -4.0:0.1:10.0),
    ((2, 4), -4.0:0.1:10.0),
])

rmathcomp_tests("pois", [
    ((0.5,), 0:20),
    ((1.0,), 0:20),
    ((2.0,), 0:20),
    ((10.0,), 0:20),
])

rmathcomp_tests("tdist", [
    ((1,), -5.0:0.1:5.0),
    ((2,), -5.0:0.1:5.0),
    ((5,), -5.0:0.1:5.0),
])


text:
model_ridge_id = Sem(
    specification = spec,
    data = dat,
    loss = (SemML, SemRidge),
    α_ridge = .001,
    which_ridge = [:x16, :x17, :x18, :x19, :x20],
    diff = semdiff
)


text:
export TestCase, runtestcase

"""
    TestCase(args[, kwargs[, tolerance]])
    TestCase(args...; kwargs...[, _tolerance])

Construct a test case.

The `args` and `kwargs` are passed to the function that will be tested.
The `tolerance` is the allowable difference between the result of the testcase
and the desired result (_e.g._ the result of the instructor function).

Note that the `tolerance` is only used when the result of the test case is
something that has a natural real valued difference between two instances; for
example the difference between two real vectors is the infinity (max) norm.
See the `compare` function for more details.

# Examples
```julia-repl
julia> TestCase([1, 2, 3])
TestCase(Any[1, 2, 3], Dict{Symbol,Any}(), 1.0e-5)

julia> TestCase([1, 2, 3], Dict(:foo => "bar"))
TestCase(Any[1, 2, 3], Dict{Symbol,Any}(:foo=>"bar"), 1.0e-5)

julia> TestCase(1, 2, 3; foo="bar", _tolerance=1e-9)
TestCase(Any[1, 2, 3], Dict{Symbol,Any}(:foo=>"bar"), 1.0e-9)
```
"""
struct TestCase
    args::Array{Any}
    kwargs::Dict{Symbol, Any}

    # Maximum numeric error tolerance.
    # This is calculated according to calculate_error.
    tolerance::Float64

    # This needs to be an inner constructor to account for when the varargs
    # constructor is used with three arguments.
    TestCase(args::Array, kwargs::Dict=Dict(), tolerance::Float64=DEFAULT_TOLERANCE) = new(
        convert(Array{Any}, args),
        convert(Dict{Symbol, Any}, kwargs),
        tolerance,
    )
end

TestCase(args...; tolerance=DEFAULT_TOLERANCE, kwargs...) = TestCase(
    Array{Any}([args...]),
    convert(Dict{Symbol, Any}, kwargs),
    tolerance,
)

struct TestCaseResult
    result::Union{Any, Nothing}
    # Note: exception is a string to avoid breaking (de)serialize if the
    # exeption happens to reference context-specific data (eg. something in the
    # _Submission module which wouldn't be present in the parent process)
    exception::Union{String, Nothing}
    backtrace::Union{String, Nothing}
    time::Float64
    stdout::String
    stderr::String
end

"""
    runtestcase(f::Function, testcase::TestCase)::TestCaseResult

Run the given test case on the specified function.
"""
function runtestcase(f::Function, testcase::TestCase)::TestCaseResult
    @debug("Running testcase for function $(repr(f)).")
    result, exception, backtrace = nothing, nothing, nothing

    start_ns = Base.time_ns()
    sub_stdout, sub_stderr = @capture_stdstreams begin
        try
            result = Base.invokelatest(f, testcase.args...; testcase.kwargs...)
        catch e
            exception = repr(e)
            backtrace = sprint(showerror, e)
        end
    end
    end_ns = Base.time_ns()

    # Get the elapsed time in seconds (1e-9 converts from nanoseconds)
    elapsedtime = (end_ns - start_ns) * 1e-9
    @debug "Finished running testcase." result exception elapsedtime

    return TestCaseResult(
        result, exception, backtrace,
        elapsedtime, sub_stdout, sub_stderr,
    )
end


text:
using Gurobi
using Mosek

# train adversarial method with kernel using constraint generation
function train_adv_kernel_cg(X::Matrix, y::Vector, C::Real=1.0,
  kernel::Symbol=:linear, kernel_params::Vector=[];
  perturb::Real=0.0, tol::Real=1e-6, psdtol::Real=1e-6, obj_reltol::Real=0.0,
  log::Real=0, n_thread::Int=0, solver::Symbol=:gurobi, verbose::Bool=true)

  n = length(y)
  # add one
  X1 = [ones(n) X]'   # transpose
  m = size(X1, 1)

  # number of class
  n_c = maximum(y)
  n_f = n_c * m   # number of features

  alpha = zeros(0)

  # array of tuple
  constraints = Tuple{Integer, Vector}[]

  # prepare saved vars
  idmi = map(i -> idi(m, i), collect(1:n_c))

  # kernel
  kernel_func = linear_kernel
  if kernel == :gaussian
    kernel_func = gaussian_kernel
  elseif kernel == :polynomial
    kernel_func = polynomial_kernel
  end

  # precompute kernel
  K = [( i >= j ? kernel_func(X1[:,i], X1[:,j], kernel_params...) : 0.0 )::Float64 for i=1:n, j=1:n]
  K = [( i >= j ? K[i,j] : K[j,i] )::Float64 for i=1:n, j=1:n]

  if solver == :gurobi
    # gurobi solver
    # gurobi environtment
    env = Gurobi.Env()
    # Method : 0=primal simplex, 1=dual simplex, 2=barrier  ; default for QP: barrier
    # Threads : default = 0 (use all threads)
    setparams!(env, PSDTol=psdtol, LogToConsole=log, Method=2, Threads=n_thread)
  elseif solver == :mosek
    # mosek environtment
    env = makeenv()
  end

  # params for Gurobi
  Q = zeros(0,0)
  nu = zeros(0)
  A = spzeros(n, 0)        # sparse matrix
  b = ones(n) * C

  Q_prev = zeros(0,0)
  nu_prev = zeros(0)
  A_prev = spzeros(n, 0)  # sparse matrix

  # save lambda dot psi_i
  n_lps = n * n_c
  LPsi = zeros(n_lps, 0)
  LPsi_prev = zeros(n_lps, 0)

  iter = 0
  dual_obj = 0.0
  dual_obj_prev = 0.0

  while true
    iter += 1

    if verbose
      println("Iteration : ", iter)
      tic();
    end

    # previous constraints
    const_prev = copy(constraints)
    # previous alpha
    alpha_prev = copy(alpha)

    ## add to constraints
    const_added = Tuple{Integer, Vector}[]

    # find constraint for each sample
    for i=1:n
      psis = psi_list_dual(alpha, LPsi, i, n_c)
      psis_id, val = best_psis(psis)      # most violated constraints

      # current xi_i
      id_i = find(x -> x[1] == i, constraints)
      xi_i_list =  map(x -> x[2], constraints[id_i])
      max_xi_i = 0
      for j = 1:length(xi_i_list)
        a = calc_const(psis::Vector, xi_i_list[j])
        if a > max_xi_i
          max_xi_i = a
        end
      end

      if val > max_xi_i
        cs = (i, sort!(psis_id))
        if findfirst(constraints .== cs) == 0
          push!(constraints, cs)
          push!(const_added, cs)
        end
      end
    end

    # if no constraints added
    if length(const_added) == 0
      break
    end

    n_const = length(constraints)

    #### Start QP ###

    if verbose
      println(">> Start QP")
      toc();
      tic();
    end

    n_prev = length(const_prev)
    n_added = length(const_added)

    # init alpha with previous iteration
    alpha = zeros(n_const)
    alpha[1:n_prev] = alpha_prev

    Q_aug = [ ( calc_dot((const_prev[i][1], const_prev[i][2], const_added[j][1], const_added[j][2]), K, y) )::Float64
          for i=1:n_prev, j=1:n_added]

    Q_aug_diag = [
          ( i >= j ? calc_dot((const_added[i][1], const_added[i][2], const_added[j][1], const_added[j][2]), K, y) : 0.0 )::Float64
          for i=1:n_added, j=1:n_added]
    Q_aug_diag = [ (i >= j ? Q_aug_diag[i,j] : Q_aug_diag[j,i] )::Float64 for i=1:n_added, j=1:n_added ]

    Q = [
          ( (i <= n_prev && j <= n_prev) ? Q_prev[i, j] : ( i <= n_prev ? Q_aug[i, j-n_prev] :
          ( j <= n_prev ? Q_aug[j, i-n_prev] : Q_aug_diag[i-n_prev, j-n_prev] ) ) )::Float64
          for i=1:n_const, j=1:n_const]

    nu_aug = [ ( calc_cconst(const_added[i][2]) )::Float64 for i=1:n_added ]
    nu = [ ( i <= n_prev ? nu_prev[i] : nu_aug[i-n_prev] )::Float64 for i=1:n_const]

    A_aug = spzeros(n, n_added)
    for j = 1:n_added
      A_aug[const_added[j][1], j] = 1.0
    end
    A = [A_prev A_aug]

    # update LPsi
    LPsi_aug = [ ( calc_dot((ceil(Int64, i/n_c), [(i%n_c==0)? n_c:i%n_c], const_added[j][1], const_added[j][2]), K, y) )::Float64
          for i=1:n_lps, j=1:n_added]
    LPsi = [LPsi_prev LPsi_aug]

    ## add perturbation
    for i=1:n_const
      Q[i,i] = Q[i,i] + perturb
    end

    if verbose
      toc();
      tic();
    end

    if solver == :gurobi

      if verbose println(">> Optim :: Gurobi") end

      ## init model
      model = gurobi_model(env,
                 sense = :minimize,
                 H = Q,
                 f = -nu,
                 A = A,
                 b = b,
                 lb = zeros(n_const)
                 )
      # Print the model to check correctness
      # print(model)

      # Solve with Gurobi
      Gurobi.optimize(model)

      if verbose
        toc();
        println("<< End QP")
      end

      dual_obj = -get_objval(model)
      # Solution
      if verbose println("Objective value: ", dual_obj) end

      # get alpha
      alpha = get_solution(model)
      # println("alpha = ", alpha)
      #
      # println("constraints = ", constraints)

      if verbose println("n constraints = ", length(constraints)) end

      ### end QP ###

    elseif solver == :mosek

      if verbose println(">> Optim :: Mosek") end

      task = maketask(env)

      # set params
      putintparam(task, Mosek.MSK_IPAR_LOG, 1)
      putintparam(task, Mosek.MSK_IPAR_LOG_CHECK_CONVEXITY, 1)
      putdouparam(task, Mosek.MSK_DPAR_CHECK_CONVEXITY_REL_TOL, psdtol)

      # variables
      appendvars(task, n_const)
      # bound on var
      for i::Int32 = 1:n_const
        putbound(task, Mosek.MSK_ACC_VAR, i, Mosek.MSK_BK_RA, 0.0, C)
      end

      # objective
      putobjsense(task, Mosek.MSK_OBJECTIVE_SENSE_MINIMIZE)
      qi = zeros(Int32, (n_const * (n_const+1)) ÷ 2 )
      qj = zeros(Int32, (n_const * (n_const+1)) ÷ 2 )
      qv = zeros(Float64, (n_const * (n_const+1)) ÷ 2 )
      ix = 1
      for j::Int32 = 1:n_const
        for i::Int32 = j:n_const
          qi[ix] = i
          qj[ix] = j
          qv[ix] = Q[i,j]
          ix += 1
        end
      end
      putqobj(task, qi, qj, qv)

      putclist(task, collect(1:n_const), -nu)

      # constraints
      ## sparse array
      appendcons(task, n)
      for i::Int32 = 1:n_const
        id_nz = A[:,i].nzind
        putacol(task, i, id_nz, ones(length(id_nz)))
      end

      for i::Int32 = 1:n
        putbound(task, Mosek.MSK_ACC_CON, i, Mosek.MSK_BK_RA, 0.0, C)
      end

      if verbose
        toc(); tic();
      end

      Mosek.optimize(task)

      if verbose
        toc();
        println("<< End QP")
      end

      # Solution
      dual_obj, _ = getsolutioninfo(task, Mosek.MSK_SOL_ITR)
      if verbose println("Objective value: ", -dual_obj) end

      # get alpha
      alpha = getxx(task, Mosek.MSK_SOL_ITR)

      if verbose println("n constraints = ", length(constraints)) end
    end

    if obj_reltol > 0.0
      if (dual_obj - dual_obj_prev) / dual_obj_prev < obj_reltol && iter > 1
        if verbose
          println((dual_obj - dual_obj_prev) / dual_obj_prev)
          println(">> Iteration STOPPED | Objective relative tolerance : ", obj_reltol)
        end
        break
      end
    end

    if verbose println() end

    Q_prev = Q
    nu_prev = nu
    A_prev = A
    LPsi_prev = LPsi
    dual_obj_prev = dual_obj
  end

  # finalizing losses
  gv_aug = zeros(n)
  gv_01 = zeros(n)
  l_adv = zeros(n)
  l_01 = zeros(n)
  for i=1:n
    psis = psi_list_dual(alpha, LPsi, i, n_c)
    psis_id, val = best_psis(psis)      # most violated constraints
    n_ps = length(psis_id)

    gv_aug[i] = val

    # compute probs
    p_hat = zeros(n_c)
    p_check = zeros(n_c)
    for j=1:n_c
      if j in psis_id
        p_hat[j] = ( (n_ps-1.0)*psis[j] - sum(psis[psis_id[psis_id .!= j]]) + 1.0 ) / n_ps
        p_check[j] = 1.0 / n_ps
      else
        p_hat[j] = 0.0
        p_check[j] = 0.0
      end
    end

    C01 = 1 - eye(n_c)    # 01 loss matrix
    v = p_hat' * C01 * p_check   # the result is vector size 1 not a number
    gv_01[i] = v[1]     #

    ## training loss
    l_adv[i] = 1.0 - p_hat[y[i]]
    l_01[i] = 1.0 - round(Int, indmax(p_hat) == y[i])
  end

  game_value_01 = mean(gv_01)
  game_value_augmented = mean(gv_aug)

  # create model
  adv_model = KernelMultiAdversarialModel(kernel, kernel_params, alpha, constraints, n_c, game_value_01, game_value_augmented, mean(l_adv), mean(l_01))

  return adv_model::KernelMultiAdversarialModel
end

function predict_adv_kernel(model::KernelMultiAdversarialModel, X_test::Matrix, X_train::Matrix, y_train::Vector)

  alpha = model.alpha
  n_c = model.n_class
  n = size(X_test, 1)

  constraints = model.constraints
  n_const = length(alpha)

  X1 = [ones(n) X_test]'   # transpose
  m = size(X1, 1)

  # training data
  n_tr = size(X_train, 1)
  X1_tr = [ones(n_tr) X_train]'   # transpose

  # kernel
  kernel = model.kernel
  kernel_params = model.kernel_params
  # kernel function
  kernel_func = linear_kernel
  if kernel == :gaussian
    kernel_func = gaussian_kernel
  elseif kernel == :polynomial
    kernel_func = polynomial_kernel
  end

  # compute Kernel
  K = [ kernel_func(X1_tr[:,i], X1[:,j], kernel_params...)::Float64 for i=1:n_tr, j=1:n]

  prob = zeros(n, n_c)
  pred = zeros(n)
  for i=1:n

    psis = zeros(n_c)
    for j = 1:n_const
      psis -= alpha[j] * calc_dotlphi( (constraints[j][1], constraints[j][2], i), K, y_train, n_c)
    end

    psis_id, val = best_psis(psis)      # most violated constraints
    n_ps = length(psis_id)

    for j=1:n_c
      if j in psis_id
        prob[i,j] = ( (n_ps-1.0)*psis[j] - sum(psis[psis_id[psis_id .!= j]]) + 1.0 ) / n_ps
      else
        prob[i,j] = 0.0
      end
    end

    pred[i] = indmax(psis)

  end

  return prob::Matrix{Float64}, pred::Vector{Float64}
end

function test_adv_kernel(model::KernelMultiAdversarialModel, X_test::Matrix, y_test::Vector, X_train::Matrix, y_train::Vector)
  n = size(X_test, 1)

  y_prob, y_pred = predict_adv_kernel(model, X_test, X_train, y_train)

  # calculate testing loss
  losses = zeros(n)
  losses01 = zeros(n)
  for i=1:n
    losses[i] = 1.0 - y_prob[i, y_test[i]]
    losses01[i] = 1.0 - round(Int, y_pred[i] == y_test[i])
  end

  loss = sum(losses) / n
  loss01 = sum(losses01) / n

  return loss::Float64, losses::Vector{Float64}, loss01::Float64, losses01::Vector{Float64},
   y_prob::Matrix{Float64}, y_pred::Vector{Float64}

end

max_stars_repo_name: rizalzaf/adversarial-multiclass

max_issues_repo_name: rizalzaf/adversarial-multiclass

max_forks_repo_name: rizalzaf/adversarial-multiclass


text:
module Tar

# TODO: add some version of this method to Base
function Base.skip(io::Union{Base.Process, Base.ProcessChain}, n::Integer)
    n < 0 && throw(ArgumentError("cannot skip backwards when reading from a process"))
    for _ = 1:n
        read(io, UInt8)
    end
end

include("header.jl")
include("create.jl")
include("extract.jl")

## official API: create, list, extract

"""
    create([ predicate, ] dir, [ tarball, ]) -> tarball

        predicate :: Function
        dir       :: AbstractString
        tarball   :: Union{AbstractString, IO}

Create a tar archive ("tarball") of the directory `dir`. The resulting archive
is written to the path `tarball` or if no path is specified, a temporary path is
created and returned by the function call. If `tarball` is an IO object then the
tarball content is written to that handle instead (the handle is left open).

If a `predicate` function is passed, it is called on each system path that is
encountered while recursively searching `dir` and `path` is only included in the
tarball if `predicate(path)` is true.
"""
function create(predicate::Function, dir::AbstractString, tarball::AbstractString)
    create_dir_check(dir)
    try open(tarball, write=true) do out
            write_tarball(predicate, out, dir)
        end
    catch
        rm(tarball, force=true)
        rethrow()
    end
    return tarball
end

function create(predicate::Function, dir::AbstractString, tarball::IO)
    create_dir_check(dir)
    write_tarball(predicate, dir, out)
    return tarball
end

function create(predicate::Function, dir::AbstractString)
    create_dir_check(dir)
    tarball, out = mktemp()
    try write_tarball(predicate, out, dir)
    catch
        rm(tarball, force=true)
        rethrow()
    finally
        close(out)
    end
    return tarball
end

create(dir::AbstractString, tarball::Union{AbstractString, IO}) =
    create(p->true, dir, tarball)
create(dir::AbstractString) = create(p->true, dir)

"""
    list(tarball; [ strict = true ]) -> Vector{Header}

        tarball   :: Union{AbstractString, IO}
        strict    :: Bool

List the contents of a tar archive ("tarball") located at the path `tarball`.
If `tarball` is an IO handle, read the tar contents from that stream. Returns
a vector of `Header` structs. See [`Header`](@ref) for details.

By default `list` will error if it encounters any tarball contents which the
`extract` function would refuse to extract. With `strict=false` it will skip
these checks and list all the the contents of the tar file whether `extract`
would extract them or not. Beware that malicious tarballs can do all sorts of
crafty and unexpected things to try to trick you into doing something bad.
"""
function list(tarball::AbstractString; strict::Bool=true)
    list_tarball_check(tarball)
    open(tarball) do io
        list_tarball(io, strict=strict)
    end
end

list(tarball::IO; strict::Bool=true) = list_tarball(tarball, strict=strict)

"""
    extract(tarball, [ dir, ]; [ force=false ]) -> dir

        tarball   :: Union{AbstractString, IO}
        dir       :: AbstractString
        force     :: Bool

Extract a tar archive ("tarball") located at the path `tarball` into the
directory `dir`. If `tarball` is an IO object instead of a path, then the
archive contents will be read from that IO stream. The archive is extracted to
`dir` which must either be an existing writeable directory or a path which can
be created as a directory with `mkdir`. By default, `extract` will refuse to
extract into a non-empty directory, but the `force=true` option allows this. If
`dir` is not specified, the archive is extracted into a temporary directory
which is returned by the `extract` function call.
"""
function extract(
    tarball::Union{AbstractString, IO},
    dir::Union{AbstractString,Nothing} = nothing;
    force::Bool = false,
)
    if tarball isa AbstractString
        extract_tarball_check(tarball)
    end
    if dir !== nothing
        extract_dir_check(dir, force=force)
        cleanup = !ispath(dir)
    else
        dir = mktempdir()
        cleanup = true
    end
    try if tarball isa IO
            extract_tarball(tarball, dir)
        else
            open(tarball) do io
                extract_tarball(io, dir)
            end
        end
    catch
        cleanup && rm(dir, force=true, recursive=true)
        rethrow()
    end
    return dir
end

## error checking utility functions

create_dir_check(dir::AbstractString) = isdir(dir) ||
    error("not a directory: $dir\nUSAGE: create([predicate,] dir, [tarball])")

list_tarball_check(tarball::AbstractString) = isfile(tarball) ||
    error("not a file: $tarball\nUSAGE: list(tarball)")

extract_tarball_check(tarball::AbstractString) = isfile(tarball) ||
    error("not a file: $tarball\nUSAGE: extract(tarball, [dir])")

function extract_dir_check(dir::AbstractString; force::Bool=false)
    st = stat(dir)
    if !isdir(st)
        ispath(st) &&
            error("not a directory: $dir\nUSAGE: extract(tarball, [dir])")
        mkdir(dir)
    else
        force || isempty(readdir(dir)) ||
            error("directory not empty: $dir\n USAGE: extract(tarball, [dir])")
    end
end

end # module


text:
@diffusion_process SIRAux begin
    :parameters
    (α, β, σ1, σ2, t, T) --> Float64
    (u, v) --> ℝ{2}

    :additional
    linear --> true
    constdiff --> true
    statespace -->BoundedStateSpace(((1,2), (0.0,0.0)), ((1,2), (1.0,1.0)))
end

function B(t, P::SIRAux)
    @SMatrix[
        (P.α*(1 - P.v[1] - P.v[2]) - P.β)  0.0;
        P.β   0.0
    ]
end

β(t, P::SIRAux) = ℝ{2}(0.0, 0.0)

function σ(t, P::SIRAux)
    @SMatrix Float64[
        (-P.σ1*(1 - v[1] - v[2])*v[1])  -P.σ2*v[1];
        0.0   P.σ2*v[1]
    ]
end


text:
using DifferentialEquations
using CUDA
using LinearAlgebra

function cpu_example(n)
    tf = 10
    u0 = float_to_32_cpu(complex(ones(2^n)))
    u0 = u0 / norm(u0)
    # generate random 2^n x 2^n positive matrix (with complex entires)
    a0 = float_to_32_cpu(rand(2^n, 2^n) + 1.0im * rand(2^n, 2^n))
    H0 = conj(transpose(a0)) * a0
    a1 = float_to_32(rand(2^n, 2^n) + 1.0im * rand(2^n, 2^n))
    H1 = conj(transpose(a0)) * a0

    function interpolate_H(t)
        """
        Returns (1-s)*H0 + s*H1 where s ∈ [0, 1].
        """
        s = t / tf
        return (1-s)*H0 + s*H1
    end

    function update_func!(A, u, p, t)
        """
        Updator for DiffEqArrayOperator
        """
        fill!(A, 0.0)
        current_H = interpolate_H(t)
        @inbounds axpy!(-1.0im, current_H, A)
    end

    diff_op = DiffEqArrayOperator(-1.0im*copy(H0), update_func = update_func!)
    @show typeof(diff_op)
    jac_cache = similar(H0)
    jac_op = DiffEqArrayOperator(jac_cache, update_func = update_func!)
    ff = ODEFunction(diff_op, jac_prototype = jac_op)

    prob = ODEProblem{true}(ff, u0, (0.0, tf))
    sol = solve(prob; alg_hints = [:nonstiff])

    return sol
end

function gpu_example(n)
    tf = 10
    u0 = float_to_32_gpu(complex(ones(2^n)))
    u0 = cu(u0 / norm(u0))
    # generate random 2^n x 2^n positive matrix (with complex entires)
    a0 = float_to_32_gpu(rand(2^n, 2^n) + 1.0im * rand(2^n, 2^n))
    H0 = cu(conj(transpose(a0)) * a0)
    a1 = float_to_32_gpu(rand(2^n, 2^n) + 1.0im * rand(2^n, 2^n))
    H1 = cu(conj(transpose(a0)) * a0)

    function interpolate_H(t)
        """
        Returns (1-s)*H0 + s*H1 where s ∈ [0, 1].
        """
        s = t / tf
        return (1-s)*H0 + s*H1
    end

    function update_func!(A, u, p, t)
        """
        Updator for DiffEqArrayOperator
        """
        fill!(A, 0.0)
        current_H = interpolate_H(t)
        @inbounds axpy!(-1.0im, current_H, A)
    end

    diff_op = cu(DiffEqArrayOperator(-1.0im*copy(H0), update_func = update_func!))
    jac_cache = similar(H0)
    jac_op = cu(DiffEqArrayOperator(jac_cache, update_func = update_func!))
    ff = ODEFunction(diff_op, jac_prototype = jac_op)

    prob = ODEProblem{true}(ff, u0, (0.0, tf))
    sol = solve(prob; alg_hints = [:nonstiff])

    return sol
end

max_stars_repo_name: naezzell/accelqat

max_issues_repo_name: naezzell/accelqat

max_forks_repo_name: naezzell/accelqat


text:
Random.seed!(1234)

p     = 2
P1    = randn(p)
P1   *= P1'
mu    = zeros(p)
mvg   = MvGaussianCanon(mu,P1)
chain = chaingraph([ Factor( (x,v)->nextevent_bps(mvg, x, v),
                              x   ->gradloglik(mvg, x),
                              i) for i in 1:3])

Random.seed!(53)

# test_local_definegraph defines a chain.
nvars = chain.structure.nvars
x0    = randn(nvars)
v0    = randn(nvars)
v0   /= norm(v0)

all_evlist = AllEventList(Float64, nvars)

for i in 1:nvars
    evi = Event(x0[i], v0[i], 0.0)
    pushevent!(all_evlist.evl[i], evi)
end

i = rand(1:nvars)
t = rand()

@test isapprox(samplelocalpath(all_evlist.evl[i], t), x0[i]+t*v0[i])

t1 = t+rand()
ev1 = Event(randn(), randn(), t1)
pushevent!(all_evlist.evl[i], ev1)
t2 = t1+rand()

@test isapprox( samplelocalpath(all_evlist.evl[i], t), x0[i]+t*v0[i])
@test isapprox( samplelocalpath(all_evlist.evl[i], t2), ev1.x+(t2-ev1.t)*ev1.v)

t3  = t2+rand()
ev2 = Event(randn(), randn(), t2)
pushevent!(all_evlist.evl[i], ev2)

ts  = sort(rand(100))
ts /= maximum(ts)
ts *= t3
ss  = samplelocalpath(all_evlist.evl[i],ts)
ses = [ ts[j] < t1 ? (x0[i]+ts[j]*v0[i]) :
            (ts[j] < t2 ? (ev1.x+(ts[j]-ev1.t)*ev1.v) :
                (ev2.x+(ts[j]-ev2.t)*ev2.v)) for j in 1:length(ts)]

# sample localpath (multiple between)
@test isapprox( norm(ss-ses)/length(ss) , 0.0 )

################
all_evlist = AllEventList(Vector{Float64}, nvars)

xs = [[0.,0.] [1.,1.] [2.,1.5] [3.,5.0] [3.4,-2.0]]
ts = [0.0, 1.0, 2.0, 2.5, 2.7]
i  = rand(1:nvars)
for k in 1:length(ts)-1
    x = xs[:,k]
    v = (xs[:,k+1]-x) / (ts[k+1]-ts[k])
    evk = Event(x, v, ts[k])
    pushevent!(all_evlist.evl[i], evk)
end

Ns = 10000
ss = range(0.0, stop=ts[end], length=Ns)
# compare with classical
ss1 = samplelocalpath(all_evlist.evl[i], ss)
ss2 = samplepath(Path(xs,ts), ss)
@test norm(sum(ss1)/Ns - sum(ss2, dims=2)/Ns) <= 1e-10
# compare quadrature
@test norm(pathmean(Path(xs,ts)) - pathmean(all_evlist.evl[i],ts[end]))<=1e-10

############
# testing pathmean (and implicitly quadpathpoly) on entire alleventlist
all_evlist = AllEventList(Float64, nvars)
for k in 1:nvars
    t = 0.0
    for e in 1:50
        pushevent!(all_evlist.evl[k], Event(randn(),randn(),t))
        t+=rand()
    end
end
i = rand(1:nvars)
@test abs.(
        pathmean(Path(all_evlist.evl[i].xs', all_evlist.evl[i].ts)) .-
        pathmean(all_evlist, all_evlist.evl[i].ts[end])[i]
        )[1] <= 1e-10