text:
module RungeKuttaToolKit

export RKOCEvaluator, evaluate_residual!, evaluate_jacobian!,
    evaluate_error_coefficients!, evaluate_error_jacobian!,
    constrain!, compute_order!, compute_stages,
    RKOCBackpropEvaluator, evaluate_residual2, evaluate_gradient!,
    populate_explicit!, populate_implicit!,
    RKOCExplicitBackpropObjectiveFunctor,
    RKOCExplicitBackpropGradientFunctor,
    RKOCImplicitBackpropObjectiveFunctor,
    RKOCImplicitBackpropGradientFunctor,
    rkoc_explicit_backprop_functors,
    rkoc_implicit_backprop_functors

using Base.Threads: @threads, nthreads, threadid
using LinearAlgebra: mul!, ldiv!, qrfactUnblocked!

include("./ExampleMethods.jl")

##################################################################### PARTITIONS

function _int_parts_impl(n::Int, max::Int)::Vector{Vector{Pair{Int,Int}}}
    if n == 0
        [Pair{Int,Int}[]]
    elseif max == 0
        Vector{Pair{Int,Int}}[]
    else
        result = Vector{Pair{Int,Int}}[]
        for m = max : -1 : 1, q = div(n, m) : -1 : 1
            for p in _int_parts_impl(n - q * m, m - 1)
                push!(result, push!(p, m => q))
            end
        end
        result
    end
end

function _int_parts_impl(n::Int, max::Int, len::Int)::Vector{Vector{Int}}
    if n == 0
        [zeros(Int, len)]
    elseif max * len < n
        Vector{Int}[]
    else
        result = Vector{Int}[]
        for m = min(n, max) : -1 : div(n + len - 1, len)
            for p in _int_parts_impl(n - m, m, len - 1)
                push!(result, push!(p, m))
            end
        end
        result
    end
end

integer_partitions(n::Int)::Vector{Vector{Pair{Int,Int}}} =
    reverse!.(_int_parts_impl(n, n))
integer_partitions(n::Int, len::Int)::Vector{Vector{Int}} =
    reverse!.(_int_parts_impl(n, n, len))

################################################################### COMBINATIONS

function previous_permutation!(items::Vector{T})::Bool where {T}
    num_items = length(items)
    if num_items == 0; return false; end
    current_item = items[num_items]
    pivot_index = num_items - 1
    while pivot_index != 0
        next_item = items[pivot_index]
        if next_item <= current_item
            pivot_index -= 1
            current_item = next_item
        else; break; end
    end
    if pivot_index == 0; return false; end
    pivot = items[pivot_index]
    successor_index = num_items
    while items[successor_index] >= pivot; successor_index -= 1; end
    items[pivot_index], items[successor_index] =
        items[successor_index], items[pivot_index]
    reverse!(view(items, pivot_index + 1 : num_items))
    return true
end

function combinations_with_replacement(
        items::Vector{T}, n::Int)::Vector{Vector{Pair{T,Int}}} where {T}
    combinations = Vector{Pair{T,Int}}[]
    for p in integer_partitions(n, length(items))
        while true
            comb = Pair{T,Int}[]
            for (item, k) in zip(items, p)
                if k > 0; push!(comb, item => k); end
            end
            push!(combinations, comb)
            if !previous_permutation!(p); break; end
        end
    end
    combinations
end

################################################################### ROOTED TREES

struct RootedTree
    order::Int
    children::Vector{Pair{RootedTree,Int}}
end

function rooted_trees(n::Int)::Vector{Vector{RootedTree}}
    result = Vector{RootedTree}[]
    for k = 1 : n
        trees = RootedTree[]
        for partition in integer_partitions(k - 1)
            combination_candidates = [
                combinations_with_replacement(result[order], multiplicity)
                for (order, multiplicity) in partition]
            for combination_sequence in Base.product(combination_candidates...)
                push!(trees, RootedTree(k, vcat(combination_sequence...)))
            end
        end
        push!(result, trees)
    end
    result
end

function rooted_tree_count(n::Int)::BigInt
    if n < 0; error("rooted_tree_count requires non-negative argument"); end
    if n == 0; return 0; end
    counts = zeros(Rational{BigInt}, n)
    counts[1] = 1
    for i = 2 : n, k = 1 : i - 1, m = 1 : div(i - 1, k)
        @inbounds counts[i] += k * counts[k] * counts[i - k * m] // (i - 1)
    end
    @assert all(isone(denominator(c)) for c in counts)
    sum(numerator.(counts))
end

function butcher_density(tree::RootedTree)::Int
    result = tree.order
    for (subtree, multiplicity) in tree.children
        result *= butcher_density(subtree)^multiplicity
    end
    result
end

function butcher_symmetry(tree::RootedTree)::Int
    result = 1
    for (subtree, multiplicity) in tree.children
        result *= butcher_symmetry(subtree) * factorial(multiplicity)
    end
    result
end

####################################################### ROOTED TREE MANIPULATION

function Base.show(io::IO, tree::RootedTree)::Nothing
    print(io, '[')
    for (subtree, multiplicity) in tree.children
        print(io, subtree)
        if multiplicity != 1
            print(io, '^', multiplicity)
        end
    end
    print(io, ']')
end

function tree_index_table(
        trees_of_order::Vector{Vector{RootedTree}})::Dict{String,Int}
    order = length(trees_of_order)
    tree_index = Dict{String,Int}()
    i = 0
    for p = 1 : order - 1, tree in trees_of_order[p]
        tree_index[string(tree)] = (i += 1)
    end
    tree_index
end

function dependency_table(
        trees_of_order::Vector{Vector{RootedTree}})::Vector{Vector{Int}}
    order = length(trees_of_order)
    tree_index = tree_index_table(trees_of_order)
    dependencies = Vector{Int}[]
    for p = 1 : order, tree in trees_of_order[p]
        num_children = length(tree.children)
        if num_children == 0
            push!(dependencies, Int[])
        elseif num_children == 1
            child, m = tree.children[1]
            if m == 1
                a = tree_index[string(child)]
                push!(dependencies, [a])
            elseif m % 2 == 0
                half_child = RootedTree(0, [child => div(m, 2)])
                a = tree_index[string(half_child)]
                push!(dependencies, [a, a])
            else
                half_child = RootedTree(0, [child => div(m, 2)])
                half_plus_one_child = RootedTree(0, [child => div(m, 2) + 1])
                a = tree_index[string(half_child)]
                b = tree_index[string(half_plus_one_child)]
                push!(dependencies, [a, b])
            end
        else
            head_factor = RootedTree(0, tree.children[1 : 1])
            tail_factor = RootedTree(0, tree.children[2 : end])
            a = tree_index[string(head_factor)]
            b = tree_index[string(tail_factor)]
            push!(dependencies, [a, b])
        end
    end
    dependencies
end

function generation_table(dependencies::Vector{Vector{Int}})::Vector{Int}
    generation = Vector{Int}(undef, length(dependencies))
    for (i, deps) in enumerate(dependencies)
        if length(deps) == 0
            generation[i] = 0
        else
            generation[i] = maximum(generation[j] for j in deps) + 1
        end
    end
    generation
end

function deficit_table(dependencies::Vector{Vector{Int}})::Vector{Int}
    deficit = Vector{Int}(undef, length(dependencies))
    for (i, deps) in enumerate(dependencies)
        if length(deps) == 0
            deficit[i] = 0
        elseif length(deps) == 1
            deficit[i] = deficit[deps[1]] + 1
        else
            deficit[i] = maximum(deficit[j] for j in deps)
        end
    end
    deficit
end

function range_table(dependencies::Vector{Vector{Int}},
        num_stages::Int)::Vector{Tuple{Int,Int}}
    deficit = deficit_table(dependencies)
    lengths = num_stages .- deficit
    clamp!(lengths, 0, typemax(Int))
    lengths[1] = 0
    cumsum!(lengths, lengths)
    table_size = lengths[end]
    ranges = vcat([(1, 0)], collect(
        zip(lengths[1 : end - 1] .+ 1, lengths[2 : end])))
end

function instruction_table(dependencies::Vector{Vector{Int}},
        ranges::Vector{Tuple{Int,Int}})::Vector{NTuple{5,Int}}
    instructions = Vector{NTuple{5,Int}}(undef, length(dependencies))
    for (i, deps) in enumerate(dependencies)
        num_deps = length(deps)
        a, b = ranges[i]
        if num_deps == 0
            instructions[i] = (i, a, b, -1, -1)
        elseif num_deps == 1
            c, d = ranges[deps[1]]
            instructions[i] = (i, a, b, c, -1)
        else
            c, d = ranges[deps[1]]
            e, f = ranges[deps[2]]
            instructions[i] = (i, a, b, d - (b - a), f - (b - a))
        end
    end
    instructions
end

############################################################# INDEX MANIPULATION

function require!(required::Vector{Bool},
        dependencies::Vector{Vector{Int}}, index::Int)::Nothing
    if !required[index]
        required[index] = true
        for dep in dependencies[index]
            require!(required, dependencies, dep)
        end
    end
end

function assert_requirements_satisfied(required::Vector{Bool},
        dependencies::Vector{Vector{Int}}, indices::Vector{Int})::Nothing
    @assert(all(required[indices]))
    for index = 1 : length(dependencies)
        if required[index]
            @assert(all(required[dep] for dep in dependencies[index]))
        end
    end
end

function included_indices(mask::Vector{Bool})::Vector{Int}
    result = Vector{Int}(undef, sum(mask))
    count = 0
    for index = 1 : length(mask)
        if mask[index]
            result[count += 1] = index
        end
    end
    result
end

function restricted_indices(mask::Vector{Bool})::Vector{Int}
    result = Vector{Int}(undef, length(mask))
    count = 0
    for index = 1 : length(mask)
        if mask[index]
            result[index] = (count += 1)
        else
            result[index] = 0
        end
    end
    result
end

function inverse_index_table(indices::Vector{Int}, n::Int)
    result = [0 for _ = 1 : n]
    for (i, j) in enumerate(indices)
        result[j] = i
    end
    result
end

function restricted_trees_dependencies(indices::Vector{Int})
    max_index = maximum(indices)
    order = 1
    while max_index > rooted_tree_count(order)
        order += 1
    end
    trees_of_order = rooted_trees(order)
    full_dependencies = dependency_table(trees_of_order)
    required = [false for _ = 1 : length(full_dependencies)]
    for index in indices
        require!(required, full_dependencies, index)
    end
    assert_requirements_satisfied(required, full_dependencies, indices)
    included = included_indices(required)
    restrict = restricted_indices(required)
    restricted_trees = vcat(trees_of_order...)[included]
    restricted_dependencies = Vector{Int}[]
    for index in included
        push!(restricted_dependencies, restrict[full_dependencies[index]])
    end
    (restricted_trees, restricted_dependencies,
        inverse_index_table(restrict[indices], length(included)))
end

##################################################### RKOC EVALUATOR CONSTRUCTOR

struct RKOCEvaluator{T <: Real}
    order::Int
    num_stages::Int
    num_vars::Int
    num_constrs::Int
    output_indices::Vector{Int}
    rounds::Vector{Vector{NTuple{5,Int}}}
    inv_density::Vector{T}
    inv_symmetry::Vector{T}
    u::Vector{T}
    vs::Vector{Vector{T}}
end

function RKOCEvaluator{T}(order::Int, num_stages::Int) where {T <: Real}
    # TODO: This assertion should not be necessary.
    @assert(order >= 2)
    num_vars = div(num_stages * (num_stages + 1), 2)
    trees_of_order = rooted_trees(order)
    num_constrs = sum(length.(trees_of_order))
    dependencies = dependency_table(trees_of_order)
    generation = generation_table(dependencies)
    ranges = range_table(dependencies, num_stages)
    table_size = ranges[end][2]
    instructions = instruction_table(dependencies, ranges)
    rounds = [NTuple{5,Int}[] for _ = 1 : maximum(generation) - 1]
    for (i, g) in enumerate(generation)
        if g > 1; push!(rounds[g - 1], instructions[i]); end
    end
    RKOCEvaluator{T}(order, num_stages, num_vars, num_constrs, Int[], rounds,
        inv.(T.(butcher_density.(vcat(trees_of_order...)))),
        inv.(T.(butcher_symmetry.(vcat(trees_of_order...)))),
        Vector{T}(undef, table_size),
        [Vector{T}(undef, table_size) for _ = 1 : nthreads()])
end

function RKOCEvaluator{T}(indices::Vector{Int},
        num_stages::Int) where {T <: Real}
    num_vars = div(num_stages * (num_stages + 1), 2)
    trees, dependencies, output_indices =
        restricted_trees_dependencies(indices)
    num_constrs = length(indices)
    generation = generation_table(dependencies)
    ranges = range_table(dependencies, num_stages)
    table_size = ranges[end][2]
    instructions = instruction_table(dependencies, ranges)
    rounds = [NTuple{5,Int}[] for _ = 1 : maximum(generation) - 1]
    for (i, g) in enumerate(generation)
        if g > 1; push!(rounds[g - 1], instructions[i]); end
    end
    RKOCEvaluator{T}(0, num_stages, num_vars, num_constrs,
        output_indices, rounds,
        inv.(T.(butcher_density.(trees))),
        inv.(T.(butcher_symmetry.(trees))),
        Vector{T}(undef, table_size),
        [Vector{T}(undef, table_size) for _ = 1 : nthreads()])
end

######################################################## RKOC EVALUATION HELPERS

function populate_u_init!(evaluator::RKOCEvaluator{T},
        x::Vector{T})::Nothing where {T <: Real}
    u = evaluator.u
    k = 1
    for i = 1 : evaluator.num_stages - 1
        @inbounds result = x[k]
        @simd for j = 1 : i - 1
            @inbounds result += x[k + j]
        end
        @inbounds u[i] = result
        k += i
    end
end

function lvm_u!(evaluator::RKOCEvaluator{T}, dst_begin::Int, dst_end::Int,
                x::Vector{T}, src_begin::Int)::Nothing where {T <: Real}
    if dst_begin > dst_end; return; end
    u = evaluator.u
    dst_size = dst_end - dst_begin + 1
    skip = evaluator.num_stages - dst_size
    index = div(skip * (skip + 1), 2)
    @inbounds u[dst_begin] = x[index] * u[src_begin]
    for i = 1 : dst_size - 1
        index += skip
        skip += 1
        @inbounds result = x[index] * u[src_begin]
        @simd for j = 1 : i
            @inbounds result += x[index + j] * u[src_begin + j]
        end
        @inbounds u[dst_begin + i] = result
    end
end

function populate_u!(evaluator::RKOCEvaluator{T},
        x::Vector{T})::Nothing where {T <: Real}
    u = evaluator.u
    populate_u_init!(evaluator, x)
    for round in evaluator.rounds
        @threads for (_, dst_begin, dst_end, src1, src2) in round
            if src2 == -1
                lvm_u!(evaluator, dst_begin, dst_end, x, src1)
            elseif src1 == src2
                @simd ivdep for i = 0 : dst_end - dst_begin
                    @inbounds u[dst_begin + i] = u[src1 + i]^2
                end
            else
                @simd ivdep for i = 0 : dst_end - dst_begin
                    @inbounds u[dst_begin + i] = u[src1 + i] * u[src2 + i]
                end
            end
        end
    end
end

function populate_v_init!(evaluator::RKOCEvaluator{T},
        v::Vector{T}, var_index::Int)::Nothing where {T <: Real}
    j = 0
    for i = 1 : evaluator.num_stages - 1
        result = (j == var_index)
        j += 1
        for _ = 1 : i - 1
            result |= (j == var_index)
            j += 1
        end
        @inbounds v[i] = T(result)
    end
end

function lvm_v!(evaluator::RKOCEvaluator{T}, v::Vector{T}, var_index::Int,
                dst_begin::Int, dst_end::Int,
                x::Vector{T}, src_begin::Int)::Nothing where {T <: Real}
    if dst_begin > dst_end; return; end
    u = evaluator.u
    dst_size = dst_end - dst_begin + 1
    skip = evaluator.num_stages - dst_size
    index = div(skip * (skip + 1), 2)
    if index == var_index + 1
        @inbounds v[dst_begin] = (x[index - 1 + 1] * v[src_begin] +
            u[src_begin + var_index - index + 1])
    else
        @inbounds v[dst_begin] = (x[index - 1 + 1] * v[src_begin])
    end
    for i = 1 : dst_size - 1
        index += skip
        skip += 1
        @inbounds result = x[index] * v[src_begin]
        @simd for j = 1 : i
            @inbounds result += x[index + j] * v[src_begin + j]
        end
        if (index <= var_index + 1) && (var_index + 1 <= index + i)
            @inbounds result += u[src_begin + var_index - index + 1]
        end
        @inbounds v[dst_begin + i] = result
    end
end

function populate_v!(evaluator::RKOCEvaluator{T}, v::Vector{T},
                     x::Vector{T}, var_index::Int)::Nothing where {T <: Real}
    u = evaluator.u
    populate_v_init!(evaluator, v, var_index)
    for round in evaluator.rounds
        for (_, dst_begin, dst_end, src1, src2) in round
            if src2 == -1
                lvm_v!(evaluator, v, var_index, dst_begin, dst_end, x, src1)
            elseif src1 == src2
                @simd ivdep for i = 0 : dst_end - dst_begin
                    temp = u[src1 + i] * v[src1 + i]
                    @inbounds v[dst_begin + i] = temp + temp
                end
            else
                @simd ivdep for i = 0 : dst_end - dst_begin
                    @inbounds v[dst_begin + i] = (u[src1 + i] * v[src2 + i] +
                                                  u[src2 + i] * v[src1 + i])
                end
            end
        end
    end
end

@inline function dot_inplace(n::Int,
        v::Vector{T}, v_offset::Int,
        w::Vector{T}, w_offset::Int)::T where {T <: Real}
    result = zero(T)
    @simd for i = 1 : n
        @inbounds result += v[v_offset + i] * w[w_offset + i]
    end
    result
end

################################################################ RKOC EVALUATION

function evaluate_residual!(res::Vector{T}, x::Vector{T},
        evaluator::RKOCEvaluator{T})::Nothing where {T <: Real}
    @assert(length(x) == evaluator.num_vars)
    u = evaluator.u
    num_stages = evaluator.num_stages
    num_vars = evaluator.num_vars
    output_indices = evaluator.output_indices
    inv_density = evaluator.inv_density
    populate_u!(evaluator, x)
    if length(output_indices) == 0
        let
            first = -one(T)
            b_idx = num_vars - num_stages + 1
            @simd for i = b_idx : num_vars
                @inbounds first += x[i]
            end
            @inbounds res[1] = first
        end
        @inbounds res[2] = dot_inplace(num_stages - 1, u, 0,
                                       x, num_vars - num_stages + 1) - T(0.5)
        for round in evaluator.rounds
            @threads for (res_index, dst_begin, dst_end, _, _) in round
                j = dst_begin - 1
                n = dst_end - j
                @inbounds res[res_index] =
                    dot_inplace(n, u, j, x, num_vars - n) -
                    inv_density[res_index]
            end
        end
    else
        let
            output_index = output_indices[1]
            if output_index > 0
                first = -one(T)
                b_idx = num_vars - num_stages + 1
                @simd for i = b_idx : num_vars
                    @inbounds first += x[i]
                end
                @inbounds res[output_index] = first
            end
        end
        let
            output_index = output_indices[2]
            if output_index > 0
                b_idx = num_vars - num_stages + 1
                @inbounds res[output_index] =
                    dot_inplace(num_stages - 1, u, 0, x, b_idx) - T(0.5)
            end
        end
        for round in evaluator.rounds
            @threads for (res_index, dst_begin, dst_end, _, _) in round
                output_index = output_indices[res_index]
                if output_index > 0
                    j = dst_begin - 1
                    n = dst_end - j
                    @inbounds res[output_index] =
                        dot_inplace(n, u, j, x, num_vars - n) -
                        inv_density[res_index]
                end
            end
        end
    end
end

function evaluate_error_coefficients!(res::Vector{T}, x::Vector{T},
        evaluator::RKOCEvaluator{T})::Nothing where {T <: Real}
    @assert(length(x) == evaluator.num_vars)
    u = evaluator.u
    num_stages = evaluator.num_stages
    num_vars = evaluator.num_vars
    output_indices = evaluator.output_indices
    inv_density = evaluator.inv_density
    inv_symmetry = evaluator.inv_symmetry
    populate_u!(evaluator, x)
    if length(output_indices) == 0
        let
            first = -one(T)
            b_idx = num_vars - num_stages + 1
            @simd for i = b_idx : num_vars
                @inbounds first += x[i]
            end
            @inbounds res[1] = first
        end
        @inbounds res[2] = dot_inplace(num_stages - 1, u, 0,
                                       x, num_vars - num_stages + 1) - T(0.5)
        for round in evaluator.rounds
            @threads for (res_index, dst_begin, dst_end, _, _) in round
                j = dst_begin - 1
                n = dst_end - j
                @inbounds res[res_index] = inv_symmetry[res_index] * (
                    dot_inplace(n, u, j, x, num_vars - n) -
                    inv_density[res_index])
            end
        end
    else
        let
            output_index = output_indices[1]
            if output_index > 0
                first = -one(T)
                b_idx = num_vars - num_stages + 1
                @simd for i = b_idx : num_vars
                    @inbounds first += x[i]
                end
                @inbounds res[output_index] = first
            end
        end
        let
            output_index = output_indices[2]
            if output_index > 0
                b_idx = num_vars - num_stages + 1
                @inbounds res[output_index] =
                    dot_inplace(num_stages - 1, u, 0, x, b_idx) - T(0.5)
            end
        end
        for round in evaluator.rounds
            @threads for (res_index, dst_begin, dst_end, _, _) in round
                output_index = output_indices[res_index]
                if output_index > 0
                    j = dst_begin - 1
                    n = dst_end - j
                    @inbounds res[output_index] = inv_symmetry[res_index] * (
                        dot_inplace(n, u, j, x, num_vars - n) -
                        inv_density[res_index])
                end
            end
        end
    end
end

########################################################### RKOC DIFFERENTIATION

function evaluate_jacobian!(jac::Matrix{T}, x::Vector{T},
        evaluator::RKOCEvaluator{T})::Nothing where {T <: Real}
    @assert(length(x) == evaluator.num_vars)
    u = evaluator.u
    num_stages = evaluator.num_stages
    num_vars = evaluator.num_vars
    output_indices = evaluator.output_indices
    populate_u!(evaluator, x)
    @threads for var_idx = 1 : num_vars
        @inbounds v = evaluator.vs[threadid()]
        populate_v!(evaluator, v, x, var_idx - 1)
        if length(output_indices) == 0
            @inbounds jac[1, var_idx] = T(var_idx + num_stages > num_vars)
            let
                n = num_stages - 1
                m = num_vars - n
                result = dot_inplace(n, v, 0, x, m)
                if var_idx + n > num_vars
                    @inbounds result += u[var_idx - m]
                end
                @inbounds jac[2, var_idx] = result
            end
            for round in evaluator.rounds
                for (res_index, dst_begin, dst_end, _, _) in round
                    n = dst_end - dst_begin + 1
                    m = num_vars - n
                    result = dot_inplace(n, v, dst_begin - 1, x, m)
                    if var_idx + n > num_vars
                        @inbounds result += u[dst_begin - 1 + var_idx - m]
                    end
                    @inbounds jac[res_index, var_idx] = result
                end
            end
        else
            let
                output_index = output_indices[1]
                if output_index > 0
                    @inbounds jac[output_index, var_idx] =
                        T(var_idx + num_stages > num_vars)
                end
            end
            let
                output_index = output_indices[2]
                if output_index > 0
                    n = num_stages - 1
                    m = num_vars - n
                    result = dot_inplace(n, v, 0, x, m)
                    if var_idx + n > num_vars
                        @inbounds result += u[var_idx - m]
                    end
                    @inbounds jac[output_index, var_idx] = result
                end
            end
            for round in evaluator.rounds
                for (res_index, dst_begin, dst_end, _, _) in round
                    output_index = output_indices[res_index]
                    if output_index > 0
                        n = dst_end - dst_begin + 1
                        m = num_vars - n
                        result = dot_inplace(n, v, dst_begin - 1, x, m)
                        if var_idx + n > num_vars
                            @inbounds result += u[dst_begin - 1 + var_idx - m]
                        end
                        @inbounds jac[output_index, var_idx] = result
                    end
                end
            end
        end
    end
end

function evaluate_error_jacobian!(jac::Matrix{T}, x::Vector{T},
        evaluator::RKOCEvaluator{T})::Nothing where {T <: Real}
    @assert(length(x) == evaluator.num_vars)
    u = evaluator.u
    num_stages = evaluator.num_stages
    num_vars = evaluator.num_vars
    output_indices = evaluator.output_indices
    inv_symmetry = evaluator.inv_symmetry
    populate_u!(evaluator, x)
    @threads for var_idx = 1 : num_vars
        @inbounds v = evaluator.vs[threadid()]
        populate_v!(evaluator, v, x, var_idx - 1)
        if length(output_indices) == 0
            @inbounds jac[1, var_idx] = T(var_idx + num_stages > num_vars)
            let
                n = num_stages - 1
                m = num_vars - n
                result = dot_inplace(n, v, 0, x, m)
                if var_idx + n > num_vars
                    @inbounds result += u[var_idx - m]
                end
                @inbounds jac[2, var_idx] = result
            end
            for round in evaluator.rounds
                for (res_index, dst_begin, dst_end, _, _) in round
                    n = dst_end - dst_begin + 1
                    m = num_vars - n
                    result = dot_inplace(n, v, dst_begin - 1, x, m)
                    if var_idx + n > num_vars
                        @inbounds result += u[dst_begin - 1 + var_idx - m]
                    end
                    result *= inv_symmetry[res_index]
                    @inbounds jac[res_index, var_idx] = result
                end
            end
        else
            let
                output_index = output_indices[1]
                if output_index > 0
                    @inbounds jac[output_index, var_idx] =
                        T(var_idx + num_stages > num_vars)
                end
            end
            let
                output_index = output_indices[2]
                if output_index > 0
                    n = num_stages - 1
                    m = num_vars - n
                    result = dot_inplace(n, v, 0, x, m)
                    if var_idx + n > num_vars
                        @inbounds result += u[var_idx - m]
                    end
                    @inbounds jac[output_index, var_idx] = result
                end
            end
            for round in evaluator.rounds
                for (res_index, dst_begin, dst_end, _, _) in round
                    output_index = output_indices[res_index]
                    if output_index > 0
                        n = dst_end - dst_begin + 1
                        m = num_vars - n
                        result = dot_inplace(n, v, dst_begin - 1, x, m)
                        if var_idx + n > num_vars
                            @inbounds result += u[dst_begin - 1 + var_idx - m]
                        end
                        result *= inv_symmetry[res_index]
                        @inbounds jac[output_index, var_idx] = result
                    end
                end
            end
        end
    end
end

############################################################## FUNCTOR INTERFACE

struct RKOCEvaluatorAdjointProxy{T <: Real}
    evaluator::RKOCEvaluator{T}
end

function Base.adjoint(evaluator::RKOCEvaluator{T}) where {T <: Real}
    RKOCEvaluatorAdjointProxy{T}(evaluator)
end

function (evaluator::RKOCEvaluator{T})(x::Vector{T}) where {T <: Real}
    residual = Vector{T}(undef, evaluator.num_constrs)
    evaluate_error_coefficients!(residual, x, evaluator)
    residual
end

function (proxy::RKOCEvaluatorAdjointProxy{T})(x::Vector{T}) where {T <: Real}
    jacobian = Matrix{T}(undef,
                         proxy.evaluator.num_constrs, proxy.evaluator.num_vars)
    evaluate_error_jacobian!(jacobian, x, proxy.evaluator)
    jacobian
end

############################################################## METHOD INSPECTION

@inline function norm2(x::Vector{T}, n::Int)::T where {T <: Number}
    result = zero(float(real(T)))
    @simd for i = 1 : n
        @inbounds result += abs2(x[i])
    end
    result
end

function constrain!(x::Vector{T},
                    evaluator::RKOCEvaluator{T})::T where {T <: Real}
    num_vars, num_constrs = evaluator.num_vars, evaluator.num_constrs
    x_new = Vector{T}(undef, num_vars)
    residual = Vector{T}(undef, num_constrs)
    jacobian = Matrix{T}(undef, num_constrs, num_vars)
    direction = Vector{T}(undef, num_vars)
    evaluate_error_coefficients!(residual, x, evaluator)
    obj_old = norm2(residual, num_constrs)
    while true
        evaluate_error_jacobian!(jacobian, x, evaluator)
        ldiv!(direction, qrfactUnblocked!(jacobian), residual)
        @simd ivdep for i = 1 : num_vars
            @inbounds x_new[i] = x[i] - direction[i]
        end
        evaluate_error_coefficients!(residual, x_new, evaluator)
        obj_new = norm2(residual, num_constrs)
        if obj_new < obj_old
            @simd ivdep for i = 1 : num_vars
                @inbounds x[i] = x_new[i]
            end
            obj_old = obj_new
        else
            return sqrt(obj_old)
        end
    end
end

function compute_order!(x::Vector{T}, threshold::T; verbose::Bool=false,
                        prefix::String="")::Int where {T <: Real}
    num_stages = compute_stages(x)
    order = 2
    while true
        # rmk(prefix, "Deriving conditions for order ", order, "...";
        #     verbose=verbose)
        evaluator = RKOCEvaluator{T}(order, num_stages)
        # rmk(prefix, "Testing constraints for order ", order, "...";
        #     verbose=verbose)
        obj_new = constrain!(x, evaluator)
        # say(prefix, "Residual for order ", lpad(string(order), 2, ' '),
        #     ": ", obj_new; verbose=verbose)
        if obj_new <= threshold
            order += 1
        else
            return order - 1
        end
    end
end

function compute_stages(x::Vector{T})::Int where {T <: Real}
    num_stages = div(isqrt(8 * length(x) + 1) - 1, 2)
    @assert(length(x) == div(num_stages * (num_stages + 1), 2))
    num_stages
end

##################################################################### ODE SOLVER

struct RKSolver{T <: Real}
    num_stages::Int
    coeffs::Vector{T}
    dimension::Int
    y_temp::Vector{T}
    k_temp::Matrix{T}
end

function RKSolver{T}(coeffs::Vector{T}, dimension::Int) where {T <: Real}
    num_stages = compute_stages(coeffs)
    RKSolver{T}(num_stages, coeffs, dimension, Vector{T}(undef, dimension),
        Matrix{T}(undef, dimension, num_stages))
end

function runge_kutta_step!(f!, y::Vector{T}, step_size::T,
        solver::RKSolver{T}) where {T <: Real}
    s, x, dim, y_temp, k = solver.num_stages, solver.coeffs,
        solver.dimension, solver.y_temp, solver.k_temp
    @inbounds f!(view(k, :, 1), y)
    @simd ivdep for d = 1 : dim
        @inbounds k[d, 1] *= step_size
    end
    n = 1
    for i = 2 : s
        @simd ivdep for d = 1 : dim
            @inbounds y_temp[d] = y[d]
        end
        for j = 1 : i - 1
            @simd ivdep for d = 1 : dim
                @inbounds y_temp[d] += x[n] * k[d, j]
            end
            n += 1
        end
        @inbounds f!(view(k, :, i), y_temp)
        @simd ivdep for d = 1 : dim
            @inbounds k[d, i] *= step_size
        end
    end
    for i = 1 : s
        @simd ivdep for d = 1 : dim
            @inbounds y[d] += x[n] * k[d, i]
        end
        n += 1
    end
end

########################################################## RKOC BACKPROP HELPERS

function compute_butcher_weights!(m::Matrix{T}, A::Matrix{T},
        dependencies::Vector{Pair{Int,Int}}) where {T <: Number}
    num_stages, num_constrs = size(m)
    @inbounds for i = 1 : num_constrs
        d, e = dependencies[i]
        if d == 0
            @simd ivdep for j = 1 : num_stages
                m[j,i] = one(T)
            end
        elseif e == 0
            for j = 1 : num_stages
                temp = zero(T)
                @simd for k = 1 : num_stages
                    temp += A[j,k] * m[k,d]
                end
                m[j,i] = temp
            end
        else
            @simd ivdep for j = 1 : num_stages
                m[j,i] = m[j,d] * m[j,e]
            end
        end
    end
    return m
end

function backprop_butcher_weights!(u::Matrix{T}, A::Matrix{T}, b::Vector{T},
        m::Matrix{T}, p::Vector{T}, children::Vector{Int},
        joined_siblings::Vector{Tuple{Int,Int}},
        sibling_indices::Vector{Pair{Int,Int}}) where {T <: Number}
    num_stages, num_constrs = size(u, 1), size(u, 2)
    @inbounds for r = 0 : num_constrs - 1
        i = num_constrs - r
        x = p[i]
        @simd ivdep for j = 1 : num_stages
            u[j,i] = x * b[j]
        end
        c = children[i]
        if c > 0
            for j = 1 : num_stages
                @simd for k = 1 : num_stages
                    u[j,i] += A[k,j] * u[k,c]
                end
            end
        end
        k1, k2  = sibling_indices[i]
        for k = k1 : k2
            (s, t) = joined_siblings[k]
            @simd ivdep for j = 1 : num_stages
                u[j,i] += m[j,s] .* u[j,t]
            end
        end
    end
    return u
end

function find_children_siblings(dependencies::Vector{Vector{Int}})
    children = [0 for _ in dependencies]
    siblings = [Tuple{Int,Int}[] for _ in dependencies]
    for (i, dep) in enumerate(dependencies)
        if length(dep) == 1
            children[dep[1]] = i
        elseif length(dep) == 2
            push!(siblings[dep[1]], (dep[2], i))
            push!(siblings[dep[2]], (dep[1], i))
        end
    end
    children, siblings
end

################################################################## RKOC BACKPROP

struct RKOCBackpropEvaluator{T <: Real}
    order::Int
    num_stages::Int
    num_constrs::Int
    dependencies::Vector{Pair{Int,Int}}
    children::Vector{Int}
    joined_siblings::Vector{Tuple{Int,Int}}
    sibling_indices::Vector{Pair{Int,Int}}
    inv_density::Vector{T}
    m::Matrix{T} # Matrix of Butcher weights
    u::Matrix{T} # Gradients of Butcher weights
    p::Vector{T} # Vector of doubled residuals
    q::Vector{T} # Dot products of Butcher weights
end

function pair_deps(dependencies::Vector{Vector{Int}})::Vector{Pair{Int,Int}}
    result = Pair{Int,Int}[]
    @inbounds for dep in dependencies
        if length(dep) == 0
            push!(result, 0 => 0)
        elseif length(dep) == 1
            push!(result, dep[1] => 0)
        else
            push!(result, dep[1] => dep[2])
        end
    end
    return result
end

function RKOCBackpropEvaluator{T}(order::Int, num_stages::Int) where {T <: Real}
    trees = rooted_trees(order)
    num_constrs = sum(length.(trees))
    dependencies = dependency_table(trees)
    children, siblings = find_children_siblings(dependencies)
    ends = cumsum(length.(siblings))
    starts = vcat([0], ends[1:end-1]) .+ 1
    inv_density = inv.(T.(butcher_density.(vcat(trees...))))
    RKOCBackpropEvaluator(order, num_stages, num_constrs,
        pair_deps(dependencies), children,
        vcat(siblings...), Pair{Int,Int}.(starts, ends),
        inv_density,
        Matrix{T}(undef, num_stages, num_constrs),
        Matrix{T}(undef, num_stages, num_constrs),
        Vector{T}(undef, num_constrs),
        Vector{T}(undef, num_constrs))
end

function evaluate_residual2(A::Matrix{T}, b::Vector{T},
        evaluator::RKOCBackpropEvaluator{T})::T where {T <: Real}
    num_constrs, inv_density = evaluator.num_constrs, evaluator.inv_density
    m, q = evaluator.m, evaluator.q
    compute_butcher_weights!(m, A, evaluator.dependencies)
    mul!(q, transpose(m), b)
    residual = zero(T)
    @inbounds @simd for j = 1 : num_constrs
        x = q[j] - inv_density[j]
        residual += x * x
    end
    residual
end

function evaluate_gradient!(gA::Matrix{T}, gb::Vector{T}, A::Matrix{T},
        b::Vector{T}, evaluator::RKOCBackpropEvaluator{T})::T where {T <: Real}
    num_stages, num_constrs = evaluator.num_stages, evaluator.num_constrs
    inv_density, children = evaluator.inv_density, evaluator.children
    m, u, p, q = evaluator.m, evaluator.u, evaluator.p, evaluator.q
    compute_butcher_weights!(m, A, evaluator.dependencies)
    mul!(q, transpose(m), b)
    residual = zero(T)
    @inbounds @simd for j = 1 : num_constrs
        x = q[j] - inv_density[j]
        residual += x * x
        p[j] = x + x
    end
    backprop_butcher_weights!(u, A, b, m, p, children,
        evaluator.joined_siblings, evaluator.sibling_indices)
    @inbounds for t = 1 : num_stages
        @simd ivdep for s = 1 : num_stages
            gA[s,t] = zero(T)
        end
    end
    @inbounds for i = 1 : num_constrs
        j = children[i]
        if j > 0
            for t = 1 : num_stages
                @simd ivdep for s = 1 : num_stages
                    gA[s,t] += u[s,j] * m[t,i]
                end
            end
        end
    end
    mul!(gb, m, p)
    residual
end

################################################### IMPLICIT/EXPLICIT CONVERSION

function populate_explicit!(A::Matrix{T}, b::Vector{T}, x::Vector{T},
        n::Int)::Nothing where {T <: Number}
    k = 0
    for i = 1 : n
        @simd ivdep for j = 1 : i - 1
            @inbounds A[i,j] = x[k + j]
        end
        k += i - 1
        @simd ivdep for j = i : n
            @inbounds A[i,j] = zero(T)
        end
    end
    @simd ivdep for i = 1 : n
        @inbounds b[i] = x[k + i]
    end
end

function populate_explicit!(x::Vector{T}, A::Matrix{T}, b::Vector{T},
        n::Int)::Nothing where {T <: Number}
    k = 0
    for i = 2 : n
        @simd ivdep for j = 1 : i - 1
            @inbounds x[k + j] = A[i,j]
        end
        k += i - 1
    end
    @simd ivdep for i = 1 : n
        @inbounds x[k + i] = b[i]
    end
end

function populate_implicit!(A::Matrix{T}, b::Vector{T}, x::Vector{T},
        n::Int)::Nothing where {T <: Number}
    n2 = n * n
    @simd ivdep for i = 1 : n2
        @inbounds A[i] = x[i]
    end
    @simd ivdep for i = 1 : n
        @inbounds b[i] = x[n2 + i]
    end
end

function populate_implicit!(x::Vector{T}, A::Matrix{T}, b::Vector{T},
        n::Int)::Nothing where {T <: Number}
    n2 = n * n
    @simd ivdep for i = 1 : n2
        @inbounds x[i] = A[i]
    end
    @simd ivdep for i = 1 : n
        @inbounds x[n2 + i] = b[i]
    end
end

##################################################### BACKPROP FUNCTOR INTERFACE

struct RKOCExplicitBackpropObjectiveFunctor{T <: Real}
    evaluator::RKOCBackpropEvaluator{T}
    A::Matrix{T}
    b::Vector{T}
end

struct RKOCImplicitBackpropObjectiveFunctor{T <: Real}
    evaluator::RKOCBackpropEvaluator{T}
    A::Matrix{T}
    b::Vector{T}
end

struct RKOCExplicitBackpropGradientFunctor{T <: Real}
    evaluator::RKOCBackpropEvaluator{T}
    A::Matrix{T}
    b::Vector{T}
    gA::Matrix{T}
    gb::Vector{T}
end

struct RKOCImplicitBackpropGradientFunctor{T <: Real}
    evaluator::RKOCBackpropEvaluator{T}
    A::Matrix{T}
    b::Vector{T}
    gA::Matrix{T}
    gb::Vector{T}
end

@inline function (of::RKOCExplicitBackpropObjectiveFunctor{T})(
        x::Vector{T})::T where {T <: Real}
    A, b, evaluator = of.A, of.b, of.evaluator
    populate_explicit!(A, b, x, evaluator.num_stages)
    return evaluate_residual2(A, b, evaluator)
end

@inline function (of::RKOCImplicitBackpropObjectiveFunctor{T})(
        x::Vector{T})::T where {T <: Real}
    A, b, evaluator = of.A, of.b, of.evaluator
    populate_implicit!(A, b, x, evaluator.num_stages)
    return evaluate_residual2(A, b, evaluator)
end

@inline function (gf::RKOCExplicitBackpropGradientFunctor{T})(
        gx::Vector{T}, x::Vector{T})::T where {T <: Real}
    A, b, gA, gb, evaluator = gf.A, gf.b, gf.gA, gf.gb, gf.evaluator
    populate_explicit!(A, b, x, evaluator.num_stages)
    result = evaluate_gradient!(gA, gb, A, b, evaluator)
    populate_explicit!(gx, gA, gb, evaluator.num_stages)
    return result
end

@inline function (gf::RKOCImplicitBackpropGradientFunctor{T})(
        gx::Vector{T}, x::Vector{T})::T where {T <: Real}
    A, b, gA, gb, evaluator = gf.A, gf.b, gf.gA, gf.gb, gf.evaluator
    populate_implicit!(A, b, x, evaluator.num_stages)
    result = evaluate_gradient!(gA, gb, A, b, evaluator)
    populate_implicit!(gx, gA, gb, evaluator.num_stages)
    return result
end

function rkoc_explicit_backprop_functors(::Type{T}, order::Int,
                                         num_stages::Int) where {T <: Real}
    evaluator = RKOCBackpropEvaluator{T}(order, num_stages)
    A = Matrix{T}(undef, num_stages, num_stages)
    b = Vector{T}(undef, num_stages)
    gA = Matrix{T}(undef, num_stages, num_stages)
    gb = Vector{T}(undef, num_stages)
    return RKOCExplicitBackpropObjectiveFunctor{T}(evaluator, A, b),
           RKOCExplicitBackpropGradientFunctor{T}(evaluator, A, b, gA, gb)
end

function rkoc_implicit_backprop_functors(::Type{T}, order::Int,
                                         num_stages::Int) where {T <: Real}
    evaluator = RKOCBackpropEvaluator{T}(order, num_stages)
    A = Matrix{T}(undef, num_stages, num_stages)
    b = Vector{T}(undef, num_stages)
    gA = Matrix{T}(undef, num_stages, num_stages)
    gb = Vector{T}(undef, num_stages)
    return RKOCImplicitBackpropObjectiveFunctor{T}(evaluator, A, b),
           RKOCImplicitBackpropGradientFunctor{T}(evaluator, A, b, gA, gb)
end

end # module RungeKuttaToolKit


text:
# Model
include_model("acrobot")

# Horizon
T = 101

# Time step
tf = 5.0
h = tf / (T - 1)

# ul <= u <= uu
ul, uu = control_bounds(model, T, -10.0, 10.0)

# Initial and final states
x1 = [0.0; 0.0; 0.0; 0.0]
xT = [π; 0.0; 0.0; 0.0]
xl, xu = state_bounds(model, T, x1 = x1, xT = xT)

# Objective
obj = quadratic_tracking_objective(
        [t < T ? Diagonal(1.0 * ones(model.n)) : Diagonal(ones(model.n)) for t = 1:T],
        [Diagonal(1.0e-1 * ones(model.m)) for t = 1:T-1],
        [xT for t = 1:T],
        [zeros(model.m) for t = 1:T])

# Problem
prob = trajectory_optimization_problem(model,
           obj,
           T,
           h = h,
           xl = xl,
           xu = xu,
           ul = ul,
           uu = uu)

# Trajectory initialization
x0 = linear_interpolation(x1, xT, T) # linear interpolation on state
u0 = random_controls(model, T, 0.001) # random controls

# Pack trajectories into vector
z0 = pack(x0, u0, prob)

# Solve
@time z , info = solve(prob, copy(z0))

# Visualize
using Plots
x, u = unpack(z, prob)
plot(hcat(x...)', width = 2.0)
plot(hcat(u...)', width = 2.0, linetype = :steppost)

include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)
visualize!(vis, model, x, Δt = h)

ext: jl


text:
using Turing: VarEstimator, add_sample!, get_var
using Test

D = 10
ve = VarEstimator{Float64}(0, zeros(D), zeros(D))

for _ = 1:10000
    s = randn(D)
    add_sample!(ve, s)
end

var = get_var(ve)

@test var ≈ ones(D) atol=0.5


text:
rand_net()    = uppercase(randstring(2))
rand_sta()    = randstring('A':'Z', rand(3:5))
rand_loc()    = rand() < 0.5 ? "" : lpad(rand(0:9), 2, "0")
rand_tmax()   = rand(1200:7200)
rand_ts()     = round(Int64, sμ*(time() - 86400.0 + randn()))
rand_irr_id() = string(rand_net(), ".", rand_sta(), ".", rand_loc(), ".OY", rand('A':'Z'))
rand_hc()     = (rand() > 0.8) ? h_crit : 1.0f0

function mk_fc(fs::Float64)
  f_min = 0.5*fs
  fc = rand(fc_vals)
  while f_min ≤ fc
    fc = rand(fc_vals)
  end
  return fc
end

function rand_reg_id(fs::Float64, s::Bool)
  bc = getbandcode(fs)
  cha, units = iccodes_and_units(bc, s)

  # faster and less memory than join(), but looks clumsy
  id = string(rand_net(), ".", randstring('A':'Z', rand(3:5)), ".", rand_loc(), ".", cha)
  return id, units
end

function rand_loc(randomize_loctype::Bool)
  y = randomize_loctype == true ? rand() : 1.0
  datum = rand_datum()
  loc = [ rand_lat(),
          rand_lon(),
          1000.0 * rand(),
          1000.0 * rand(),
          rand_lat(),
          90.0 * rand() ]
  if y > 0.25
    return GeoLoc(datum, loc...)
  else
    g = GenLoc(loc)
    g.datum = datum
    return g
  end
end

function rand_resp(fc::Float64, n::Int64)
  if n > 0
    T = Float64
    i = max(1, div(n,2))
  else
    T = rand() < 0.5 ? Float32 : Float64
    i = rand(2:2:8)
  end
  pstub = T(10.0).*rand(T, i)
  z = zeros(complex(T), 2*i)
  p = vcat(pstub .+ pstub.*im, pstub .- pstub*im)
  if T == Float32
    resp = PZResp(1.0f0, Float32(fc), p, z)
  else
    resp = PZResp64(1.0, fc, p, z)
  end
  return resp
end

# function rand_t(fs::Float64, nx::Int64)
function rand_t(fs::Float64, nx::Int64, n::Int64, gs::Int64)
  ts = rand_ts()

  ngaps = n < 0 ? rand(0:9) : n
  L = ngaps + 2
  Δ = sμ/fs
  δ = 0.5*Δ + 1.0
  t = zeros(Int64, L, 2)

  # first row is always start time; no gap in last row for now
  t[1,1] = 1
  t[1,2] = ts
  t[L,1] = nx

  # rest are random-length time gaps
  if ngaps > 0
    r = ((gs > 1) && ((nx-1)/gs > ngaps)) ? range(gs, step=gs, stop=nx-1) : range(2, step=1, stop=nx-1)
    gi = rand(r, ngaps)
    sort!(gi)
    unique!(gi)
    lg = length(gi)
    while lg < ngaps
      append!(gi, rand(r))
      sort!(gi)
      unique!(gi)
      lg = length(gi)
    end

    # Generate Gaussian-distributed gap lengths
    gl = ceil.(Int64, max.(δ, Δ .* min.(1.0e5, 10.0 .^ abs.(randn(ngaps)))))
    # .* ((-1).^rand(Bool, ngaps))

    t[2:ngaps+1,1] .= gi
    t[2:ngaps+1,2] .= gl
  end
  return t
end

function pop_chan_tail!(C::GphysChannel, c::Bool, nx::Int64)
  src_str = string("randSeisChannel(c=", c, ", nx=",  nx, ")")
  setfield!(C, :name, randstring(rand(12:64)))        # :name
  setfield!(C, :gain, rand()*10^rand(0:10))           # :gain
  setfield!(C, :loc, rand_loc(true))                  # :loc
  setfield!(C, :src, src_str)                         # :src
  setfield!(C, :misc, rand_misc(rand(4:24)))          # :misc
  note!(C, "+source ¦ " * src_str)
  return nothing
end

# Populate a channel with irregularly-sampled (campaign-style) data
function populate_irr!(C::SeisChannel, nx::Int64)
  # number of samples
  L = nx < 1 ? 2^rand(6:12) : nx

  # start time
  ts = rand_ts()

  # generate time matrix that always starts at ts with no duplicate sample times
  t_max = round(Int64, sμ*rand_tmax())
  ti = rand(1:t_max, L)
  ti[1] = zero(Int64)
  unique!(ti)
  sort!(ti)
  broadcast!(+, ti, ti, ts)
  L = length(ti)
  t  = zeros(Int64, L, 2)
  t[:,1] .= 1:L
  t[:,2] .= ti

  # eh, let's use a uniform distribution for irregularly-sampled time series
  x = (rand() > 0.5) ? rand(Float64, L) .- 0.5 : rand(Float64, L)

  # fill fields of C
  C.id      = rand_irr_id()
  C.fs      = zero(Float64)
  C.units   = rand(irregular_units)
  C.t       = t
  C.x       = x
  pop_chan_tail!(C, true, nx)

  return nothing
end

# Populate a channel with regularly-sampled (time-series) data
function populate_chan!(C::SeisChannel, s::Bool, nx::Int64, fs_min::Float64, f0::Float64)
  # determine length, fs, fc
  fs = rand(fs_vals)
  while fs < fs_min
    fs = rand(fs_vals)
  end
  Lx = nx < 1 ? ceil(Int64, fs*rand_tmax()) : nx
  gs = (fs_min ≤ 0.0) ? round(Int64, 2*fs_min) : 1
  fc = mk_fc(fs)
  resp = (f0 > 0.0) ? fctoresp(Float32(f0), rand_hc()) : rand_resp(fc, 0)

  # populate channel
  C.id, C.units = rand_reg_id(fs, s)
  C.fs     = fs                                                     # fs
  C.resp   = resp                                                   # resp
  C.t      = rand_t(fs, Lx, -1, gs)                                 # t
  C.x      = randn(rand() < 0.5 ? Float32 : Float64, Lx)            # x
  pop_chan_tail!(C, false, nx)

  return nothing
end

"""
    randSeisChannel()

Generate a random channel of geophysical time-series data as a SeisChannel.

### Keywords
| KW      | Default | Type      | Meaning                                 |
|-------- |:--------|:----------|:-------------------------------         |
| s       | false   | Bool      | force channel to have seismic data?     |
| c       | false   | Bool      | force channel to have irregular data?   |
| nx      | 0       | Int64     | number of samples in channel [^1]       |
| fs_min  | 0.0     | Float64   | channels will have fs ≥ fs_min          |
| fc      | 0.0     | Float64   | rolloff frequency [^2]                  |

[^1]: if `nx ≤ 0`, the number of samples is determined randomly
[^2]: specifing `fc` with `c=false` returns a geophone instrument response

See also: `randSeisData`, `fctoresp`
"""
function randSeisChannel(; c::Bool=false, s::Bool=false, nx::Int64=0, fs_min::Float64=0.0, fc::Float64=0.0)
  Ch = SeisChannel()
  if c == true
    populate_irr!(Ch, nx)
  else
    populate_chan!(Ch, s, nx, fs_min, fc)
  end
  return Ch
end


text:
include("idl_types.jl")
include("common-funcs.jl")
include("common-macros.jl")

if is_apple()
    cd(idl_lib_dir) do
        Libdl.dlopen("libidl")
    end
end

function init()
    ecode = ccall((:IDL_Init, idlcall), Cint, (Cint, Ptr{Cint}, Ptr{Ptr{UInt8}}),
                  0, C_NULL, C_NULL)
    if ecode == 0
        error("IDL.init: IDL init failed")
    end
    global output_cb
    ccall((:IDL_ToutPush, idlcall), Void, (Ptr{Void},), output_cb)
end

# function execute{T<:AbstractString}(strarr::Array{T,1})
#     println("Strarray")
#     strarr =  ASCIIString[string(s) for s in strarr]
#     ecode = ccall((:IDL_Execute, idlcall), Cint, (Cint, Ptr{Ptr{UInt8}},), length(strarr), strarr)
#     if ecode != 0
#         # since error get printed by IDL, we just reset error state
#         ecode = ccall((:IDL_ExecuteStr, idlcall), Cint, (Ptr{UInt8},), "message, /RESET")
#     end
#     return nothing
# end

function execute_converted(str::AbstractString)
    # does no conversion of interpolated vars, continuation chars, or newlines
    ecode = ccall((:IDL_ExecuteStr, idlcall), Cint, (Ptr{UInt8},), str)
    if ecode != 0
        # since error get printed by IDL, we just reset error state
        ecode = ccall((:IDL_ExecuteStr, idlcall), Cint, (Ptr{UInt8},), "message, /RESET")
    end
    return true
end

function get_output(flags::Cint, buf::Ptr{UInt8}, n::Cint)
    line = unsafe_string(buf, n)
    stderr = (flags & IDL_TOUT_F_STDERR) != 0
    newline = (flags & IDL_TOUT_F_NLPOST) != 0
    if newline line = line*"\n" end
    print(line)
    return
end

output_cb = cfunction(get_output, Void, (Cint, Ptr{UInt8},Cint))

# function exit()
#     # probably better to do a .full_reset instead
#     ecode = ccall((:IDL_Cleanup, idlcall), Cint, (Cint,), IDL_TRUE)
#     return nothing
# end

# hold a ref to imported variable so they don't get gc'ed
const var_refs = Dict{Ptr{UInt8}, Any}()

function done_with_var(p::Ptr{UInt8})
    if !haskey(var_refs, p)
        error("IDL.done_with_var: ptr not found: "*string(p))
    end
    delete!(var_refs, p)
    return
end

free_cb = cfunction(done_with_var, Void, (Ptr{UInt8},))

function put_var{T,N}(arr::Array{T,N}, name::AbstractString)
    if !isbits(eltype(arr)) || (idl_type(arr) < 0)
        error("IDL.put_var: only works with some vars containing bits types")
    end
    dim = zeros(Int, IDL_MAX_ARRAY_DIM)
    dim[1:N] = [size(arr)...]
    vptr = ccall((:IDL_ImportNamedArray, idlcall), Ptr{IDL_Variable},
                 (Ptr{UInt8}, Cint, IDL_ARRAY_DIM, Cint, Ptr{UInt8}, IDL_ARRAY_FREE_CB , Ptr{Void}),
                 name, N, dim, idl_type(arr), pointer(arr), free_cb, C_NULL)
    if vptr == C_NULL
        error("IDL.put_var: failed")
    end
    var_refs[pointer(arr)] = (name, vptr, arr)
    return
end

function put_var(x, name::AbstractString)
    # Sort of a HACK: import as one-element array and then truncate to scalar
    # IDL_ImportArray(int n_dim, IDL_MEMINT dim[], int type,
    #                 UCHAR *data, IDL_ARRAY_FREE_CB free_cb, void *s)
    if !isbits(x) || (idl_type(x) < 0)
        error("IDL.put_var: only works with some vars containing bits types")
    end
    dim = zeros(Int, IDL_MAX_ARRAY_DIM)
    dim[1] = 1
    ccall((:IDL_ImportNamedArray, idlcall), Ptr{Void},
          (Ptr{UInt8}, Cint, Ptr{IDL_MEMINT}, Cint, Ptr{UInt8}, Ptr{Void}, Ptr{Void}),
          name, 1, dim, idl_type(x), pointer([x]), C_NULL, C_NULL)
    execute("$name = $name[0]")
    return
end

function put_var(str::AbstractString, name::AbstractString)
    # Sort of a HACK: do direcly since ImportNamedArray doesn't work
    execute("$name = '$str'")
    return
end

function get_name(vptr::Ptr{IDL_Variable})
    str = ccall((:IDL_VarName, idlcall), Ptr{Cchar}, (Ptr{IDL_Variable},), vptr)
    return unsafe_string(str)
end

function get_vptr(name::AbstractString)
    # returns C_NULL if name not in scope
    name = uppercase(name)
    vptr = ccall((:IDL_GetVarAddr, idlcall), Ptr{IDL_Variable}, (Ptr{UInt8},), name)
    vptr
end

function get_var(name::AbstractString)
    name = uppercase(name)
    vptr = ccall((:IDL_GetVarAddr, idlcall), Ptr{IDL_Variable}, (Ptr{UInt8},), name)
    if vptr == C_NULL
        error("IDL.get_var: variable $name does not exist")
    end
    get_var(vptr)
end


text:
using Test
using Dictionaries
using Indexing

include("Indices.jl")
include("Dictionary.jl")
include("ArrayIndices.jl")
include("ArrayDictionary.jl")
include("PairDictionary.jl")
include("FillDictionary.jl")
include("indexing.jl")
include("foreach.jl")
include("map.jl")
include("broadcast.jl")
include("filter.jl")
include("find.jl")
include("reverse.jl")
include("show.jl")

# Run the following test without julia --check-bounds=yes mode
cmd = deepcopy(Base.julia_cmd())
filter!(a->!startswith(a, "--check-bounds="), cmd.exec)
@test process_exited(run(`$cmd auto_boundscheck.jl`))


text:
include("helpers.jl")
include("algo_util.jl")

abstract type ZerothOrder end

function solve(M::ZerothOrder, f, x, max_iters; num_eval_termination=true)
    init!(M, x)
    x_hist = [x]
    y, terminate = f(x), false
    while !terminate
        x, y, terminate = step!(M, f, x, y)
        push!(x_hist, x)
        if num_eval_termination && (count(f) >= max_iters - M.evals_per_iter)
            break
        end
    end
    return x, x_hist
end

# Coordinate Descent
Base.@kwdef mutable struct CoordinateDescentAcceleration <: ZerothOrder
    ϵ = 1e-4
    approx_line_search = backtracking_line_search
    n = nothing
end
function init!(M::CoordinateDescentAcceleration, x)
    M.n = length(x)
end
function step!(M::CoordinateDescentAcceleration, f, ∇f, x)
    ϵ, n, approx_line_search = M.ϵ, M.n, M.approx_line_search
    terminate = false

    x_old = copy(x)
    for i in 1:n
        d = basis(i,n)
        α = backtracking_line_search(f, ∇f, x, d, 1e-2)
        x += α * d
    end
    α = backtracking_line_search(f, ∇f, x, x - x_old, 1e-2)
    x += α * (x - x_old)

    if abs(norm(x - x_old)) <= ϵ
        terminate = true
    end
    return x, terminate
end

function solve(M::CoordinateDescentAcceleration, f, ∇f, x, max_iters; num_eval_termination=true)
    init!(M, x)
    x_hist = [x]
    terminate = false
    while !terminate
        x, terminate = step!(M, f, ∇f, x)
        push!(x_hist, x)
        if num_eval_termination && (count(f) >= max_iters - M.evals_per_iter)
            break
        end
    end
    return x, x_hist
end

# Hooke Jeeves
Base.@kwdef mutable struct HookeJeeves <: ZerothOrder
    α = 1e-2
    ϵ = 1e-4
    γ = 0.5
    n = nothing
    evals_per_iter = nothing
end
function init!(M::HookeJeeves, x)
    M.n = length(x)
    M.evals_per_iter = 2 * M.n
end
function step!(M::HookeJeeves, f, x, y)
    α, ϵ, γ, n = M.α, M.ϵ, M.γ, M.n
    improved, terminate = false, false

    x_best, y_best = x, y
    xs_new = [x + sgn * α * basis(i,n) for i in 1:n for sgn in (-1, +1)]
    ys_new = [f(x_new) for x_new in xs_new]
    y_best_iter, y_best_iter_idx = findmin(ys_new)

    if y_best_iter < y_best
        x_best = xs_new[y_best_iter_idx]
        y_best = y_best_iter
        improved = true
    end

    M.α *= (!improved ? γ : 1)
    terminate = (M.α <= ϵ ? true : false)
    return x_best, y_best, terminate
end

# Hooke Jeeves Dynamic w/ Eager execution
Base.@kwdef mutable struct HookeJeevesDynamic <: ZerothOrder
    α = 1e-2
    ϵ = 1e-4
    γ = 0.5
    n = nothing
    evals_per_iter = nothing
    D = nothing # directions to search in, need to store it for permutating order
end
function init!(M::HookeJeevesDynamic, x)
    M.n = length(x)
    M.evals_per_iter = 2 * M.n
    M.D = [sgn * basis(i, M.n) for i in 1:M.n for sgn in (-1, +1)]
end
function step!(M::HookeJeevesDynamic, f, x, y, idx_best_prev)
    α, ϵ, γ = M.α, M.ϵ, M.γ
    improved, terminate, idx_best = false, false, 1

    x_best, y_best = x, y
    d_best_prev = M.D[idx_best_prev]
    M.D = pushfirst!(deleteat!(M.D, idx_best_prev), d_best_prev)
    xs_new = [x + α * d for d in M.D]

    for (idx, x_new) in enumerate(xs_new)
        y_new = f(x_new)
        if y_new < y_best
            x_best, y_best = x_new, y_new
            improved = true
            idx_best = idx
            break
        end
    end

    M.α *= (!improved ? γ : 1)
    terminate = (M.α <= ϵ ? true : false)
    return x_best, y_best, terminate, idx_best
end

function solve(M::HookeJeevesDynamic, f, x, max_iters; num_eval_termination=true)
    init!(M, x)
    x_hist = [x]
    y, terminate, idx_best = f(x), false, 1

    while !terminate
        x, y, terminate, idx_best = step!(M, f, x, y, idx_best)
        push!(x_hist, x)

        if num_eval_termination && (count(f) >= max_iters - M.evals_per_iter)
            break
        end
    end
    return x, x_hist
end

# one-off function to generate cool plot. do not actually use it.
function solve_diag(M::HookeJeevesDynamic, f, x, max_iters; num_eval_termination=true)
    init!(M, x)
    M.D = [[√2, √2], [-√2, -√2], [-√2, √2], [√2, -√2]]
    x_hist = [x]
    y, terminate, idx_best = f(x), false, 1

    while !terminate
        x, y, terminate, idx_best = step!(M, f, x, y, idx_best)
        push!(x_hist, x)

        if num_eval_termination && (count(f) >= max_iters - M.evals_per_iter)
            break
        end
    end
    return x, x_hist
end


text:
if !(VERSION.major==0&&VERSION.minor<=6)
	using Test
else
	using Base.Test
end
using FinancialToolbox

#Test Parameters
testToll=1e-14;

print_colored("Starting Standard Test\n",:green)
spot=10;K=10;r=0.02;T=2.0;sigma=0.2;d=0.01;
assert_(value,toll)=@test abs(value)<toll
#EuropeanCall Option
PriceCall=blsprice(spot,K,r,T,sigma,d);
PriceCallBlack=blkprice(spot,K,r,T,sigma);
DeltaCall=blsdelta(spot,K,r,T,sigma,d);
ThetaCall=blstheta(spot,K,r,T,sigma,d);
RhoCall=blsrho(spot,K,r,T,sigma,d);
VannaCall=blsvanna(spot,K,r,T,sigma,d);
PsiCall=blspsi(spot,K,r,T,sigma,d);
LambdaCall=blslambda(spot,K,r,T,sigma,d);
SigmaCall=blsimpv(spot, K, r, T, PriceCall, d);
SigmaCallBlack=blkimpv(spot, K, r, T, PriceCallBlack);

#EuropeanPut Option
PricePut=blsprice(spot,K,r,T,sigma,d,false);
PricePutBlack=blkprice(spot,K,r,T,sigma,false);
DeltaPut=blsdelta(spot,K,r,T,sigma,d,false);
ThetaPut=blstheta(spot,K,r,T,sigma,d,false);
VannaPut=blsvanna(spot,K,r,T,sigma,d,false);
PsiPut=blspsi(spot,K,r,T,sigma,d,false);
LambdaPut=blslambda(spot,K,r,T,sigma,d,false);
RhoPut=blsrho(spot,K,r,T,sigma,d,false);
SigmaPut=blsimpv(spot,K,r,T,PricePut,d,false);
SigmaPutBlack=blkimpv(spot,K,r,T,PricePutBlack,false);

#Equals for both Options
Gamma=blsgamma(spot,K,r,T,sigma,d);
Vega=blsvega(spot,K,r,T,sigma,d);

print_colored("---  European Call: Price and Sensitivities\n",:yellow)
#Standard Test European Call Option
print_colored("-----Testing Price\n",:blue);
assert_(PriceCall-1.191201316999582,testToll)
print_colored("-----Testing Black Price\n",:blue);
assert_(PriceCallBlack-1.080531820066428,testToll)
print_colored("-----Testing Delta\n",:blue);
assert_(DeltaCall-0.572434050810368,testToll)
print_colored("-----Testing Theta\n",:blue);
assert_(ThetaCall+0.303776337550247,testToll)
print_colored("-----Testing Rho\n",:blue);
assert_(RhoCall-9.066278382208203,testToll)
print_colored("-----Testing Lambda\n",:blue);
assert_(LambdaCall-4.805518955034612,testToll)
print_colored("-----Testing Implied Volatility\n",:blue);
assert_(SigmaCall-0.2,testToll)
print_colored("-----Testing Implied Volatility Black\n",:blue);
assert_(SigmaCallBlack-0.2,testToll)
#Low Xtol for blsimpv
tol_low=1e-16;
tol_high=1e-3;
SigmaLowXTol=blsimpv(spot, K, r, T, PriceCall, d,true,tol_low,tol_high);
print_colored("-----Testing Implied Volatility Low X tol, High Y tol\n",:blue);
@test((abs(blsprice(spot,K,r,T,SigmaLowXTol,d,true)-PriceCall)<tol_high)||abs(SigmaLowXTol-sigma)<tol_low)
#Low Ytol for blsimpv
SigmaLowYTol=blsimpv(spot, K, r, T, PriceCall, d,true,tol_high,tol_low);
print_colored("-----Testing Implied Volatility Low Y tol, High X tol\n",:blue);
@test((abs(blsprice(spot,K,r,T,SigmaLowYTol,d,true)-PriceCall)<tol_low)||abs(SigmaLowYTol-sigma)<tol_high)

print_colored("---  European Put: Price and Sensitivities\n",:yellow)
#Standard Test European Put Option
print_colored("-----Testing Price\n",:blue);
assert_(PricePut-0.997108975455260,testToll)
print_colored("-----Testing Price Black\n",:blue);
assert_(PricePutBlack-1.080531820066428,testToll)
print_colored("-----Testing Delta\n",:blue);
assert_(DeltaPut+0.407764622496387,testToll)
print_colored("-----Testing Theta\n",:blue);
assert_(ThetaPut+0.209638317050458,testToll)
print_colored("-----Testing Rho\n",:blue);
assert_(RhoPut+10.149510400838260,testToll)
print_colored("-----Testing Lambda\n",:blue);
assert_(LambdaPut+4.089468980160465,testToll)
print_colored("-----Testing Implied Volatility\n",:blue);
assert_(SigmaPut-0.2,testToll)
print_colored("-----Testing Implied Volatility Black\n",:blue);
assert_(SigmaPutBlack-0.2,testToll)

#Standard Test for Common Sensitivities
print_colored("-----Testing Gamma\n",:blue);
assert_(Gamma-0.135178479404601,testToll)
print_colored("-----Testing Vega\n",:blue);
assert_(Vega-5.407139176184034,testToll)
print_colored("Standard Test Passed\n",:green)
println("")

#TEST OF INPUT VALIDATION
print_colored("Starting Input Validation Test Real\n",:magenta)
print_colored("----Testing Negative Spot Price\n",:cyan)
@test_throws(ErrorException, blsprice(-spot,K,r,T,sigma,d))
@test_throws(ErrorException, blkprice(-spot,K,r,T,sigma))
@test_throws(ErrorException, blsdelta(-spot,K,r,T,sigma,d))
@test_throws(ErrorException, blsgamma(-spot,K,r,T,sigma,d))
@test_throws(ErrorException, blstheta(-spot,K,r,T,sigma,d))
@test_throws(ErrorException, blsrho(-spot,K,r,T,sigma,d))
@test_throws(ErrorException, blsvega(-spot,K,r,T,sigma,d));
@test_throws(ErrorException, blspsi(-spot,K,r,T,sigma,d))
@test_throws(ErrorException, blslambda(-spot,K,r,T,sigma,d))
@test_throws(ErrorException, blsvanna(-spot,K,r,T,sigma,d))
@test_throws(ErrorException, blsimpv(-spot,K,r,T,PriceCall,d))
@test_throws(ErrorException, blkimpv(-spot,K,r,T,PriceCall))

print_colored("----Testing Negative Strike Price\n",:cyan)
@test_throws(ErrorException, blsprice(spot,-K,r,T,sigma,d))
@test_throws(ErrorException, blkprice(spot,-K,r,T,sigma))
@test_throws(ErrorException, blsdelta(spot,-K,r,T,sigma,d))
@test_throws(ErrorException, blsgamma(spot,-K,r,T,sigma,d))
@test_throws(ErrorException, blstheta(spot,-K,r,T,sigma,d))
@test_throws(ErrorException, blsrho(spot,-K,r,T,sigma,d))
@test_throws(ErrorException, blsvega(spot,-K,r,T,sigma,d))
@test_throws(ErrorException, blspsi(spot,-K,r,T,sigma,d))
@test_throws(ErrorException, blslambda(spot,-K,r,T,sigma,d))
@test_throws(ErrorException, blsvanna(spot,-K,r,T,sigma,d))
@test_throws(ErrorException, blsimpv(spot,-K,r,T,PriceCall,d))
@test_throws(ErrorException, blkimpv(spot,-K,r,T,PriceCall))

print_colored("----Testing Negative Time to Maturity\n",:cyan)
@test_throws(ErrorException, blsprice(spot,K,r,-T,sigma,d))
@test_throws(ErrorException, blkprice(spot,K,r,-T,sigma))
@test_throws(ErrorException, blsdelta(spot,K,r,-T,sigma,d))
@test_throws(ErrorException, blsgamma(spot,K,r,-T,sigma,d))
@test_throws(ErrorException, blstheta(spot,K,r,-T,sigma,d))
@test_throws(ErrorException, blsrho(spot,K,r,-T,sigma,d))
@test_throws(ErrorException, blsvega(spot,K,r,-T,sigma,d))
@test_throws(ErrorException, blspsi(spot,K,r,-T,sigma,d))
@test_throws(ErrorException, blslambda(spot,K,r,-T,sigma,d))
@test_throws(ErrorException, blsvanna(spot,K,r,-T,sigma,d))
@test_throws(ErrorException, blsimpv(spot,K,r,-T,PriceCall,d))
@test_throws(ErrorException, blkimpv(spot,K,r,-T,PriceCall))

print_colored("----Testing Negative Volatility\n",:cyan)
@test_throws(ErrorException, blsprice(spot,K,r,T,-sigma,d))
@test_throws(ErrorException, blkprice(spot,K,r,T,-sigma))
@test_throws(ErrorException, blsdelta(spot,K,r,T,-sigma,d))
@test_throws(ErrorException, blsgamma(spot,K,r,T,-sigma,d))
@test_throws(ErrorException, blstheta(spot,K,r,T,-sigma,d))
@test_throws(ErrorException, blsrho(spot,K,r,T,-sigma,d))
@test_throws(ErrorException, blspsi(spot,K,r,T,-sigma,d))
@test_throws(ErrorException, blslambda(spot,K,r,T,-sigma,d))
@test_throws(ErrorException, blsvanna(spot,K,r,T,-sigma,d))
@test_throws(ErrorException, blsvega(spot,K,r,T,-sigma,d))
print_colored("----Testing Negative Option Price\n",:cyan)
@test_throws(ErrorException, blsimpv(spot,K,r,T,-PriceCall,d))
@test_throws(ErrorException, blkimpv(spot,K,r,T,-PriceCall))

print_colored("----Testing Negative Tollerance\n",:cyan)
@test_throws(ErrorException, blsimpv(spot,K,r,T,PriceCall,d,true,-1e-12,1e-12))
@test_throws(ErrorException, blkimpv(spot,K,r,T,PriceCall,true,-1e-12,1e-12))
@test_throws(ErrorException, blsimpv(spot,K,r,T,PriceCall,d,true,1e-12,-1e-12))
@test_throws(ErrorException, blkimpv(spot,K,r,T,PriceCall,true,1e-12,-1e-12))

#Too low tollerance
@test_throws(ErrorException, blsimpv(spot,K,r,T,PriceCall,d,true,0.0,0.0))
@test_throws(ErrorException, blkimpv(spot,K,r,T,PriceCall,true,0.0,0.0))

print_colored("Real Input Validation Test Passed\n",:magenta)

#End of the Test


text:
function calc_error_norms(func, u, t, analyzer,
                          mesh::ParallelTreeMesh{2}, equations, initial_condition,
                          dg::DGSEM, cache, cache_analysis)
  # call the method accepting a general `mesh::TreeMesh{2}`
  # TODO: MPI, we should improve this; maybe we should dispatch on `u`
  #       and create some MPI array type, overloading broadcasting and mapreduce etc.
  #       Then, this specific array type should also work well with DiffEq etc.
  l2_error, linf_error = invoke(calc_error_norms,
    Tuple{typeof(func), typeof(u), typeof(t), typeof(analyzer), TreeMesh{2},
          typeof(equations), typeof(initial_condition), typeof(dg), typeof(cache),
          typeof(cache_analysis)},
    func, u, t, analyzer, mesh, equations, initial_condition, dg, cache, cache_analysis)

  # Since the local L2 norm is already normalized and square-rooted, we need to undo this first
  total_volume = mesh.tree.length_level_0^ndims(mesh)
  global_l2_error = Vector(l2_error.^2 .* total_volume)
  global_linf_error = Vector(linf_error)
  MPI.Reduce!(global_l2_error, +, mpi_root(), mpi_comm())
  MPI.Reduce!(global_linf_error, max, mpi_root(), mpi_comm())
  if mpi_isroot()
    l2_error   = convert(typeof(l2_error),   global_l2_error)
    linf_error = convert(typeof(linf_error), global_linf_error)
  else
    l2_error   = convert(typeof(l2_error),   NaN * global_l2_error)
    linf_error = convert(typeof(linf_error), NaN * global_linf_error)
  end

  l2_error = @. sqrt(l2_error / total_volume)

  return l2_error, linf_error
end

function integrate_via_indices(func::Func, u,
                               mesh::ParallelTreeMesh{2}, equations, dg::DGSEM, cache,
                               args...; normalize=true) where {Func}
  # call the method accepting a general `mesh::TreeMesh{2}`
  # TODO: MPI, we should improve this; maybe we should dispatch on `u`
  #       and create some MPI array type, overloading broadcasting and mapreduce etc.
  #       Then, this specific array type should also work well with DiffEq etc.
  local_integral = invoke(integrate_via_indices,
    Tuple{typeof(func), typeof(u), TreeMesh{2}, typeof(equations),
          typeof(dg), typeof(cache), map(typeof, args)...},
    func, u, mesh, equations, dg, cache, args..., normalize=normalize)

  # OBS! Global results are only calculated on MPI root, all other domains receive `nothing`
  global_integral = MPI.Reduce!(Ref(local_integral), +, mpi_root(), mpi_comm())
  if mpi_isroot()
    integral = convert(typeof(local_integral), global_integral[])
  else
    integral = convert(typeof(local_integral), NaN * local_integral)
  end

  return integral
end


text:
module SCEDC

using Dates, DelimitedFiles, Distributed, Random
using AWS, AWSS3, CSV, DataFrames, SeisIO
using SeisIO.Quake
@service S3
@service Athena

# include modules
include("catalog.jl")
include("phases.jl")
include("events.jl")
include("s3.jl")
include("athena.jl")
include("ec2.jl")

end


text:
const DEFAULT_WORKER_CPU = 1
const DEFAULT_WORKER_MEMORY = "4Gi"

# Notifies tasks that the abnormal worker deregistration warning has been emitted
const DEREGISTER_ALERT = Condition()

struct K8sClusterManager <: ClusterManager
    np::Int
    pod_name::String
    image::String
    cpu::String
    memory::String

    pending_timeout::Int
    configure::Function
end

"""
    K8sClusterManager(np::Integer; kwargs...)

A cluster manager using Kubernetes (k8s) which spawns additional pods as workers. Attempts
to spawn `np` workers but may launch with less workers if the cluster has less resources
available.

## Arguments

- `np`: Desired number of worker pods to be launched.

## Keywords

- `namespace`: the Kubernetes namespace to launch worker pods within. Defaults to
  `current_namespace()`.
- `manager_pod_name`: the name of the manager pod. Defaults to `ENV["HOSTNAME"]` which is
  the name of the pod when executed inside of a Kubernetes pod.
- `image`: Docker image to use for the workers. Defaults to using the image of the Julia
  caller if running within a pod using a single container otherwise is a required argument.
- `cpu`: [CPU resources requested](https://kubernetes.io/docs/concepts/configuration/manage-resources-containers/#meaning-of-cpu)
  for each worker. Defaults to `$(repr(DEFAULT_WORKER_CPU))`,
- `memory`: [Memory resource requested](https://kubernetes.io/docs/concepts/configuration/manage-resources-containers/#meaning-of-memory)
  for each worker in bytes. Requests may provide a unit suffix (e.g. "G" for Gigabytes and
  "GiB" for Gibibytes). Defaults to `$(repr(DEFAULT_WORKER_MEMORY))`.
- `pending_timeout`: The maximum number of seconds to wait for a "Pending" worker pod to
  enter the "Running" phase. Once the timeout has been reached the manager will continue
  with the number of workers available (`<= np`). Defaults to `180` (3 minutes).
- `configure`: A function which allows modification of the worker pod specification before
  their creation. Defaults to `identity`.
"""
function K8sClusterManager(np::Integer;
                           namespace::String=current_namespace(),
                           manager_pod_name::String=get(ENV, "HOSTNAME", "localhost"),
                           image=nothing,
                           cpu=DEFAULT_WORKER_CPU,
                           memory=DEFAULT_WORKER_MEMORY,
                           pending_timeout::Real=180,
                           configure=identity)

    # Default to using the image of the pod if possible
    if image === nothing
        pod = get_pod(manager_pod_name)
        images = map(c -> c["image"], pod["spec"]["containers"])

        if length(images) == 1
            image = first(images)
        elseif length(images) > 0
            error("Unable to determine image from pod \"$manager_pod_name\" which uses multiple containers")
        else
            error("Unable to find any images for pod \"$manager_pod_name\"")
        end
    end

    return K8sClusterManager(np, manager_pod_name, image, string(cpu), string(memory), pending_timeout, configure)
end

struct TimeoutException <: Exception
    msg::String
end

Base.showerror(io::IO, e::TimeoutException) = print(io, "TimeoutException: ", e.msg)

function worker_pod_spec(manager::K8sClusterManager; kwargs...)
    pod = worker_pod_spec(; manager_name=manager.pod_name,
                          image=manager.image,
                          cpu=manager.cpu,
                          memory=manager.memory,
                          kwargs...)

    return pod
end

function Distributed.launch(manager::K8sClusterManager, params::Dict, launched::Array, c::Condition)
    exename = params[:exename]
    exeflags = params[:exeflags]

    cmd = `$exename $exeflags --worker=$(cluster_cookie())`

    worker_manifest = worker_pod_spec(manager; cmd)

    # Note: User-defined `configure` function may or may-not be mutating
    worker_manifest = manager.configure(worker_manifest)

    @sync for i in 1:manager.np
        @async begin
            pod_name = create_pod(worker_manifest)

            pod = try
                wait_for_running_pod(pod_name; timeout=manager.pending_timeout)
            catch e
                e isa TimeoutException || rethrow()

                delete_pod(pod_name; wait=false)
                @warn sprint(showerror, e)
                nothing
            end

            if pod !== nothing
                @info "$pod_name is up"

                # Redirect any stdout/stderr from the worker to be displayed on the manager.
                # Note: `start_worker` (via `--worker`) automatically redirects stderr to
                # stdout.
                p = open(detach(`$(kubectl()) logs -f pod/$pod_name`), "r+")

                config = WorkerConfig()
                config.io = p.out
                config.userdata = (; pod_name=pod_name)

                push!(launched, config)
                notify(c)
            end
        end
    end
end

function Distributed.manage(manager::K8sClusterManager, id::Integer, config::WorkerConfig, op::Symbol)
    pod_name = config.userdata.pod_name

    if op === :register
        # Note: Labelling the pod with the worker ID is only a nice-to-have. We may want to
        # make this fail gracefully if "patch" access is unavailable.
        label_pod(pod_name, "worker-id" => id)

    elseif op === :interrupt
        os_pid = config.ospid
        if os_pid !== nothing
            try
                exec_pod(pod_name, `bash -c "kill -2 $os_pid"`)
            catch e
                @error "Error sending a Ctrl-C to julia worker $id on pod $pod_name:\n" * sprint(showerror, e)
            end
        else
            # This state can happen immediately after an addprocs
            @error "Worker $id cannot be presently interrupted."
        end

    elseif op === :deregister
        # As the deregister `manage` call occurs before remote workers are told to
        # deregister we should avoid unnecessarily blocking.
        @async begin
            # In the event of a worker pod failure Julia may notice the worker socket
            # close before Kubernetes has the pod's final status available.
            #
            # Note: The wait duration of 30 seconds was picked as it's the default
            # "grace period" used for pod termination.
            # https://kubernetes.io/docs/concepts/workloads/pods/pod-lifecycle/#pod-termination
            state = reason = nothing
            timedwait(30.0, pollint=1.0) do
                state, reason = pod_status(pod_name)
                state == "terminated"
            end

            # Report any abnormal terminations of worker pods
            if state == "terminated" && reason != "Completed"
                @warn "Worker $id on pod $pod_name was terminated due to: $reason"
            end

            notify(DEREGISTER_ALERT; all=true)
        end
    end
end


text:
using Test
using CodecZlib
using GraphIO.LGCompressed

@testset "LGCompressed" begin
    for g in values(allgraphs)
        readback_test(LGCompressedFormat(), g; testfail=true)
    end
end


text:
immutable NoncentralBeta <: ContinuousUnivariateDistribution
    alpha::Float64
    beta::Float64
    ncp::Float64
    function NonCentralBeta(a::Real, b::Real, nc::Real)
    	if a > 0.0 && b > 0.0 && nc >= 0.0
    		new(float64(a), float64(b), float64(nc))
    	else
    		error("alpha and beta must be > 0 and ncp >= 0")
    	end
    end
end

@_jl_dist_3p NoncentralBeta nbeta


text:
using Logic
using Base.Test

# number
@test zero == identity(Number, +)
@test zero == identity(Number, -)
@test one == identity(Number, *)
@test one == identity(Number, /)
@test one == identity(Number, \)

@test (-) == inverse(Number, +)
@test (+) == inverse(Number, -)
@test :(1/n) == inverse(Number, *)
@test :(1\n) == inverse(Number, /)
@test :(n/1) == inverse(Number, \)

a = Number(5)
@test a == a + identity(a, +)
@test a == a - identity(a, -)
@test a == a * identity(a, *)
@test a == a / identity(a, /)
@test a == identity(a, \) \ a

@test 0 == a + inverse(a, +)
@test 0 == a - inverse(a, -)
@test 1 == a * inverse(a, *)
@test 1 == a / inverse(a, /)
@test 1 == inverse(a, \) \ a

# matrix
@test eye == identity(Matrix, *)
@test inv == inverse(Matrix, *)

A = [4 7; 2 6]
@test A == A * identity(A, *)
@test [1 0;0 1] ≈ A * inverse(A, *)
@test [1 0 0; 0 1 0; 0 0 1] == identity([3 2 1; 4 3 2; 4 4 3], *)

# logic
@test zero == logic(Number, +, identity)
@test one == logic(Number, *, identity)
@test (-) == logic(Number, +, inverse)
@test :(1/n) == logic(Number, *, inverse)
@test eye == logic(Matrix, *, identity)
@test inv == logic(Matrix, *, inverse)
@test 0 == logic(5, +, identity)
@test 1 == logic(5, *, identity)
@test 0 == 5 + logic(5, +, inverse)
@test 1 == 5 * logic(5, *, inverse)
@test [1 0;0 1] == logic([4 7;2 6], *, identity)
@test [0.6 -0.7;-0.2 0.4] ≈ logic([4 7;2 6], *, inverse)


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

using GLM
using DataFrames
using StatsBase
using StatsModels
#import StatsBase: residuals

const INNERMOD = Union{GLM.GeneralizedLinearModel, GLM.LinearModel}
#const UNIONC = Union{Type{Nothing}, Type{LinearAlgebra.Cholesky}, Type{PositiveFactorizations.Positive}}
#============
General GLM methods
=============#
StatsModels.modelmatrix(m::StatsModels.DataFrameRegressionModel{T}) where T<:INNERMOD = m.mm.m
StatsModels.modelmatrix(m::T) where T<:INNERMOD = m.pp.X

residuals(m::StatsModels.DataFrameRegressionModel{T}) where T<:INNERMOD = m.model.rr.wrkresid
residuals(m::T) where T<:INNERMOD = m.rr.y .- m.rr.mu

modelweights(m::StatsModels.DataFrameRegressionModel{T}) where T<:INNERMOD = m.model.rr.wrkwt
modelweights(m::GLM.GeneralizedLinearModel) = m.rr.wrkwt
modelweights(m::GLM.LinearModel) = m.rr.wts

smplweights(m::StatsModels.DataFrameRegressionModel{T}) where T<:INNERMOD = m.model.rr.wts
smplweights(m::GLM.GeneralizedLinearModel) = m.rr.wts
smplweights(m::GLM.LinearModel) = m.rr.wts

choleskyfactor(m::StatsModels.DataFrameRegressionModel{T}) where T<:INNERMOD = m.model.pp.chol.UL
choleskyfactor(m::T) where T<:INNERMOD = m.pp.chol.UL

unweighted_nobs(m::StatsModels.DataFrameRegressionModel{T}) where T<:INNERMOD = size(modelmatrix(m), 1)
unweighted_nobs(m::T) where T<:INNERMOD = size(modelmatrix(m), 1)

unweighted_dof_residual(m::StatsModels.DataFrameRegressionModel{T}) where T<:INNERMOD = unweighted_nobs(m) - length(coef(m))
unweighted_dof_residual(m::T) where T<:INNERMOD = unweighted_nobs(m) - length(coef(m))

function installxuw!(cache, m::T) where T<:INNERMOD
    copyto!(cache.X, copy(modelmatrix(m)))
    copyto!(cache.u, copy(residuals(m)))
    if !isempty(m.rr.wts)
        broadcast!(*, cache.u, cache.u, smplweights(m))
    end
end

function esteq!(cache, m::StatsModels.RegressionModel, k::T) where T<:Union{HC, CRHC}
    broadcast!(*, cache.q, cache.X, cache.u)
    return cache.q
end

function esteq!(cache, m::RegressionModel, k::HAC)
    X = copy(modelmatrix(m))
    u = copy(residuals(m))
    if !isempty(smplweights(m))
        broadcast!(*, u, u, smplweights(m))
    end
    broadcast!(*, X, X, u)
end

pseudohessian(m::StatsModels.DataFrameRegressionModel{T}) where T<:INNERMOD = GLM.invchol(m.model.pp)
pseudohessian(m::T) where T<:INNERMOD = GLM.invchol(m.pp)

#==============
HAC GLM Methods
===============#

function HACCache(m::StatsModels.DataFrameRegressionModel{T}; prewhiten = false) where T<:INNERMOD
    HACCache(m.model, prewhiten = prewhiten)
end

function HACCache(m::T; prewhiten = false) where T<:INNERMOD
    HACCache(copy(modelmatrix(m)), prewhiten = prewhiten)
end

function StatsBase.vcov(m::StatsModels.DataFrameRegressionModel{T}, k::HAC, args...; kwargs...) where T<:INNERMOD
    set_bw_weights!(k, m)
    vcov(m.model, k, args...; kwargs...)
end

function StatsBase.vcov(m::T, k::K, returntype::Type{T1} = CovarianceMatrix, factortype::Type{T2} = Cholesky; kwargs...) where {T<:INNERMOD, K<:HAC, T1<:Union{CovarianceMatrix, Matrix}, T2<:Union{Nothing, Factorization}}
    cache = HACCache(m, prewhiten = k.prewhiten)
    vcov(m, k, cache, returntype, factortype; kwargs...)
end

function StatsBase.vcov(m::T, k::K, cache::HACCache, returntype::Type{T1} = CovarianceMatrix, factortype::Type{T2} = Cholesky; demean::Bool = false, dof_adjustment::Bool = true) where {T<:INNERMOD, K<:HAC, T1<:Union{CovarianceMatrix, Matrix}, T2<:Union{Nothing, Factorization}}
    mf = esteq!(cache, m, k)
    br = pseudohessian(m)
    set_bw_weights!(k, m)
    Ω = covariance(mf, k, cache, Matrix, demean = demean, scale = 1)
    V = br*Ω*br'
    rmul!(V, size(mf, 1)^2)
    scale = dof_adjustment ? unweighted_dof_residual(m) : unweighted_nobs(m)
    finalize(V, T2, T1, k, scale)
end

# setcholesky!(cache, cholesky::Type{Nothing}, V) = nothing

# function setcholesky!(cache, cholesky::Type{Cholesky}, V)
#     chol = LinearAlgebra.cholesky(Symmetric(V))
#     copyto!(cache.chol.UL.data, chol.UL.data)
#     copyto!(cache.chol.U.data, chol.U.data)
#     copyto!(cache.chol.L.data, chol.L.data)
# end

# function setcholesky!(cache, cholesky::Type{PositiveFactorizations.Positive}, V)
#     chol = LinearAlgebra.cholesky(Positive, V)
#     copyto!(cache.chol.UL.data, chol.UL.data)
#     copyto!(cache.chol.U.data, chol.U.data)
#     copyto!(cache.chol.L.data, chol.L.data)
# end

function set_bw_weights!(k::HAC, m::StatsModels.DataFrameRegressionModel{T}) where T<:INNERMOD
    if isempty(k.weights)
        for j in eachindex(coef(m))
            push!(k.weights, 1.0)
        end
    end
    if length(k.weights) == length(coef(m))
        "(Intercept)" ∈ coefnames(m) ? (k.weights .= 1.0; k.weights[1] = 0.0) : k.weights .= 1.0
    else
        error("Bandwidth weights have wrong dimension for the problem")
    end
end

function set_bw_weights!(k::HAC, m::T) where T<:INNERMOD
    if isempty(k.weights)
        for j in eachindex(coef(m))
            push!(k.weights, 1.0)
        end
    end
end

#==============
HC GLM Methods
===============#
function HCCache(m::StatsModels.DataFrameRegressionModel{T}) where T<:INNERMOD
    HCCache(m.model)
end

function HCCache(m::T) where {T<:INNERMOD}
    n, p = unweighted_nobs(m), length(coef(m))
    HCCache(similar(Array{eltype(modelmatrix(m)), 2}(undef, n, p)))
end

function StatsBase.vcov(m::StatsModels.DataFrameRegressionModel{T},
                        k::K, returntype::Type{T1} = CovarianceMatrix,
                        factortype::Type{T2} = Cholesky;
                        kwargs...) where {T<:INNERMOD,
                                          K<:HC,
                                          T1<:Union{CovarianceMatrix, Matrix},
                                          T2<:Union{Nothing, Factorization}}
    vcov(m.model, k, HCCache(m.model), returntype, factortype; kwargs...)
end

function StatsBase.vcov(m::T,
    k::K, returntype::Type{T1} = CovarianceMatrix,
    factortype::Type{T2} = Cholesky;
    kwargs...) where {T<:INNERMOD,
                      K<:HC,
                      T1<:Union{CovarianceMatrix, Matrix},
                      T2<:Union{Nothing, Factorization}}
    vcov(m, k, HCCache(m), returntype, factortype; kwargs...)
end

# function StatsBase.vcov(m::T,
#                         k::K,
#                         cache::HCCache,
#                         returntype::Type{T1} = CovarianceMatrix,
#                         factortype::Type{T2} = Cholesky; kwargs...) where {T<:INNERMOD,
#                                                                            K<:HC,
#                                                                            T1<:Union{CovarianceMatrix, Matrix},
#                                                                            T2<:Union{Nothing, Factorization}}

#     vcov(m, k, HCCache(m), returntype, factortype; kwargs...)
# end

function StatsBase.vcov(m::T,
                        k::K,
                        cache::HCCache,
                        returntype::Type{T1} = CovarianceMatrix,
                        factortype::Type{T2} = Cholesky;
                        demean::Bool = false) where {T<:INNERMOD,
                                                            K<:HC,
                                                            T1<:Union{CovarianceMatrix, Matrix},
                                                            T2<:Union{Nothing, Factorization}}
    CovarianceMatrices.installxuw!(cache, m)
    mf = CovarianceMatrices.esteq!(cache, m, k)
    br = CovarianceMatrices.pseudohessian(m)
    CovarianceMatrices.adjfactor!(cache, m, k)
    mul!(cache.x, mf', broadcast!(*, cache.X, mf, cache.η))
    V = br*cache.x*br'
    finalize(V, T2, T1, k, 1.0)
end

## -----
## DataFramesRegressionModel/AbstractGLM methods
## -----

function hatmatrix(cache, m::T) where T<:INNERMOD
      z = cache.X  ## THIS ASSUME THAT X IS WEIGHTED BY SQRT(W)
      if !isempty(modelweights(m))
          z .= z.*sqrt.(modelweights(m))
      end
      cf = choleskyfactor(m)::UpperTriangular
      rdiv!(z, cf)
      sum!(cache.v, z.^2)
 end

  adjfactor!(cache, m::StatsModels.RegressionModel, k::HC0) = cache.η .= one(eltype(cache.u))
  adjfactor!(cache, m::StatsModels.RegressionModel, k::HC1) = cache.η .= unweighted_nobs(m)./unweighted_dof_residual(m)
  adjfactor!(cache, m::StatsModels.RegressionModel, k::HC2) = cache.η .= one(eltype(cache.u))./(one(eltype(cache.u)).-hatmatrix(cache, m))
  adjfactor!(cache, m::StatsModels.RegressionModel, k::HC3) = cache.η .= one(eltype(cache.u))./(one(eltype(cache.u)).-hatmatrix(cache, m)).^2

  function adjfactor!(cache, m::StatsModels.RegressionModel, k::HC4)
      n, p = size(cache.X)
      tone = one(eltype(cache.u))
      h = hatmatrix(cache, m)
      @inbounds for j in eachindex(h)
          delta = min(4, n*h[j]/p)
          cache.η[j] = tone/(tone-h[j])^delta
      end
      cache.η
  end

  function adjfactor!(cache, m::StatsModels.RegressionModel, k::HC4m)
      n, p = size(cache.X)
      tone = one(eltype(cache.u))
      h = hatmatrix(cache, m)
      @inbounds for j in eachindex(h)
          delta = min(tone, n*h[j]/p) + min(1.5, n*h[j]/p)
          cache.η[j] = tone/(tone-h[j])^delta
      end
      cache.η
  end

  function adjfactor!(cache, m::StatsModels.RegressionModel, k::HC5)
      n, p = size(cache.X)
      tone = one(eltype(cache.u))
      h     = hatmatrix(cache, m)
      mx    = max(n*0.7*maximum(h)/p, 4)
      @inbounds for j in eachindex(h)
          alpha =  min(n*h[j]/p, mx)
          cache.η[j] = sqrt(tone/(tone-h[j])^alpha)
      end
      cache.η
  end

#==========
Cluster GLM
=========#

function CRHCCache(m::StatsModels.DataFrameRegressionModel{T}, cl::AbstractArray{F}) where {T<:INNERMOD, F}
    CRHCCache(m.model, cl, returntype=returntype)
end

function CRHCCache(m::T, cl::AbstractVector{F}) where {T<:INNERMOD, F}
    n, p = unweighted_nobs(m), length(coef(m))
    CRHCCache(similar(modelmatrix(m)), cl)
end

function StatsBase.vcov(m::StatsModels.DataFrameRegressionModel{T},
                        k::K,
                        returntype::Type{T1} = CovarianceMatrix,
                        factortype::Type{T2} = Cholesky; kwargs...) where {T<:INNERMOD,
                                                                           K<:CRHC,
                                                                           T1<:Union{CovarianceMatrix, Matrix},
                                                                           T2<:Union{Nothing, Factorization}}
    vcov(m.model, k, CRHCCache(m.model, k.cl), returntype, factortype; kwargs...)
end

function StatsBase.vcov(m::T,
                        k::K,
                        returntype::Type{T1} = CovarianceMatrix,
                        factortype::Type{T2} = Cholesky;
                        kwargs...) where {T<:INNERMOD,
                                          K<:CRHC,
                                          T1<:Union{CovarianceMatrix, Matrix},
                                          T2<:Union{Nothing, Factorization}}
    vcov(m, k, CRHCCache(m, k.cl), returntype, factortype; kwargs...)
end

function StatsBase.vcov(m::T,
                        k::CRHC,
                        cache::CRHCCache,
                        returntype::Type{T1} = CovarianceMatrix,
                        factortype::Type{T2} = Cholesky;
                        sorted::Bool = false) where {T<:INNERMOD,
                                                     K<:HC,
                                                     T1<:Union{CovarianceMatrix, Matrix},
                                                     T2<:Union{Nothing, Factorization}}
    B = CovarianceMatrices.pseudohessian(m)
    CovarianceMatrices.installsortedxuw!(cache, m, k, Val{sorted})
    bstarts = (searchsorted(cache.clus, j[2]) for j in enumerate(unique(cache.clus)))
    CovarianceMatrices.adjresid!(k, cache, B, bstarts)
    CovarianceMatrices.esteq!(cache, m, k)
    CovarianceMatrices.clusterize!(cache, bstarts)
    V = dof_adjustment(cache, k, bstarts).*(B*cache.M*B)
    finalize(V, T2, T1, k, 1.0)
end

StatsBase.stderror(m::StatsModels.DataFrameRegressionModel{T}, k::RobustVariance; kwargs...) where T<:INNERMOD = sqrt.(diag(vcov(m, k, Matrix; kwargs...)))


text:
function DSSTSolarRadiationPressure(arg0::ExtendedPVCoordinatesProvider, arg1::jdouble, arg2::RadiationSensitive, arg3::jdouble)
    return DSSTSolarRadiationPressure((ExtendedPVCoordinatesProvider, jdouble, RadiationSensitive, jdouble), arg0, arg1, arg2, arg3)
end

function DSSTSolarRadiationPressure(arg0::jdouble, arg1::jdouble, arg2::ExtendedPVCoordinatesProvider, arg3::jdouble, arg4::RadiationSensitive, arg5::jdouble)
    return DSSTSolarRadiationPressure((jdouble, jdouble, ExtendedPVCoordinatesProvider, jdouble, RadiationSensitive, jdouble), arg0, arg1, arg2, arg3, arg4, arg5)
end

function DSSTSolarRadiationPressure(arg0::jdouble, arg1::jdouble, arg2::ExtendedPVCoordinatesProvider, arg3::jdouble, arg4::jdouble)
    return DSSTSolarRadiationPressure((jdouble, jdouble, ExtendedPVCoordinatesProvider, jdouble, jdouble), arg0, arg1, arg2, arg3, arg4)
end

function DSSTSolarRadiationPressure(arg0::jdouble, arg1::jdouble, arg2::jdouble, arg3::jdouble, arg4::ExtendedPVCoordinatesProvider, arg5::jdouble, arg6::jdouble)
    return DSSTSolarRadiationPressure((jdouble, jdouble, jdouble, jdouble, ExtendedPVCoordinatesProvider, jdouble, jdouble), arg0, arg1, arg2, arg3, arg4, arg5, arg6)
end

function get_equatorial_radius(obj::DSSTSolarRadiationPressure)
    return jcall(obj, "getEquatorialRadius", jdouble, ())
end

function get_events_detectors(obj::DSSTSolarRadiationPressure)
    return jcall(obj, "getEventsDetectors", Vector{EventDetector}, ())
end

function get_field_events_detectors(obj::DSSTSolarRadiationPressure, arg0::Field)
    return jcall(obj, "getFieldEventsDetectors", Vector{FieldEventDetector}, (Field,), arg0)
end

function get_spacecraft(obj::DSSTSolarRadiationPressure)
    return jcall(obj, "getSpacecraft", RadiationSensitive, ())
end


text:
#    Copyright (c) 2021 J.A. Duffek
#    Copyright (c) 2000 D.M. Spink
#
#    This program is free software; you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation; either version 2 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program; if not, see <http://www.gnu.org/licenses/>.

"""
    vecrotz(angle::F)::Matrix{F} where {F<:AbstractFloat}

Transformation matrix for a rotation around the z axis.

Return the (4x4) Transformation matrix for a rotation about the z axis by the
defined angle.

The matrix is:

    [  cos(angle)  -sin(angle)       0          0]
    [ -sin(angle)   cos(angle)       0          0]
    [      0            0            1          0]
    [      0            0            0          1]

# Arguments:
- `angle`: rotation angle defined in radians

# Output:
- `rz`: (4x4) Transformation matrix.

# Examples:
```julia
julia> rz = vecrotz(angle)
```

Rotate the NURBS line `[0.0;0.0;0.0]` - `[3.0;3.0;3.0]` by 45 degrees around
the z-axis

```julia
julia> line = nrbline(vec([0.0 0.0 0.0]),vec([3.0 3.0 3.0]));
julia> rotate = vecrotz(pi/4);
julia> rline = nrbtform(line, rotate)
```

# Reference
- see also: [`nrbtform`](@ref)
"""
function vecrotz(angle::F)::Matrix{F} where {F<:AbstractFloat}
sn = sin(angle);
cn = cos(angle);
rz = [ cn -sn   0   0;
       sn  cn   0   0;
        0   0   1   0;
        0   0   0   1];
return rz;
end # vecrotz


text:
using Test

include("adder.jl")
include("divider.jl")

@testset "All the tests" begin
    @test adder(2, 5) == 7
    @test divider(10, 3) ≈ 3.333333333
end

max_stars_repo_name: wigging/julia-computing

max_issues_repo_name: wigging/julia-computing

max_forks_repo_name: wigging/julia-computing


text:
SRC_DIR = joinpath(Base.source_dir(), "..", "..", "src")
OUT_DIR = joinpath(Base.source_dir(), "..", "..", "output")
mkpath(OUT_DIR)
#NOTE additional packages that need installing:
#=
    CSV, DataFrames, RollingFunctions

=#

#include(joinpath(SRC_DIR, "BridgeSDEInference.jl"))
#using Main.BridgeSDEInference
include(joinpath(SRC_DIR, "BridgeSDEInference_for_tests.jl"))
#using BridgeSDEInference
#const BSI = BridgeSDEInference
using Distributions # to define priors
using Random        # to seed the random number generator
using DataFrames
using CSV
using LinearAlgebra

DIR = "auxiliary"
include(joinpath(SRC_DIR, DIR, "read_and_write_data.jl"))
include(joinpath(SRC_DIR, DIR, "transforms.jl"))
include(joinpath(SRC_DIR, DIR, "utility_functions.jl"))

# decide if first passage time observations or partially observed diffusion
fptObsFlag = false

# simulate dataset
include(joinpath(SRC_DIR, "..", "scripts", "data_generation",
                 "simulate_repeated_fhn.jl"))

param = :complexConjug
θ_init = [10.0, -8.0, 15.0, 0.0, 3.0]
P˟ = FitzhughDiffusion(param, θ_init...)
P̃s = [std_aux_laws(FitzhughDiffusionAux, param, θ_init, map(λ->λ[1], xx.yy), xx.tt, 1) for xx in XX]
L = @SMatrix [1. 0.]
Σ = @SMatrix [10^(-10)]

model_setups = [DiffusionSetup(P˟, P̃, PartObs()) for P̃ in P̃s]
for (k,(model_setup, P̃, xx))  in enumerate(zip(model_setups, P̃s, XX))
      obs, obs_times = map(λ->λ[1], xx.yy), xx.tt
      set_observations!(model_setup, [L for _ in P̃], [Σ for _ in P̃], obs, obs_times)
      set_imputation_grid!(model_setup, 1/1000)
      set_x0_prior!(model_setup, KnownStartingPt(x0[k]))
      initialise!(eltype(x0[k]), model_setup, Vern7(), false, NoChangePt(100))
      set_auxiliary!(model_setup; skip_for_save=10^0, adaptive_prop=NoAdaptation())
end

mcmc_setup = MCMCSetup(
      Imputation(NoBlocking(), 0.975, Vern7()),
      ParamUpdate(MetropolisHastingsUpdt(), 5, θ_init,
                  UniformRandomWalk(0.5, true), ImproperPosPrior(),
                  UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃s[1], 5))
                  ),
      ParamUpdate(ConjugateUpdt(), [1,2,3], θ_init, nothing,
                  MvNormal(fill(0.0, 3), diagm(0=>fill(1000.0, 3))),
                  UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃s[1], [1,2,3]))
                  ))

schedule = MCMCSchedule(5*10^4, [[1,2,3]],
                        (save=3*10^2, verbose=10^2, warm_up=100,
                         readjust=(x->x%100==0), fuse=(x->false)))

Random.seed!(4)
out, elapsed = @timeit mcmc(mcmc_setup, schedule, model_setups)
#out = mcmc(mcmc_setup, schedule, model_setup)

include(joinpath(SRC_DIR, DIR, "plotting_fns.jl"))
plot_chains(out[2]; truth=[10.0, -8.0, 15.0, 0.0, 3.0])
plot_paths(out[1], out[2], schedule; transf=[(x,θ)->x, (x,θ)->conjugToRegular(x, θ[1], 0)],
           obs=(times=obs_time[2:end], vals=obs[2:end], indices=1))


text:
# ---
# title: 1343. Number of Sub-arrays of Size K and Average Greater than or Equal to Threshold
# id: problem1343
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array
# link: <https://leetcode.com/problems/number-of-sub-arrays-of-size-k-and-average-greater-than-or-equal-to-threshold/description/>
# hidden: true
# ---
#
# Given an array of integers `arr` and two integers `k` and `threshold`.
#
# Return _the number of sub-arrays_ of size `k` and average greater than or
# equal to `threshold`.
#
#
#
# **Example 1:**
#
#
#
#     Input: arr = [2,2,2,2,5,5,5,8], k = 3, threshold = 4
#     Output: 3
#     Explanation: Sub-arrays [2,5,5],[5,5,5] and [5,5,8] have averages 4, 5 and 6 respectively. All other sub-arrays of size 3 have averages less than 4 (the threshold).
#
#
# **Example 2:**
#
#
#
#     Input: arr = [1,1,1,1,1], k = 1, threshold = 0
#     Output: 5
#
#
# **Example 3:**
#
#
#
#     Input: arr = [11,13,17,23,29,31,7,5,2,3], k = 3, threshold = 5
#     Output: 6
#     Explanation: The first 6 sub-arrays of size 3 have averages greater than 5. Note that averages are not integers.
#
#
# **Example 4:**
#
#
#
#     Input: arr = [7,7,7,7,7,7,7], k = 7, threshold = 7
#     Output: 1
#
#
# **Example 5:**
#
#
#
#     Input: arr = [4,4,4,4], k = 4, threshold = 1
#     Output: 1
#
#
#
#
# **Constraints:**
#
#   * `1 <= arr.length <= 10^5`
#   * `1 <= arr[i] <= 10^4`
#   * `1 <= k <= arr.length`
#   * `0 <= threshold <= 10^4`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end