text:
using LinearAlgebra
using Roots
using Distances

# Shortcuts (already defined)
# MDI = DynamicalSystemsBase.MinimalDiscreteIntegrator
# DEI = DynamicalSystemsBase.DiffEqBase.DEIntegrator

# TODO: Add example with rectangle, perhaps not clear

export transit_time_statistics, transit_return

"""
    transit_time_statistics(ds::DynamicalSystem, u₀, εs, T; diffeq...) → exits, entries
Collect transit time statistics for a ball/box centered at `u₀` with radii `εs` (see below),
in the state space of the given dynamical system.
Return the exit and (re-)entry return times to the set(s), where each of these is a vector
containing all collected times for the respective `ε`-radius set, for `ε ∈ εs`.

Use `transit_return(exits, entries)` to transform the result to transit and return time
statistics instead.

## Keywords
* `diffeq...`: All keywords are propagated to the integrators of DifferentialEquations.jl
  for continuous systems (discrete have no keywords).

## Description
Transit time statistics are important for the transport properties of dynamical systems[^Meiss1997]
and can even be connected with the fractal dimension of chaotic sets[^Boev2014].

The current algorithm collects exit and re-entry times to given sets in the state space,
which are centered at `u₀` (**algorithm always starts at `u₀`** and the initial state of
`ds` is irrelevant). `εs` is always a `Vector`.

The sets around `u₀` are nested hyper-spheres of radius `ε ∈ εs`, if each entry of
`εs` is a real number. The sets can also be
hyper-rectangles (boxes), if each entry of `εs` is a vector itself.
Then, the `i`-th box is defined by the space covered by `u0 .± εs[i]` (thus the actual
box size is `2εs[i]`!).

The reason to input multiple `εs` at once is purely for performance.

For discrete systems, exit time is recorded immediatelly after exitting of the set, and
re-entry is recorded immediatelly on re-entry. This means that if an orbit needs
1 step to leave the set and then it re-enters immediatelly on the next step, the return time
is 1. For continuous systems high-order
interpolation is done to accurately record the time of exactly crossing the `ε`-ball/box.

[^Meiss1997]: Meiss, J. D. *Average exit time for volume-preserving maps*, Chaos (1997)](https://doi.org/10.1063/1.166245)

[^Boev2014]: Boev, Vadivasova, & Anishchenko, *Poincaré recurrence statistics as an indicator of chaos synchronization*, Chaos (2014)](https://doi.org/10.1063/1.4873721)
"""
function transit_time_statistics end

"""
    transit_return(exits, entries) → transit, return
Convert the output of [`return_time_statistics`](@ref) to the transit and return times.
"""
function transit_return(exits, entries)
    # the main reason this function exists is because entry times can be one less
    # than entry times. (otherwise you could just directly subtract)
    returns = [en .- view(ex, 1:length(en)) for (en, ex) in zip(entries, exits)]
    transits = similar(entries)
    for (j, (en, ex)) in enumerate(zip(entries, exits))
        M, N = length(en), length(ex)
        enr, exr = M == N ? (1:N-1, 2:N) : (1:M, 2:N)
        transits[j] = view(ex, exr) .- view(en, enr)
    end
    return transits, returns
end

##########################################################################################
# ε-distances
##########################################################################################
function check_εs_sorting(εs, L)
    correct = if εs[1] isa Real
        issorted(εs; rev = true)
    elseif εs[1] isa AbstractVector
        @assert all(e -> length(e) == L, εs) "Boxes must have same dimension as state space!"
        for j in 1:L
            if !issorted([εs[i][j] for i in 1:length(εs)]; rev = true)
                return false
            end
        end
        true
    end
    if !correct
        throw(ArgumentError("`εs` must be sorted from largest to smallest ball/box size."))
    end
    return correct
end

# Support both types of sets: balls and boxes
"Return `true` if state is outside ε-ball"
function isoutside(u, u0, ε::AbstractVector)
    @inbounds for i in 1:length(u)
        abs(u[i] - u0[i]) > ε[i] && return true
    end
    return false
end
isoutside(u, u0, ε::Real) = euclidean(u, u0) > ε

"Return the **signed** distance of state to ε-ball (negative means inside ball)"
function εdistance(u, u0, ε::AbstractVector)
    m = eltype(u)(-Inf)
    @inbounds for i in 1:length(u)
        m2 = abs(u[i] - u0[i]) - ε[i]
        m2 > m && (m = m2)
    end
    return m
end
εdistance(u, u0, ε::Real) = euclidean(u, u0) - ε

##########################################################################################
# Discrete systems
##########################################################################################
function transit_time_statistics(ds::DiscreteDynamicalSystem, u0, εs, T; diffeq...)
    check_εs_sorting(εs, length(u0))
    integ = integrator(ds, u0)
    transit_time_statistics(integ, u0, εs, T)
end

function transit_time_statistics(integ::MDI, u0, εs, T)
    E = length(εs)
    pre_outside = fill(false, E) # `true` if outside the ball. Previous step
    cur_outside = copy(pre_outside)       # current step.
    exits   = [Int[] for _ in 1:E]
    entries = [Int[] for _ in 1:E]

    while integ.t < T
        step!(integ)

        # here i gives the index of the largest ε-ball that the trajectory is out of.
        # It is guaranteed that the trajectory is thus outside all other boxes
        i = first_outside_index(integ, u0, εs, E) # TODO: Continuous version
        cur_outside[i:end] .= true
        cur_outside[1:i-1] .= false

        update_exit_times!(exits, i, pre_outside, cur_outside, integ)
        update_entry_times!(entries, i, pre_outside, cur_outside, integ)
        pre_outside .= cur_outside
    end
    return exits, entries
end

function first_outside_index(integ::MDI, u0, εs, E)::Int
    i = findfirst(e -> isoutside(integ.u, u0, e), εs)
    return isnothing(i) ? E+1 : i
end

function update_exit_times!(exits, i, pre_outside, cur_outside, integ::MDI)
    @inbounds for j in i:length(pre_outside)
        cur_outside[j] && !pre_outside[j] && push!(exits[j], integ.t)
    end
end

function update_entry_times!(entries, i, pre_outside, cur_outside, integ::MDI)
    # TODO: Can I use `i` here?
    @inbounds for j in 1:length(pre_outside)
        pre_outside[j] && !cur_outside[j] && push!(entries[j], integ.t)
    end
end

##########################################################################################
# Continuous
##########################################################################################
# using OrdinaryDiffEq: Tsit5
using DynamicalSystemsBase.DiffEqBase: ODEProblem, solve
using DynamicalSystemsBase.DiffEqBase: ContinuousCallback, CallbackSet

function transit_time_statistics(ds::ContinuousDynamicalSystem, u0, εs, T;
        diffeq...
    )
    # alg = Tsit5()

    error("Continuous system version is not yet ready.")

    eT = eltype(ds.t0)
    check_εs_sorting(εs, length(u0))
    exits = [eT[] for _ in 1:length(εs)]
    entries = [eT[] for _ in 1:length(εs)]

    # Make the magic callbacks:
    callbacks = ContinuousCallback[]
    for i in eachindex(εs)
        crossing(u, t, integ) = ChaosTools.εdistance(u, u0, εs[i])
        negative_affect!(integ) = push!(entries[i], integ.t)
        positive_affect!(integ) = push!(exits[i], integ.t)
        cb = ContinuousCallback(crossing, positive_affect!, negative_affect!;
            save_positions = (false, false)
        )
        push!(callbacks, cb)
    end
    cb = CallbackSet(callbacks...)

    prob = ODEProblem(ds, (eT(0), eT(T)); u0 = u0)
    sol = solve(prob, alg;
        callback=cb, save_everystep = false, dense = false,
        save_start=false, save_end = false, diffeq...
    )
    return exits, entries
end


text:
for (i, f) ∈ enumerate((
  :second,
  :third,
  :fourth,
  :fifth,
  :sixth,
  :seventh,
  :eighth,
  :ninth,
  :tenth,
  :eleventh,
  :twelfth,
  :thirteenth,
))
  @eval @inline $f(x) = @inbounds getindex(x, $(i + 1))
end

ext: jl


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

mutable struct ContainerServiceSshConfiguration <: SwaggerModel
    publicKeys::Any # spec type: Union{ Nothing, Vector{ContainerServiceSshPublicKey} } # spec name: publicKeys

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

const _property_map_ContainerServiceSshConfiguration = Dict{Symbol,Symbol}(Symbol("publicKeys")=>Symbol("publicKeys"))
const _property_types_ContainerServiceSshConfiguration = Dict{Symbol,String}(Symbol("publicKeys")=>"Vector{ContainerServiceSshPublicKey}")
Base.propertynames(::Type{ ContainerServiceSshConfiguration }) = collect(keys(_property_map_ContainerServiceSshConfiguration))
Swagger.property_type(::Type{ ContainerServiceSshConfiguration }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ContainerServiceSshConfiguration[name]))}
Swagger.field_name(::Type{ ContainerServiceSshConfiguration }, property_name::Symbol) =  _property_map_ContainerServiceSshConfiguration[property_name]

function check_required(o::ContainerServiceSshConfiguration)
    (getproperty(o, Symbol("publicKeys")) === nothing) && (return false)
    true
end

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


text:
"""
    T1A_SS!(ψₒ, ψᵢ, dx, ops)

Compute `ψₒ`, i.e. `ψᵢ` advanced a step `dx` forward using a symplectic first order integrator of type A for the Sasa-Satsuma equation. The structure `ops::Operators` contains the FFT plans and the kinetic energy operators.

See also: [`solve!`](@ref), [`Operators`](@ref)
"""
function T1A_SS!(ψₒ, ψᵢ, dx, ops)
    # Nonlinear
    @. ψₒ = cis(-dx * (-1*abs2(ψᵢ)))*ψᵢ

    # Dispersion
    ops.F̂*ψₒ
    ψₒ .= ops.K̂(dx) .* ψₒ
    ops.F̃̂*ψₒ

    # Burger
    set_u!(ops.B̂, ψₒ)
    step!(ops.B̂)
    ψₒ .= ops.B̂.u

    set_u!(ops.B̂B, ψₒ)
    step!(ops.B̂B)
    ψₒ .= ops.B̂B.u
end #T₁ʰ


text:
@testset "8.10 Formal derivatives" begin
    (F, a, b, m, n, x, ) = @variables F a b m n x

    #= ::Package:: =#

    #= ::Title:: =#
    #=Differentiation*Function*Integration*Problems=#

    @test_int [f'(x), x, 1, f(x)]
    @test_int [f''(x), x, 1, f'(x)]
    @test_int [f'''(x), x, 1, f''(x)]
    @test_int [Derivative(n)(f)(x), x, 1, Derivative(n-1)(f)(x)]
    @test_int [Derivative(1)(u)(x)*Derivative(2)(u)(x), x, 2, Derivative(1)(u)(x)^2/2]

    @test_int [f'(x)/f(x), x, 2, log(f(x))]
    @test_int [f'(x)/(a+b*f(x)), x, 2, log(a+b*f(x))/b]

    @test_int [f(x)*f'(x), x, 2, f(x)^2/2]
    @test_int [(a+b*f(x))*f'(x), x, 2, a*f(x) + (1/2)*b*f(x)^2]

    @test_int [f'(x)/sqrt(f(x)), x, 2, 2*sqrt(f(x))]
    @test_int [f'(x)/sqrt(a+b*f(x)), x, 2, (2*sqrt(a+b*f(x)))/b]

    @test_int [f(x)^n*f'(x), x, 2, f(x)^(n+1)/(n+1)]
    @test_int [(a+b*f(x))^n*f'(x), x, 2, (a+b*f(x))^(n+1)/(b*(n+1))]

    @test_int [f''(x)/f'(x), x, 2, log(f'(x))]
    @test_int [f''(x)/(a+b*f'(x)), x, 2, log(a+b*f'(x))/b]

    @test_int [f'(x)*f''(x), x, 2, f'(x)^2/2]
    @test_int [(a+b*f'(x))*f''(x), x, 2, a*f'(x) + (1/2)*b*f'(x)^2]

    @test_int [f''(x)/sqrt(f'(x)), x, 2, 2*sqrt(f'(x))]
    @test_int [f''(x)/sqrt(a+b*f'(x)), x, 2, (2*sqrt(a+b*f'(x)))/b]

    @test_int [f'(x)^n*f''(x), x, 2, f'(x)^(n+1)/(n+1)]
    @test_int [(a+b*f'(x))^n*f''(x), x, 2, (a+b*f'(x))^(n+1)/(b*(n+1))]

    @test_int [f(g(x))*g'(x), x, 1, CannotIntegrate(f(g(x))*g'(x), x)]
    @test_int [f(g'(x))*g''(x), x, 1, CannotIntegrate(f(g'(x))*g''(x), x)]

    @test_int [f'(x)*g(x) + f(x)*g'(x), x, -1, f(x)*g(x)]
    @test_int [(f'(x)*g(x) - f(x)*g'(x))/g(x)^2, x, 1, f(x)/g(x)]
    @test_int [(f'(x)*g(x) - f(x)*g'(x))/(f(x)*g(x)), x, 1, log(f(x)/g(x))]

    @test_int [(g(x)*Derivative(1)(f)(x) + f(x)*g'(x))/(1 + f(x)^2*g(x)^2), x, 2, atan(f(x)*g(x))]
    @test_int [(g(x)*Derivative(1)(f)(x) - f(x)*g'(x))/(f(x)^2 + g(x)^2), x, 2, atan(f(x)/g(x))]
    @test_int [-((g(x)*Derivative(1)(f)(x) + f(x)*g'(x))/(1 + f(x)^2*g(x)^2)), x, 2, -atan(f(x)*g(x))]

    @test_int [(g(x)*Derivative(1)(f)(x) + f(x)*g'(x))/(1 - f(x)^2*g(x)^2), x, 2, atanh(f(x)*g(x))]
    @test_int [((-g(x))*Derivative(1)(f)(x) + f(x)*g'(x))/(f(x)^2 - g(x)^2), x, 2, atanh(f(x)/g(x))]

    @test_int [(f(x)^(-1 + m)*g(x)^(-1 + n)*(m*g(x)*Derivative(1)(f)(x) + n*f(x)*g'(x)))/(1 - f(x)^(2*m)*g(x)^(2*n)), x, 2, atanh(f(x)^m*g(x)^n)]
    @test_int [(f(x)^(-1 + m)*g(x)^(-1 + n)*((-m)*g(x)*Derivative(1)(f)(x) + n*f(x)*g'(x)))/(f(x)^(2*m) - g(x)^(2*n)), x, 3, atanh(g(x)^n/f(x)^m)]
    @test_int [(f(x)^(-1 + m)*g(x)^(-1 - n)*((-m)*g(x)*Derivative(1)(f)(x) - n*f(x)*g'(x)))/(f(x)^(2*m) - g(x)^(-2*n)), x, 3, atanh(1/(f(x)^m*g(x)^n))]

    @test_int [(g(x)*Derivative(1)(f)(x) + f(x)*g'(x))/(a + b*f(x)^1*g(x)^1), x, 2, log(a + b*f(x)*g(x))/b]
    @test_int [(g(x)*Derivative(1)(f)(x) + f(x)*g'(x))/(a + b*f(x)^2*g(x)^2), x, 2, atan((sqrt(b)*f(x)*g(x))/sqrt(a))/(sqrt(a)*sqrt(b))]
    @test_int [(g(x)*Derivative(1)(f)(x) + f(x)*g'(x))/(a + b*f(x)^3*g(x)^3), x, 7, -(atan((a^(1/3) - 2*b^(1/3)*f(x)*g(x))/(sqrt(3)*a^(1/3)))/(sqrt(3)*a^(2/3)*b^(1/3))) + log(a^(1/3) + b^(1/3)*f(x)*g(x))/(3*a^(2/3)*b^(1/3)) - log(a^(2/3) - a^(1/3)*b^(1/3)*f(x)*g(x) + b^(2/3)*f(x)^2*g(x)^2)/(6*a^(2/3)*b^(1/3))]

    @test_int [(g(x)*Derivative(1)(f)(x) + f(x)*g'(x))/(a + b*(f(x)*g(x))^(1/2)), x, 4, (2*sqrt(f(x)*g(x)))/b - (2*a*log(a + b*sqrt(f(x)*g(x))))/b^2]
    @test_int [(g(x)*Derivative(1)(f)(x) + f(x)*g'(x))/(a + b*(f(x)*g(x))^(3/2)), x, 8, -((2*atan((a^(1/3) - 2*b^(1/3)*sqrt(f(x)*g(x)))/(sqrt(3)*a^(1/3))))/(sqrt(3)*a^(1/3)*b^(2/3))) - (2*log(a^(1/3) + b^(1/3)*sqrt(f(x)*g(x))))/(3*a^(1/3)*b^(2/3)) + log(a^(2/3) + b^(2/3)*f(x)*g(x) - a^(1/3)*b^(1/3)*sqrt(f(x)*g(x)))/(3*a^(1/3)*b^(2/3))]
    @test_int [(g(x)*Derivative(1)(f)(x) + f(x)*g'(x))/(a + b*(f(x)*g(x))^(5/2)), x, 8, -((sqrt(2*(5 - sqrt(5)))*atan(sqrt((1/5)*(5 - 2*sqrt(5))) + (2*sqrt(2/(5 + sqrt(5)))*b^(1/5)*sqrt(f(x)*g(x)))/a^(1/5)))/(5*a^(3/5)*b^(2/5))) - (sqrt(2*(5 + sqrt(5)))*atan(sqrt((1/5)*(5 + 2*sqrt(5))) - (sqrt((2/5)*(5 + sqrt(5)))*b^(1/5)*sqrt(f(x)*g(x)))/a^(1/5)))/(5*a^(3/5)*b^(2/5)) - (2*log(a^(1/5) + b^(1/5)*sqrt(f(x)*g(x))))/(5*a^(3/5)*b^(2/5)) + ((1 - sqrt(5))*log(2*a^(2/5) + 2*b^(2/5)*f(x)*g(x) - a^(1/5)*b^(1/5)*sqrt(f(x)*g(x)) - sqrt(5)*a^(1/5)*b^(1/5)*sqrt(f(x)*g(x))))/(10*a^(3/5)*b^(2/5)) + ((1 + sqrt(5))*log(2*a^(2/5) + 2*b^(2/5)*f(x)*g(x) - a^(1/5)*b^(1/5)*sqrt(f(x)*g(x)) + sqrt(5)*a^(1/5)*b^(1/5)*sqrt(f(x)*g(x))))/(10*a^(3/5)*b^(2/5))]

    @test_int [(g(x)*Derivative(1)(f)(x) + f(x)*g'(x))/(a + b*(f(x)*g(x))^n), x, 2, (f(x)*g(x)*HypergeometricFunctions._₂F₁(1, 1/n, 1 + 1/n, -((b*(f(x)*g(x))^n)/a)))/a]
    @test_int [(g(x)*Derivative(1)(f)(x) + f(x)*g'(x))/(a + b*f(x)^n*g(x)^n), x, 2, CannotIntegrate((g(x)*Derivative(1)(f)(x))/(a + b*f(x)^n*g(x)^n), x) + CannotIntegrate((f(x)*g'(x))/(a + b*f(x)^n*g(x)^n), x)]

    @test_int [cos(x)*g(ℯ^x)*Derivative(1)(f)(sin(x)) + ℯ^x*f(sin(x))*Derivative(1)(g)(ℯ^x), x, -1, f(sin(x))*g(ℯ^x)]

    @test_int [F^(a + b*x)*Derivative(3)(f)(x), x, 3, b^2*F^(a + b*x)*f(x)*log(F)^2 - b^3*CannotIntegrate(F^(a + b*x)*f(x), x)*log(F)^3 - b*F^(a + b*x)*log(F)*Derivative(1)(f)(x) + F^(a + b*x)*Derivative(2)(f)(x)]
    @test_int [F^(a + b*x)*Derivative(2)(f)(x), x, 2, (-b)*F^(a + b*x)*f(x)*log(F) + b^2*CannotIntegrate(F^(a + b*x)*f(x), x)*log(F)^2 + F^(a + b*x)*Derivative(1)(f)(x)]
    @test_int [F^(a + b*x)*Derivative(1)(f)(x), x, 1, F^(a + b*x)*f(x) - b*CannotIntegrate(F^(a + b*x)*f(x), x)*log(F)]
    @test_int [F^(a + b*x)*Derivative(0)(f)(x), x, 0, CannotIntegrate(F^(a + b*x)*f(x), x)]
    @test_int [F^(a + b*x)*Derivative(-1)(f)(x), x, 1, -(CannotIntegrate(F^(a + b*x)*f(x), x)/(b*log(F))) + (F^(a + b*x)*Derivative(-1)(f)(x))/(b*log(F))]
    @test_int [F^(a + b*x)*Derivative(-2)(f)(x), x, 2, CannotIntegrate(F^(a + b*x)*f(x), x)/(b^2*log(F)^2) + (F^(a + b*x)*Derivative(-2)(f)(x))/(b*log(F)) - (F^(a + b*x)*Derivative(-1)(f)(x))/(b^2*log(F)^2)]
    @test_int [F^(a + b*x)*Derivative(-3)(f)(x), x, 3, -(CannotIntegrate(F^(a + b*x)*f(x), x)/(b^3*log(F)^3)) + (F^(a + b*x)*Derivative(-3)(f)(x))/(b*log(F)) - (F^(a + b*x)*Derivative(-2)(f)(x))/(b^2*log(F)^2) + (F^(a + b*x)*Derivative(-1)(f)(x))/(b^3*log(F)^3)]

    @test_int [F^(a + b*x)*Derivative(3)(f)(x) + b^3*F^(a + b*x)*f(x)*log(F)^3, x, 4, b^2*F^(a + b*x)*f(x)*log(F)^2 - b*F^(a + b*x)*log(F)*Derivative(1)(f)(x) + F^(a + b*x)*Derivative(2)(f)(x)]

    @test_int [sin(a + b*x)*Derivative(3)(f)(x), x, 3, b^3*CannotIntegrate(cos(a + b*x)*f(x), x) - b^2*f(x)*sin(a + b*x) - b*cos(a + b*x)*Derivative(1)(f)(x) + sin(a + b*x)*Derivative(2)(f)(x)]
    @test_int [sin(a + b*x)*Derivative(2)(f)(x), x, 2, (-b)*cos(a + b*x)*f(x) - b^2*CannotIntegrate(f(x)*sin(a + b*x), x) + sin(a + b*x)*Derivative(1)(f)(x)]
    @test_int [sin(a + b*x)*Derivative(1)(f)(x), x, 1, (-b)*CannotIntegrate(cos(a + b*x)*f(x), x) + f(x)*sin(a + b*x)]
    @test_int [sin(a + b*x)*Derivative(0)(f)(x), x, 0, CannotIntegrate(f(x)*sin(a + b*x), x)]
    @test_int [sin(a + b*x)*Derivative(-1)(f)(x), x, 1, CannotIntegrate(cos(a + b*x)*f(x), x)/b - (cos(a + b*x)*Derivative(-1)(f)(x))/b]
    @test_int [sin(a + b*x)*Derivative(-2)(f)(x), x, 2, -(CannotIntegrate(f(x)*sin(a + b*x), x)/b^2) - (cos(a + b*x)*Derivative(-2)(f)(x))/b + (sin(a + b*x)*Derivative(-1)(f)(x))/b^2]
    @test_int [sin(a + b*x)*Derivative(-3)(f)(x), x, 3, -(CannotIntegrate(cos(a + b*x)*f(x), x)/b^3) - (cos(a + b*x)*Derivative(-3)(f)(x))/b + (sin(a + b*x)*Derivative(-2)(f)(x))/b^2 + (cos(a + b*x)*Derivative(-1)(f)(x))/b^3]

    @test_int [sin(a + b*x)*Derivative(3)(f)(x) - b^3*cos(a + b*x)*f(x), x, 4, (-b^2)*f(x)*sin(a + b*x) - b*cos(a + b*x)*Derivative(1)(f)(x) + sin(a + b*x)*Derivative(2)(f)(x)]

    @test_int [cos(a + b*x)*Derivative(3)(f)(x), x, 3, (-b^2)*cos(a + b*x)*f(x) - b^3*CannotIntegrate(f(x)*sin(a + b*x), x) + b*sin(a + b*x)*Derivative(1)(f)(x) + cos(a + b*x)*Derivative(2)(f)(x)]
    @test_int [cos(a + b*x)*Derivative(2)(f)(x), x, 2, (-b^2)*CannotIntegrate(cos(a + b*x)*f(x), x) + b*f(x)*sin(a + b*x) + cos(a + b*x)*Derivative(1)(f)(x)]
    @test_int [cos(a + b*x)*Derivative(1)(f)(x), x, 1, cos(a + b*x)*f(x) + b*CannotIntegrate(f(x)*sin(a + b*x), x)]
    @test_int [cos(a + b*x)*Derivative(0)(f)(x), x, 0, CannotIntegrate(cos(a + b*x)*f(x), x)]
    @test_int [cos(a + b*x)*Derivative(-1)(f)(x), x, 1, -(CannotIntegrate(f(x)*sin(a + b*x), x)/b) + (sin(a + b*x)*Derivative(-1)(f)(x))/b]
    @test_int [cos(a + b*x)*Derivative(-2)(f)(x), x, 2, -(CannotIntegrate(cos(a + b*x)*f(x), x)/b^2) + (sin(a + b*x)*Derivative(-2)(f)(x))/b + (cos(a + b*x)*Derivative(-1)(f)(x))/b^2]
    @test_int [cos(a + b*x)*Derivative(-3)(f)(x), x, 3, CannotIntegrate(f(x)*sin(a + b*x), x)/b^3 + (sin(a + b*x)*Derivative(-3)(f)(x))/b + (cos(a + b*x)*Derivative(-2)(f)(x))/b^2 - (sin(a + b*x)*Derivative(-1)(f)(x))/b^3]

    @test_int [cos(a + b*x)*Derivative(-2)(f)(x) + cos(a + b*x)*(f(x)/b^2), x, 3, (sin(a + b*x)*Derivative(-2)(f)(x))/b + (cos(a + b*x)*Derivative(-1)(f)(x))/b^2]

    @test_int [cos(f(x)*g(x))*(g(x)*Derivative(1)(f)(x) + f(x)*g'(x)), x, 2, sin(f(x)*g(x))]
    @test_int [cos(g(x)*Derivative(m)(f)(x))*(g'(x)*Derivative(m)(f)(x) + g(x)*Derivative(1 + m)(f)(x)), x, 2, sin(g(x)*Derivative(m)(f)(x))]
    @test_int [cos(Derivative(-1 + m)(f)(x)*Derivative(-1 + n)(g)(x))*(Derivative(m)(f)(x)*Derivative(-1 + n)(g)(x) + Derivative(-1 + m)(f)(x)*Derivative(n)(g)(x)), x, 2, sin(Derivative(-1 + m)(f)(x)*Derivative(-1 + n)(g)(x))]

    @test_int [(g(x)*Derivative(1)(f)(x) + f(x)*g'(x))/(a + b*f(x)^2*g(x)^2), x, 2, atan((sqrt(b)*f(x)*g(x))/sqrt(a))/(sqrt(a)*sqrt(b))]
    @test_int [(g'(x)*Derivative(m)(f)(x) + g(x)*Derivative(1 + m)(f)(x))/(a + b*g(x)^2*Derivative(m)(f)(x)^2), x, 2, atan((sqrt(b)*g(x)*Derivative(m)(f)(x))/sqrt(a))/(sqrt(a)*sqrt(b))]
    @test_int [(Derivative(1 + m)(f)(x)*Derivative(n)(g)(x) + Derivative(m)(f)(x)*Derivative(1 + n)(g)(x))/(a + b*Derivative(m)(f)(x)^2*Derivative(n)(g)(x)^2), x, 2, atan((sqrt(b)*Derivative(m)(f)(x)*Derivative(n)(g)(x))/sqrt(a))/(sqrt(a)*sqrt(b))]

    @test_int [Derivative(1)(F)(f(x)*g(x))*(g(x)*Derivative(1)(f)(x) + f(x)*g'(x)), x, 2, F(f(x)*g(x))]
    @test_int [Derivative(1)(F)(g(x)*Derivative(m)(f)(x))*(g'(x)*Derivative(m)(f)(x) + g(x)*Derivative(1 + m)(f)(x)), x, 2, F(g(x)*Derivative(m)(f)(x))]
    @test_int [Derivative(1)(F)(Derivative(-1 + m)(f)(x)*Derivative(-1 + n)(g)(x))*(Derivative(m)(f)(x)*Derivative(-1 + n)(g)(x) + Derivative(-1 + m)(f)(x)*Derivative(n)(g)(x)), x, 2, F(Derivative(-1 + m)(f)(x)*Derivative(-1 + n)(g)(x))]

    @test_int [cos(f(x)^2*g(x))*f(x)*(2*g(x)*Derivative(1)(f)(x) + f(x)*g'(x)), x, 2, sin(f(x)^2*g(x))]
    @test_int [cos(g(x)^2*Derivative(m)(f)(x))*g(x)*(2*g'(x)*Derivative(m)(f)(x) + g(x)*Derivative(1 + m)(f)(x)), x, 2, sin(g(x)^2*Derivative(m)(f)(x))]
    @test_int [cos(g(x)*Derivative(m)(f)(x)^2)*Derivative(m)(f)(x)*(g'(x)*Derivative(m)(f)(x) + 2*g(x)*Derivative(1 + m)(f)(x)), x, 2, sin(g(x)*Derivative(m)(f)(x)^2)]
    @test_int [cos(Derivative(-1 + m)(f)(x)^2*Derivative(-1 + n)(g)(x))*Derivative(-1 + m)(f)(x)*(2*Derivative(m)(f)(x)*Derivative(-1 + n)(g)(x) + Derivative(-1 + m)(f)(x)*Derivative(n)(g)(x)), x, 2, sin(Derivative(-1 + m)(f)(x)^2*Derivative(-1 + n)(g)(x))]

    @test_int [(f(x)*(2*g(x)*Derivative(1)(f)(x) + f(x)*g'(x)))/(a + b*f(x)^4*g(x)^2), x, 2, atan((sqrt(b)*f(x)^2*g(x))/sqrt(a))/(sqrt(a)*sqrt(b))]
    @test_int [(g(x)*(2*g'(x)*Derivative(m)(f)(x) + g(x)*Derivative(1 + m)(f)(x)))/(a + b*g(x)^4*Derivative(m)(f)(x)^2), x, 2, atan((sqrt(b)*g(x)^2*Derivative(m)(f)(x))/sqrt(a))/(sqrt(a)*sqrt(b))]
    @test_int [(Derivative(m)(f)(x)*(g'(x)*Derivative(m)(f)(x) + 2*g(x)*Derivative(1 + m)(f)(x)))/(a + b*g(x)^2*Derivative(m)(f)(x)^4), x, 2, atan((sqrt(b)*g(x)*Derivative(m)(f)(x)^2)/sqrt(a))/(sqrt(a)*sqrt(b))]
    @test_int [(Derivative(m)(f)(x)*(2*Derivative(1 + m)(f)(x)*Derivative(n)(g)(x) + Derivative(m)(f)(x)*Derivative(1 + n)(g)(x)))/(a + b*Derivative(m)(f)(x)^4*Derivative(n)(g)(x)^2), x, 2, atan((sqrt(b)*Derivative(m)(f)(x)^2*Derivative(n)(g)(x))/sqrt(a))/(sqrt(a)*sqrt(b))]

    @test_int [f(x)*Derivative(1)(F)(f(x)^2*g(x))*(2*g(x)*Derivative(1)(f)(x) + f(x)*g'(x)), x, 2, F(f(x)^2*g(x))]
    @test_int [g(x)*Derivative(1)(F)(g(x)^2*Derivative(m)(f)(x))*(2*g'(x)*Derivative(m)(f)(x) + g(x)*Derivative(1 + m)(f)(x)), x, 2, F(g(x)^2*Derivative(m)(f)(x))]
    @test_int [Derivative(1)(F)(g(x)*Derivative(m)(f)(x)^2)*Derivative(m)(f)(x)*(g'(x)*Derivative(m)(f)(x) + 2*g(x)*Derivative(1 + m)(f)(x)), x, 2, F(g(x)*Derivative(m)(f)(x)^2)]
    @test_int [Derivative(1)(F)(Derivative(-1 + m)(f)(x)^2*Derivative(-1 + n)(g)(x))*Derivative(-1 + m)(f)(x)*(2*Derivative(m)(f)(x)*Derivative(-1 + n)(g)(x) + Derivative(-1 + m)(f)(x)*Derivative(n)(g)(x)), x, 2, F(Derivative(-1 + m)(f)(x)^2*Derivative(-1 + n)(g)(x))]

    @test_int [cos(f(x)^2*g(x)^3)*f(x)*g(x)^2*(2*g(x)*Derivative(1)(f)(x) + 3*f(x)*g'(x)), x, 2, sin(f(x)^2*g(x)^3)]
    @test_int [cos(g(x)^3*Derivative(m)(f)(x)^2)*g(x)^2*Derivative(m)(f)(x)*(3*g'(x)*Derivative(m)(f)(x) + 2*g(x)*Derivative(1 + m)(f)(x)), x, 2, sin(g(x)^3*Derivative(m)(f)(x)^2)]
    @test_int [cos(Derivative(m)(f)(x)^2*Derivative(n)(g)(x)^3)*Derivative(m)(f)(x)*Derivative(n)(g)(x)^2*(2*Derivative(1 + m)(f)(x)*Derivative(n)(g)(x) + 3*Derivative(m)(f)(x)*Derivative(1 + n)(g)(x)), x, 2, sin(Derivative(m)(f)(x)^2*Derivative(n)(g)(x)^3)]

    @test_int [(f(x)*g(x)^2*(2*g(x)*Derivative(1)(f)(x) + 3*f(x)*g'(x)))/(a + b*f(x)^4*g(x)^6), x, 2, atan((sqrt(b)*f(x)^2*g(x)^3)/sqrt(a))/(sqrt(a)*sqrt(b))]
    @test_int [(g(x)^2*Derivative(m)(f)(x)*(3*g'(x)*Derivative(m)(f)(x) + 2*g(x)*Derivative(1 + m)(f)(x)))/(a + b*g(x)^6*Derivative(m)(f)(x)^4), x, 2, atan((sqrt(b)*g(x)^3*Derivative(m)(f)(x)^2)/sqrt(a))/(sqrt(a)*sqrt(b))]
    @test_int [Derivative(m)(f)(x)*Derivative(n)(g)(x)^2*(2*Derivative(1 + m)(f)(x)*Derivative(n)(g)(x) + 3*Derivative(m)(f)(x)*Derivative(1 + n)(g)(x))/(a + b*Derivative(m)(f)(x)^4*Derivative(n)(g)(x)^6), x, 2, atan((sqrt(b)*Derivative(m)(f)(x)^2*Derivative(n)(g)(x)^3)/sqrt(a))/(sqrt(a)*sqrt(b))]

    @test_int [f(x)*g(x)^2*Derivative(1)(F)(f(x)^2*g(x)^3)*(2*g(x)*Derivative(1)(f)(x) + 3*f(x)*g'(x)), x, 2, F(f(x)^2*g(x)^3)]
    @test_int [g(x)^2*Derivative(1)(F)(g(x)^3*Derivative(m)(f)(x)^2)*Derivative(m)(f)(x)*(3*g'(x)*Derivative(m)(f)(x) + 2*g(x)*Derivative(1 + m)(f)(x)), x, 2, F(g(x)^3*Derivative(m)(f)(x)^2)]
    @test_int [Derivative(1)(F)(Derivative(m)(f)(x)^2*Derivative(n)(g)(x)^3)*Derivative(m)(f)(x)*Derivative(n)(g)(x)^2*(2*Derivative(1 + m)(f)(x)*Derivative(n)(g)(x) + 3*Derivative(m)(f)(x)*Derivative(1 + n)(g)(x)), x, 2, F(Derivative(m)(f)(x)^2*Derivative(n)(g)(x)^3)]
    end

max_stars_repo_path: test/8 Special functions/8.10 Formal derivatives.jl


text:
"""
Penalties that are applied to the entire parameter array only
"""
abstract type ArrayPenalty <: Penalty end
name(p::ArrayPenalty) = replace(string(typeof(p)), "PenaltyFunctions." => "")

#------------------------------------------------------------------# abstract methods
value(p::ArrayPenalty, A::AA{<:Number}, λ::Number) = λ * value(p, A)

#----------------------------------------------------------------# NuclearNormPenalty
struct NuclearNormPenalty <: ArrayPenalty end
function value(p::NuclearNormPenalty, A::AbstractMatrix{<:Number})
    >(size(A)...) ? tr(sqrt(A'A)) : tr(sqrt(A * A'))
end
function prox!(p::NuclearNormPenalty, A::AbstractMatrix{<:Number}, λ::Number)
    svdecomp = svd!(A)
    soft_thresh!(svdecomp.S, λ)
    copyto!(A, Matrix(svdecomp))
end

#-----------------------------------------------------------------# GroupLassoPenalty
"Group Lasso Penalty.  Able to set the entire vector (group) to 0."
struct GroupLassoPenalty <: ArrayPenalty end
value(p::GroupLassoPenalty, A::AA{<:Number}) = norm(A)
function prox!(p::GroupLassoPenalty, A::AA{T}, λ::Number) where {T <: Number}
    denom = norm(A)
    if denom <= λ
        fill!(A, zero(T))
    else
        scaling = λ / denom
        for i in eachindex(A)
            @inbounds A[i] = (1 - scaling) * A[i]
        end
    end
    A
end

#-----------------------------------------------------------------# MahalanobisPenalty
"""
    MahalanobisPenalty(C)

Supports a Mahalanobis distance penalty (`xᵀCᵀCx` for a vector `x`).
"""
mutable struct MahalanobisPenalty{T <: Number, S <: AA{T,2}} <: ArrayPenalty
    C::S
    CtC::S
    CtC_Iλ::LU{T, Matrix{T}} # LU factorization of C'C + I/λ
    λ::T
end
function MahalanobisPenalty(C::AbstractMatrix{T}, λ::T = one(T)) where {T<:Number}
    MahalanobisPenalty(C, C'C, lu(C'C + I), λ)
end
function value(p::MahalanobisPenalty{T}, x::AbstractVector{T}) where {T <: Number}
    inv(T(2)) * T(sum(abs2, p.C * x))
end
function prox!(p::MahalanobisPenalty{T}, A::AA{T, 1}, λ::Number) where {T <: Number}
    if λ != p.λ
        p.λ = λ
        p.CtC_Iλ = lu(p.CtC + I / λ)
    end
    rmul!(A, one(T) / λ)
    ldiv!(p.CtC_Iλ, A) # overwrites result in A
end

#--------------------------------------------------------------------------------# scaled
struct ScaledArrayPenalty{T, P <: ArrayPenalty} <: ArrayPenalty
    penalty::P
    λ::T
end
scaled(p::ArrayPenalty, λ::Number) = (_scale_check(λ); ScaledArrayPenalty(p, λ))

Base.show(io::IO, sp::ScaledArrayPenalty) = print(io, "$(sp.λ) * ($(sp.penalty))")

value(p::ScaledArrayPenalty, θ::AA{<:Number}) = p.λ * value(p.penalty, θ)
prox!(p::ScaledArrayPenalty, θ::AA{<:Number}) = prox!(p.penalty, θ, p.λ)


text:
"""
    QuadratureRule

Integration rule for approximating the continuous integrals for the equations of motion.
Currently divided into two classes of integrators:
* [`Explicit`](@ref)
* [`Implicit`](@ref)

## Interface
All integrators need to have a constructor of the form:

    Q(n, m)

where `n` and `m` are the state and control dimensions of the dynamical system.
"""
abstract type QuadratureRule end

# General member-wise copy for all integrators
@generated function Base.copy(int::T) where T <: QuadratureRule
    copyexprs = [:(copy(int.$field)) for field in fieldnames(T)]
    return :(T($(copyexprs...)))
end

"""
    Explicit <: QuadratureRule

Integration rules of the form
```math
x_{k+1} = f(x_k,u_k)
```.

## Interface
All explicit integrators need to define the following methods:

    integrate(::Q, model, x, u, t, h)
    integrate!(::Q, model, xn, x, u, t, h)
    jacobian!(::Q, sig, model, J, xn, x, u, t, h)

where `h` is the time step (i.e. `dt`). The Jacobian `J` should have size `(n, n + m)`.
The `sig` is a [`FunctionSignature`](@ref) that can be used to dispatch the Jacobian on
the desired function signature. For example, when using `StaticReturn` intermediate
cache variables aren't usually needed, but are when using `InPlace` methods. These cache
variables can be stored in the integrator itself (see the [`ADVector`](@ref) type provided
by this package).

!!! note
    The `Q` in these methods should be replaced by your explicit integrator,
    e.g. `MyRungeKutta`.

"""
abstract type Explicit <: QuadratureRule end

"""
    Implicit <: QuadratureRule

Integration rules of the form
```math
f(x_{k+1},u_{k+1},x_k,u_k) = 0
```.

## Interface
All implicit integrators need to define the following methods:

    dynamics_error(::Q, model, z2, z1)
    dynamics_error!(::Q, model, y2, y1, z2, z1)
    dynamics_error_jacobian!(::Q, sig, model, J2, J1, y2, y1, z2, z1)

where `model` is a [`ContinuousDynamics`](@ref) model. For the in place method,
the output should be stored in the `y2` vector. For the Jacobian method, `J2`
holds the Jacobian with respect to the state and control of `z2` (the knot point
at the next time step) and `J1` holds the Jacobian with respect to the state
and control of `z1` (the knot point at the current time step).

## Treating an Implicit method as an Explicit method
Implicit methods can be treated as an explicit integrator. In addition to the the
methods above that evaluate the dynamics error residual for the integrator, they
overload `discrete_dynamics`, `discrete_dynamics!`, and `jacobian`. The forward simulation
functions `discrete_dynamics` and `discrete_dynamics!` use Newton's method. The dynamics
Jacobians are calculated using the Implicit Function Theorem.

If the state and control vector passed to `jacobian!` are the same as the last call
to `jacobian!`, `discrete_dynamics` or `discrete_dynamics!` the Jacobians and matrix
factorization from the Newton solve are used to speed up the computation. Note that
the factorization is only cached for `InPlace` methods, not `StaticReturn`.

To call `jacobian!` on `DiscretizedDynamics` with an implicit dynamics model, it must be
called using the `ImplicitFunctionTheorem` method, which also specified the
[`DiffMethod`](@ref) that should be used to evaluate `dynamics_error_jacobian!`, e.g.

    jacobian!(InPlace(), ImplicitFunctionTheorem(ForwardDiff), implicitmodel, J, y, z)

## Defining an Implicit Integrator
An implicit integrator should define the signatures above for `dynamics_error` and
`dynamics_error_jacobian!`. If the method wants to provide the functionality of an
explicit integrator (as described in the previous section), it should store internally
a [`ImplicitNewtonCache`](@ref), which must be returned using the getter function
`getnewtoncache`.
"""
abstract type Implicit <: QuadratureRule end
setnewtontolerance(integrator::Implicit, tol) = getnewtoncache(integrator).newton_tol = tol
setnewtoniters(integrator::Implicit, iter) = getnewtoncache(integrator).newton_iters = iter

# Create an integrator from a dynamics model
(::Type{Q})(model::AbstractModel) where {Q<:QuadratureRule} = Q(state_dim(model), control_dim(model))

"""
    DiscretizedDynamics

Represents a [`DiscreteDynamics`](@ref) model formed by integrating a continuous
dynamics model. It is essentially a [`ContinuousDynamics`](@ref) paired with a
[`QuadratureRule`](@ref) that defines how to use the [`dynamics!`](@ref) function
to get a [`discrete_dynamics!`](@ref) function.

## Constructor
A `DiscretedDynamics` type can be created using either of the following signatures:

    DiscretizedDynamics(dynamics::ContinuousDynamics, Q::QuadratureRule)
    DiscretizedDynamics{Q}(dynamics::ContinuousDynamics) where Q <: QuadratureRule

In the second case, the integrator is constructed by calling
`Q(state_dim(dynamics), control_dim(dynamics))`.

## Usage
A `DiscretizedDynamics` model is used just like any other [`DiscreteDynamics`](@ref) model.
The state, control, and error state dimensions are all taken from the
continuous time dynamics model, and the [`StateVectorType`](@ref) and
[`FunctionInputs`](@ref) traits are inherited from the continuous time dynamics.

The `default_diffmethod`, however, is set to `ForwardAD`, since our benchmarks show that
it is usually faster to make a single call to `ForwardDiff.jacobian!` than using the chain
rule and multiple calls to query the Jacobian of the continuous dynamics at multiple points.
If the continuous model has a `UserDefined` Jacobian method, calling `jacobian!` on
a `DiscretizedDynamics` model will use the chain rule with analytical Jacobians of
the integrator. If the user wants to use a combination of ForwardDiff (or any other
differentiation method, for that matter) with the integrator, they are free to define their
own [`DiffMethod`](@ref) to dispatch on.

## Implicit Dynamics
If an [`Implicit`](@ref) integrator is used, the interface changes slightly. The dynamics
residual is calculated and it's Jacobian are calculated using

    dynamics_error(model, z2, z1)
    dynamics_error!(model, y2, y1, z2, z1)
    dynamics_error_jacobian!(sig, diff, model, J2, J1, y2, y1, z2, z1)

Where the dynamics residual is stored in `y2` for the in-place version, and the Jacobians
with respect to the state and control at the next and current knot points are stored in `J2`
and `J1`, respectively.

Implicit dynamics integrators can also support the normal API
(see documentation for [`Implicit`](@ref)). Note that calls to `discrete_dynamics` will
use a Newton's solve to calculate the next state and calls to `jacobian!` will use the
Implicit Function Theorem. In order to use the implicit function theorem to get the
Jacobians, `jacobian!` must be called using the `ImplicitFunctionTheorem`
[`DiffMethod`](@ref), which also specifies the `DiffMethod` to use to on
`dynamics_error_jacobian!`.

When using Newton's method during `discrete_dynamics`, the `DiffMethod` returned by
`default_diffmethod` for the continuous dynamics model is used to evaluate the
`dynamics_error_jacobian!`. There is no option to change this at run time, since
`discrete_dynamics` and `discrete_dynamics!` do not take a `DiffMethod` as an argument.

The convergence tolerance and max number of iterations of the Newton's method can be set
using

    setnewtontolerance(model::ImplicitDynamicsModel, tol)
    setnewtoniters(model::ImplicitDynamicsModel, iters)
"""
@autodiff struct DiscretizedDynamics{L,Q} <: DiscreteDynamics
    continuous_dynamics::L
    integrator::Q
    function DiscretizedDynamics(
        continuous_dynamics::L,
        integrator::Q,
    ) where {L<:ContinuousDynamics,Q<:QuadratureRule}
        new{L,Q}(continuous_dynamics, integrator)
    end
end
function DiscretizedDynamics{Q}(
    continuous_dynamics::L,
) where {L<:ContinuousDynamics} where {Q}
    n, m = dims(continuous_dynamics)
    DiscretizedDynamics(continuous_dynamics, Q(n, m))
end

for method in (:state_dim, :control_dim, :output_dim, :errstate_dim, :statevectortype,
               :functioninputs)
    @eval $method(model::DiscretizedDynamics) = $method(model.continuous_dynamics)
end
default_diffmethod(::DiscretizedDynamics) = ForwardAD()
default_signature(model::DiscretizedDynamics) = default_signature(model.continuous_dynamics)

# Shallow copy
function Base.copy(model::DiscretizedDynamics)
    DiscretizedDynamics(model.continuous_dynamics, model.integrator)
end

function Base.deepcopy(model::DiscreteDynamics)
    DiscreteDynamics(copy(model.continuous_dynamics), copy(model.integrator))
end

@inline integration(model::DiscretizedDynamics) = model.integrator
discrete_dynamics(model::DiscretizedDynamics, x, u, t, dt) =
    integrate(integration(model), model.continuous_dynamics, x, u, t, dt)
discrete_dynamics!(model::DiscretizedDynamics, xn, x, u, t, dt) =
    integrate!(integration(model), model.continuous_dynamics, xn, x, u, t, dt)

jacobian!(sig::FunctionSignature, ::UserDefined, model::DiscretizedDynamics, J, xn, z) =
    jacobian!(
        integration(model),
        sig,
        model.continuous_dynamics,
        J,
        xn,
        state(z),
        control(z),
        time(z),
        timestep(z),
    )

########################################
# Implicit Dynamics
########################################
const ImplicitDynamicsModel{L,Q} = DiscretizedDynamics{L,Q} where {L,Q<:Implicit}

@inline setnewtontolerance(model::ImplicitDynamicsModel, tol) =
    setnewtontolerance(integration(model), tol)
@inline setnewtoniters(model::ImplicitDynamicsModel, iter) =
    setnewtoniters(integration(model), iter)

discrete_dynamics(model::ImplicitDynamicsModel, z::AbstractKnotPoint) =
    integrate(integration(model), model, z)
# discrete_dynamics!(model::ImplicitDynamicsModel, xn, z::AbstractKnotPoint) =
#     integrate!(integration(model), model, xn, z)

discrete_dynamics!(model::ImplicitDynamicsModel, xn, x, u, t, dt) =
    integrate!(integration(model), model, xn, x, u, t, dt)

function jacobian!(sig::FunctionSignature, ::ImplicitFunctionTheorem{D}, model::ImplicitDynamicsModel, J, xn, z) where D
    jacobian!(integration(model), sig, D(), model, J, xn, z)
end

for sig in (:InPlace, :StaticReturn), diff in (:ForwardAD, :FiniteDifference, :UserDefined)
    @eval begin
        function jacobian!(::$sig, ::$diff, ::ImplicitDynamicsModel, J, xn, z)
            throw(ArgumentError("Must call jacobian! on an ImplicitDynamicsModel with the ImplicitFunctionTheorem diff method. Got $($diff)."))
        end
    end
end

default_diffmethod(model::ImplicitDynamicsModel) =
    ImplicitFunctionTheorem(default_diffmethod(model.continuous_dynamics))

# TODO: [#19] overwrite `evaluate` to solve for the next state using Newton
# TODO: [#19] overwrite `jacobian` to provide A,B using implicit function theorem

dynamics_error(
    model::ImplicitDynamicsModel,
    z2::AbstractKnotPoint,
    z1::AbstractKnotPoint,
) = dynamics_error(integration(model), model.continuous_dynamics, z2, z1)

dynamics_error!(
    model::ImplicitDynamicsModel,
    y2,
    y1,
    z2::AbstractKnotPoint,
    z1::AbstractKnotPoint,
) = dynamics_error!(integration(model), model.continuous_dynamics, y2, y1, z2, z1)

dynamics_error_jacobian!(
    sig::FunctionSignature,
    ::UserDefined,
    model::ImplicitDynamicsModel,
    J2,
    J1,
    y2,
    y1,
    z2::AbstractKnotPoint,
    z1::AbstractKnotPoint,
) = dynamics_error_jacobian!(
    integration(model),
    sig,
    model.continuous_dynamics,
    J2,
    J1,
    y2,
    y1,
    z2,
    z1,
)


text:
using Mocking

create_ipset_1_patch = @patch function create_ipset(a...; b...)
    return Dict{String, Any}("ChangeToken" => "abcd12f2-46da-4fdb-b8d5-fbd4c466928f", "IPSet" => Dict{String, Any}("IPSetId" => "example1ds3t-46da-4fdb-b8d5-abc321j569j5", "Name" => "MyIPSetFriendlyName", "IPSetDescriptors" => Any[Dict{String, Any}("Value" => "192.0.2.44/32", "Type" => "IPV4")]))
end

create_rule_1_patch = @patch function create_rule(a...; b...)
    return Dict{String, Any}("ChangeToken" => "abcd12f2-46da-4fdb-b8d5-fbd4c466928f", "Rule" => Dict{String, Any}("MetricName" => "WAFByteHeaderRule", "Name" => "WAFByteHeaderRule", "RuleId" => "WAFRule-1-Example", "Predicates" => Any[Dict{String, Any}("Negated" => false, "DataId" => "MyByteMatchSetID", "Type" => "ByteMatch")]))
end

create_size_constraint_set_1_patch = @patch function create_size_constraint_set(a...; b...)
    return Dict{String, Any}("ChangeToken" => "abcd12f2-46da-4fdb-b8d5-fbd4c466928f", "SizeConstraintSet" => Dict{String, Any}("SizeConstraintSetId" => "example1ds3t-46da-4fdb-b8d5-abc321j569j5", "SizeConstraints" => Any[Dict{String, Any}("FieldToMatch" => Dict{String, Any}("Type" => "QUERY_STRING"), "TextTransformation" => "NONE", "Size" => 0, "ComparisonOperator" => "GT")], "Name" => "MySampleSizeConstraintSet"))
end

create_sql_injection_match_set_1_patch = @patch function create_sql_injection_match_set(a...; b...)
    return Dict{String, Any}("ChangeToken" => "abcd12f2-46da-4fdb-b8d5-fbd4c466928f", "SqlInjectionMatchSet" => Dict{String, Any}("SqlInjectionMatchSetId" => "example1ds3t-46da-4fdb-b8d5-abc321j569j5", "Name" => "MySQLInjectionMatchSet", "SqlInjectionMatchTuples" => Any[Dict{String, Any}("FieldToMatch" => Dict{String, Any}("Type" => "QUERY_STRING"), "TextTransformation" => "URL_DECODE")]))
end

create_web_acl_1_patch = @patch function create_web_acl(a...; b...)
    return Dict{String, Any}("WebACL" => Dict{String, Any}("DefaultAction" => Dict{String, Any}("Type" => "ALLOW"), "Rules" => Any[Dict{String, Any}("Priority" => 1, "Action" => Dict{String, Any}("Type" => "ALLOW"), "RuleId" => "WAFRule-1-Example")], "MetricName" => "CreateExample", "WebACLId" => "example-46da-4444-5555-example", "Name" => "CreateExample"), "ChangeToken" => "abcd12f2-46da-4fdb-b8d5-fbd4c466928f")
end

create_xss_match_set_1_patch = @patch function create_xss_match_set(a...; b...)
    return Dict{String, Any}("XssMatchSet" => Dict{String, Any}("XssMatchSetId" => "example1ds3t-46da-4fdb-b8d5-abc321j569j5", "Name" => "MySampleXssMatchSet", "XssMatchTuples" => Any[Dict{String, Any}("FieldToMatch" => Dict{String, Any}("Type" => "QUERY_STRING"), "TextTransformation" => "URL_DECODE")]), "ChangeToken" => "abcd12f2-46da-4fdb-b8d5-fbd4c466928f")
end

delete_byte_match_set_1_patch = @patch function delete_byte_match_set(a...; b...)
    return Dict{String, Any}("ChangeToken" => "abcd12f2-46da-4fdb-b8d5-fbd4c466928f")
end

delete_ipset_1_patch = @patch function delete_ipset(a...; b...)
    return Dict{String, Any}("ChangeToken" => "abcd12f2-46da-4fdb-b8d5-fbd4c466928f")
end

delete_rule_1_patch = @patch function delete_rule(a...; b...)
    return Dict{String, Any}("ChangeToken" => "abcd12f2-46da-4fdb-b8d5-fbd4c466928f")
end

delete_size_constraint_set_1_patch = @patch function delete_size_constraint_set(a...; b...)
    return Dict{String, Any}("ChangeToken" => "abcd12f2-46da-4fdb-b8d5-fbd4c466928f")
end

delete_sql_injection_match_set_1_patch = @patch function delete_sql_injection_match_set(a...; b...)
    return Dict{String, Any}("ChangeToken" => "abcd12f2-46da-4fdb-b8d5-fbd4c466928f")
end

delete_web_acl_1_patch = @patch function delete_web_acl(a...; b...)
    return Dict{String, Any}("ChangeToken" => "abcd12f2-46da-4fdb-b8d5-fbd4c466928f")
end

delete_xss_match_set_1_patch = @patch function delete_xss_match_set(a...; b...)
    return Dict{String, Any}("ChangeToken" => "abcd12f2-46da-4fdb-b8d5-fbd4c466928f")
end

get_byte_match_set_1_patch = @patch function get_byte_match_set(a...; b...)
    return Dict{String, Any}("ByteMatchSet" => Dict{String, Any}("ByteMatchSetId" => "exampleIDs3t-46da-4fdb-b8d5-abc321j569j5", "ByteMatchTuples" => Any[Dict{String, Any}("FieldToMatch" => Dict{String, Any}("Data" => "referer", "Type" => "HEADER"), "TextTransformation" => "NONE", "TargetString" => "badrefer1", "PositionalConstraint" => "CONTAINS")], "Name" => "ByteMatchNameExample"))
end

get_change_token_1_patch = @patch function get_change_token(a...; b...)
    return Dict{String, Any}("ChangeToken" => "abcd12f2-46da-4fdb-b8d5-fbd4c466928f")
end

get_change_token_status_1_patch = @patch function get_change_token_status(a...; b...)
    return Dict{String, Any}("ChangeTokenStatus" => "PENDING")
end

get_ipset_1_patch = @patch function get_ipset(a...; b...)
    return Dict{String, Any}("IPSet" => Dict{String, Any}("IPSetId" => "example1ds3t-46da-4fdb-b8d5-abc321j569j5", "Name" => "MyIPSetFriendlyName", "IPSetDescriptors" => Any[Dict{String, Any}("Value" => "192.0.2.44/32", "Type" => "IPV4")]))
end

get_rule_1_patch = @patch function get_rule(a...; b...)
    return Dict{String, Any}("Rule" => Dict{String, Any}("MetricName" => "WAFByteHeaderRule", "Name" => "WAFByteHeaderRule", "RuleId" => "example1ds3t-46da-4fdb-b8d5-abc321j569j5", "Predicates" => Any[Dict{String, Any}("Negated" => false, "DataId" => "MyByteMatchSetID", "Type" => "ByteMatch")]))
end

get_sampled_requests_1_patch = @patch function get_sampled_requests(a...; b...)
    return Dict{String, Any}("PopulationSize" => 50, "SampledRequests" => Any[Dict{String, Any}("Request" => Dict{String, Any}("Headers" => Any[Dict{String, Any}("Value" => "BadBot ", "Name" => "User-Agent")], "Method" => "HEAD", "HTTPVersion" => "HTTP/1.1", "Country" => "US", "ClientIP" => "192.0.2.44"), "Weight" => 1, "Action" => "BLOCK", "Timestamp" => "2016-09-27T14:55Z")], "TimeWindow" => Dict{String, Any}("StartTime" => "2016-09-27T14:50Z", "EndTime" => "2016-09-27T15:50Z"))
end

get_size_constraint_set_1_patch = @patch function get_size_constraint_set(a...; b...)
    return Dict{String, Any}("SizeConstraintSet" => Dict{String, Any}("SizeConstraintSetId" => "example1ds3t-46da-4fdb-b8d5-abc321j569j5", "SizeConstraints" => Any[Dict{String, Any}("FieldToMatch" => Dict{String, Any}("Type" => "QUERY_STRING"), "TextTransformation" => "NONE", "Size" => 0, "ComparisonOperator" => "GT")], "Name" => "MySampleSizeConstraintSet"))
end

get_sql_injection_match_set_1_patch = @patch function get_sql_injection_match_set(a...; b...)
    return Dict{String, Any}("SqlInjectionMatchSet" => Dict{String, Any}("SqlInjectionMatchSetId" => "example1ds3t-46da-4fdb-b8d5-abc321j569j5", "Name" => "MySQLInjectionMatchSet", "SqlInjectionMatchTuples" => Any[Dict{String, Any}("FieldToMatch" => Dict{String, Any}("Type" => "QUERY_STRING"), "TextTransformation" => "URL_DECODE")]))
end

get_web_acl_1_patch = @patch function get_web_acl(a...; b...)
    return Dict{String, Any}("WebACL" => Dict{String, Any}("DefaultAction" => Dict{String, Any}("Type" => "ALLOW"), "Rules" => Any[Dict{String, Any}("Priority" => 1, "Action" => Dict{String, Any}("Type" => "ALLOW"), "RuleId" => "WAFRule-1-Example")], "MetricName" => "CreateExample", "WebACLId" => "createwebacl-1472061481310", "Name" => "CreateExample"))
end

get_xss_match_set_1_patch = @patch function get_xss_match_set(a...; b...)
    return Dict{String, Any}("XssMatchSet" => Dict{String, Any}("XssMatchSetId" => "example1ds3t-46da-4fdb-b8d5-abc321j569j5", "Name" => "MySampleXssMatchSet", "XssMatchTuples" => Any[Dict{String, Any}("FieldToMatch" => Dict{String, Any}("Type" => "QUERY_STRING"), "TextTransformation" => "URL_DECODE")]))
end

list_ipsets_1_patch = @patch function list_ipsets(a...; b...)
    return Dict{String, Any}("IPSets" => Any[Dict{String, Any}("IPSetId" => "abcd12f2-46da-4fdb-b8d5-fbd4c466928f", "Name" => "MyIPSetFriendlyName")])
end

list_rules_1_patch = @patch function list_rules(a...; b...)
    return Dict{String, Any}("Rules" => Any[Dict{String, Any}("Name" => "WAFByteHeaderRule", "RuleId" => "WAFRule-1-Example")])
end

list_size_constraint_sets_1_patch = @patch function list_size_constraint_sets(a...; b...)
    return Dict{String, Any}("SizeConstraintSets" => Any[Dict{String, Any}("SizeConstraintSetId" => "example1ds3t-46da-4fdb-b8d5-abc321j569j5", "Name" => "MySampleSizeConstraintSet")])
end

list_sql_injection_match_sets_1_patch = @patch function list_sql_injection_match_sets(a...; b...)
    return Dict{String, Any}("SqlInjectionMatchSets" => Any[Dict{String, Any}("SqlInjectionMatchSetId" => "example1ds3t-46da-4fdb-b8d5-abc321j569j5", "Name" => "MySQLInjectionMatchSet")])
end

list_web_acls_1_patch = @patch function list_web_acls(a...; b...)
    return Dict{String, Any}("WebACLs" => Any[Dict{String, Any}("WebACLId" => "webacl-1472061481310", "Name" => "WebACLexample")])
end

list_xss_match_sets_1_patch = @patch function list_xss_match_sets(a...; b...)
    return Dict{String, Any}("XssMatchSets" => Any[Dict{String, Any}("XssMatchSetId" => "example1ds3t-46da-4fdb-b8d5-abc321j569j5", "Name" => "MySampleXssMatchSet")])
end

update_byte_match_set_1_patch = @patch function update_byte_match_set(a...; b...)
    return Dict{String, Any}("ChangeToken" => "abcd12f2-46da-4fdb-b8d5-fbd4c466928f")
end

update_ipset_1_patch = @patch function update_ipset(a...; b...)
    return Dict{String, Any}("ChangeToken" => "abcd12f2-46da-4fdb-b8d5-fbd4c466928f")
end

update_rule_1_patch = @patch function update_rule(a...; b...)
    return Dict{String, Any}("ChangeToken" => "abcd12f2-46da-4fdb-b8d5-fbd4c466928f")
end

update_size_constraint_set_1_patch = @patch function update_size_constraint_set(a...; b...)
    return Dict{String, Any}("ChangeToken" => "abcd12f2-46da-4fdb-b8d5-fbd4c466928f")
end

update_sql_injection_match_set_1_patch = @patch function update_sql_injection_match_set(a...; b...)
    return Dict{String, Any}("ChangeToken" => "abcd12f2-46da-4fdb-b8d5-fbd4c466928f")
end

update_web_acl_1_patch = @patch function update_web_acl(a...; b...)
    return Dict{String, Any}("ChangeToken" => "abcd12f2-46da-4fdb-b8d5-fbd4c466928f")
end

update_xss_match_set_1_patch = @patch function update_xss_match_set(a...; b...)
    return Dict{String, Any}("ChangeToken" => "abcd12f2-46da-4fdb-b8d5-fbd4c466928f")
end


text:
using HDF5
const test_path = splitdir(@__FILE__)[1]

# Create a new file
fn = joinpath(tempdir(),"test.h5")
f = h5open(fn, "w")
# Write scalars
f["Float64"] = 3.2
f["Int16"] = int16(4)
# Create arrays of different types
A = randn(3,5)
write(f, "Afloat64", float64(A))
write(f, "Afloat32", float32(A))
Ai = rand(1:20, 2, 4)
write(f, "Aint8", int8(Ai))
f["Aint16"] = int16(Ai)
write(f, "Aint32", int32(Ai))
write(f, "Aint64", int64(Ai))
write(f, "Auint8", uint8(Ai))
write(f, "Auint16", uint16(Ai))
write(f, "Auint32", uint32(Ai))
write(f, "Auint64", uint64(Ai))
# Test strings
salut = "Hi there"
ucode = "uniçº∂e"
write(f, "salut", salut)
write(f, "ucode", ucode)
# Arrays of strings
salut_split = ["Hi", "there"]
write(f, "salut_split", salut_split)
# Empty arrays
empty = Array(Uint32, 0)
write(f, "empty", empty)
# Empty strings
empty_string = ""
write(f, "empty_string", empty_string)
# Empty array of strings
empty_string_array = ASCIIString[]
write(f, "empty_string_array", empty_string_array)
# Array of empty string
empty_array_of_strings = ASCIIString[""]
write(f, "empty_array_of_strings", empty_array_of_strings)
# Attributes
species = [["N", "C"], ["A", "B"]]
attrs(f)["species"] = species
C∞ = 42
attrs(f)["C∞"] = C∞
dset = f["salut"]
label = "This is a string"
attrs(dset)["typeinfo"] = label
close(dset)
# Group
g = g_create(f, "mygroup")
# Test dataset with compression
R = rand(1:20, 20, 40);
g["CompressedA", "chunk", (5,6), "compress", 9] = R
close(g)
# Copy group containing dataset
o_copy(f, "mygroup", f, "mygroup2")
# Copy dataset
g = g_create(f, "mygroup3")
o_copy(f["mygroup/CompressedA"], g, "CompressedA")
close(g)
# Writing hyperslabs
dset = d_create(f,"slab",datatype(Float64),dataspace(20,20,5),"chunk",(5,5,1))
Xslab = randn(20,20,5)
for i = 1:5
    dset[:,:,i] = Xslab[:,:,i]
end
# More complex hyperslab and assignment with "incorrect" types (issue #34)
d = d_create(f, "slab2", datatype(Float64), ((10,20),(100,200)), "chunk", (1,1))
d[:,:] = 5
d[1,1] = 4
# 1d indexing
d = d_create(f, "slab3", datatype(Int), ((10,),(-1,)), "chunk", (5,))
@assert d[:] == zeros(Int, 10)
d[3:5] = 3:5
# Create a dataset designed to be deleted
f["deleteme"] = 17.2
close(f)
# Test the h5read/write interface
W = reshape(1:120, 15, 8)
h5write(fn, "newgroup/W", W)

# Read the file back in
fr = h5open(fn)
x = read(fr, "Float64")
@assert x == 3.2 && isa(x, Float64)
y = read(fr, "Int16")
@assert y == 4 && isa(y, Int16)
Af32 = read(fr, "Afloat32")
@assert float32(A) == Af32
@assert eltype(Af32) == Float32
Af64 = read(fr, "Afloat64")
@assert float64(A) == Af64
@assert eltype(Af64) == Float64
@assert eltype(fr["Afloat64"]) == Float64  # issue 167
Ai8 = read(fr, "Aint8")
@assert Ai == Ai8
@assert eltype(Ai8) == Int8
Ai16 = read(fr, "Aint16")
@assert Ai == Ai16
@assert eltype(Ai16) == Int16
Ai32 = read(fr, "Aint32")
@assert Ai == Ai32
@assert eltype(Ai32) == Int32
Ai64 = read(fr, "Aint64")
@assert Ai == Ai64
@assert eltype(Ai64) == Int64
Ai8 = read(fr, "Auint8")
@assert Ai == Ai8
@assert eltype(Ai8) == Uint8
Ai16 = read(fr, "Auint16")
@assert Ai == Ai16
@assert eltype(Ai16) == Uint16
Ai32 = read(fr, "Auint32")
@assert Ai == Ai32
@assert eltype(Ai32) == Uint32
Ai64 = read(fr, "Auint64")
@assert Ai == Ai64
@assert eltype(Ai64) == Uint64
salutr = read(fr, "salut")
@assert salut == salutr
ucoder = read(fr, "ucode")
@assert ucode == ucoder
salut_splitr = read(fr, "salut_split")
@assert salut_splitr == salut_split
Rr = read(fr, "mygroup/CompressedA")
@assert Rr == R
Rr2 = read(fr, "mygroup2/CompressedA")
@assert Rr2 == R
Rr3 = read(fr, "mygroup3/CompressedA")
@assert Rr3 == R
dset = fr["mygroup/CompressedA"]
@assert get_chunk(dset) == (5,6)
@assert name(dset) == "/mygroup/CompressedA"
Xslabr = read(fr, "slab")
@assert Xslabr == Xslab
Xslabr = h5read(fn, "slab", (:, :, :))  # issue #87
@assert Xslabr == Xslab
Xslab2r = read(fr, "slab2")
target = fill(5, 10, 20)
target[1] = 4
@assert Xslab2r == target
dset = fr["slab3"]
@assert dset[3:5] == [3:5]
emptyr = read(fr, "empty")
@assert isempty(emptyr)
empty_stringr = read(fr, "empty_string")
@assert empty_stringr == empty_string
empty_string_arrayr = read(fr, "empty_string_array")
@assert empty_string_arrayr == empty_string_array
empty_array_of_stringsr = read(fr, "empty_array_of_strings")
@assert empty_array_of_stringsr == empty_array_of_strings
@assert a_read(fr, "species") == species
@assert a_read(fr, "C∞") == C∞
dset = fr["salut"]
@assert a_read(dset, "typeinfo") == label
close(dset)
# Test ref-based reading
Aref = fr["Afloat64"]
sel = (2:3, 1:2:5)
Asub = Aref[sel...]
@assert Asub == A[sel...]
close(Aref)
# Test iteration, name, and parent
for obj in fr
    @assert filename(obj) == fn
    n = name(obj)
    p = parent(obj)
end
# Test reading multiple vars at once
z = read(fr, "Float64", "Int16")
@assert z == (3.2, 4)
@assert typeof(z) == (Float64, Int16)
# Test function syntax
read(fr, "Float64") do x
	@assert x == 3.2
end
read(fr, "Float64", "Int16") do x, y
	@assert x == 3.2
	@assert y == 4
end
# Test reading entire file at once
z = read(fr)
@assert z["Float64"] == 3.2
close(fr)

# Test object deletion
fr = h5open(fn, "r+")
@assert exists(fr, "deleteme")
o_delete(fr, "deleteme")
@assert !exists(fr, "deleteme")
close(fr)

# Test the h5read interface
Wr = h5read(fn, "newgroup/W")
@assert Wr == W
rng = (2:3:15, 3:5)
Wr = h5read(fn, "newgroup/W", rng)
@assert Wr == W[rng...]

# more do syntax
h5open(fn, "w") do fid
    g_create(fid, "mygroup") do g
        write(g, "x", 3.2)
    end
end
fid = h5open(fn, "r")
@assert names(fid) == ASCIIString["mygroup"]
g = fid["mygroup"]
@assert names(g) == ASCIIString["x"]
close(g)
close(fid)

d = h5read(joinpath(test_path, "compound.h5"), "/data")
@assert typeof(d) == HDF5.HDF5Compound
@assert typeof(d.data) == Array{Uint8,1}
@assert length(d.data) == 128
@assert d.membertype == Type[Float64, HDF5.FixedArray{Float64,(HDF5.DimSize{3},)}, HDF5.FixedArray{Float64,(HDF5.DimSize{3},)}, Float64]
@assert d.membername == ASCIIString["wgt", "xyz", "uvw", "E"]
@assert d.memberoffset == Uint64[0x00, 0x08, 0x20, 0x38]


text:
promote_rule(::Type{UncertainDataset}, ::Type{UncertainValueDataset}) = UncertainValueDataset
promote_rule(::Type{UncertainValueDataset}, ::Type{UncertainDataset}) = UncertainValueDataset

promote_rule(::Type{UncertainDataset}, ::Type{UncertainIndexDataset}) = UncertainValueDataset
promote_rule(::Type{UncertainIndexDataset}, ::Type{UncertainDataset}) = UncertainValueDataset

promote_rule(::Type{UncertainValueDataset}, ::Type{UncertainIndexDataset}) = UncertainValueDataset
promote_rule(::Type{UncertainIndexDataset}, ::Type{UncertainValueDataset}) = UncertainValueDataset

convert(::Type{UncertainValueDataset}, udata::T) where {T <: AbstractUncertainValueDataset} =
    UncertainValueDataset(udata.values)
convert(::Type{UncertainValueDataset}, udata::UncertainIndexDataset) =
    UncertainValueDataset(udata.indices)

convert(::Type{UncertainIndexDataset}, udata::T) where {T <: AbstractUncertainValueDataset} =
    UncertainIndexDataset(udata.values)
convert(::Type{UncertainIndexDataset}, udata::UncertainIndexDataset) =
    UncertainIndexDataset(udata.indices)

convert(::Type{UncertainDataset}, udata::T) where {T <: AbstractUncertainValueDataset} =
    UncertainDataset(udata.values)
convert(::Type{UncertainDataset}, udata::UncertainIndexDataset) =
    UncertainDataset(udata.indices)

# Converting vectors of any input to uncertain datasets
function convert(::Type{UncertainValueDataset}, uvec::AbstractVector)
    uvals = [UncertainValue(x) for x in uvec]
    UncertainValueDataset(uvals)
end

function convert(::Type{UncertainValueDataset},
        uvec::AbstractVector{T}) where {T<:AbstractUncertainValue}
    UncertainValueDataset(uvec)
end

UncertainValueDataset(udata::UncertainDataset) = convert(UncertainValueDataset, udata)
UncertainValueDataset(udata::UncertainIndexDataset) = convert(UncertainValueDataset, udata)
UncertainValueDataset(udata::UncertainValueDataset) = udata

UncertainIndexDataset(udata::UncertainDataset) = convert(UncertainIndexDataset, udata)
UncertainIndexDataset(udata::UncertainIndexDataset) = udata
UncertainIndexDataset(udata::UncertainValueDataset) = convert(UncertainIndexDataset, udata)

UncertainDataset(udata::UncertainDataset) = udata
UncertainDataset(udata::UncertainIndexDataset) = convert(UncertainDataset, udata)
UncertainDataset(udata::UncertainValueDataset) = convert(UncertainDataset, udata)


text:
module ConfigReader
    include("./FileReader.jl")
    import .FileReader

    include("./Chronometer.jl")
    import .Chronometer

    export read, get_configs

    """Функция получения конфигурации сервера из YAML или JSON файла.

    config_path: строка с полным путём до файла конфигурации
    """
    function read(config_path::String)::Union{Dict{String, Any}, Nothing}
        is_yaml::Bool = findfirst(".yml", config_path) === nothing ? false : true
        is_json::Bool = findfirst(".json", config_path) === nothing ? false : true
        config::Union{Dict{String, Any}, Nothing} = nothing
        try
            if is_yaml === true
                config = FileReader.from_yaml(config_path)
            elseif  is_json === true
                config = FileReader.from_json(config_path)
            else
                config = nothing
            end
        catch error
            @error Chronometer.message_with_time("Произошла ошибка при запуске сервера") error
            throw(error)
        end

        return config
    end

    function get_configs(path::String)
        common_config::Union{Dict{String, Any}, Nothing} = read(path)

        rabbit_config::Union{Dict{String, Any}, Nothing} = nothing
        channels_config::Union{Dict{String, Any}, Nothing} = nothing
        broker_config::Union{Dict{String, Any}, Nothing} = nothing

        broker_field::String = "broker"
        channels_field::String = "channels"

        default_value::Nothing = nothing

        if common_config !== nothing
            broker_config = get(common_config, broker_field, default_value)
            channels_config = get(common_config, channels_field, default_value)
        end

        return broker_config, channels_config
    end
end


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

const max_ccall_threads = parse(Int, get(ENV, "UV_THREADPOOL_SIZE", "4"))
const thread_notifiers = [Nullable{Condition}() for i in 1:max_ccall_threads]
const threadcall_restrictor = Semaphore(max_ccall_threads)

function notify_fun(idx)
    global thread_notifiers
    notify(get(thread_notifiers[idx]))
    return
end

function init_threadcall()
    global c_notify_fun = cfunction(notify_fun, Void, Tuple{Cint})
end

"""
    @threadcall((cfunc, clib), rettype, (argtypes...), argvals...)

The `@threadcall` macro is called in the same way as [`ccall`](@ref) but does the work
in a different thread. This is useful when you want to call a blocking C
function without causing the main `julia` thread to become blocked. Concurrency
is limited by size of the libuv thread pool, which defaults to 4 threads but
can be increased by setting the `UV_THREADPOOL_SIZE` environment variable and
restarting the `julia` process.

Note that the called function should never call back into Julia.
"""
macro threadcall(f, rettype, argtypes, argvals...)
    # check for usage errors
    isa(argtypes,Expr) && argtypes.head == :tuple ||
        error("threadcall: argument types must be a tuple")
    length(argtypes.args) == length(argvals) ||
        error("threadcall: wrong number of arguments to C function")

    # hygiene escape arguments
    f = esc(f)
    rettype = esc(rettype)
    argtypes = map(esc, argtypes.args)
    argvals = map(esc, argvals)

    # construct non-allocating wrapper to call C function
    wrapper = :(function wrapper(args_ptr::Ptr{Void}, retval_ptr::Ptr{Void})
        p = args_ptr
    end)
    body = wrapper.args[2].args
    args = Symbol[]
    for (i,T) in enumerate(argtypes)
        arg = Symbol("arg", i)
        push!(body, :($arg = unsafe_load(convert(Ptr{$T}, p))))
        push!(body, :(p += sizeof($T)))
        push!(args, arg)
    end
    push!(body, :(ret = ccall($f, $rettype, ($(argtypes...),), $(args...))))
    push!(body, :(unsafe_store!(convert(Ptr{$rettype}, retval_ptr), ret)))
    push!(body, :(return sizeof($rettype)))

    # return code to generate wrapper function and send work request thread queue
    :(let
        $wrapper
        do_threadcall(wrapper, $rettype, Any[$(argtypes...)], Any[$(argvals...)])
    end)
end

function do_threadcall(wrapper::Function, rettype::Type, argtypes::Vector, argvals::Vector)
    # generate function pointer
    fun_ptr = cfunction(wrapper, Int, Tuple{Ptr{Void}, Ptr{Void}})

    # cconvert, root and unsafe_convert arguments
    roots = Any[]
    args_size = isempty(argtypes) ? 0 : sum(sizeof, argtypes)
    args_arr = Vector{UInt8}(args_size)
    ptr = pointer(args_arr)
    for (T, x) in zip(argtypes, argvals)
        y = cconvert(T, x)
        push!(roots, y)
        unsafe_store!(convert(Ptr{T}, ptr), unsafe_convert(T, y))
        ptr += sizeof(T)
    end

    # create return buffer
    ret_arr = Vector{UInt8}(sizeof(rettype))

    # wait for a worker thread to be available
    acquire(threadcall_restrictor)
    idx = findfirst(isnull, thread_notifiers)
    thread_notifiers[idx] = Nullable{Condition}(Condition())

    # queue up the work to be done
    ccall(:jl_queue_work, Void,
        (Ptr{Void}, Ptr{UInt8}, Ptr{UInt8}, Ptr{Void}, Cint),
        fun_ptr, args_arr, ret_arr, c_notify_fun, idx)

    # wait for a result & return it
    wait(get(thread_notifiers[idx]))
    thread_notifiers[idx] = Nullable{Condition}()
    release(threadcall_restrictor)

    unsafe_load(convert(Ptr{rettype}, pointer(ret_arr)))
end

max_stars_repo_name: Mikewl/julia

max_issues_repo_name: Mikewl/julia

max_forks_repo_name: Mikewl/julia


text:
#!/usr/bin/env julia
# Tim Sterne-Weiler 2015

using Pkg

const dir = abspath( splitdir(@__FILE__)[1] )
const ver = readline(open(dir * "/../VERSION"))

start = time_ns()
println( stderr, "Whippet $ver loading and compiling... " )

Pkg.activate(dir * "/..")

using ArgParse
using Libz

using Whippet

function parse_cmd()
   s = ArgParseSettings()
   @add_arg_table s begin
    "--tpm"
      help     = "File containing two columns, isoform names, and transcripts-per-million (TpM)"
      arg_type = String
    "--index", "-x"
      help     = "Output prefix for saving index 'dir/prefix' (default Whippet/index/graph[.jls])"
      arg_type = String
      default  = fixpath( "$dir/../../index/graph" )
    "--out", "-o"
      help     = "Where should the gzipped output go 'dir/prefix'?"
      arg_type = String
      default  = fixpath( "./output" )
   end
   return parse_args(s)
end

function main()

   args = parse_cmd()

   println(stderr, " $( round( (time_ns()-start)/1e9, digits=6 ) ) seconds" )

   indexpath = fixpath( args["index"] )
   println(stderr, "Loading splice graph index... $( indexpath ).jls")
   @timer const lib = open(deserialize, "$( indexpath ).jls")

   data = load_tpm_file( fixpath(args["tpm"]) )

   tpms = Vector{Vector{Float64}}(length(lib.graphs))
   allocate_tpms!( tpms, lib )
   fill_tpms!( tpms, lib, data )

   io = open( args["out"], "w" )
   stream = ZlibDeflateOutputStream( io )

   output_tpms( stream, tpms, lib )

   close(stream)
   close(io)

   println(stderr, "Whippet $ver done." )
end

function load_tpm_file( filename::String, header=true )
   tpmfile = Dict{String,Float64}()
   fh = open( filename )
   for l in eachline(fh)
      if header
         header = false
         continue
      end
      s = split(l)
      tpmfile[String(s[1])] = parse(Float64, s[2])
   end
   tpmfile
end

function allocate_tpms!( tpms, lib )
   for i in 1:length(lib.graphs)
      tpms[i] = zeros(length(lib.graphs[i].annoname))
   end
   tpms
end

function fill_tpms!( tpms, lib, data )
   for i in 1:length(lib.graphs)
      for j in 1:length(lib.graphs[i].annoname)
         if haskey( data, lib.graphs[i].annoname[j] )
            tpms[i][j] = data[ lib.graphs[i].annoname[j] ]
         end
      end
   end
   tpms
end

function flanking_spliced( edges, node, name )
   upstream   = false
   downstream = false
   for i in node:-1:1
      if edges[i] in (EDGETYPE_RR, EDGETYPE_LR, EDGETYPE_SR)
         upstream = true
         #print("$name $(edges[i]):$i for n:$node")
         break
      elseif edges[i] in (EDGETYPE_LS, EDGETYPE_SL, EDGETYPE_RS)
         return false
      end
   end
   for i in node+1:length(edges)
      if edges[i] in (EDGETYPE_LL, EDGETYPE_LR, EDGETYPE_LS)
         downstream = true
         #println(" to $(edges[i]):$i")
         break
      elseif edges[i] in (EDGETYPE_SR, EDGETYPE_RS, EDGETYPE_SL)
         return false
      end
   end
   return (upstream && downstream)
end

function output_tpms( stream, tpms, lib )
   for i in 1:length(lib.graphs)
      for n in 1:length(lib.graphs[i].nodelen)
         if lib.graphs[i].edgetype[n]   in (EDGETYPE_RR, EDGETYPE_LR, EDGETYPE_LL) &&
            lib.graphs[i].edgetype[n+1] in (EDGETYPE_LL, EDGETYPE_LR, EDGETYPE_RR) &&
            flanking_spliced(lib.graphs[i].edgetype, n, lib.info[i].gene)

             tpm_total = 0.0
             tpm_node  = 0.0
             for j in 1:length(lib.graphs[i].annopath)
                curpath = lib.graphs[i].annopath[j]
                if n in first(curpath):last(curpath) # is within bounds
                   tpm_total += tpms[i][j]
                   if n in curpath
                      tpm_node += tpms[i][j]
                   end
                end
             end
             if tpm_total > 0.0
                psi = tpm_node / tpm_total
                write( stream, lib.info[i].gene * "\t" )
                write( stream, lib.info[i].name * ":" * string(lib.graphs[i].nodecoord[n]) *
                                                  "-" * string(lib.graphs[i].nodecoord[n]+lib.graphs[i].nodelen[n]-1) * "\t" )
                write( stream, string(n) * "\t" )
                write( stream, string(psi) * "\n" )
             end
         end
      end
   end
end

@timer main()


text:
__precompile__(true)

module Rows

export rows, columns

rows(M::Matrix) = map(x->reshape(getindex(M, x, :), :, size(M)[2]), 1:size(M)[1])
columns(M::Matrix) = map(x->getindex(M, :, x), 1:size(M)[2])

end # module


text:
using Funcle
using Test

@testset "Funcle.jl" begin
    # Write your tests here.
end


text:
function CC.OptimizationState(frame::InferenceState, params::OptimizationParams, interp::XInterpreter)
    ret = @invoke OptimizationState(frame::InferenceState, params::OptimizationParams, interp::AbstractInterpreter)

    return ret
end

# NOTE this one won't be called in the `XInterpreter` pipeline
function CC.OptimizationState(linfo::MethodInstance, src::CodeInfo, params::OptimizationParams, interp::XInterpreter)
    ret = @invoke OptimizationState(linfo::MethodInstance, src::CodeInfo, params::OptimizationParams, interp::AbstractInterpreter)

    return ret
end

# NOTE this one won't be called in the `XInterpreter` pipeline
function CC.OptimizationState(linfo::MethodInstance, params::OptimizationParams, interp::XInterpreter)
    ret = @invoke OptimizationState(linfo::MethodInstance, params::OptimizationParams, interp::AbstractInterpreter)

    return ret
end

function CC.finish(interp::XInterpreter, opt::OptimizationState, params::OptimizationParams, ir::IRCode, @nospecialize(result))
    ret = @invoke finish(interp::AbstractInterpreter, opt::OptimizationState, params::OptimizationParams, ir::IRCode, @nospecialize(result))

    return ret
end

function CC.optimize(interp::XInterpreter, opt::OptimizationState, params::OptimizationParams, @nospecialize(result))
    ret = @invoke optimize(interp::AbstractInterpreter, opt::OptimizationState, params::OptimizationParams, @nospecialize(result))

    return ret
end


text:
using Waveform
using JOLI
using Base.Test

# write your own tests here
include("test_3D.jl")


text:
using VarianceComponentModels, Distributions
"""
VarianceComponent
this VarianceComponent type stores A, B , CholA and CholB so we don't have to compute the cholesky decomposition inside the loop.
"""
struct VarianceComponent
	Σ::Matrix{Float64} # n_traits by n_traits
	V::Matrix{Float64} # n_people by n_people
	CholΣ::Array{Float64,2} # cholesky decomposition of A
	CholV::Array{Float64,2} # cholesky decomposition of B
	function VarianceComponent(Σ, V) #inner constructor given A, B
		CholΣ = cholesky(Symmetric(Σ)).factors
		CholV = cholesky(Symmetric(V)).factors
		return new(Σ, V, CholΣ, CholV) # stores these values (this is helpful so we don't have it inside the loop)
	end
end

function Base.show(io::IO, x::VarianceComponent)
    print(io, "Variance Component\n")
    print(io, "  * number of traits: $(ntraits(x))\n")
    print(io, "  * sample size: $(nsamplesize(x))")
end

#make our new type implement the interface defined above
nsamplesize(vc::VarianceComponent) = size(vc.V, 1)
ntraits(vc::VarianceComponent) = size(vc.Σ, 1)

# For a single Variance Component, algorithm that will transform standard normal distribution
function simulate_matrix_normal!(Z::Matrix, vc::VarianceComponent)
	randn!(Z)
	BLAS.trmm!('L', 'U', 'T', 'N', 1.0, vc.CholV, Z)
	BLAS.trmm!('R', 'U', 'N', 'N', 1.0, vc.CholΣ, Z)
	return Z #adds onto Z the effects of each variance component
end

"""
	VCM_trait_simulation(mu::Matrix{Float64}, vc::Vector{VarianceComponent})
For an evaluated mean matrix and vector of VarianceComponent objects, simulate from VCM.
"""
function VCM_trait_simulation(Y::Matrix, Z::Matrix, mu::Matrix{Float64}, vc::Vector{VarianceComponent}) # for an evaluated matrix
	for i in eachindex(vc)
		TraitSimulation.simulate_matrix_normal!(Z, vc[i]) # this step aggregates the variance components by
		axpy!(1.0, Z, Y) # summing the independent matrix normals to Y, rewriting over Z for each variance component
	end
	axpy!(1.0, mu, Y) # add the mean back to shift the matrix normal
	return Y
end


text:
using Documenter, DNC

makedocs(;
    modules=[DNC],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/sjrodahl/DNC.jl/blob/{commit}{path}#L{line}",
    sitename="DNC.jl",
    authors="Sondre Rodahl",
    assets=String[],
)

deploydocs(;
    repo="github.com/sjrodahl/DNC.jl",
)


text:
using DimensionalityReduction

# TODO: approximate principal components rather than principal components could be more efficient for N > T
# TODO: time series structure of factors: how many lags of the factors to include?
type DynamicFactorModel
    coefficients::Array{Float64, 1}
    coefficient_covariance::Array{Float64, 2}
    y::Array{Float64, 1}  # regressor
    w::Array{Float64, 2}  # regressands (e.g. lags of y, constant and variables known to affect y directly)
    x::Array{Float64, 2}  # variables to calculate factors from
    design_matrix::Array{Float64, 2}  # w and the factors
    targeted_predictors::BitArray
    number_of_factors::Int64  # columns of factors we use (which capture a certain percentage of the variation e.g.)
    number_of_factor_lags::Int64  # columns of factors we use (which capture a certain percentage of the variation e.g.)
    break_indices::Array{Int, 1}  # gives indices of breakpoints which are to be taken into account in estimation
    #rotation::Array{Float64, 2}  # rotation matrix to calculate factors from x (inverse of transpose of factor loadings lambda? Yes and also equal to factor loadings due to ortogonality?!)
    factors::Array{Array{Float64, 2}, 1} #Array{Float64, 2}
    loadings::Array{Array{Float64, 2}, 1} #Array{Float64, 2}
    t_stats::Array{Float64, 1}
    residuals::Array{Float64, 1}  # residuals of the regression of y on w and the factors
    factor_residuals::Array{Float64, 2}  # residuals from the factor estimation  x = factors * lambda
    factor_type::String
    number_of_factors_criterion::String
    number_of_factors_criterion_value::Float64

    # workhorse method, other methods exist which e.g. determine some of the arguments and then call this function
    function DynamicFactorModel(y::Array{Float64,1}, w::Matrix{Float64}, x::Matrix{Float64}, number_of_factors::Int64=ceil(minimum(size(x))/2), number_of_factors_criterion::String="", factor_type::String="principal components", targeted_predictors::BitArray=trues(size(x, 2)), number_of_factor_lags::Int64=0, break_indices::Array{Int64, 1}=Array(Int64, 0))
        # TODO: so far we only have a static factor model. factors need to be defined as in Stock, Watson (2010) page 3
        # TODO: include lagged factors into the regression (how many?)
        factors, loadings, number_of_factors = calculate_factors(x, factor_type, targeted_predictors, number_of_factors, break_indices)
        # TODO: check if correct factors are used (p 1136 on Bai Ng 2006)
        factor_residuals = x .- apply(vcat, [factors[i][:, 1:number_of_factors] * loadings[i][:, 1:number_of_factors]' for i in 1:length(factors)])  # we take account of possible breaks -> we calculate residuals per subperiod
        #factor_residuals = x - x*loadings[:, 1:number_of_factors]*inv(loadings[:, 1:number_of_factors]'loadings[:, 1:number_of_factors])*loadings[:, 1:number_of_factors]' # rescaled version from french dude
        # TODO: for lags we have to regress residuals on its lags and use info criterio to choose lag length p. this should be offloaded into another method
        # TODO: solve chicken and egg problem between number_of_factors and number_of_factor_lags. Both require residuals. --> maybe select one after the other repeatedly until convergence?
        # TODO: do the above e.g. factors, loadings, number_of_factors = augment_factors_by_lags(factors, number_of_factor_lags)  # TODO: do we have to
        # for structural breaks we augment the design matrix with factors multiplied with index which is 0 from t in [1, ..., t*-1] 1 from T in [t*, ..., T]

        design_matrix = make_factor_model_design_matrix(w, factors, number_of_factors, number_of_factor_lags, break_indices)
        coefficients = inv(design_matrix'design_matrix)*design_matrix'y
        residuals = y - design_matrix*coefficients
        hat_matrix = design_matrix*inv(design_matrix'design_matrix)*design_matrix'
        residual_variance = (residuals.^2)./(1.-diag(hat_matrix))  # HC 2
        #residual_variance = residuals.^2  # HC 0
        coefficient_covariance = inv(design_matrix'design_matrix)*(design_matrix'diagm(residual_variance)design_matrix)*inv(design_matrix'design_matrix)
        #coefficient_covariance = inv(design_matrix'design_matrix)*(1/length(residual_variance)*sum(residual_variance))
        t_stats = coefficients./(sqrt(diag(coefficient_covariance)))

        return(calculate_criterion(new(coefficients, coefficient_covariance, y, w, x, design_matrix, targeted_predictors, number_of_factors, number_of_factor_lags, break_indices, factors, loadings, t_stats, residuals, factor_residuals, factor_type, number_of_factors_criterion)))
    end

    function DynamicFactorModel(y::Array{Float64,1}, w::Matrix{Float64}, x::Matrix{Float64}, number_of_factors_criterion::String, factor_type::String="principal components", targeted_predictors::BitArray=trues(size(x, 2)), number_of_lags::Int=0, break_indices::Array{Int, 1}=Array(Int, 0))
        max_factors = int(ceil(minimum(size(x))/2))
        println("finding optimal DFM (with maximum number of factors=", max_factors, ")")
        # number of factors to include is not given but a criterion --> we use the criterion to determine the number of factors
        # if number of factors are to be calculated according to a criterion we need to estimate the model until criterion is optimal
        # we can probably stop when the criterion starts growing, I think the criterion functions are convex
        # TODO: what to do with chicken and egg problem: next line tries to
        # calculate number of factors using among other things the residuals
        # but number_of_factor lags depends on residuals as well
        models = [apply(DynamicFactorModel, (y, w, x, number_of_factors, number_of_factors_criterion, factor_type, targeted_predictors, number_of_lags, break_indices)) for number_of_factors in 1:max_factors]  # we keep all the models in memory which can be a problem depending on the dimensions of x. TODO: will refactor later when debugging and testing is done
        criteria = [model.number_of_factors_criterion_value for model in models]
        println("criteria_values: ", criteria)
        return models[indmin(criteria[1:max_factors])]  # keep the model with the best information criterion
    end
end

# calculate the factors and the rotation matrix to transform data into the space spanned by the factors
function calculate_factors(x::Matrix, factor_type::String="principal components", targeted_predictors=1:size(x, 2), number_of_factors=ceil(minimum(size(x))/2), break_indices=[])
    T, N = size(x)  # T: time dimension, N cross-sectional dimension
    break_indices = vcat(1, break_indices, T+1)

    function principal_components(x)  # calculate factors and loadings for a given x matrix
        # see Stock, Watson (1998)
        if T >= N
            eigen_values, eigen_vectors = eig(x'x)
            eigen_values, eigen_vectors = reverse(eigen_values), eigen_vectors[:, size(eigen_vectors, 2):-1:1]  # reverse the order from highest to lowest eigenvalue
            loadings = sqrt(N) * eigen_vectors  # we may rescale the eigenvectors say Bai, Ng 2002
            factors = x*loadings/N  # this is from Bai, Ng 2002 except that for me the inverse of the loadings matrix primed is not the same as the loadings matrix
            #factors = factors*(factors'factors/T)^(1/2)  # TODO: not sure it is correct to rescale like this (see Bai, Ng 2002 p.198)
            #factors = 1/N*x*loadings'  # this is from Stock, Watson 2010
            #loadings = (x'x)*loadings/(N*T)  # and this is from C. Hurlin from University of Orléans... rescales factor loadings to have the same variance as x
        end
        if N > T
            eigen_values, eigen_vectors = eig(x*x')
            eigen_values, eigen_vectors = reverse(eigen_values), eigen_vectors[:, size(eigen_vectors, 2):-1:1]  # reverse the order from highest to lowest eigenvalue
            factors = sqrt(T) * eigen_vectors  # sqrt(T) comes from normalization F'F/T = eye(r) where r is number of factors (see Bai 2003), factors are Txr
            loadings = x'factors/T  # see e.g. Bai 2003 p.6 or Breitung, Eickmeier 2011 p. 80. Dimension of loadings is Nxr where r here is also N
            # betahat = loadings * chol(loadings'loadings/N)  # C. Hurlin from University of Orléans rescales like this
        end
            # resultingly factors*loadings' estimates x
        return factors, loadings
    end

    if factor_type == "principal components"
        factors_loadings = [principal_components(x[break_indices[i-1]:break_indices[i]-1, :]) for i in 2:length(break_indices)]
        factors = convert(Array{Array{Float64, 2}, 1}, [factors_loadings[i][1] for i in 1:length(factors_loadings)])
        loadings = convert(Array{Array{Float64, 2}, 1}, [factors_loadings[i][2] for i in 1:length(factors_loadings)])
        max_factor_number = ceil(minimum(size(x))/2)
    elseif factor_type == "squared principal components"  # include squares of X
        pca_res = pca([x[:, targeted_predictors] x[:, targeted_predictors].^2]; center=false, scale=false)
        max_factor_number = minimum([size(x, 1), size(x,2)*2])
    elseif factor_type == "quadratic principal components"  # include squares of X and interaction terms - better only use in combination with targeted_predictors
        pca_cols = x[:, targeted_predictors]  # columns to use for principal components
        for i in 1:size(x[:, targeted_predictors], 2)
            for j in 1:size(x[:, targeted_predictors], 2)
                pca_cols = hcat(pca_cols, x[:, i].*x[:, j])
            end
        end
        max_factor_number = minimum(size(pca_cols))
        pca_res = pca(pca_cols; center=false, scale=false)
        # TODO: or replace with Bai Ng 2002: Determining the number of factors in approximate factor models
    end
    if number_of_factors > max_factor_number
        number_of_factors = max_factor_number
        warn("can not estimate more than `minimum(size(x))` factors with $factor_type. Number of factors set to $number_of_factors")
    end
    return factors, loadings, number_of_factors
end

# transforms x to the space spanned by the factors and optionally only selects active factors
#   type="active" returns only the active factors (which explain enough of the variance)
function get_factors(dfm::DynamicFactorModel, x::Matrix, factors="active")  # TODO: do I have to divide by N if T>N? (see Bai, Ng 2002 p. 198)
    rotation = dfm.loadings * inv(dfm.loadings'dfm.loadings) # simplifies to dfm.loadings if T>N and loadings = inv(loadings')  which makes sense given x = F * loadings'
    (normalize(x[:, dfm.targeted_predictors], (mean(dfm.x), std(dfm.x)))*dfm.rotation)[:, factors=="active" ? (1:dfm.number_of_factors) : (1:end)]
end

function make_factor_model_design_matrix(w, factors, number_of_factors, number_of_factor_lags, break_indices)
    return(hcat(w, apply(vcat, factors)[:, 1:number_of_factors]))  # we could assume breaks in the forecasting equation as well
    # TODO: unfinished business: factor lags to get dynamic factor models (how many? see page
end

function calculate_criterion(dfm::DynamicFactorModel)
    if dfm.number_of_factors_criterion != ""
        number_of_factors_criterion = dfm.number_of_factors_criterion
        dfm.number_of_factors_criterion_value = eval(symbol("criterion_$number_of_factors_criterion"))(dfm)
    end
    return(dfm)
end

function Base.show(io::IO, dfm::DynamicFactorModel)
    @printf io "Dynamic Factor Model\n"
    @printf io "Dimensions of X: %s\n" size(dfm.x)
    @printf io "Number of factors used: %s\n" dfm.number_of_factors
    @printf io "Factors calculated by: %s\n" dfm.factor_type
    @printf io "Factor model residual variance: %s\n" sum(dfm.factor_residuals.^2)/apply(*, size(dfm.x))
end

# prediction needs w and (normalized) x
function predict(dfm::DynamicFactorModel, w, x)
    design_matrix = hcat(w, get_factors(dfm, x))
    return design_matrix*dfm.coefficients
end

include("criteria.jl")  # defines the criteria in Bai and Ng 2002


text:
@inline function update_batch_edge(mp::T, adj, E::Fill{S,2,Axes}, X::CuMatrix, u) where {T<:MessagePassing,S,Axes}
    E = fill(E.value, E.axes)
    update_batch_edge(mp, adj, E, X, u)
end

@inline function update_batch_edge(mp::T, adj, E::CuMatrix, X::Fill{S,2,Axes}, u) where {T<:MessagePassing,S,Axes}
    X = fill(X.value, X.axes)
    update_batch_edge(mp, adj, E, X, u)
end

@inline function update_batch_edge(mp::T, adj, E::AbstractMatrix, X::CuMatrix, u) where {T<:MessagePassing}
    E = convert(typeof(X), E)
    update_batch_edge(mp, adj, E, X, u)
end

@inline function update_batch_edge(mp::T, adj, E::CuMatrix, X::AbstractMatrix, u) where {T<:MessagePassing}
    X = convert(typeof(E), X)
    update_batch_edge(mp, adj, E, X, u)
end

@inline function update_batch_edge(mp::T, adj, E::CuMatrix, X::CuMatrix, u) where {T<:MessagePassing}
    n = size(adj, 1)
    edge_idx = edge_index_table(adj)
    hcat([apply_batch_message(mp, i, adj[i], edge_idx, E, X, u) for i in 1:n]...)
end

@inline function apply_batch_message(mp::T, i, js, edge_idx, E::CuMatrix, X::CuMatrix, u) where {T<:MessagePassing}
    hcat([message(mp, get_feature(X, i), get_feature(X, j), get_feature(E, edge_idx[(i,j)])) for j = js]...)
end

@inline function update_batch_vertex(mp::T, M::AbstractMatrix, X::CuMatrix, u) where {T<:MessagePassing}
    M = convert(typeof(X), M)
    update_batch_vertex(mp, M, X, u)
end

@inline function update_batch_vertex(mp::T, M::CuMatrix, X::AbstractMatrix, u) where {T<:MessagePassing}
    X = convert(typeof(M), X)
    update_batch_vertex(mp, M, X, u)
end

@inline function update_batch_vertex(mp::T, M::CuMatrix, X::CuMatrix, u) where {T<:MessagePassing}
    hcat([update(mp, get_feature(M, i), get_feature(X, i)) for i = 1:size(X,2)]...)
end


text:
# syntax sugar
const U = Union

# TODO UN{t...} = U{Nothing, t...}
UN{t} = U{Nothing, t}


text:
module DependencyParsing
  using ..DependencyParser.Core

  include("settings.jl")
  include("parsing_system.jl")
  include("transition.jl")
  include("tree_node.jl")
  include("dependency_tree.jl")
  include("connlu.jl")
  include("configuration.jl")
  include("arc_eager/arc_eager.jl")
  include("training_context.jl")
  include("model.jl")
  include("dep_parser.jl")
end

max_stars_repo_name: LarsWl/DependencyParser

max_issues_repo_name: LarsWl/DependencyParser

max_forks_repo_name: LarsWl/DependencyParser


text:
function benchmark_correlate()
    results_min = DataFrame(
        Samples = SAMPLES,
        sCPU_time_1ant = zeros(Float64,length(SAMPLES)),
        sCPU_time_4ant = zeros(Float64,length(SAMPLES)),
        sCPU_time_16ant = zeros(Float64,length(SAMPLES)),
        GPU_time_1ant = zeros(Float64,length(SAMPLES)),
        GPU_time_4ant = zeros(Float64,length(SAMPLES)),
        GPU_time_16ant = zeros(Float64,length(SAMPLES)),
        sGPU_time_1ant = zeros(Float64,length(SAMPLES)),
        sGPU_time_4ant = zeros(Float64,length(SAMPLES)),
        sGPU_time_16ant = zeros(Float64,length(SAMPLES))
    )
    results_med = DataFrame(
        Samples = SAMPLES,
        sCPU_time_1ant = zeros(Float64,length(SAMPLES)),
        sCPU_time_4ant = zeros(Float64,length(SAMPLES)),
        sCPU_time_16ant = zeros(Float64,length(SAMPLES)),
        GPU_time_1ant = zeros(Float64,length(SAMPLES)),
        GPU_time_4ant = zeros(Float64,length(SAMPLES)),
        GPU_time_16ant = zeros(Float64,length(SAMPLES)),
        sGPU_time_1ant = zeros(Float64,length(SAMPLES)),
        sGPU_time_4ant = zeros(Float64,length(SAMPLES)),
        sGPU_time_16ant = zeros(Float64,length(SAMPLES))
    )
    results_mean = DataFrame(
        Samples = SAMPLES,
        sCPU_time_1ant = zeros(Float64,length(SAMPLES)),
        sCPU_time_4ant = zeros(Float64,length(SAMPLES)),
        sCPU_time_16ant = zeros(Float64,length(SAMPLES)),
        GPU_time_1ant = zeros(Float64,length(SAMPLES)),
        GPU_time_4ant = zeros(Float64,length(SAMPLES)),
        GPU_time_16ant = zeros(Float64,length(SAMPLES)),
        sGPU_time_1ant = zeros(Float64,length(SAMPLES)),
        sGPU_time_4ant = zeros(Float64,length(SAMPLES)),
        sGPU_time_16ant = zeros(Float64,length(SAMPLES))
    )
    rowpos = Int32(1)

    for N in SAMPLES
        for M in ANTENNA
            # init csv column write position
            M==1 ? columnpos=2 : M==4 ? columnpos=3 : M==16 ? columnpos=4 : throw("Error in antenna num")

            #init signals
            correlator = EarlyPromptLateCorrelator(0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im)
            gpsl1 = GNSSSignals.GPSL1()
            early_late_sample_shift = 1
            gpucode = CuArray(GNSSSignals.get_code(
                gpsl1,
                (1 - early_late_sample_shift:N + early_late_sample_shift) * 1023e3 / 2.5e6,
                1
            ))
            cpucode = Array{ComplexF32}(gpucode)
            cpudwnsignal = Array{ComplexF32}(cpucode .+ zeros(N+2) .* im)
            scpudwnsignal = StructArray{ComplexF32}((real(cpudwnsignal),imag(cpudwnsignal)))
            gpudwnsignal = CuArray{ComplexF32}(cpudwnsignal)
            sgpudwnsignal = StructArray{ComplexF32}((real(gpudwnsignal),imag(gpudwnsignal)))

            # init signals as matrices
            if M > 1
                correlator = EarlyPromptLateCorrelator(NumAnts(M))
                cpudwnsignal = zeros(ComplexF32, N, M)
                gpudwnsignal = CUDA.zeros(ComplexF32, N, M)
                scpudwnsignal = StructArray{ComplexF32}(
                    (real(zeros(Float32,N,M)),
                    imag(zeros(Float32,N,M))))
                sgpudwnsignal = StructArray{ComplexF32}(
                    (real(CUDA.zeros(Float32,N,M)),
                    imag(CUDA.zeros(Float32,N,M))))
            end

            # do the benchmark
            println("Benchmarking correlate on CPU: StructArray{ComplexF32} ", N," samples, ", M, " antenna...")
            result = @benchmark gpu_correlate(
                $correlator,
                $scpudwnsignal,
                $cpucode,
                $early_late_sample_shift,
                1,
                $N,
                1.0,
                2,
                Val(7)
            )
            println(minimum(result).time)
            results_min[rowpos, columnpos] = minimum(result).time
            results_med[rowpos, columnpos] = median(result).time
            results_mean[rowpos, columnpos] = mean(result).time

            println("Benchmarking the correlator on GPU: CuArray{ComplexF32} ", N," samples, ", M, " antenna...")
            result = @benchmark CUDA.@sync gpu_correlate(
                $correlator,
                $gpudwnsignal,
                $gpucode,
                $early_late_sample_shift,
                1,
                $N,
                1.0,
                2,
                Val(7),
            )
            println(minimum(result).time)
            results_min[rowpos, columnpos+3] = minimum(result).time
            results_med[rowpos, columnpos+3] = median(result).time
            results_mean[rowpos, columnpos+3] = mean(result).time

            println("Benchmarking downconvert on GPU: StructArray{ComplexF32}(CuArray,CuArray) ", N, " samples, ", M, " antenna...")
            result = @benchmark CUDA.@sync sgpu_correlate(
                $correlator,
                $sgpudwnsignal.re,
                $sgpudwnsignal.im,
                $gpucode,
                $early_late_sample_shift,
                1,
                $N,
                1.0,
                2,
                Val(7),
            )
            println(minimum(result).time)
            results_min[rowpos, columnpos+6] = minimum(result).time
            results_med[rowpos, columnpos+6] = median(result).time
            results_mean[rowpos, columnpos+6] = mean(result).time

        end
        rowpos += 1
    end
    CSV.write("data/correlate_min.csv", results_min)
    CSV.write("data/correlate_med.csv", results_med)
    CSV.write("data/correlate_mean.csv", results_mean)
end