text:
################################################################################
#                                 Quaternions
################################################################################

export dquat, eye, norm, quat_to_angle, quat_to_angleaxis, quat_to_dcm, vect

################################################################################
#                                 Initializers
################################################################################

"""
    function Quaternion(q0::T0, q1::T1, q2::T2, q3::T3) where {T0,T1,T2,T3}

Create the following quaternion:

    q0 + q1.i + q2.j + q3.k

in which:

* `q0` is the real part of the quaternion.
* `q1` is the X component of the quaternion vectorial part.
* `q2` is the Y component of the quaternion vectorial part.
* `q3` is the Z component of the quaternion vectorial part.

"""
function Quaternion(q0::T0, q1::T1, q2::T2, q3::T3) where {T0,T1,T2,T3}
    T = promote_type(T0,T1,T2,T3)
    Quaternion{T}(q0,q1,q2,q3)
end

"""
    function Quaternion(v::AbstractVector)

If the vector `v` has 3 components, then create a quaternion in which the real
part is `0` and the vectorial or imaginary part has the same components of the
vector `v`. In other words:

    q = 0 + v[1].i + v[2].j + v[3].k

Otherwise, if the vector `v` has 4 components, then create a quaternion in which
the elements match those of the input vector:

    q = v[1] + v[2].i + v[3].j + v[4].k

"""
function Quaternion(v::AbstractVector)
    # The vector must have 3 or 4 components.
    if length(v) != 3 && length(v) != 4
        throw(ArgumentError("The input vector must have 3 or 4 components."))
    end

    if length(v) == 3
        return Quaternion(0, v[1], v[2], v[3])
    else
        return Quaternion(v[1], v[2], v[3], v[4])
    end
end

"""
    function Quaternion(r::Number, v::AbstractVector)

Create a quaternion with real part `r` and vectorial or imaginary part `v`:

    r + v[1].i + v[2].j + v[3].k

"""
Quaternion(r::Number, v::AbstractVector) = Quaternion(r, v[1], v[2], v[3])

"""
    function Quaternion(u::UniformScaling{T}) where T
    function Quaternion{T}(u::UniformScaling) where T

Create the quaternion `u.λ + 0.i + 0.j + 0.k`.

"""
Quaternion(u::UniformScaling{T}) where T = Quaternion{T}(T(u.λ), T(0), T(0), T(0))
Quaternion{T}(u::UniformScaling) where T = Quaternion{T}(T(u.λ), T(0), T(0), T(0))

"""
    function Quaternion(::UniformScaling,::Quaternion{T}) where T

Create an identity quaternion of type `T`:

    T(1) + T(0).i + T(0).j + T(0).k

"""
Quaternion(::UniformScaling,::Quaternion{T}) where T = Quaternion{T}(I)

################################################################################
#                                  Operations
################################################################################

# Operation: +
# ==============================================================================

"""
    @inline function +(qa::Quaternion, qb::Quaternion)

Compute `qa + qb`.

"""
@inline +(qa::Quaternion, qb::Quaternion) =
    Quaternion(qa.q0 + qb.q0, qa.q1 + qb.q1, qa.q2 + qb.q2, qa.q3 + qb.q3)

"""
    @inline function +(u::UniformScaling, q::Quaternion)
    @inline function +(q::Quaternion, u::UniformScaling)

Compute `qu + q` or `q + qu`, in which `qu` is the scaled identity quaternion
`qu = u.λ * I`.

"""
@inline +(u::UniformScaling, q::Quaternion) =
    Quaternion(u.λ + q.q0, q.q1, q.q2, q.q3)

@inline +(q::Quaternion, u::UniformScaling) = u+q

# Operation: -
# ==============================================================================

"""
    @inline function -(qa::Quaternion, qb::Quaternion)

Compute `qa - qb`.

"""
@inline -(qa::Quaternion, qb::Quaternion) =
    Quaternion(qa.q0 - qb.q0, qa.q1 - qb.q1, qa.q2 - qb.q2, qa.q3 - qb.q3)

"""
    @inline function -(u::UniformScaling, q::Quaternion)
    @inline function -(q::Quaternion, u::UniformScaling)

Compute `qu - q` or `q - qu`, in which `qu` is the scaled identity quaternion
`qu = u.λ * I`.

"""
@inline -(u::UniformScaling, q::Quaternion) =
    Quaternion{T}(u.λ - q.q0, -q.q1, -q.q2, -q.q3)

@inline -(q::Quaternion, u::UniformScaling) = (-u)+q

# Operation: *
# ==============================================================================

"""
    @inline function *(λ::Number, q::Quaternion)
    @inline function *(q::Quaternion, λ::Number)

Compute `λ*q` or `q*λ`, in which `λ` is a scalar.

"""
@inline *(λ::Number, q::Quaternion) = Quaternion(λ*q.q0, λ*q.q1, λ*q.q2, λ*q.q3)
@inline *(q::Quaternion, λ::Number) = Quaternion(λ*q.q0, λ*q.q1, λ*q.q2, λ*q.q3)

"""
    @inline function *(q1::Quaternion, q2::Quaternion)

Compute the quaternion multiplication `q1*q2` (Hamilton product).

"""
@inline *(q1::Quaternion, q2::Quaternion) =
    Quaternion(q1.q0*q2.q0 - q1.q1*q2.q1 - q1.q2*q2.q2 - q1.q3*q2.q3,
               q1.q0*q2.q1 + q1.q1*q2.q0 + q1.q2*q2.q3 - q1.q3*q2.q2,
               q1.q0*q2.q2 - q1.q1*q2.q3 + q1.q2*q2.q0 + q1.q3*q2.q1,
               q1.q0*q2.q3 + q1.q1*q2.q2 - q1.q2*q2.q1 + q1.q3*q2.q0)

"""
    @inline function *(v::AbstractVector, q::Quaternion)
    @inline function *(q::Quaternion, v::AbstractVector)

Compute the multiplication `qv*q` or `q*qv` in which `qv` is a quaternion with
real part `0` and vectorial/imaginary part `v` (Hamilton product).

"""
@inline *(v::AbstractVector, q::Quaternion) =
    Quaternion(-v[1]*q.q1 - v[2]*q.q2 - v[3]*q.q3,
               +v[1]*q.q0 + v[2]*q.q3 - v[3]*q.q2,
               -v[1]*q.q3 + v[2]*q.q0 + v[3]*q.q1,
               +v[1]*q.q2 - v[2]*q.q1 + v[3]*q.q0)

@inline *(q::Quaternion, v::AbstractVector) =
    Quaternion(           - q.q1*v[1] - q.q2*v[2] - q.q3*v[3],
               q.q0*v[1]              + q.q2*v[3] - q.q3*v[2],
               q.q0*v[2] - q.q1*v[3]              + q.q3*v[1],
               q.q0*v[3] + q.q1*v[2] - q.q2*v[1]             )

"""
    @inline function *(u::UniformScaling, q::Quaternion)
    @inline function *(q::Quaternion, u::UniformScaling)

Compute `qu*q` or `q*qu` (Hamilton product), in which `qu` is the scaled
identity quaternion `qu = u.λ * I`.

"""
@inline *(u::UniformScaling, q::Quaternion) = Quaternion(u)*q
@inline *(q::Quaternion, u::UniformScaling) = q*Quaternion(u)

# Operation: /
# ==============================================================================

"""
    @inline function /(λ::Number, q::Quaternion)
    @inline function /(q::Quaternion, λ::Number)

Compute the division `λ/q` or `q/λ`, in which `λ` is a scalar.

"""
@inline function /(λ::Number, q::Quaternion)
    # Compute `λ*(1/q)`.
    norm_q² = q.q0*q.q0 + q.q1*q.q1 + q.q2*q.q2 + q.q3*q.q3
    Quaternion(λ*q.q0/norm_q², -λ*q.q1/norm_q², -λ*q.q2/norm_q², -λ*q.q3/norm_q²)
end

@inline /(q::Quaternion, λ::Number) = q*(1/λ)

"""
    @inline /(q1::Quaternion, q2::Quaternion) = q1*inv(q2)

Compute `q1*inv(q2)` (Hamilton product).

"""
@inline /(q1::Quaternion, q2::Quaternion) = q1*inv(q2)

"""
    @inline function /(u::UniformScaling, q::Quaternion)
    @inline function /(q::Quaternion, u::UniformScaling)

Compute `qu/q` or `q/qu` (Hamilton product), in which `qu` is the scaled
identity quaternion `qu = u.λ * I`.

"""
@inline /(u::UniformScaling, q::Quaternion) = Quaternion(u)/q
@inline /(q::Quaternion, u::UniformScaling) = q/Quaternion(u)

# Operation: \
# ==============================================================================

"""
    @inline \\(q1::Quaternion, q2::Quaternion) = inv(q1)*q2

Compute `inv(q1)*q2`.

"""
@inline \(q1::Quaternion, q2::Quaternion) = inv(q1)*q2

"""
    @inline \\(q::Quaternion, v::AbstractVector)
    @inline \\(v::AbstractVector, q::Quaternion)

Compute `inv(q)*qv` or `inv(qv)*q` in which `qv` is a quaternion with real part
`0` and vectorial/imaginary part `v` (Hamilton product).

"""
@inline \(q::Quaternion,     v::AbstractVector) = inv(q)*v
@inline \(v::AbstractVector, q::Quaternion)     = inv(Quaternion(v))*q

"""
    @inline function \\(u::UniformScaling, q::Quaternion)
    @inline function \\(q::Quaternion, u::UniformScaling)

Compute `inv(qu)*q` or `inv(q)*qu` (Hamilton product), in which `qu` is the
scaled identity quaternion `qu = u.λ * I`.

"""
@inline \(u::UniformScaling, q::Quaternion) = inv(Quaternion(u))*q
@inline \(q::Quaternion, u::UniformScaling) = inv(q)*Quaternion(u)

# Operation: [:]
# ==============================================================================

"""
    @inline function getindex(q::Quaternion, ::Colon)

Transform the quaternion into a 4x1 vector of type `T`.

"""
@inline getindex(q::Quaternion, ::Colon) = [q.q0;q.q1;q.q2;q.q3]

################################################################################
#                                  Functions
################################################################################

"""
    @inline function conj(q::Quaternion)

Compute the complex conjugate of the quaternion `q`:

    q0 - q1.i - q2.j - q3.k

"""
@inline conj(q::Quaternion) = Quaternion(q.q0, -q.q1, -q.q2, -q.q3)

"""
    @inline function copy(q::Quaternion{T}) where T

Create a copy of the quaternion `q`.

"""
@inline copy(q::Quaternion{T}) where T = Quaternion{T}(q.q0, q.q1, q.q2, q.q3)

"""
    @inline function imag(q::Quaternion)

Return the vectorial or imaginary part of the quaternion `q` represented by a
3 × 1 vector of type `SVector{3}`.

"""
@inline imag(q::Quaternion) = SVector{3}(q.q1, q.q2, q.q3)

"""
    @inline function inv(q::Quaternion)

Compute the inverse of the quaternion `q`:

    conj(q)
    -------
      |q|²

"""
@inline function inv(q::Quaternion)
    # Compute the inverse of the quaternion.
    norm_q² = q.q0*q.q0 + q.q1*q.q1 + q.q2*q.q2 + q.q3*q.q3
    Quaternion(q.q0/norm_q², -q.q1/norm_q², -q.q2/norm_q², -q.q3/norm_q²)
end

"""
    @inline function norm(q::Quaternion)

Compute the Euclidean norm of the quaternion `q`:

    sqrt(q0² + q1² + q2² + q3²)

"""
@inline norm(q::Quaternion) = sqrt(q.q0*q.q0 + q.q1*q.q1 + q.q2*q.q2 + q.q3*q.q3)

"""
    @inline function real(q::Quaternion)

Return the real part of the quaternion `q`: `q0`.

"""
@inline real(q::Quaternion) = q.q0

"""
    @inline function vect(q::Quaternion)

Return the vectorial or imaginary part of the quaternion `q` represented by a
3 × 1 vector of type `SVector{3}`.

"""
@inline vect(q::Quaternion) = SVector{3}(q.q1, q.q2, q.q3)

"""
    @inline function zeros(::Type{Quaternion{T}}) where T

Create the null quaternion of type `T`:

    T(0) + T(0).i + T(0).j + T(0).k

If the type `T` is omitted, then it defaults to `Float64`.

# Example

```julia-repl
julia> zeros(Quaternion{Float32})
Quaternion{Float32}:
  + 0.0 + 0.0.i + 0.0.j + 0.0.k

julia> zeros(Quaternion)
Quaternion{Float64}:
  + 0.0 + 0.0.i + 0.0.j + 0.0.k
```

"""
@inline zeros(::Type{Quaternion{T}}) where T =
    Quaternion{T}(T(0), T(0), T(0), T(0))

@inline zeros(::Type{Quaternion}) = Quaternion{Float64}(0.0,0.0,0.0,0.0)

"""
    @inline function zeros(q::Quaternion{T}) where T

Create the null quaternion with the same type `T` of another quaternion `q`:

    T(0) + T(0).i + T(0).j + T(0).k

# Example

```julia-repl
julia> q1 = Quaternion{Float32}(cosd(45/2),sind(45/2),0,0);

julia> zeros(q1)
Quaternion{Float32}:
  + 0.0 + 0.0.i + 0.0.j + 0.0.k
```

"""
@inline zeros(q::Quaternion{T}) where T = zeros(Quaternion{T})

################################################################################
#                                      IO
################################################################################

"""
    function show(io::IO, q::Quaternion{T}) where T
    function show(io::IO, mime::MIME"text/plain", q::Quaternion{T}) where T

Print the quaternion `q` to the stream `io`.

"""
function show(io::IO, q::Quaternion{T}) where T
    # Get the absolute values.
    aq0 = abs(q.q0)
    aq1 = abs(q.q1)
    aq2 = abs(q.q2)
    aq3 = abs(q.q3)

    # Get the signs.
    sq0 = (q.q0 >= 0) ? "+" : "-"
    sq1 = (q.q1 >= 0) ? "+" : "-"
    sq2 = (q.q2 >= 0) ? "+" : "-"
    sq3 = (q.q3 >= 0) ? "+" : "-"

    print(io, "Quaternion{$(T)}:")
    print(io, " $(sq0) $(aq0) $(sq1) $(aq1).i $(sq2) $(aq2).j $(sq3) $(aq3).k")

    nothing
end

function show(io::IO, mime::MIME"text/plain", q::Quaternion{T}) where T
    # Check if the user wants colors.
    color = get(io, :color, false)

    b = (color) ? "\x1b[1m"         : ""
    d = (color) ? "\x1b[0m"         : ""
    g = (color) ? "\x1b[1m\x1b[32m" : ""
    y = (color) ? "\x1b[1m\x1b[33m" : ""

    # Get the absolute values.
    aq0 = abs(q.q0)
    aq1 = abs(q.q1)
    aq2 = abs(q.q2)
    aq3 = abs(q.q3)

    # Get the signs.
    sq0 = (q.q0 >= 0) ? "+" : "-"
    sq1 = (q.q1 >= 0) ? "+" : "-"
    sq2 = (q.q2 >= 0) ? "+" : "-"
    sq3 = (q.q3 >= 0) ? "+" : "-"

    # Unitary vectors.
    i = "$(b)i$d"
    j = "$(b)j$d"
    k = "$(b)k$d"

    println(io, "Quaternion{$(T)}:")
    print(io, "  $(sq0) $(aq0) $(sq1) $(aq1).$i $(sq2) $(aq2).$j $(sq3) $(aq3).$k")

    nothing
end

################################################################################
#                                 Conversions
################################################################################

# Direction Cosine Matrix
# ==============================================================================

"""
    function quat_to_dcm(q::Quaternion)

Convert the quaternion `q` to a Direction Cosine Matrix (DCM).

# Example

```julia-repl
julia> q = Quaternion(cosd(45/2), sind(45/2), 0, 0);

julia> quat_to_dcm(q)
3×3 StaticArrays.SArray{Tuple{3,3},Float64,2,9}:
 1.0   0.0       0.0
 0.0   0.707107  0.707107
 0.0  -0.707107  0.707107
```

"""
function quat_to_dcm(q::Quaternion)
    # Auxiliary variables.
    q0 = q.q0
    q1 = q.q1
    q2 = q.q2
    q3 = q.q3

    DCM(q0^2+q1^2-q2^2-q3^2,   2(q1*q2+q0*q3)   ,   2(q1*q3-q0*q2),
          2(q1*q2-q0*q3)   , q0^2-q1^2+q2^2-q3^2,   2(q2*q3+q0*q1),
          2(q1*q3+q0*q2)   ,   2(q2*q3-q0*q1)   , q0^2-q1^2-q2^2+q3^2)'
end

# Euler Angle and Axis
# ==============================================================================

"""
    function quat_to_angleaxis(q::Quaternion{T}) where T

Convert the quaternion `q` to a Euler angle and axis representation (see
`EulerAngleAxis`). By convention, the Euler angle will be kept between [0, π]
rad.

# Remarks

This function will not fail if the quaternion norm is not 1. However, the
meaning of the results will not be defined, because the input quaternion does
not represent a 3D rotation. The user must handle such situations.

# Example

```julia-repl
julia> q = Quaternion(cosd(45/2), sind(45/2), 0, 0);

julia> quat_to_angleaxis(q)
EulerAngleAxis{Float64}(0.7853981633974484, [1.0, 0.0, 0.0])
```

"""
function quat_to_angleaxis(q::Quaternion{T}) where T

    # If `q0` is 1 or -1, then we have an identity rotation.
    if abs(q.q0) >= 1 - eps()
        return EulerAngleAxis( T(0), SVector{3,T}(0, 0, 0) )
    else
        # Compute sin(θ/2).
        sθo2 = sqrt( q.q1*q.q1 + q.q2*q.q2 + q.q3*q.q3 )

        # Compute θ in range [0, 2π].
        θ = 2acos(q.q0)

        # Keep θ between [0, π].
        s = +1
        if θ > π
            θ = 2π - θ
            s = -1
        end

        return EulerAngleAxis( θ, s*[q.q1, q.q2, q.q3]/sθo2 )
    end

end

# Euler Angles
# ==============================================================================

"""
    function quat_to_angle(q::Quaternion, rot_seq::Symbol = :ZYX)

Convert the quaternion `q` to Euler Angles (see `EulerAngles`) given a rotation
sequence `rot_seq`.

The rotation sequence is defined by a `:Symbol`. The possible values are:
`:XYX`, `XYZ`, `:XZX`, `:XZY`, `:YXY`, `:YXZ`, `:YZX`, `:YZY`, `:ZXY`, `:ZXZ`,
`:ZYX`, and `:ZYZ`. If no value is specified, then it defaults to `:ZYX`.

# Example

```julia-repl
julia> q = Quaternion(cosd(45/2), sind(45/2), 0, 0);

julia> quat_to_angle(q,:XYZ)
EulerAngles{Float64}(0.7853981633974484, 0.0, -0.0, :XYZ)
```

"""
function quat_to_angle(q::Quaternion, rot_seq::Symbol=:ZYX)
    # TODO: This function calls uses `angle_to_dcm` to convert the quaternion to
    # a Direction Cosine Matrix. It **must** be rewritten to avoid this
    # intermediate step to increase the performance.

    # Convert the quaternion to DCM.
    dcm = quat_to_dcm(q)

    # Convert the DCM to the Euler Angles.
    dcm_to_angle(dcm, rot_seq)
end

################################################################################
#                                  Kinematics
################################################################################

"""
    function dquat(qba::Quaternion, wba_b::AbstractVector)

Compute the time-derivative of the quaternion `qba` that rotates a reference
frame `a` into alignment to the reference frame `b` in which the angular
velocity of `b` with respect to `a`, and represented in `b`, is `wba_b`.

# Example

```julia-repl
julia> q = Quaternion(1.0I);

julia> dquat(q,[1;0;0])
Quaternion{Float64}:
  + 0.0 + 0.5.i + 0.0.j + 0.0.k
```

"""
function dquat(qba::Quaternion, wba_b::AbstractVector)
    # Auxiliary variable.
    w = wba_b

    # Check the dimensions.
    if length(wba_b) != 3
        throw(ArgumentError("The angular velocity vector must have three components."))
    end

    # Return the time-derivative.
    #         1         x
    #   dq = --- (wba_b) . q
    #         2

    Quaternion(               -w[1]/2*qba.q1 -w[2]/2*qba.q2 -w[3]/2*qba.q3,
               +w[1]/2*qba.q0                +w[3]/2*qba.q2 -w[2]/2*qba.q3,
               +w[2]/2*qba.q0 -w[3]/2*qba.q1                +w[1]/2*qba.q3,
               +w[3]/2*qba.q0 +w[2]/2*qba.q1 -w[1]/2*qba.q2                )
end


text:
@testset "Determinant by minor expansion" begin
    for d = 1:5
        for testcase = 1:10
            mat_space = Nemo.MatrixSpace(Nemo.QQ, d, d)
            rnd_matrix = mat_space([mod(rand(Int), 1000) for i = 1:d, j = 1:d])
            @test det(rnd_matrix) == det_minor_expansion(rnd_matrix)
        end
    end
end

ext: jl


text:
using DiffEqSensitivity
using QuasiMonteCarlo
using CairoMakie

Nsamples = 10000
Nparams = length(parameter_set.names)
total_samples = Nsamples * (Nparams + 2)
@show total_samples

# Directory to which to save the files generated in this script
dir = joinpath(directory, "sensitivity_analysis")
isdir(dir) || mkdir(dir)

# Generate design matrices where rows are samples
@assert parameter_set.names == training.free_parameters.names
bds = [bounds(name, parameter_set) for name in parameter_set.names]
lb = getindex.(bds, 1)
ub = getindex.(bds, 2)
sampler = SobolSample()
A,B = QuasiMonteCarlo.generate_design_matrices(Nsamples,lb,ub,sampler,2)

big_training = inverse_problem(four_day_suite_path_2m, Nsamples, training_times)
big_eki = EnsembleKalmanInversion(big_training; noise_covariance, pseudo_stepping, resampler, tikhonov = true)

using ParameterEstimocean.EnsembleKalmanInversions: eki_objective

# Given `Nsamples`, the algorithm will compute sensitivity indices
# based on `Nsamples * (Nparams + 2)` samples of the forward map.
function f(p)
    @show size(p)

    Φs = []
    for i = 1:(Nparams+2)
        last = i*Nsamples
        X_sobol = p[:, (last-Nsamples+1):(last)] # Nparam x Nsamples
        G_sobol = forward_map(big_training, X_sobol) # Noutput x Nsamples
        Φs = vcat(Φs, [sum(eki_objective(eki, X_sobol[:,j], G_sobol[:,j]; constrained = true)) for j in 1:Nsamples])
        @show Φs
    end

    return Φs
end

sobol_result = gsa(f, Sobol(), A, B; batch=true)

# Figure and Axis
fig = Figure(fontsize=24, resolution=(1200,400))
ax = Axis(fig[1,1], xticks = (1:Nparams, string.(collect(parameter_set.names))),
        ylabel = "Sensitivity Indices",
        xlabel = "Parameter",
        xticklabelrotation = pi/2)

# Plot
colors = Makie.wong_colors()
first_order_indices = sobol_result.S1
total_order_indices = sobol_result.ST
x = vcat([1:Nparams...], [1:Nparams...])
height = vcat(first_order_indices, total_order_indices)

grp = vcat([1 for _ in first_order_indices], [2 for _ in total_order_indices])
barplot!(ax, x, height, dodge = grp, color = colors[grp])

# Legend
labels = ["First Order", "Total Order"]
elements = [PolyElement(polycolor = colors[i]) for i in eachindex(labels)]
title = "Index Type"
Legend(fig[1,2], elements, labels, title)

save(joinpath(dir, "sensitivity_indices.png"), fig)

# using ParameterEstimocean.Parameters: unconstrained_prior
# free_parameters = training.free_parameters
# unconstrained_priors = NamedTuple(name => unconstrained_prior(free_parameters.priors[name]) for name in free_parameters.names)
# unconstrained_parameters = [rand(unconstrained_priors[i]) for i=1:Nθ, k=1:Nsamples]
# y = observation_map(big_training)

# using FileIO
# file = joinpath(dir, "G.jld2")
# G_sobol = forward_map(big_training, params)
# save(file, G)
# G_sobol = load(file)["G"]


text:
using DataFrames
using CSV
using Printf
using Random
using Dates

using POMDPs
using ParticleFilters
using MCTS
using POMCPOW
using QMDP
using DiscreteValueIteration
using POMDPSimulators
using POMDPPolicies

using SunbergTypes
using IPFT
using Scenarios

using Plots
pyplot()

seed = Random.GLOBAL_RNG.seed[1]

@show max_time = 1.0
@show max_depth = 30
@show initial_state = -5.0

pomdp = ContinuousLightDark(big_step=3, transition_std=0.1, observation_std_base=0.5)
@show pomdp

solvers = Dict{String, Union{Solver,Policy}}(

    "pomcpow" => begin
        rng = MersenneTwister(seed)
        solver = POMCPOWSolver(tree_queries=10_000_000,
                               criterion=MaxUCB(90.0),
                               max_depth=max_depth,
                               max_time=max_time,
                               enable_action_pw=false,
                               k_observation=5.,
                               alpha_observation=1/15.0,
                               estimate_value=RolloutEstimator(RandomSolver(rng)),
                               check_repeat_obs=false,
                               tree_in_info=false,
                               rng=rng
                              )
        solve(solver, pomdp)
    end,

    "ipft" => begin
        rng = MersenneTwister(seed)
        m = 20
        ifm = GMMEntropy(SilvermansRule())
        discount_information_gain = false
        solver = IPFTSolver(n_iterations=typemax(Int),
                        criterion=MultiObjectiveUCB([1.0,60.0],100.0),
                        depth=max_depth,
                        max_time=max_time,
                        check_repeat_action=false,
                        reuse_states=false,
                        k_observation = 5.0,
                        alpha_observation = 1/20,
                        estimate_value=RolloutEstimator(RandomSolver(rng)),
                        enable_action_pw=false,
                        tree_in_info=true,
                        rng=rng
                        )
        belief_mdp = InformationRewardBeliefMDP(pomdp,
                                         LowVarianceResampler(m),
                                         ifm,
                                         discount_information_gain,
                                         0.05
                                        )
        solve(solver, belief_mdp)
    end,
)

histories = []

state_histories=Dict()

for (k, planner) in solvers
    @show k
    state_histories[k] = []

    be_updater = planner isa InformationMCTSPlanner

    Random.seed!(planner, seed+50_000)

    N_p = 10000
    up = ObsAdaptiveParticleFilter(pomdp,
                                       LowVarianceResampler(N_p),
                                       0.05, MersenneTwister(seed+90_000))
    if be_updater
        up = InformationFilter(up, planner.mdp.im)
    end

    # hr = HistoryRecorder(max_steps=100)
    # h = simulate(hr, pomdp, planner, up, initialstate_distribution(pomdp), -5.0)
    # push!(histories, h)
    # display(plot(h.state_hist))
    for (s,a,r,sp,o) in stepthrough(pomdp, planner, up, initialstate_distribution(pomdp), initial_state, "sarspo", max_steps=50)
        @show s, a, r, sp, o
        push!(state_histories[k], s)
    end
    display(plot(state_histories[k]))
end


text:
# implementation
function start_time(file)
  VideoIO.FFMPEG.ffprobe() do exe
    txt = read(`$exe  -v error -show_entries format=start_time -of default=noprint_wrappers=1:nokey=1 $file`, String)
    parse(Float64, txt)
  end
end

struct Segment
  file::String
  video::VideoIO.StreamContext
  start_time::Float64
  Segment(file::String, opts...) = new(file, VideoIO.openvideo(file, opts...), start_time(file))
end

mutable struct SegmentedVideo
  segments::Vector{Segment}
  current::Int
  SegmentedVideo(files::Vector{AbstractString}, opts...) = new(Segment.(files, opts...), 1)
end

VideoIO.openvideo(files::Vector{AbstractString}, opts...) = SegmentedVideo(files, opts...)

get_current(v::SegmentedVideo, f::Symbol) = getfield(v.segments[v.current], f)

VideoIO.aspect_ratio(v::SegmentedVideo) = VideoIO.aspect_ratio(get_current(v, :video))

function VideoIO.seek(v::SegmentedVideo, t::Number)
  i = findlast(s -> s.start_time < t, v.segments)
  v.current = isnothing(i) ? 1 : i
  seek(get_current(v, :video), t)
end

function VideoIO.seek(v::SegmentedVideo, file::String, t::Number)
  v.current = findfirst(s -> s.file == file, v.segments)
  seek(get_current(v, :video), get_current(v, :start_time) + t)
end

function VideoIO.read!(v::SegmentedVideo, img)
  if eof(get_current(v, :video)) && v.current < length(v.segments)
    v.current += 1
  end
  read!(get_current(v, :video), img)
end

function VideoIO.read(v::SegmentedVideo)
  if eof(get_current(v, :video)) && v.current < length(v.segments)
    v.current += 1
  end
  read(get_current(v, :video))
end

VideoIO.gettime(v::SegmentedVideo) = gettime(get_current(v, :video))


text:
@testset "Day7.jl" begin
    # Write your tests here.
    lines = read_infile(Day7_1, "$(@__DIR__)/assets/day7_test1.txt")
    @test solve(Day7_1, lines) == 4

    @test solve(Day7_2, lines) == 32

    lines = read_infile(Day7_1, "$(@__DIR__)/assets/day7_test2.txt")
    @test solve(Day7_2, lines) == 126

end


text:
"""
The R NAMED property, represented by 2 bits in the info field. This can take
values 0,1 or 2, corresponding to whether it is bound to 0,1 or 2 or more
symbols. See
http://cran.r-project.org/doc/manuals/r-patched/R-exts.html#Named-objects-and-copying
"""
function bound{S<:Sxp}(s::Ptr{S})
    u = unsafe_load(convert(UnknownSxpPtr,s))
    (u.info >>> 6) & 0x03
end

"""
Sxp methods for `length` return the R length.

`Rf_xlength` handles Sxps that are not vector-like and R's
"long vectors", which have a negative value for the `length` member.
"""
length{S<:Sxp}(s::Ptr{S}) =
    Int(ccall((:Rf_xlength,libR),Cptrdiff_t,(Ptr{S},),s))
length(r::RObject) = length(r.p)

## Predicates applied to an Sxp
##
## Many of these are unneeded but a few extra definitions is not a big deal
for sym in (:isArray,:isComplex,:isEnvironment,:isExpression,:isFactor,
            :isFrame,:isFree,:isFunction,:isInteger,:isLanguage,:isList,
            :isLogical,:isSymbol,:isMatrix,:isNewList,:isNull,:isNumeric,
            :isNumber,:isObject,:isOrdered,:isPairListSxp,:isPrimitiveSxp,
            :isReal,:isS4,:isString,:isTs,:isUnordered,:isUnsorted,
            :isUserBinop,:isValidString,:isValidStringF,:isVector,
            :isVectorAtomicSxp,:isVectorizable,:isVectorListSxp)
    @eval begin
        $sym{S<:Sxp}(s::Ptr{S}) = ccall(($(string("Rf_",sym)),libR),Bool,(Ptr{SxpPtrInfo},),s)
        $sym(r::RObject) = $sym(r.p)
    end
end

const voffset = Ref{UInt}()

"""
Pointer to start of the data array in a SEXPREC. Corresponds to DATAPTR C macro.
"""
dataptr{S<:VectorSxp}(s::Ptr{S}) = convert(Ptr{eltype(S)}, s+voffset[])

"""
Represent the contents of a VectorSxp type as a `Vector`.

This does __not__ copy the contents.  If the argument is not named (in R) or
otherwise protected from R's garbage collection (e.g. by keeping the
containing RObject in scope) the contents of this vector can be modified or
could cause a memory error when accessed.

The contents are as stored in R.  Missing values (NA's) are represented
in R by sentinels.  Missing data values in RealSxp and CplxSxp show
up as `NaN` and `NaN + NaNim`, respectively.  Missing data in IntSxp show up
as `-2147483648`, the minimum 32-bit integer value.  Internally a `LglSxp` is
represented as `Vector{Int32}`.  The convention is that `0` is `false`,
`-2147483648` is `NA` and all other values represent `true`.
"""
unsafe_vec{S<:VectorSxp}(s::Ptr{S}) = Compat.unsafe_wrap(Array, dataptr(s), length(s))
unsafe_vec{S<:VectorSxp}(r::RObject{S}) = unsafe_vec(r.p)

"""
The same as `unsafe_vec`, except returns an appropriately sized array.
"""
unsafe_array{S<:VectorSxp}(s::Ptr{S}) =  Compat.unsafe_wrap(Array, dataptr(s), size(s))
unsafe_array{S<:VectorSxp}(r::RObject{S}) = unsafe_array(r.p)

"""
Indexing into `VectorSxp` types uses Julia indexing into the `vec` result,
except for `StrSxp` and the `VectorListSxp` types, which must apply `sexp`
to the `Ptr{Void}` obtained by indexing into the `vec` result.
"""
getindex{S<:VectorAtomicSxp}(s::Ptr{S}, I::Real) = getindex(unsafe_vec(s),I)
getindex{S<:VectorAtomicSxp}(s::Ptr{S}, I::AbstractVector) = getindex(unsafe_vec(s),I)
getindex{S<:VectorAtomicSxp}(s::Ptr{S}, I::Real...) = getindex(unsafe_array(s),I...)

getindex{S<:VectorListSxp}(s::Ptr{S}, I::Real) = sexp(getindex(unsafe_vec(s),I))
getindex{S<:VectorListSxp}(s::Ptr{S}, I::AbstractVector) = sexp(getindex(unsafe_vec(s),I))
getindex{S<:VectorListSxp}(s::Ptr{S}, I::Real...) = sexp(getindex(unsafe_array(s),I...))

"""
String indexing finds the first element with the matching name
"""
function getindex{S<:VectorSxp}(s::Ptr{S}, label::AbstractString)
    ls = unsafe_vec(getnames(s))
    for (i,l) in enumerate(ls)
        if rcopy(l) == label
            return s[i]
        end
    end
    throw(BoundsError())
end
getindex{S<:VectorSxp}(s::Ptr{S}, label::Symbol) = getindex(s,string(label))

getindex{S<:VectorAtomicSxp}(r::RObject{S}, I...) = getindex(sexp(r), I...)
getindex{S<:VectorAtomicSxp}(r::RObject{S}, I::AbstractArray) = getindex(sexp(r), I)

getindex{S<:VectorListSxp}(r::RObject{S}, I...) = RObject(getindex(sexp(r), I...))
getindex{S<:VectorListSxp}(r::RObject{S}, I::AbstractArray) = map(RObject,getindex(sexp(r),I))

function setindex!{S<:VectorAtomicSxp}(s::Ptr{S}, value, I...)
    setindex!(unsafe_array(s), value, I...)
end
function setindex!{S<:VectorAtomicSxp}(s::Ptr{S}, value, I)
    setindex!(unsafe_vec(s), value, I)
end
function setindex!(s::Ptr{StrSxp}, value::CharSxpPtr, key::Integer)
    1 <= key <= length(s) || throw(BoundsError())
    ccall((:SET_STRING_ELT,libR), Void,
          (Ptr{StrSxp},Cptrdiff_t, CharSxpPtr),
          s, key-1, value)
    value
end
function setindex!(s::Ptr{StrSxp}, value::AbstractString, key::Integer)
    setindex!(s,sexp(CharSxp,value),key)
end

function setindex!{S<:Union{VecSxp,ExprSxp},T<:Sxp}(s::Ptr{S}, value::Ptr{T}, key::Integer)
    1 <= key <= length(s) || throw(BoundsError())
    ccall((:SET_VECTOR_ELT,libR), Ptr{T},
          (Ptr{S},Cptrdiff_t, Ptr{T}),
          s, key-1, value)
end
function setindex!{S<:Union{VecSxp,ExprSxp}}(s::Ptr{S}, value, key::Integer)
    setindex!(s,sexp(value),key)
end

setindex!(r::RObject, value, keys...) = setindex!(sexp(r), value, keys...)

start{S<:VectorSxp}(s::Ptr{S}) = 0
next{S<:VectorSxp}(s::Ptr{S},state) = (state += 1;(s[state],state))
done{S<:VectorSxp}(s::Ptr{S},state) = state ≥ length(s)

# PairListSxps

cdr{S<:PairListSxp}(s::Ptr{S}) = sexp(ccall((:CDR,libR),UnknownSxpPtr,(Ptr{S},),s))
car{S<:PairListSxp}(s::Ptr{S}) = sexp(ccall((:CAR,libR),UnknownSxpPtr,(Ptr{S},),s))
tag{S<:PairListSxp}(s::Ptr{S}) = sexp(ccall((:TAG,libR),UnknownSxpPtr,(Ptr{S},),s))

function setcar!{S<:PairListSxp,T<:Sxp}(s::Ptr{S}, c::Ptr{T})
    ccall((:SETCAR,libR),Ptr{Void},(Ptr{S},Ptr{T}),s,c)
    nothing
end
setcar!{S<:PairListSxp,T<:Sxp}(s::Ptr{S}, c::RObject{T}) = setcar!(s,sexp(c))

function settag!{S<:PairListSxp,T<:Sxp}(s::Ptr{S}, c::Ptr{T})
    ccall((:SET_TAG,libR),Void,(Ptr{S},Ptr{T}),s,c)
    nothing
end
settag!{S<:PairListSxp,T<:Sxp}(s::Ptr{S}, c::RObject{T}) = settag!(s,sexp(c))

function setcdr!{S<:PairListSxp,T<:Sxp}(s::Ptr{S}, c::Ptr{T})
    ccall((:SETCDR,libR),Ptr{Void},(Ptr{S},Ptr{T}),s,c)
    nothing
end
setcdr!{S<:PairListSxp,T<:Sxp}(s::Ptr{S}, c::RObject{T}) = setcdr!(s,sexp(c))

start{S<:PairListSxp}(s::Ptr{S}) = s
function next{S<:PairListSxp,T<:PairListSxp}(s::Ptr{S},state::Ptr{T})
    t = tag(state)
    c = car(state)
    (t,c), cdr(state)
end
done{S<:PairListSxp,T<:PairListSxp}(s::Ptr{S},state::Ptr{T}) = state == sexp(Const.NilValue)

"extract the i-th element of LangSxp l"
function getindex{S<:PairListSxp}(l::Ptr{S},I::Integer)
    1 ≤ I ≤ length(l) || throw(BoundsError())
    for i in 2:I
        l = cdr(l)
    end
    car(l)
end

getindex{S<:PairListSxp}(r::RObject{S},I::Integer) = RObject(getindex(sexp(r),I))

"assign value v to the i-th element of LangSxp l"
function setindex!{S<:PairListSxp,T<:Sxp}(l::Ptr{S},v::Ptr{T},I::Integer)
    1 ≤ I ≤ length(l) || throw(BoundsError())
    for i in 2:I
        l = cdr(l)
    end
    setcar!(l,v)
end
function setindex!{S<:PairListSxp}(s::Ptr{S}, value, key::Integer)
    setindex!(s,sexp(value),key)
end

"Return a particular attribute of an RObject"
function getattrib{S<:Sxp}(s::Ptr{S}, sym::Ptr{SymSxp})
    sexp(ccall((:Rf_getAttrib,libR),UnknownSxpPtr,(Ptr{S},Ptr{SymSxp}),s,sym))
end
getattrib{S<:Sxp}(s::Ptr{S}, sym::RObject{SymSxp}) = getattrib(s,sexp(sym))
getattrib{S<:Sxp}(s::Ptr{S}, sym::Symbol) = getattrib(s,sexp(SymSxp,sym))
getattrib{S<:Sxp}(s::Ptr{S}, sym::AbstractString) = getattrib(s,sexp(SymSxp,sym))

getattrib(r::RObject, sym) = RObject(getattrib(r.p,sym))

"Set a particular attribute of an RObject"
function setattrib!{S<:Sxp,T<:Sxp}(s::Ptr{S},sym::Ptr{SymSxp},t::Ptr{T})
    ccall((:Rf_setAttrib,libR),Ptr{Void},(Ptr{S},Ptr{SymSxp},Ptr{T}),s,sym,t)
    return nothing
end
setattrib!{S<:Sxp,T<:Sxp}(s::Ptr{S},sym::RObject{SymSxp},t::Ptr{T}) = setattrib!(s,sexp(sym),t)
setattrib!{S<:Sxp,T<:Sxp}(s::Ptr{S},sym::Symbol,t::Ptr{T}) = setattrib!(s,sexp(SymSxp,sym),t)
setattrib!{S<:Sxp,T<:Sxp}(s::Ptr{S},sym::AbstractString,t::Ptr{T}) = setattrib!(s,sexp(SymSxp,sym),t)
setattrib!{S<:Sxp}(s::Ptr{S},sym,t) = setattrib!(s,sym,sexp(t))

setattrib!(r::RObject, sym, t) = setattrib!(r.p, sym, t)

attributes(s::SxpHead) = sexp(s.attrib)
attributes(s::Sxp) = attributes(s.head)
attributes{S<:Sxp}(s::Ptr{S}) = attributes(unsafe_load(s))

function size{S<:Sxp}(s::Ptr{S})
    isArray(s) || return (length(s),)
    tuple(convert(Array{Int},unsafe_vec(getattrib(s,Const.DimSymbol)))...)
end
size(r::RObject) = size(sexp(r))

"""
Returns the names of an R vector.
"""
getnames{S<:Sxp}(s::Ptr{S}) = getattrib(s,Const.NamesSymbol)
getnames(r::RObject) = RObject(getnames(sexp(r)))

"""
Returns the names of an R vector, the result is converted to a Julia symbol array.
"""
names(r::RObject) = rcopy(Array{Symbol}, getnames(sexp(r)))

"""
Set the names of an R vector.
"""
setnames!{S<:Sxp}(s::Ptr{S}, n::Ptr{StrSxp}) = setattrib!(s,Const.NamesSymbol,n)
setnames!(r::RObject,n) = RObject(setnames!(sexp(r),sexp(StrSxp,n)))

"""
Returns the class of an R object.
"""
getclass{S<:Sxp}(s::Ptr{S}) = getattrib(s,Const.ClassSymbol)
getclass(r::RObject) = RObject(getclass(sexp(r)))

"""
Set the class of an R object.
"""
setclass!{S<:Sxp}(s::Ptr{S},c::Ptr{StrSxp}) = setattrib!(s,Const.ClassSymbol,c)
setclass!(r::RObject,c) = RObject(setclass!(sexp(r)),sexp(StrSxp,c))

allocList(n::Int) = ccall((:Rf_allocList,libR),Ptr{ListSxp},(Cint,),n)
allocArray{S<:Sxp}(::Type{S}, n::Integer) =
    ccall((:Rf_allocVector,libR),Ptr{S},(Cint,Cptrdiff_t),sexpnum(S),n)

allocArray{S<:Sxp}(::Type{S}, n1::Integer, n2::Integer) =
    ccall((:Rf_allocMatrix,libR),Ptr{S},(Cint,Cint,Cint),sexpnum(S),n1,n2)

allocArray{S<:Sxp}(::Type{S}, n1::Integer, n2::Integer, n3::Integer) =
    ccall((:Rf_alloc3DArray,libR),Ptr{S},(Cint,Cint,Cint,Cint),sexpnum(S),n1,n2,n3)

function allocArray{S<:Sxp}(::Type{S}, dims::Integer...)
    sdims = sexp(IntSxp,[dims...])
    ccall((:Rf_allocArray,libR),Ptr{S},(Cint,Ptr{IntSxp}),sexpnum(S),sdims)
end

"""
Check if values correspond to R's NULL object.
"""
isnull{S<:Sxp}(s::Ptr{S}) = isNull(s)
isnull{S<:Sxp}(r::RObject{S}) = isnull(r.p)

"""
NA element for each type
"""
naeltype(::Type{LglSxp}) = Const.NaInt
naeltype(::Type{IntSxp}) = Const.NaInt
naeltype(::Type{RealSxp}) = Const.NaReal
naeltype(::Type{CplxSxp}) = complex(Const.NaReal,Const.NaReal)
naeltype(::Type{StrSxp}) = sexp(Const.NaString)
naeltype(::Type{VecSxp}) = sexp(LglSxp,Const.NaInt) # used for setting

"""
Check if values correspond to R's sentinel NA values.
"""
isna(x::Complex128) = real(x) === Const.NaReal && imag(x) === Const.NaReal
isna(x::Float64) = x === Const.NaReal
isna(x::Int32) = x == Const.NaInt
isna(a::AbstractArray) = reshape(bitpack([isna(aa) for aa in a]),size(a))
isna(s::CharSxpPtr) = s === sexp(Const.NaString)

# this doesn't allow us to check VecSxp s
function isna{S<:VectorSxp}(s::Ptr{S})
    b = BitArray(size(s)...)
    for (i,e) in enumerate(s)
        b[i] = isna(e)
    end
    b
end

isna(r::RObject) = isna(r.p)

"""
Check if there are any NA values in the vector.
"""
function anyna{S<:VectorSxp}(s::Ptr{S})
    for i in s
        if isna(i)
            return true
        end
    end
    return false
end
anyna{S<:VectorSxp}(r::RObject{S}) = anyna(r.p)

# StrSxp
"""
Determines the encoding of the CharSxp. This is determined by the 'gp' part of the sxpinfo (this is the middle 16 bits).
 * 0x00_0002_00 (bit 1): set of bytes (no known encoding)
 * 0x00_0004_00 (bit 2): Latin-1
 * 0x00_0008_00 (bit 3): UTF-8
 * 0x00_0040_00 (bit 6): ASCII

We only support ASCII and UTF-8.
"""
function isascii(s::CharSxp)
    if s.head.info & 0x00_0040_00 != 0
        return true
    elseif s.head.info & 0x00_0008_00 != 0
        return false
    else
        error("Unsupported string type.")
    end
end
isascii(s::CharSxpPtr) = isascii(unsafe_load(s))
isascii(r::RObject{CharSxp}) = isascii(sexp(r))

function isascii(s::StrSxpPtr)
    ind = true
    for c in s
        ind &= isna(c) || isascii(c)
    end
    return ind
end
isascii(r::RObject{StrSxp}) = isascii(sexp(r))

# EnvSxp

"extract the value of symbol s in the environment e"
function getindex(e::Ptr{EnvSxp},s::Ptr{SymSxp})
    v = ccall((:Rf_findVarInFrame,libR),UnknownSxpPtr,(Ptr{EnvSxp},Ptr{SymSxp}),e,s)
    v == sexp(Const.UnboundValue) && error("$s is not defined in the environment")
    sexp(v)
end
getindex(e::Ptr{EnvSxp},s) = getindex(e,sexp(SymSxp,s))
getindex(e::RObject{EnvSxp},s) = RObject(getindex(sexp(e),s))

"assign value v to symbol s in the environment e"
function setindex!{S<:Sxp}(e::Ptr{EnvSxp},v::Ptr{S},s::Ptr{SymSxp})
    # This should be done more carefully.  First check for the symbol in the
    # frame.  If it is defined call Rf_setVar, otherwise call Rf_defineVar.
    # As it stands this segfaults if the symbol is bound in, say, the base
    # environment.
    ccall((:Rf_defineVar,libR),Void,(Ptr{SymSxp},Ptr{S},Ptr{EnvSxp}),s,v,e)
end
function setindex!(e::Ptr{EnvSxp},v,s)
    nprotect = 0
    try
        sv = protect(sexp(v))
        nprotect += 1
        ss = protect(sexp(SymSxp,s))
        nprotect += 1
        setindex!(e,sv,ss)
    finally
        unprotect(nprotect)
    end
end
setindex!(e::RObject{EnvSxp},v,s) = setindex!(sexp(e),v,s)

"""
    newEnvironment([env])

Create a new environment which extends environment `env` (`globalEnv` by default).
"""
function newEnvironment(env::Ptr{EnvSxp})
    ccall((:Rf_NewEnvironment,libR),Ptr{EnvSxp},
            (Ptr{NilSxp},Ptr{NilSxp},Ptr{EnvSxp}),sexp(Const.NilValue),sexp(Const.Const.NilValue),env)
end
newEnvironment(env::RObject{EnvSxp}) = newEnvironment(sexp(env))
newEnvironment() = newEnvironment(globalEnv)

"find namespace by name of the namespace, it is not error tolerant."
function findNamespace(str::Compat.String)
    ccall((:R_FindNamespace,libR),Ptr{EnvSxp}, (Ptr{StrSxp},), sexp(str))
end

"get namespace by name of the namespace. It is safer to be used than findNamespace as it checks bound."
getNamespace(str::Compat.String) = reval(rlang_p(RCall.Const.BaseNamespace["getNamespace"], str))


text:
@testset "Dynamics" begin
    T = 3
    num_state = 2
    num_action = 1
    num_parameter = 0
    parameter_dimensions = [num_parameter for t = 1:T]

    function pendulum(z, u, w)
        mass = 1.0
        lc = 1.0
        gravity = 9.81
        damping = 0.1
        [z[2], (u[1] / ((mass * lc * lc)) - gravity * sin(z[1]) / lc - damping * z[2] / (mass * lc * lc))]
    end

    function euler_implicit(y, x, u, w)
        h = 0.1
        y - (x + h * pendulum(y, u, w))
    end

    dt = Dynamics(euler_implicit, num_state, num_state, num_action,
        num_parameter=num_parameter);
    dynamics = [dt for t = 1:T-1]

    x1 = ones(num_state)
    u1 = ones(num_action)
    w1 = zeros(num_parameter)
    X = [x1 for t = 1:T]
    U = [u1 for t = 1:T]
    W = [w1 for t = 1:T]
    idx_dyn = DTO.constraint_indices(dynamics)
    idx_jac = DTO.jacobian_indices(dynamics)

    d = zeros(DTO.num_constraint(dynamics))
    j = zeros(DTO.num_jacobian(dynamics))

    dt.evaluate(dt.evaluate_cache, x1, x1, u1, w1)
    # @benchmark $dt.evaluate($dt.evaluate_cache, $x1, $x1, $u1, $w1)
    @test norm(dt.evaluate_cache - euler_implicit(x1, x1, u1, w1)) < 1.0e-8
    dt.jacobian(dt.jacobian_cache, x1, x1, u1, w1)
    jac_dense = zeros(dt.num_next_state, dt.num_state + dt.num_action + dt.num_next_state)
    for (i, ji) in enumerate(dt.jacobian_cache)
        jac_dense[dt.jacobian_sparsity[1][i], dt.jacobian_sparsity[2][i]] = ji
    end
    jac_fd = ForwardDiff.jacobian(a -> euler_implicit(a[num_state + num_action .+ (1:num_state)], a[1:num_state], a[num_state .+ (1:num_action)], w1), [x1; u1; x1])
    @test norm(jac_dense - jac_fd) < 1.0e-8

    DTO.constraints!(d, idx_dyn, dynamics, X, U, W)
    @test norm(vcat(d...) - vcat([euler_implicit(X[t+1], X[t], U[t], W[t]) for t = 1:T-1]...)) < 1.0e-8
    # info = @benchmark DTO.constraints!($d, $idx_dyn, $dynamics, $X, $U, $W)

    DTO.jacobian!(j, idx_jac, dynamics, X, U, W)
    s = DTO.sparsity_jacobian(dynamics, DTO.dimensions(dynamics)[1:2]...)
    jac_dense = zeros(DTO.num_constraint(dynamics), DTO.num_state_action_next_state(dynamics))
    for (i, ji) in enumerate(j)
        jac_dense[s[i][1], s[i][2]] = ji
    end

    @test norm(jac_dense - [jac_fd zeros(dynamics[2].num_state, dynamics[2].num_action + dynamics[2].num_next_state); zeros(dynamics[2].num_next_state, dynamics[1].num_state + dynamics[1].num_action) jac_fd]) < 1.0e-8
    # info = @benchmark DTO.jacobian!($j, $idx_jac, $dynamics, $X, $U, $W)

    x_idx = DTO.state_indices(dynamics)
    u_idx = DTO.action_indices(dynamics)
    xu_idx = DTO.state_action_indices(dynamics)
    xuy_idx = DTO.state_action_next_state_indices(dynamics)

    nz = sum([t < T ? dynamics[t].num_state : dynamics[t-1].num_next_state for t = 1:T]) + sum([dynamics[t].num_action for t = 1:T-1])
    z = rand(nz)
    x = [zero(z[x_idx[t]]) for t = 1:T]
    u = [[zero(z[u_idx[t]]) for t = 1:T-1]..., zeros(0)]

    DTO.trajectory!(x, u, z, x_idx, u_idx)
    z̄ = zero(z)
    for (t, idx) in enumerate(x_idx)
        z̄[idx] .= x[t]
    end
    for (t, idx) in enumerate(u_idx)
        z̄[idx] .= u[t]
    end

    @test norm(z - z̄) < 1.0e-8
    # info = @benchmark DTO.trajectory!($x, $u, $z, $x_idx, $u_idx)
end


text:
import Base: convert, ndims, size

# --------------
# Abstract types
# --------------

@doc """
# Interface

* probability(P, s, t, a) -> the probability of transitioning from state s to state t given action a

""" ->
abstract AbstractTransitionProbability

@doc """
# Interface

* probability(P, a) -> probability transition matrix for action a

## Optional

* getindex(P, s, t, a) -> the probability of transitioning from state s to state t given action a

""" ->
abstract AbstractArrayTransitionProbability <: AbstractTransitionProbability

# ----------
# Array type
# ----------

immutable ArrayTransitionProbability{T<:Real} <: AbstractArrayTransitionProbability
    array::Array{T,3}

    function ArrayTransitionProbability(array::Array{T,3})
        is_square_stochastic(array) || error("Not valid transition probability.")
        new(array)
    end
end

ArrayTransitionProbability{T}(array::Array{T,3}) = ArrayTransitionProbability{T}(array)

TransitionProbability{T}(A::Array{T,3}) = ArrayTransitionProbability(A)

# Base
# ----

ndims(P::ArrayTransitionProbability) = ndims(P.array)

size(P::ArrayTransitionProbability) = size(P.array)
size(P::ArrayTransitionProbability, dims...) = size(P.array, dims...)

getindex(P::ArrayTransitionProbability, dims...) = getindex(P.array, dims...)

# num
# ---

num_actions(P::ArrayTransitionProbability) = size(P, 3)

num_states(P::ArrayTransitionProbability) = size(P, 1)

# probability
# -----------

@doc """
# Transition probability

The probability of transitioning from state `s` to
state `t` given that action `a` is taken.

""" ->
probability(P::ArrayTransitionProbability, s, t, a) = getindex(P, s, t, a)

probability(P::ArrayTransitionProbability, a) = P[:, :, a]

# -----------------
# Sparse array type
# -----------------

immutable SparseArrayTransitionProbability{Tv<:Real,Ti} <: AbstractArrayTransitionProbability
    array::Vector{SparseMatrixCSC{Tv,Ti}}
end

TransitionProbability{T<:SparseMatrixCSC}(A::Vector{T}) =
    SparseArrayTransitionProbability(A)

num_actions(P::SparseArrayTransitionProbability) = length(P.array)

num_states(P::SparseArrayTransitionProbability) = size(getindex(P.array, 1), 1)

probability(P::SparseArrayTransitionProbability, s, t, a) =
    getindex(getindex(P.array, a), s, t)

probability(P::SparseArrayTransitionProbability, a) = getindex(P.array, a)

# -------------
# Function type
# -------------

immutable FunctionTransitionProbability <: AbstractTransitionProbability
    func::Function
    states::Int
    actions::Int
end

TransitionProbability(F::Function, states, actions) = FunctionTransitionProbability(F, states, actions)

num_actions(P::FunctionTransitionProbability) = P.actions

num_states(P::FunctionTransitionProbability) = P.states

probability(P::FunctionTransitionProbability, s, t, a) = P.func(s, t, a)


text:
using
  CUDA,
  FourierFlows,
  Statistics,
  Random,
  FFTW,
  Test

import # use 'import' rather than 'using' for submodules to keep namespace clean
  GeophysicalFlows.TwoDNavierStokes,
  GeophysicalFlows.BarotropicQG,
  GeophysicalFlows.BarotropicQGQL,
  GeophysicalFlows.MultilayerQG,
  GeophysicalFlows.SurfaceQG

using FourierFlows: parsevalsum
using GeophysicalFlows: lambdipole, peakedisotropicspectrum

# the devices on which tests will run
devices = (CPU(),)
@has_cuda devices = (CPU(), GPU())

const rtol_lambdipole = 1e-2 # tolerance for lamb dipole tests
const rtol_twodnavierstokes = 1e-13 # tolerance for twodnavierstokes forcing tests
const rtol_barotropicQG = 1e-13 # tolerance for barotropicqg forcing tests
const rtol_multilayerqg = 1e-13 # tolerance for multilayerqg forcing tests
const rtol_surfaceqg = 1e-13 # tolerance for multilayerqg forcing tests

"Get the CFL number, assuming a uniform grid with `dx=dy`."
cfl(u, v, dt, dx) = maximum([maximum(abs.(u)), maximum(abs.(v))]*dt/dx)
cfl(prob) = cfl(prob.vars.u, prob.vars.v, prob.clock.dt, prob.grid.dx)

# Run tests
testtime = @elapsed begin
for dev in devices

  println("testing on "*string(typeof(dev)))

  @testset "Utils" begin
    include("test_utils.jl")

    @test testpeakedisotropicspectrum(dev)
    @test_throws ErrorException("the domain is not square") testpeakedisotropicspectrum_rectangledomain()
  end

  @testset "TwoDNavierStokes" begin
    include("test_twodnavierstokes.jl")

    @test test_twodnavierstokes_advection(0.0005, "ForwardEuler", dev)
    @test test_twodnavierstokes_lambdipole(256, 1e-3, dev)
    @test test_twodnavierstokes_deterministicforcing_energybudget(dev)
    @test test_twodnavierstokes_stochasticforcing_energybudget(dev)
    @test test_twodnavierstokes_deterministicforcing_enstrophybudget(dev)
    @test test_twodnavierstokes_stochasticforcing_enstrophybudget(dev)
    @test test_twodnavierstokes_energyenstrophy(dev)
    @test test_twodnavierstokes_problemtype(dev, Float32)
    @test TwoDNavierStokes.nothingfunction() == nothing
  end

  @testset "BarotropicQG" begin
    include("test_barotropicqg.jl")

    @test test_bqg_rossbywave("ETDRK4", 1e-2, 20, dev)
    @test test_bqg_rossbywave("FilteredETDRK4", 1e-2, 20, dev)
    @test test_bqg_rossbywave("RK4", 1e-2, 20, dev)
    @test test_bqg_rossbywave("FilteredRK4", 1e-2, 20, dev)
    @test test_bqg_rossbywave("AB3", 1e-3, 200, dev)
    @test test_bqg_rossbywave("FilteredAB3", 1e-3, 200, dev)
    @test test_bqg_rossbywave("ForwardEuler", 1e-4, 2000, dev)
    @test test_bqg_rossbywave("FilteredForwardEuler", 1e-4, 2000, dev)
    @test test_bqg_stochasticforcingbudgets(dev)
    @test test_bqg_deterministicforcingbudgets(dev)
    @test test_bqg_advection(0.0005, "ForwardEuler", dev)
    @test test_bqg_formstress(0.01, "ForwardEuler", dev)
    @test test_bqg_energyenstrophy(dev)
    @test test_bqg_meanenergyenstrophy(dev)
    @test test_bqg_problemtype(dev, Float32)
    @test BarotropicQG.nothingfunction() == nothing
  end

  @testset "BarotropicQGQL" begin
    include("test_barotropicqgql.jl")

    @test test_bqgql_rossbywave("ETDRK4", 1e-2, 20, dev)
    @test test_bqgql_rossbywave("FilteredETDRK4", 1e-2, 20, dev)
    @test test_bqgql_rossbywave("RK4", 1e-2, 20, dev)
    @test test_bqgql_rossbywave("FilteredRK4", 1e-2, 20, dev)
    @test test_bqgql_rossbywave("AB3", 1e-3, 200, dev)
    @test test_bqgql_rossbywave("FilteredAB3", 1e-3, 200, dev)
    @test test_bqgql_rossbywave("ForwardEuler", 1e-4, 2000, dev)
    @test test_bqgql_rossbywave("FilteredForwardEuler", 1e-4, 2000, dev)
    @test test_bqgql_deterministicforcingbudgets(dev)
    @test test_bqgql_stochasticforcingbudgets(dev)
    @test test_bqgql_advection(0.0005, "ForwardEuler", dev)
    @test test_bqgql_energyenstrophy(dev)
    @test test_bqgql_problemtype(dev, Float32)
    @test BarotropicQGQL.nothingfunction() == nothing
  end

  @testset "MultilayerQG" begin
    include("test_multilayerqg.jl")

    @test test_pvtofromstreamfunction_2layer(dev)
    @test test_pvtofromstreamfunction_3layer(dev)
    @test test_mqg_rossbywave("RK4", 1e-2, 20, dev)
    @test test_mqg_nonlinearadvection(0.005, "ForwardEuler", dev)
    @test test_mqg_linearadvection(0.005, "ForwardEuler", dev)
    @test test_mqg_energies(dev)
    @test test_mqg_energysinglelayer(dev)
    @test test_mqg_fluxes(dev)
    @test test_mqg_fluxessinglelayer(dev)
    @test test_mqg_setqsetψ(dev)
    @test test_mqg_paramsconstructor(dev)
    @test test_mqg_stochasticforcedproblemconstructor(dev)
    @test test_mqg_problemtype(dev, Float32)
    @test MultilayerQG.nothingfunction() == nothing
  end

  @testset "SurfaceQG" begin
    include("test_surfaceqg.jl")

    @test test_sqg_kineticenergy_buoyancyvariance(dev)
    @test test_sqg_advection(0.0005, "ForwardEuler", dev)
    @test test_sqg_deterministicforcing_buoyancy_variance_budget(dev)
    @test test_sqg_stochasticforcing_buoyancy_variance_budget(dev)
    @test test_sqg_stochasticforcedproblemconstructor(dev)
    @test test_sqg_problemtype(dev, Float32)
    @test test_sqg_paramsconstructor(dev)
    @test test_sqg_noforcing(dev)
    @test SurfaceQG.nothingfunction() == nothing
  end

end

end # time

println("Total test time: ", testtime)


text:
# export predict

function apply_tree_by_leaf(idxOrig::Vector{Int}, leaves_of_tree, t::Tree, features::DataFrame)
	return apply_tree_by_leaf(idxOrig, leaves_of_tree, t.rootnode, features)
end

################################################################################################################
# these functions are used internally (e.g. by the boosting algorithm)
################################################################################################################

"""
This is the core apply_tree_by_leaf! function which is used for boosting (and other models)
    apply_tree_by_leaf!(fit::Vector{Float64},leaf::Vector{Int},idx::Vector{Int},t::Union{Leaf,Node{UInt8},Node{UInt16}},features::DataFrame)
"""
function apply_tree_by_leaf!(fit::Vector{Float64}, leaf::Vector{Int}, idx::Vector{Int}, t::Union{Leaf,Node{UInt8},Node{UInt16}}, features::DataFrame)
    nobs = size(features, 1)
    @assert nobs == length(fit) == length(leaf)
    apply_tree_by_leaf_iteration!(idx, t, features, fit, leaf)
    return nothing
end

function apply_tree_by_leaf(idx::Vector{Int}, t::Union{Leaf,Node{UInt8},Node{UInt16}}, features::DataFrame)
  # as we always supply an index here, the function will always return two arrays of length size(features,2)
  # thus the "training" observations will be empty if this function is called on a validation index only!!
    nobs = size(features, 1)
    fit = zeros(Float64, nobs)
    leaf = zeros(Int, nobs)
    apply_tree_by_leaf_iteration!(idx, t, features, fit, leaf)
    return fit, leaf
end

function apply_tree_by_leaf_iteration!(idx::Vector{Int}, t::Node{T}, features::DataFrame, fit::Vector{Float64}, leaf::Vector{Int}) where T <: Unsigned
    if length(idx) == 0
       # @info "DTM: No data (an empty index) was provided to apply tree function"
        return nothing
    end
    subset = t.subset
    idxlInteger, idxrInteger = lrIndices(idx, features[!,t.featid_new_positive], subset)
    apply_tree_by_leaf_iteration!(idxlInteger, t.left, features, fit, leaf)
    apply_tree_by_leaf_iteration!(idxrInteger, t.right, features, fit, leaf)
    return nothing
end

function apply_tree_by_leaf_iteration!(idx, t::Leaf, features::DataFrame, fit, leaf)
	  @inbounds for i in idx
        leaf[i] = t.id
        fit[i] = t.fitted
    end
    return nothing
end

################################################################################################################
# functions for predicting estimates on unseen data
################################################################################################################

"""
Apply tree by leaf function for 'new' data
"""
function apply_tree_by_leaf(t::Union{Leaf,Node{UInt8},Node{UInt16}}, features::DataFrame)
  # fit::Vector{Float64},leaf::Vector{Int},
    nobs = size(features, 1)
    fit = zeros(Float64, nobs)
    leafnrs = zeros(Int, nobs)
    idx = collect(1:nobs)
    apply_tree_by_leaf_iteration!(idx, t, features, fit, leafnrs)
    return fit, leafnrs
end

function predictConsistentFeatures(x::Tree, f::DataFrame)
    @assert assert_consistent_features(x.featurepools, f)
    return apply_tree_by_leaf(x.rootnode, f)
end

function predictConsistentFeatures(x::BoostedTree, f::DataFrame;prroduceEstAndLeafMatrices::Bool=false)
    @assert assert_consistent_features(x.featurepools, f)

    obs = size(f, 1)
    sett = x.settings
    trn_meanobservedvalue = x.meanobserved
    iterations = sett.iterations

    estimatedRatio = ones(obs) .* trn_meanobservedvalue
    currentRelativity = ones(obs)
    indicatedRelativityForApplyTree_reused = zeros(obs)
    reused_fitted_leafnr_vector = zeros(Int, obs)
    est_SmoothedEstFromScores = zeros(obs)
    est_UnSmoothedEstFromScores = zeros(obs)

    thisidx = collect(1:obs)
    empty_validx = ones(Int, 0)

    if prroduceEstAndLeafMatrices
		est_matrix = Array{Float64}(undef, obs, iterations + 1)
		est_matrixFromScores = copy(est_matrix)
		MatrixOfLeafNumbers = Array{Int}(undef, obs, iterations + 1)
		MatrixOfLeafNumbers[:,1] = 0
		est_matrix[:,1] = copy(transpose(estimatedRatio))
		est_matrixFromScores[:,1] = copy(transpose(estimatedRatio))
	else
		est_matrix = Array{Float64}(undef, 0, 0)
		est_matrixFromScores = copy(est_matrix)
		MatrixOfLeafNumbers = Array{Int}(undef, 0, 0)
	end

    for iter = 1:iterations
        current_mdf = x.moderationvector[iter]

    # estimatedNumerator=estimatedRatio.*denominator #these are for instance the estimated losses for a LR model
    # vectorOfLeafArrays[iter+1]=create_leaves_array(res[iter])

        apply_tree_by_leaf!(indicatedRelativityForApplyTree_reused, reused_fitted_leafnr_vector, thisidx, x.trees[iter], f)
    # Moderate the estimate
        _moderate!(estimatedRatio, indicatedRelativityForApplyTree_reused, current_mdf)
        update_current_rels!(currentRelativity, estimatedRatio, trn_meanobservedvalue)

        if prroduceEstAndLeafMatrices
            write_column!(est_matrix, iter + 1, estimatedRatio)
            MatrixOfLeafNumbers[:,iter + 1] = reused_fitted_leafnr_vector
        end
    end

  # derive scores from raw rels
    scores = map(z->searchsortedfirst(x.maxRawRelativityPerScoreSorted, z), currentRelativity)
    max_score = length(x.rawObservedRatioPerScore)
    @inbounds for zz in 1:length(scores)
    # @inbounds scores[zz]=min(scores[zz],max_score)
        if scores[zz] > max_score
            scores[zz] = max_score
        end
    end
  # derive estimates
    update_mapped_estFromScores!(est_UnSmoothedEstFromScores, x.rawObservedRatioPerScore, scores)
    update_mapped_estFromScores!(est_SmoothedEstFromScores, x.ScoreToSmoothedEstimate, scores)

  # NOTE: we have three different estimates
    est_RawRelativities = currentRelativity .* trn_meanobservedvalue
    est_UnSmoothedEstFromScores
    est_SmoothedEstFromScores

    result = DataFrame()
    result[!,:scores] = scores
    result[!,:SmoothedEstimate] = est_SmoothedEstFromScores
    result[!,:UnsmoothedEstimate] = est_UnSmoothedEstFromScores
    result[!,:RawEstimate] = est_RawRelativities

    return result
end

function assert_consistent_features(fp, f::DataFrame)
    @assert length(fp) == size(f, 2)
    # the above condition could also be weakend (as could the loop below) we only need to check the variables that were actually used by the model.
    for j = 1:length(fp)
        if !(fp[j] == f[!,j].pool)
            # Note: this condition could be weakend: the features of the data which is provided only needs to be a subset of the pools used during modelling
            @warn("DTM: Features do not match model features!\r\nAlgorithm may abort.")
            return false
        end
    end
    return true
end


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

JLLWrappers.@generate_main_file_header("Qt5Base")
JLLWrappers.@generate_main_file("Qt5Base", UUID("ea2cea3b-5b76-57ae-a6ef-0a8af62496e1"))
end  # module Qt5Base_jll


text:
# lic checking file
# -----------------

# license check empty function
# ----------------------------
function emptyliccheck(lic::Vector{Cint})
    return lic
end

function touchlic(path)
    f = open(path)
    close(f)
end

function get_xpauthpath(xpauth_path = "", verbose::Bool = true)
    XPAUTH = "xpauth.xpr"

    candidates = []

    # user sent the complete path
    push!(candidates, xpauth_path)

    # user sent directory
    push!(candidates, joinpath(xpauth_path, XPAUTH))

    # default env (not metioned in manual)
    if haskey(ENV, "XPAUTH_PATH")
        xpauth_path = replace(ENV["XPAUTH_PATH"], "\"" => "")
        push!(candidates, joinpath(xpauth_path, XPAUTH))
    end

    # default lib dir
    if haskey(ENV, "XPRESSDIR")
        xpressdir = replace(ENV["XPRESSDIR"], "\"" => "")
        push!(candidates, joinpath(xpressdir, "bin", XPAUTH))
    end

    # user´s lib dir
    push!(candidates, joinpath(dirname(dirname(libxprs)), "bin", XPAUTH))

    for i in candidates
        if isfile(i)
            if verbose && !haskey(ENV, "XPRESS_JL_NO_INFO")
                @info("Xpress: Found license file $i")
            end
            return i
        end
    end

    error("Could not find xpauth.xpr license file. Check XPRESSDIR or XPAUTH_PATH environment variables.")
end
"""
    userlic(; liccheck::Function = emptyliccheck, xpauth_path::String = "" )
Performs license chhecking with `liccheck` validation function on dir `xpauth_path`
"""
function userlic(; verbose::Bool = true, liccheck::Function = emptyliccheck, xpauth_path::String = "" )

    # change directory to reach all libs
    initdir = pwd()
    if isdir(dirname(libxprs))
        cd(dirname(libxprs))
    end

    # open and free xpauth.xpr (touches the file to release it)
    path_lic = get_xpauthpath(xpauth_path, verbose)
    touchlic(path_lic)

    # pre allocate vars
    lic = Cint[1]
    slicmsg =  path_lic #xpauth_path == "dh" ? Array{Cchar}(undef, 1024*8) :

    # FIRST call do xprslicense to get BASE LIC
    license(lic, slicmsg)

    # convert BASE LIC to GIVEN LIC
    lic = liccheck(lic)

    # Send GIVEN LIC to XPRESS lib
    buffer = Array{Cchar}(undef, 1024*8)
    buffer_p = pointer(buffer)
    ierr = GC.@preserve buffer begin
        license(lic, Cstring(buffer_p))
    end

    # check LIC TYPE
    if ierr == 16
        # DEVELOPER
        if verbose && !haskey(ENV, "XPRESS_JL_NO_INFO")
            @info("Xpress: Development license detected.")
        end
    elseif ierr != 0
        # FAIL
        @info("Xpress: Failed to find working license.")
        error(getlicerrmsg())
    else
        # USER
        if verbose && !haskey(ENV, "XPRESS_JL_NO_INFO")
            @info("Xpress: User license detected.")
            @info(  unsafe_string(pointer(slicmsg))  )
        end
    end

    # go back to initial folder
    cd(initdir)
end


text:
##########################
# FunctionCell Interface #
##########################
"""
	mlkernel([σ,] κ; center=false)

Constructs an untrained function cell using the memory layout `σ` and kernel `κ`.

# Arguments
  * `σ::MLKernels.PairwiseFunctions.MemoryLayout` specifies whether to work with row or column major data
  (default `MLKernels.ColumnMajor() i.e. columns are observations)
  * `κ::MLKernels.Kernel` is the kernel object

# Keyword arguments
  * `center::Bool=false` specifies whether to center the kernel or not

A list of common used kernels in `MLKernels.jl` can be found below (definitions from the documentation of to MLKernels.jl):

# Kernels
  * `(x,y)-> a⋅xᵀy + c		MLKernels.LinearKernel(a>0, c≥0)`
  * `(x,y)-> (a⋅xᵀy+c)ᵈ		MLKernels.PolynomialKernel(a>0, c≥0, d∈Ζ₊)`
  * `(x,y)-> tanh(a⋅xᵀy+c)		MLKernels.SigmoidKernel(a>0, c≥0)`
  * `(x,y)-> exp(a⋅xᵀy)		MLKernels.ExponentiatedKernel(a>0)`
  * `(x,y)-> exp(-α*‖x-y‖)		MLKernels.ExponentialKernel(α>0)`
  * `(x,y)-> exp(-α⋅‖x-y‖²)		MLKernels.SquaredExponentialKernel(α>0)`
  * `(x,y)-> exp(-α⋅‖x-y‖ᵞ)		MLKernels.GammaExponentialKernel(α>0, 0<γ<1)`
  * `(x,y)-> (1 + α⋅‖x-y‖²)⁻ᵝ	MLKernels.RationalQuadraticKernel(α>0, β>0)`
  * `(x,y)-> (1 + α⋅‖x-y‖ᵞ)⁻ᵝ	MLKernels.GammaRationalKernel(α>0, β>0, 0<γ,1)`
  * `(x,y)-> exp{-α⋅Σᵢsin²[p⋅(xᵢ-yᵢ)]}	MLKernels.PeriodicKernel(p>0, α>0)`
  * `(x,y)-> ‖x-y‖²ᵞ		MLKernels.PowerKernel(0<γ≤1) (negative definite)`
  * `(x,y)-> log(1+α⋅‖x-y‖²ᵞ)	MLKernels.LogKernel(α>0, 0<γ≤1) (negative definite)`

Read the `MLKernels.jl` documentation for more information.

# Examples
```
julia> K = mlkernel(MLKernels.LinearKernel(1,2), center=false)
ML Kernel σ=MLKernels.PairwiseFunctions.ColumnMajor(), κ=LinearKernel(1.0,2.0), center=false, varying I/O dimensions, untrained

julia> x = [1,2,3]; y=[4,5,6,7];

julia> y |> K(x)
3×4 Array{Float64,2}:
  6.0   7.0   8.0   9.0
 10.0  12.0  14.0  16.0
 14.0  17.0  20.0  23.0
``
"""
mlkernel(κ::MLKernels.Kernel; center::Bool=false) =
	mlkernel(MLKernels.PairwiseFunctions.ColumnMajor(), κ; center=center)

mlkernel(σ::MLKernels.PairwiseFunctions.MemoryLayout, κ::MLKernels.Kernel; center::Bool=false) =
	FunctionCell(mlkernel, (σ,κ), ModelProperties(), "ML Kernel σ=$(σ), κ=$(κ), center=$center"; center=center)

############################
# DataCell/Array Interface #
############################
"""
	mlkernel(x, [σ,] κ; center=false)

Trains the function cell using the memory layout `σ`, kernel `κ` and data `x` for future kernel calculations.
"""
# Training
mlkernel(x::T where T<:CellData, κ::MLKernels.Kernel; center::Bool=false) =
	mlkernel(x, MLKernels.PairwiseFunctions.ColumnMajor(), κ; center=center)

mlkernel(x::T where T<:CellData, σ::MLKernels.PairwiseFunctions.MemoryLayout, κ::MLKernels.Kernel; center::Bool=false) =
	mlkernel(getx!(x), MLKernels.PairwiseFunctions.ColumnMajor(), κ; center=center)

mlkernel(x::T where T<:AbstractArray, κ::MLKernels.Kernel; center::Bool=false) =
	mlkernel(x, MLKernels.PairwiseFunctions.ColumnMajor(), κ; center=center)

mlkernel(x::T where T<:AbstractVector, σ::MLKernels.PairwiseFunctions.MemoryLayout, κ::MLKernels.Kernel; center::Bool=false) =
	mlkernel(mat(x, LearnBase.ObsDim.Constant{2}()), σ, κ; center=center)

mlkernel(x::T where T<:AbstractMatrix, σ::MLKernels.PairwiseFunctions.MemoryLayout, κ::MLKernels.Kernel; center::Bool=false) =
begin
	# Build model properties
	modelprops = ModelProperties(nvars(x),nobs(x))

	# Returned trained cell
	FunctionCell(mlkernel, Model((getobs(x), σ, κ, center), modelprops), "ML Kernel σ=$(σ), κ=$(κ), center=$center")
end

# Execution
mlkernel(x::T where T<:CellData, model::Model{<:Tuple{<:AbstractMatrix,<:MLKernels.PairwiseFunctions.MemoryLayout, <:MLKernels.Kernel, Bool}}) =
	datacell(mlkernel(getx!(x), model), gety(x))
mlkernel(x::T where T<:AbstractVector, model::Model{<:Tuple{<:AbstractMatrix,<:MLKernels.PairwiseFunctions.MemoryLayout, <:MLKernels.Kernel, Bool}}) =
	mlkernel(mat(x, LearnBase.ObsDim.Constant{2}()), model)
mlkernel(x::T where T<:AbstractMatrix, model::Model{<:Tuple{<:AbstractMatrix, <:MLKernels.PairwiseFunctions.MemoryLayout, <:MLKernels.Kernel, Bool}}) =
begin
	out = MLKernels.kernelmatrix(model.data[2], model.data[3], model.data[1], getobs(x))
	if model.data[4]
		MLKernels.centerkernelmatrix!(out)
	end
	return out
end


text:
#=
  Copyright (c) 2021
  https://github.com/Ismael-VC/Varvara.jl/blob/main/CONTRIBUTORS.md

  Permission to use, copy, modify, and distribute this software for any
  purpose with or without fee is hereby granted, provided that the above
  copyright notice and this permission notice appear in all copies.

  THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
  WITH REGARD TO THIS SOFTWARE.
=#

module VarvaraUtils

export low, high, concat, bool

low(x::UInt16)::UInt8 = x & 0xff
high(x::UInt16)::UInt8 = x >> 0x8
concat(low::UInt8, high::UInt8)::UInt16 = UInt16(high) << 8 | low
bool(x::Number)::Bool = !iszero(x)
bool(x::Union{AbstractArray,AbstractString}) = !isempty(x)

end  # module


text:
using DSP

function main()

	low_freqs = [697, 770, 852, 941]
	high_freqs = [1209, 1336, 1477, 1633]

	periods = Matrix{Float64}(undef, length(low_freqs), length(high_freqs))

	println("Fundametal periods for frequency combinations")

	for i in range(1, length(low_freqs))
		for j in range(1, length(high_freqs))
			periods[i,j] = 1/gcd(low_freqs[i], high_freqs[j])
			println("($(low_freqs[i]) Hz, $(high_freqs[j]) Hz) = $(periods[i,j]) s - $(1/periods[i,j]) Hz")
		end
	end
end

if PROGRAM_FILE == "FundametalPeriods.jl"
	main()
end

max_stars_repo_name: PatricioIribarneCatella/TP-Telefono

max_issues_repo_name: PatricioIribarneCatella/TP-Telefono

max_forks_repo_name: PatricioIribarneCatella/TP-Telefono


text:
using MotionPlanning

ISRR_2H{T<:AbstractFloat}(::Type{T} = Float64) = Compound2D(
    Shape2D{T}[
        Box2D([T(.0), T(.16)], [T(.36), T(.5)]),
        Box2D([T(.4), T(.5)], [T(.19), T(.35)]),
        Box2D([T(.22), T(.46)], [T(.57), T(.75)]),
        Box2D([T(.75), T(1.)], [T(.64), T(.77)]),
        Box2D([T(.22), T(.8)], [T(.34), T(.51)])
    ]
)

TRI_BALLS{T<:AbstractFloat}(::Type{T} = Float64) = Compound2D(
    Shape2D{T}[
        Polygon([(T(.3), T(.3)), (T(.7), T(.3)), (T(.5), T(.65))]),
        Circle([T(.3), T(.3)], T(.15)),
        Circle([T(.7), T(.3)], T(.15)),
        Circle([T(.5), T(.65)], T(.15))
    ]
)

ISRR_POLY{T<:AbstractFloat}(::Type{T} = Float64) = Compound2D(
    Shape2D{T}[
        Polygon([(T(.0), T(.25)), (T(.27), T(.28)), (T(.17), T(.4)), (T(.0), T(.4))]),
        Polygon([(T(.5), T(.2)), (T(.2), T(.5)), (T(.25), T(.7)), (T(.4), T(.8)), (T(.6), T(.8)), (T(.7), T(.5))]),
        Polygon([(T(.55), T(.2)), (T(.75), T(.5)), (T(.85), T(.5)), (T(.85), T(.2))]),
        # Polygon([(T(.3), T(.6)), (T(.15), T(.85)), (T(.4), T(.6))]),
        Circle([T(.9), T(.65)], T(.1))
    ]
)

ISRR_POLY_WITH_SPIKE{T<:AbstractFloat}(::Type{T} = Float64) = Compound2D(
    Shape2D{T}[
        Polygon([(T(.0), T(.25)), (T(.27), T(.28)), (T(.17), T(.4)), (T(.0), T(.4))]),
        Polygon([(T(.5), T(.2)), (T(.2), T(.5)), (T(.25), T(.7)), (T(.4), T(.8)), (T(.6), T(.8)), (T(.7), T(.5))]),
        Polygon([(T(.55), T(.2)), (T(.75), T(.5)), (T(.85), T(.5)), (T(.85), T(.2))]),
        Polygon([(T(.3), T(.6)), (T(.15), T(.85)), (T(.4), T(.6))]),
        Circle([T(.9), T(.65)], T(.1))
    ]
)

EMPTY_2D{T<:AbstractFloat}(::Type{T} = Float64) = Compound2D(Shape2D{T}[])

nothing


text:
## ArbMatrix AbstractMatrix interface
Base.size(A::arb_mat_struct) = (A.r, A.c)
Base.size(A::Union{ArbMatrix,ArbRefMatrix}) = size(A.arb_mat)

function Base.getindex(A::arb_mat_struct, i::Integer, j::Integer)
    return ccall(
        @libarb(arb_mat_entry_ptr),
        Ptr{arb_struct},
        (Ref{arb_mat_struct}, Int, Int),
        A,
        i - 1,
        j - 1,
    )
end
Base.@propagate_inbounds function Base.getindex(A::ArbMatrix, i::Integer, j::Integer)
    @boundscheck checkbounds(A, i, j)
    x = Arb(prec = precision(A))
    x[] = A.arb_mat[i, j]
    x
end
Base.@propagate_inbounds function Base.getindex(A::ArbRefMatrix, i::Integer, j::Integer)
    @boundscheck checkbounds(A, i, j)
    return ArbRef(A.arb_mat[i, j], precision(A), cstruct(A))
end

"""
    ref(A::ArbMatrix, i, j)

Similar to `A[i,j]` but instead of an `Arb` returns an `ArbRef` which still shares the
memory with the `(i,j)`-th entry of `A`.
"""
Base.@propagate_inbounds function ref(
    A::Union{ArbMatrix,ArbRefMatrix},
    i::Integer,
    j::Integer,
)
    @boundscheck checkbounds(A, i, j)
    return ArbRef(A.arb_mat[i, j], precision(A), cstruct(A))
end

Base.setindex!(A::arb_mat_struct, x, i::Integer, j::Integer) = (set!(A[i, j], x); x)
Base.@propagate_inbounds function Base.setindex!(
    A::Union{ArbMatrix,ArbRefMatrix},
    x,
    i::Integer,
    j::Integer,
)
    @boundscheck checkbounds(A, i, j)
    ref(A, i, j)[] = x
    return x
end

## AcbMatrix AbstractMatrix interface
Base.size(A::acb_mat_struct) = (A.r, A.c)
Base.size(A::Union{AcbMatrix,AcbRefMatrix}) = size(A.acb_mat)

function Base.getindex(A::acb_mat_struct, i::Integer, j::Integer)
    return ccall(
        @libarb(acb_mat_entry_ptr),
        Ptr{acb_struct},
        (Ref{acb_mat_struct}, Int, Int),
        A,
        i - 1,
        j - 1,
    )
end
Base.@propagate_inbounds function Base.getindex(A::AcbMatrix, i::Integer, j::Integer)
    @boundscheck checkbounds(A, i, j)
    x = Acb(prec = precision(A))
    x[] = A.acb_mat[i, j]
    x
end
Base.@propagate_inbounds function Base.getindex(A::AcbRefMatrix, i::Integer, j::Integer)
    @boundscheck checkbounds(A, i, j)
    return AcbRef(A.acb_mat[i, j], precision(A), cstruct(A))
end

"""
    ref(A::AcbMatrix, i, j)

Similar to `A[i,j]` but instead of an `Acb` returns an `AcbRef` which still shares the
memory with the `(i,j)`-th entry of `A`.
"""
Base.@propagate_inbounds function ref(
    A::Union{AcbMatrix,AcbRefMatrix},
    i::Integer,
    j::Integer,
)
    @boundscheck checkbounds(A, i, j)
    return AcbRef(A.acb_mat[i, j], precision(A), cstruct(A))
end

Base.setindex!(A::acb_mat_struct, x, i::Integer, j::Integer) = (set!(A[i, j], x); x)
Base.@propagate_inbounds function Base.setindex!(
    A::Union{AcbMatrix,AcbRefMatrix},
    x,
    i::Integer,
    j::Integer,
)
    @boundscheck checkbounds(A, i, j)
    ref(A, i, j)[] = x
    return x
end

## Common methods

# General constructor
for T in [:ArbMatrix, :ArbRefMatrix, :AcbMatrix, :AcbRefMatrix]
    @eval function $T(A::AbstractMatrix; prec::Integer = _precision(first(A)))
        B = $T(size(A)...; prec = prec)
        # ensure to handle all kind of indices
        ax1, ax2 = axes(A)
        for (i, i′) in enumerate(ax1), (j, j′) in enumerate(ax2)
            B[i, j] = A[i′, j′]
        end
        return B
    end

    @eval function $T(v::AbstractVector; prec::Integer = _precision(first(v)))
        A = $T(length(v), 1; prec = prec)
        for (i, vᵢ) in enumerate(v)
            A[i, 1] = vᵢ
        end
        return A
    end
end

const Matrices = Union{ArbMatrix,ArbRefMatrix,AcbMatrix,AcbRefMatrix}

Base.copy(A::T) where {T<:Matrices} = copy!(T(size(A)...; prec = precision(A)), A)
function Base.copy!(A::T, B::T) where {T<:Matrices}
    @boundscheck size(A) == size(B) || throw(DimensionMismatch())
    set!(A, B)
    A
end
function Base.copyto!(A::T, B::T) where {T<:Matrices}
    @boundscheck size(A) == size(B) || throw(DimensionMismatch())
    set!(A, B)
    A
end

# Basic arithmetic
for (jf, af) in [(:+, :add!), (:-, :sub!)]
    @eval function Base.$jf(A::T, B::T) where {T<:Matrices}
        @boundscheck (
            size(A) == size(B) ||
            throw(DimensionMismatch("Matrix sizes are not compatible."))
        )
        C = T(size(A, 1), size(B, 2); prec = max(precision(A), precision(B)))
        $af(C, A, B)
        C
    end
end

Base.:(-)(A::Matrices) = neg!(similar(A), A)

function LinearAlgebra.mul!(C::T, A::T, B::T) where {T<:Matrices}
    @boundscheck (
        (size(C) == (size(A, 1), size(B, 2)) && size(A, 2) == size(B, 1)) ||
        throw(DimensionMismatch("Matrix sizes are not compatible."))
    )
    Arblib.mul!(C, A, B)
    C
end

function Base.:(*)(A::T, B::T) where {T<:Matrices}
    C = T(size(A, 1), size(B, 2); prec = max(precision(A), precision(B)))
    LinearAlgebra.mul!(C, A, B)
end

# scalar multiplication
function Base.:(*)(c::ArbLike, A::T) where {T<:Matrices}
    C = similar(A)
    Arblib.mul!(C, A, c)
end
function Base.:(*)(c::AcbLike, A::AcbMatrixLike)
    C = similar(A)
    Arblib.mul!(C, A, c)
end
function Base.:(*)(c::AcbLike, A::ArbMatrixLike)
    C = AcbMatrix(A)
    Arblib.mul!(C, C, c)
end
Base.:(*)(A::T, c::AcbLike) where {T<:Matrices} = c * A
Base.:(*)(A::T, c::ArbLike) where {T<:Matrices} = c * A

# scalar division
function Base.:(\)(c::ArbLike, A::T) where {T<:Matrices}
    C = similar(A)
    Arblib.div!(C, A, c)
end
function Base.:(\)(c::AcbLike, A::AcbMatrixLike)
    C = similar(A)
    Arblib.div!(C, A, c)
end
function Base.:(\)(c::AcbLike, A::ArbMatrixLike)
    C = AcbMatrix(A)
    Arblib.div!(C, C, c)
end
Base.:(/)(A::T, c::AcbLike) where {T<:Matrices} = c \ A
Base.:(/)(A::T, c::ArbLike) where {T<:Matrices} = c \ A

# lu factorization
function LinearAlgebra.lu!(A::T) where {T<:Matrices}
    ipiv = zeros(Int, size(A, 2))
    retcode = lu!(ipiv, A, A; prec = precision(A))
    LinearAlgebra.LU(A, ipiv, retcode > 0 ? 0 : 1)
end
function LinearAlgebra.lu(A::Matrices)
    lu = similar(A)
    ipiv = zeros(Int, size(A, 2))
    retcode = lu!(ipiv, lu, A; prec = precision(lu))
    LinearAlgebra.LU(lu, ipiv, retcode > 0 ? 0 : 1)
end

# ldiv! and \
function LinearAlgebra.ldiv!(Y::T, A::T, B::T) where {T<:Matrices}
    @boundscheck (
        (size(Y) == size(B) && size(A, 1) == size(A, 2) && size(A, 1) == size(B, 1)) ||
        throw(DimensionMismatch("Matrix sizes are not compatible."))
    )
    LinearAlgebra.ldiv!(Y, LinearAlgebra.lu(A), B)
end

function LinearAlgebra.ldiv!(A::T, B::T) where {T<:Matrices}
    @boundscheck (
        (size(A, 1) == size(A, 2)) ||
        throw(DimensionMismatch("Expected a square matrix as the left hand side."))
    )
    LinearAlgebra.ldiv!(B, LinearAlgebra.lu(A), B)
end

function LinearAlgebra.ldiv!(Y::T, A::LinearAlgebra.LU{<:Any,T}, B::T) where {T<:Matrices}
    @boundscheck (
        (size(Y) == size(B) && size(A, 1) == size(A, 2) && size(A, 1) == size(B, 1)) ||
        throw(DimensionMismatch("Matrix sizes are not compatible."))
    )
    Arblib.solve_lu_precomp!(Y, A.ipiv, A.factors, B)
    Y
end
LinearAlgebra.ldiv!(A::LinearAlgebra.LU{<:Any,T}, B::T) where {T<:Matrices} =
    LinearAlgebra.ldiv!(B, A, B)

function Base.:(\)(A::T, B::T) where {T<:Matrices}
    Y = T(size(A, 2), size(B, 2); prec = max(precision(A), precision(B)))
    LinearAlgebra.ldiv!(Y, A, B)
end
function Base.:(\)(A::LinearAlgebra.LU{<:Any,T}, B::T) where {T<:Matrices}
    Y = T(size(A, 2), size(B, 2); prec = max(precision(A.factors), precision(B)))
    LinearAlgebra.ldiv!(Y, A, B)
end

# inverse
LinearAlgebra.inv(A::Matrices) = (B = similar(A); inv!(B, A); B)


text:
module UnitaryDecomposition

using LinearAlgebra
using Optim

include("Pauli.jl")
include("BlochRotation.jl")
include("BasicGates.jl")

export UDecompose, Pauli, BlochRotation

#This code will decompose an unitary into blocks of single qubit rotations
#and cnots between adjacent qubits

function U(n,p,nq)
Ut=1
for i=1:nq
 Ut=kron(Ut,BlochRotation(n[i,:],p[i]))
end
for j=0:nq-2
  Ut=Ut*kron(Diagonal(ones(2^j)),CNOT12,Diagonal(ones(2^(nq-j-2))))
end
return Ut
end

function Utestlist(U)
# U will have to have the form [n[,] a[] p[] nq]
# Only leave U[i] on the non fixed parameters, and fix the ones I want
  n=[U[1] U[2] U[3]; U[4] U[5] U[6]] # This depends on nq
  p=[U[7] U[8]]
  nq=2 #usually fixed, never optimize over nq
  return U(n,p,nq)
end

function UDecompose(U,N)
#U is the unitary we want to decompose
#N is the number of blocks we want
f(x)=norm(U-Ublock(N)(x))
nq=2
x0=rand(N,4*nq)
res=optimize(f,x0)
return res.minimizer, res.minimum
end

function Ublock(N)
# N is the number of layers
if N==1
return Utestlist
end
#if N==2
#Ublock2(U)=Utestlist(U[1,:])*Utestlist(U[2,:])
#return Ublock2
#end
if N>1
Ub(U)=Utestlist(U[1,:])*Ublock(N-1)(U[2:end,:])
end
end

end

max_stars_repo_name: JuliMaristany/UnitaryDecomposition

max_issues_repo_name: JuliMaristany/UnitaryDecomposition

max_forks_repo_name: JuliMaristany/UnitaryDecomposition


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

module StorageManagementClient

using Random
using Dates
using HTTP
using Swagger
import Swagger: field_name, property_type, hasproperty, validate_property, SwaggerApi, SwaggerModel
import Base: convert, propertynames

include("modelincludes.jl")

include("api_OperationsApi.jl")
include("api_SkusApi.jl")
include("api_StorageAccountsApi.jl")
include("api_UsageApi.jl")

# export models
export convert, AccountSasParameters
export convert, CheckNameAvailabilityResult
export convert, CustomDomain
export convert, Dimension
export convert, Encryption
export convert, EncryptionService
export convert, EncryptionServices
export convert, Endpoints
export convert, IPRule
export convert, Identity
export convert, KeyVaultProperties
export convert, ListAccountSasResponse
export convert, ListServiceSasResponse
export convert, MetricSpecification
export convert, NetworkRuleSet
export convert, Operation
export convert, OperationDisplay
export convert, OperationListResult
export convert, OperationProperties
export convert, Resource
export convert, Restriction
export convert, SKUCapability
export convert, ServiceSasParameters
export convert, ServiceSpecification
export convert, Sku
export convert, StorageAccountCheckNameAvailabilityParameters
export convert, StorageAccountCreateParameters
export convert, StorageAccountKey
export convert, StorageAccountListKeysResult
export convert, StorageAccountListResult
export convert, StorageAccountProperties
export convert, StorageAccountPropertiesCreateParameters
export convert, StorageAccountPropertiesUpdateParameters
export convert, StorageAccountRegenerateKeyParameters
export convert, StorageAccountUpdateParameters
export convert, StorageSkuListResult
export convert, Usage
export convert, UsageListResult
export convert, UsageName
export convert, VirtualNetworkRule
export convert, StorageAccount

# export operations
export convert, OperationsApi, SkusApi, StorageAccountsApi, UsageApi

export check_required, field_name, property_type, hasproperty, propertynames, validate_property, convert

end