text:
using Arblib: Acb, AcbRef, AcbVector, AcbRefVector, AcbMatrix, AcbRefMatrix

export certify,
    SolutionCertificate,
    CertificationResult,
    CertificationCache,
    is_certified,
    is_real,
    is_complex,
    is_positive,
    solution_candidate,
    certified_solution_interval,
    certified_solution_interval_after_krawczyk,
    certificate_index,
    solution_approximation,
    certificates,
    distinct_certificates,
    ncertified,
    nreal_certified,
    ncomplex_certified,
    ndistinct_certified,
    ndistinct_real_certified,
    ndistinct_complex_certified,
    show_straight_line_program,
    save,
    # deprecated
    initial_solution,
    certified_solution

"""
    SolutionCertificate

Result of [`certify`](@ref) for a single solution. Contains the initial solutions
and if the certification was successfull a vector of complex intervals where the true
solution is contained in.
"""
Base.@kwdef struct SolutionCertificate
    solution_candidate::AbstractVector
    certified::Bool
    real::Bool = false
    complex::Bool = false
    index::Union{Nothing,Int} = nothing
    prec::Int = 53
    # I, I′ ∈ 𝕀ℂⁿ and a certified solution has I′ ⊊ I
    I::Union{Nothing,AcbMatrix} = nothing
    I′::Union{Nothing,AcbMatrix} = nothing
    x̃::Union{Nothing,AcbMatrix} = nothing
    Y::Union{Nothing,AcbMatrix} = nothing

    # We also store a double precision representation of the midpoint of x₁
    # as the best available double precision estimate of the solution
    solution::Union{Nothing,Vector{ComplexF64}} = nothing
end

"""
    solution_candidate(certificate::SolutionCertificate)

Returns the given provided solution candidate.
"""
solution_candidate(C::SolutionCertificate) = C.solution_candidate
@deprecate initial_solution(cert::SolutionCertificate) solution_candidate(cert)

"""
    is_certified(certificate::SolutionCertificate)

Returns `true` if `certificate` is a certificate that `certified_solution_interval(certificate)`
contains a unique zero.
"""
is_certified(C::SolutionCertificate) = C.certified

"""
    is_real(certificate::SolutionCertificate)

Returns `true` if `certificate` certifies that the certified solution interval
contains a true real zero of the system.
If `false` is returned then this does not necessarily mean that the true solution is not real.
"""
is_real(C::SolutionCertificate) = C.real

"""
    is_complex(certificate::SolutionCertificate)

Returns `true` if `certificate` certifies that the certified solution interval
contains a true complex zero of the system.
"""
is_complex(C::SolutionCertificate) = C.complex

"""
    is_positive(certificate::SolutionCertificate)

Returns `true` if `is_certified(certificate)` is `true` and the unique zero contained
in `certified_solution_interval(certificate)` is real and positive.
"""
function is_positive(C::SolutionCertificate)
    if !is_certified(C) || !is_real(C)
        return false
    else
        all(i -> Arblib.is_positive(real(Arblib.ref(C.I′, i, 1))), 1:length(C.I′))
    end
end

"""
    is_positive(certificate::SolutionCertificate, i)

Returns `true` if `is_certified(certificate)` is `true` and `i`-th coordinate of the
unique zero contained in `certified_solution_interval(certificate)` is real and positive.
"""
function is_positive(C::SolutionCertificate, i::Integer)
    if !is_certified(C) || !is_real(C)
        return false
    else
        Arblib.is_positive(real(Arblib.ref(C.I′, i, 1)))
    end
end

"""
    certified_solution_interval(certificate::SolutionCertificate)

Returns an `Arblib.AcbMatrix` representing a vector of complex intervals where a unique
zero of the system is contained in.
Returns `nothing` if `is_certified(certificate)` is `false`.
"""
certified_solution_interval(certificate::SolutionCertificate) = certificate.I
@deprecate certified_solution(certificate) certified_solution_interval(certificate)

"""
    certified_solution_interval_after_krawczyk(certificate::SolutionCertificate)

Returns an `Arblib.AcbMatrix` representing a vector of complex intervals where a unique
zero of the system is contained in.
This is the result of applying the Krawczyk operator to `certified_solution_interval(certificate)`.
Returns `nothing` if `is_certified(certificate)` is `false`.
"""
certified_solution_interval_after_krawczyk(certificate::SolutionCertificate) =
    certificate.I′

"""
    certificate_index(certificate::SolutionCertificate)

Return the index of the solution certificate. Here the index refers to the index of the
provided solution candidates.
"""
certificate_index(certificate::SolutionCertificate) = C.index

"""
    precision(certificate::SolutionCertificate)

Returns the maximal precision used to produce the given `certificate`.
"""
Base.precision(certificate::SolutionCertificate) = certificate.prec

"""
    solution_approximation(certificate::SolutionCertificate)

If `is_certified(certificate)` is `true` this returns the midpoint of the
[`certified_solution_interval`](@ref) of the given `certificate` as a `Vector{ComplexF64}`.
Returns `nothing` if `is_certified(certificate)` is `false`.
"""
solution_approximation(certificate::SolutionCertificate) = certificate.solution

"""
    krawczyk_operator_parameters(cert::SolutionCertificate)

Returns a `NamedTuple` `(x, Y)` with the parameters of the Krawczyk operator following [BRT20].
"""
krawczyk_operator_parameters(cert::SolutionCertificate) = (x = cert.x̃, Y = cert.Y)

function Base.show(f::IO, cert::SolutionCertificate)
    println(f, "SolutionCertificate:")
    println(f, "solution_candidate = [")
    for z in solution_candidate(cert)
        print(f, "  ")
        print(f, z)
        println(f, ",")
    end
    println(f, "]")
    print(f, "is_certified = ", is_certified(cert))
    if !isnothing(certified_solution_interval(cert))
        println(f)
        println(f, "certified_solution_interval = [")
        for z in certified_solution_interval(cert)
            print(f, "  ")
            print(f, z)
            println(f, ",")
        end
        println(f, "]")
        println(f, "precision = ", cert.prec)
        print(f, "is_real = ", is_real(cert))
    end
    if !isnothing(cert.index)
        println(f)
        print(f, "index = ", cert.index)
    end
end

"""
    CertificationResult

The result of [`certify`](@ref) for multiple solutions.
Contains a vector of [`SolutionCertificate`](@ref) as well as a list of certificates
which correspond to the same true solution.
"""
struct CertificationResult
    certificates::Vector{SolutionCertificate}
    duplicates::Vector{Vector{Int}}
    slp::ModelKit.Interpreter
    slp_jacobian::ModelKit.Interpreter
end

"""
    certificates(R::CertificationResult)

Obtain the stored [`SolutionCertificate`](@ref)s.
"""
certificates(R::CertificationResult) = R.certificates

"""
    distinct_certificates(R::CertificationResult)

Obtain the certificates corresponding to the determined distinct solution intervals.
"""
function distinct_certificates(C::CertificationResult)
    cs = certificates(C)
    isempty(C.duplicates) && return cs
    indices = trues(length(cs))
    for d in C.duplicates, k = 2:length(d)
        indices[d[k]] = false
    end
    cs[indices]
end

"""
    ncertified(R::CertificationResult)

Returns the number of certified solutions.
"""
ncertified(R::CertificationResult) = count(is_certified, R.certificates)

"""
    nreal_certified(R::CertificationResult)

Returns the number of certified real solutions.
"""
nreal_certified(R::CertificationResult) =
    count(r -> is_certified(r) && is_real(r), R.certificates)

"""
    ncomplex_certified(R::CertificationResult)

Returns the number of certified complex solutions.
"""
ncomplex_certified(R::CertificationResult) =
    count(r -> is_certified(r) && is_complex(r), R.certificates)

"""
    ndistinct_certified(R::CertificationResult)

Returns the number of distinct certified solutions.
"""
function ndistinct_certified(R::CertificationResult)
    ncert = ncertified(R)
    if isempty(R.duplicates)
        return ncert
    else
        return ncert - sum(length, R.duplicates) + length(R.duplicates)
    end
end

"""
    ndistinct_real_certified(R::CertificationResult)

Returns the number of distinct certified real solutions.
"""
function ndistinct_real_certified(R::CertificationResult)
    ncert = nreal_certified(R)
    if isempty(R.duplicates)
        return ncert
    else
        ncert - sum(R.duplicates) do dup
            is_real(R.certificates[dup[1]]) ? length(dup) - 1 : 0
        end
    end
end

"""
    ndistinct_complex_certified(R::CertificationResult)

Returns the number of distinct certified complex solutions.
"""
function ndistinct_complex_certified(R::CertificationResult)
    ncert = ncomplex_certified(R)
    if isempty(R.duplicates)
        return ncert
    else
        ncert - sum(R.duplicates) do dup
            is_complex(R.certificates[dup[1]]) ? length(dup) - 1 : 0
        end
    end
end

"""
    show_straight_line_program(R::CertificationResult)
    show_straight_line_program(io::IO, R::CertificationResult)

Print a representation of the used straight line program.
"""
show_straight_line_program(R::CertificationResult) = show_straight_line_program(stdout, R)
show_straight_line_program(io::IO, R::CertificationResult) =
    ModelKit.show_instructions(io, R.slp)

function Base.show(io::IO, R::CertificationResult)
    println(io, "CertificationResult")
    println(io, "===================")
    println(io, "• $(length(R.certificates)) solution candidates given")
    ncert = ncertified(R)
    print(io, "• $ncert certified solution intervals")
    nreal = nreal_certified(R)
    ncomplex = ncomplex_certified(R)
    print(io, " ($nreal real, $ncomplex complex")
    if nreal + ncomplex < ncert
        println(io, ", $(ncert - (nreal + ncomplex)) undecided)")
    else
        println(io, ")")
    end
    ndist_cert = ndistinct_certified(R)
    print(io, "• $ndist_cert distinct certified solution intervals")
    ndist_real = ndistinct_real_certified(R)
    ndist_complex = ndistinct_complex_certified(R)
    print(io, " ($ndist_real real, $ndist_complex complex")
    if ndist_real + ndist_complex < ndist_cert
        print(io, ", $(ndist_cert - (ndist_real + ndist_complex)) undecided)")
    else
        print(io, ")")
    end
end

"""
    save(filename, C::CertificationResult)

Store a text representation of the certification result `C` on disk.
"""
function save(filename, R::CertificationResult)
    open(filename, "w") do f
        println(f, "## Summary")
        show(f, R)
        println(f, "\n\n## Certificates")
        for cert in certificates(R)
            show(f, cert)
            println(f)
        end
    end
    filename
end

"""
    CertificationCache(F::AbstractSystem)

Contains the necessary data structures for [`certify`](@ref).
"""
Base.@kwdef mutable struct CertificationCache{T₁,T₂}
    eval_interpreter::ModelKit.Interpreter{T₁}
    jac_interpreter::ModelKit.Interpreter{T₂}
    newton_cache::NewtonCache{MatrixWorkspace{Matrix{ComplexF64}}}
    # data for krawczyc_step
    C::Matrix{ComplexF64}
    r₀::Vector{IComplexF64}
    Δx₀::Vector{IComplexF64}
    ix̃₀::Vector{IComplexF64}
    Jx₀::Matrix{IComplexF64}
    M::Matrix{IComplexF64}
    δx::Vector{IComplexF64}
    interpreter_cache::ModelKit.InterpreterCache{Vector{IComplexF64},Nothing}

    arb_prec::Int
    arb_C::AcbRefMatrix
    arb_r₀::AcbRefMatrix # m × 1
    arb_Δx₀::AcbRefMatrix # m × 1
    arb_x̃₀::AcbRefMatrix # m × 1
    arb_x₀::AcbRefMatrix # m × 1
    arb_x₁::AcbRefMatrix # m × 1
    arb_J_x₀::AcbRefMatrix
    arb_M::AcbRefMatrix
    arb_δx::AcbRefMatrix # m × 1
    arb_mag::Arblib.Mag
    arb_interpreter_cache::ModelKit.InterpreterCache{AcbRefVector,Acb}
end
function CertificationCache(F::AbstractSystem)
    m, n = size(F)
    m == n || error("We can only certify square systems.")
    f = System(F)
    interpreter = ModelKit.interpreter(f)
    jac_interpreter = ModelKit.jacobian_interpreter(f)
    tape_length = max(
        ModelKit.cache_min_length(interpreter),
        ModelKit.cache_min_length(jac_interpreter),
    )
    arb_prec = 128
    CertificationCache(;
        eval_interpreter = interpreter,
        jac_interpreter = jac_interpreter,
        newton_cache = NewtonCache(F; optimize_data_structure = false),
        C = zeros(ComplexF64, m, m),
        r₀ = zeros(IComplexF64, m),
        Δx₀ = zeros(IComplexF64, m),
        ix̃₀ = zeros(IComplexF64, m),
        Jx₀ = zeros(IComplexF64, m, m),
        M = zeros(IComplexF64, m, m),
        δx = zeros(IComplexF64, m),
        interpreter_cache = ModelKit.InterpreterCache(IComplexF64, tape_length),
        arb_prec = arb_prec,
        arb_C = AcbRefMatrix(m, m; prec = arb_prec),
        arb_r₀ = AcbRefMatrix(m, 1; prec = arb_prec),
        arb_Δx₀ = AcbRefMatrix(m, 1; prec = arb_prec),
        arb_x̃₀ = AcbRefMatrix(m, 1; prec = arb_prec),
        arb_x₀ = AcbRefMatrix(m, 1; prec = arb_prec),
        arb_x₁ = AcbRefMatrix(m, 1; prec = arb_prec),
        arb_J_x₀ = AcbRefMatrix(m, m; prec = arb_prec),
        arb_M = AcbRefMatrix(m, m; prec = arb_prec),
        arb_δx = AcbRefMatrix(m, 1; prec = arb_prec),
        arb_mag = Arblib.Mag(),
        arb_interpreter_cache = ModelKit.InterpreterCache(
            Acb,
            tape_length;
            prec = arb_prec,
        ),
    )
end

Base.setprecision(M::AcbRefMatrix, p::Int) = AcbRefMatrix(M.acb_mat, p)
function set_arb_precision!(cache::CertificationCache, p::Int)
    cache.arb_prec == p && return cache
    cache.arb_prec = p
    cache.arb_r₀ = setprecision(cache.arb_r₀, p)
    cache.arb_Δx₀ = setprecision(cache.arb_Δx₀, p)
    cache.arb_x̃₀ = setprecision(cache.arb_x̃₀, p)
    cache.arb_x₀ = setprecision(cache.arb_x₀, p)
    cache.arb_x₁ = setprecision(cache.arb_x₁, p)
    cache.arb_J_x₀ = setprecision(cache.arb_J_x₀, p)
    cache.arb_M = setprecision(cache.arb_M, p)
    cache.arb_δx = setprecision(cache.arb_δx, p)
    ModelKit.set_arb_precision!(cache.arb_interpreter_cache, p)

    cache
end

"""
    certify(F, solutions, [p, certify_cache]; options...)
    certify(F, result, [p, certify_cache]; options...)

Attempt to certify that the given approximate `solutions` correspond to true solutions
of the polynomial system ``F(x;p)``. The system ``F`` has to be an (affine) square
polynomial system. Also attemps to certify for each solutions whether it
approximates a real solution. The certification is done using interval arithmetic and the
Krawczyk method[^Moo77]. Returns a [`CertificationResult`](@ref) which additionall returns
the number of distinct solutions. For more details of the implementation see [^BRT20].

## Options

* `show_progress = true`: If `true` shows a progress bar of the certification process.
* `max_precision = 256`: The maximal accuracy (in bits) that is used in the certification process.
* `compile = false`: See the [`solve`](@ref) documentation.

## Example
We take the [first example](https://www.juliahomotopycontinuation.org/guides/introduction/#a-first-example) from our
introduction guide.
```julia
@var x y
# define the polynomials
f₁ = (x^4 + y^4 - 1) * (x^2 + y^2 - 2) + x^5 * y
f₂ = x^2+2x*y^2 - 2y^2 - 1/2
F = System([f₁, f₂], variables = [x,y])
result = solve(F)
```
```
Result with 18 solutions
========================
• 18 paths tracked
• 18 non-singular solutions (4 real)
• random seed: 0xcaa483cd
• start_system: :polyhedral
```
We see that we obtain 18 solutions and it seems that 4 solutions are real. However,
this is based on heuristics. To be absolute certain we can certify the result

```julia
certify(F, result)
```
```
CertificationResult
===================
• 18 solution candidates given
• 18 certified solution intervals (4 real, 14 complex)
• 18 distinct certified solution intervals (4 real, 14 complex)
```

and see that there are indeed 18 solutions and that they are all distinct.

[^Moo77]: Moore, Ramon E. "A test for existence of solutions to nonlinear systems." SIAM Journal on Numerical Analysis 14.4 (1977): 611-615.
[^BRT20]: Breiding, P., Rose, K. and Timme, S. "Certifying zeros of polynomial systems using interval arithmetic." arXiv:2011.05000.
"""
function certify end

struct CertificationParameters
    params::Vector{ComplexF64}
    interval_params::Vector{IComplexF64}
    arb_interval_params::AcbRefVector
end

function CertificationParameters(p::AbstractVector; prec::Int = 256)
    arb_ip = AcbRefVector(length(p); prec = prec)
    for (i, p_i) in enumerate(p)
        x = arb_ip[i]
        x[] = p_i
    end
    CertificationParameters(
        convert(Vector{ComplexF64}, p),
        convert(Vector{IComplexF64}, p),
        arb_ip,
    )
end

certification_parameters(p::AbstractVector; prec::Int = 256) = CertificationParameters(p)
certification_parameters(::Nothing; prec::Int = 256) = nothing

complexF64_params(C::CertificationParameters) = C.params
complexF64_interval_params(C::CertificationParameters) = C.interval_params
arb_interval_params(C::CertificationParameters) = C.arb_interval_params
complexF64_params(::Nothing) = nothing
complexF64_interval_params(::Nothing) = nothing
arb_interval_params(::Nothing) = nothing

function _certify(
    F::AbstractSystem,
    solution_candidates::AbstractArray{<:AbstractArray{<:Number}},
    p::Union{Nothing,CertificationParameters},
    cache::CertificationCache;
    show_progress::Bool = true,
    max_precision::Int = 256,
    check_distinct::Bool = true,
)
    certificates = SolutionCertificate[]

    m, n = size(F)
    m == n || throw(ArgumentError("We can only certify solutions to square systems."))

    if isnothing(p) && nparameters(System(F)) > 0
        throw(ArgumentError("The given system expects parameters but none are given."))
    end

    if !show_progress
        for (i, s) in enumerate(solution_candidates)
            push!(certificates, certify_solution(F, s, p, cache, i))
        end
    else
        # Create progress meter
        n = length(solution_candidates)
        desc = "Certifying $n solutions... "
        barlen = min(ProgressMeter.tty_width(desc, stdout, false), 40)
        progress = ProgressMeter.Progress(
            n;
            dt = 0.2,
            desc = desc,
            barlen = barlen,
            color = :green,
            output = stdout,
        )
        progress.tlast += progress.dt

        # Go over all solution
        ncertified = 0
        nreal_certified = 0
        for (i, s) in enumerate(solution_candidates)
            r = certify_solution(F, s, p, cache, i)
            push!(certificates, r)
            ncertified += is_certified(r)
            nreal_certified += is_certified(r) && is_real(r)
            update_certify_progress!(progress, i, ncertified, nreal_certified)
        end
    end
    if check_distinct
        duplicates = find_duplicates(certificates)
    else
        duplicates = Vector{Vector{Int}}()
    end

    CertificationResult(
        certificates,
        duplicates,
        cache.eval_interpreter,
        cache.jac_interpreter,
    )
end

function update_certify_progress!(progress, k, ncertified, nreal_certified)
    t = time()
    if k == progress.n || t > progress.tlast + progress.dt
        showvalues = make_certify_showvalues(k, ncertified, nreal_certified)
        ProgressMeter.update!(progress, k; showvalues = showvalues)
    end
    nothing
end
@noinline function make_certify_showvalues(k, ncertified, nreal_certified)
    (
        ("# solutions candidates considered", k),
        ("# certified solution intervals (real)", "$(ncertified) ($nreal_certified)"),
    )
end

function certify_solution(
    F::AbstractSystem,
    solution_candidate::AbstractVector,
    cert_params::Union{Nothing,CertificationParameters},
    cert_cache::CertificationCache,
    index::Int,
    max_precision::Int = 256,
    refine_solution::Bool = true,
)
    @unpack C, arb_C, arb_x̃₀ = cert_cache
    n = size(C, 1)

    # refine solution to machine precicision
    res = newton(
        F,
        solution_candidate,
        complexF64_params(cert_params),
        InfNorm(),
        cert_cache.newton_cache;
        rtol = 8 * eps(),
        atol = 0.0,
        # this should already be an approximate zero
        extended_precision = true,
        max_iters = 8,
    )

    x̃₀ = solution(res)
    LA.inv!(C, cert_cache.newton_cache.J)

    certified, x₁, x₀, is_real = ε_inflation_krawczyk(x̃₀, cert_params, C, cert_cache)

    if certified
        return SolutionCertificate(
            solution_candidate = solution_candidate,
            certified = true,
            real = is_real,
            complex = any(xi -> !(0.0 in imag(xi)), x₁),
            index = index,
            prec = 53,
            I = AcbMatrix(x₀; prec = 53),
            I′ = AcbMatrix(x₁; prec = 53),
            x̃ = AcbMatrix(x̃₀; prec = 53),
            Y = AcbMatrix(C; prec = 53),
            solution = mid.(x₁),
        )
    end

    # If not certified we do the computation in higher precision using Arb
    prec = 128
    set_arb_precision!(cert_cache, prec)

    # We keep the same C matrix for now.
    for j = 1:n, i = 1:n
        arb_C[i, j][] = C[i, j]
    end
    for (i, x̃₀_i) in enumerate(x̃₀)
        arb_x̃₀[i][] = x̃₀_i
    end
    while (prec <= max_precision)
        certified, arb_x₁, arb_x₀, is_real =
            arb_ε_inflation_krawczyk(arb_x̃₀, cert_params, arb_C, cert_cache; prec = prec)

        if certified
            I = AcbMatrix(n, 1; prec = prec)
            I′ = AcbMatrix(n, 1; prec = prec)
            x̃ = AcbMatrix(n, 1; prec = prec)
            Y = AcbMatrix(n, n; prec = prec)
            Arblib.set!(I, arb_x₀)
            Arblib.set!(I′, arb_x₁)
            Arblib.set!(x̃, arb_x̃₀)
            Arblib.set!(Y, arb_C)
            return SolutionCertificate(
                solution_candidate = solution_candidate,
                certified = true,
                real = is_real,
                complex = any(xi -> !Arblib.contains_zero(Arblib.imagref(xi)), arb_x₁),
                index = index,
                prec = prec,
                I = I,
                I′ = I′,
                x̃ = x̃,
                Y = Y,
                solution = [ComplexF64(arb_x₁[i]) for i = 1:n],
            )
        end

        prec += 64
        set_arb_precision!(cert_cache, prec)
    end

    # certification failed
    return SolutionCertificate(
        solution_candidate = solution_candidate,
        certified = false,
        index = index,
    )
end

function Base.Vector{ComplexF64}(A::Arblib.AcbMatrixLike)
    @assert size(A, 2) == 1
    [ComplexF64(Arblib.ref(A, i, 1).acb_ptr) for i = 1:size(A, 1)]
end

function ε_inflation_krawczyk(x̃₀, p::Union{Nothing,CertificationParameters}, C, cert_cache)
    @unpack C, r₀, Δx₀, ix̃₀, Jx₀, M, δx, interpreter_cache = cert_cache
    J_x₀ = Jx₀

    ix̃₀ .= IComplexF64.(x̃₀)
    # r₀ = F(ix̃₀)
    ModelKit.execute!(
        r₀,
        cert_cache.eval_interpreter,
        ix̃₀,
        complexF64_interval_params(p),
        interpreter_cache,
    )
    # iΔx = C * r₀
    sqr_mul!(Δx₀, C, r₀)

    # Perform ε-inflation. We choose a different εᵢ per coordinate. This matches our strategy
    # to use a weighted norm.
    x₀ = map(x̃₀, Δx₀) do x̃₀_i, Δx₀_i
        εᵢ = 512 * max(mag(Δx₀_i), eps())
        complex(
            Interval(real(x̃₀_i) - εᵢ, real(x̃₀_i) + εᵢ),
            Interval(imag(x̃₀_i) - εᵢ, imag(x̃₀_i) + εᵢ),
        )
    end

    ModelKit.execute!(
        nothing,
        J_x₀,
        cert_cache.jac_interpreter,
        x₀,
        complexF64_interval_params(p),
        interpreter_cache,
    )

    x₁ = similar(x₀)
    # x₁ = (x̃₀ - C * F([x̃₀])) + (I - C * J(x₀)) * (x₀ - x̃₀)
    #    = (x̃₀ - C * F([x̃₀])) - (C * J(x₀) - I) * (x₀ - x̃₀)
    #    = (x̃₀ - Δx₀) - (C * J(x₀) - I) * (x₀ - x̃₀)

    # Define M = C * J(x₀) - I
    sqr_mul!(M, C, J_x₀)
    for i = 1:size(M, 1)
        M[i, i] -= 1
    end
    # Necessary condition is ||M|| < 1 / √2
    # We lower bound 1 / √2 by 0.7071
    if IntervalArithmetic.inf_norm_bound(M) < 0.7071
        δx .= x₀ .- x̃₀
        sqr_mul!(x₁, M, δx)
        x₁ .= (x̃₀ .- Δx₀) .- x₁
        certified = all2(isinterior, x₁, x₀)
        is_real =
            certified ? all2((x₁_i, x₀_i) -> isinterior(conj(x₁_i), x₀_i), x₁, x₀) : false
    else
        certified = false
        is_real = false
    end
    certified, x₁, x₀, is_real
end

function arb_ε_inflation_krawczyk(
    x̃₀::AcbRefMatrix,
    p::Union{Nothing,CertificationParameters},
    C::Arblib.AcbMatrixLike,
    cert_cache;
    prec::Int,
)
    @unpack arb_r₀, arb_Δx₀, arb_x̃₀, arb_J_x₀, arb_M, arb_δx, arb_mag, arb_x₀, arb_x₁ =
        cert_cache
    r₀, Δx₀, x̃₀, x₀, x₁, J_x₀, M, δx, m =
        arb_r₀, arb_Δx₀, arb_x̃₀, arb_x₀, arb_x₁, arb_J_x₀, arb_M, arb_δx, arb_mag
    @unpack arb_interpreter_cache = cert_cache

    m = arb_mag
    x̃₀′ = x̄₁ = Δx₀
    # Perform a simple newton refinement using arb_C until we cannot improve the
    # accuracy anymore
    acc = Inf
    max_iters = 10
    for i = 1:max_iters
        ModelKit.execute!(
            r₀,
            cert_cache.eval_interpreter,
            x̃₀,
            arb_interval_params(p),
            arb_interpreter_cache,
        )
        Arblib.mul!(Δx₀, C, r₀)

        acc_new = Arblib.get(Arblib.bound_inf_norm!(m, Δx₀))
        if acc_new > acc || i == max_iters
            break
        end
        acc = acc_new
        Arblib.sub!(x̃₀, x̃₀, Δx₀)
        Arblib.get_mid!(x̃₀, x̃₀)
    end

    # Perform ε-inflation
    n = length(x̃₀)
    Arblib.get_mid!(x₀, x̃₀)
    # We increase the radius by 2^(prec/4) to avoid hitting the precision limit.
    # We choose a dynamic increase to account for bad situations where any fixed choice
    # would be insufficient.
    incr_factor = exp2(div(prec, 4))
    mach_eps = exp2(-prec)
    for i = 1:n
        m[] = max(magF64(Δx₀[i], m), mach_eps) * incr_factor
        Arblib.add_error!(x₀[i], m)
    end

    ModelKit.execute!(
        nothing,
        J_x₀,
        cert_cache.jac_interpreter,
        x₀,
        arb_interval_params(p),
        arb_interpreter_cache,
    )

    # x₁ = (x̃₀ - C * F([x̃₀])) + (I - C * J(x₀)) * (x₀ - x̃₀)
    #    = (x̃₀ - C * F([x̃₀])) - (C * J(x₀) - I) * (x₀ - x̃₀)
    #    = (x̃₀ - Δx₀) - (C * J(x₀) - I) * (x₀ - x̃₀)

    # Define M = C * J(x₀) - I
    Arblib.mul!(M, C, J_x₀)
    for i = 1:n
        Arblib.sub!(M[i, i], M[i, i], 1)
    end

    # Necessary condition is ||M|| < 1 / √2
    # We lower bound 1 / √2 by 0.7071
    if Arblib.get(Arblib.bound_inf_norm!(m, M)) < 0.7071
        # x₀ - x̃₀
        Arblib.sub!(δx, x₀, x̃₀)
        # (C * J(x₀) - I) * (x₀ - x̃₀)
        Arblib.mul!(x₁, M, δx)
        # (x̃₀ - Δx₀)
        Arblib.sub!(x̃₀′, x̃₀, Δx₀)
        # (x̃₀ - Δx₀) - (C * J(x₀) - I) * (x₀ - x̃₀)
        Arblib.sub!(x₁, x̃₀′, x₁)

        certified = Bool(Arblib.contains(x₀, x₁))

        # check realness
        Arblib.conjugate!(x̄₁, x₁)
        is_real = Bool(Arblib.contains(x₀, x̄₁))
    else
        certified = false
        is_real = false
    end

    if !certified
        # Update the approximation of the inverse
        Arblib.get_mid!(J_x₀, J_x₀)
        Arblib.inv!(C, J_x₀)
    end

    certified, AcbMatrix(x₁), AcbMatrix(x₀), is_real
end

magF64(x, mag) = Arblib.get(Arblib.get!(mag, x))

"Simple matrix multiplication using the muladd method which is much faster four our case."
function sqr_mul!(C, A, B)
    n = size(A, 1)
    C .= zero(eltype(C))
    @inbounds for j = 1:size(B, 2)
        for k = 1:n
            bkj = B[k, j]
            for i = 1:n
                C[i, j] = muladd(A[i, k], bkj, C[i, j])
            end
        end
    end
    C
end

"""
    find_duplicates(certificates)

Find all duplicate solutions. A solution is considered to be duplicate if the unique
region overlap.
"""
function find_duplicates(certificates::AbstractVector{SolutionCertificate})
    # The strategy to indentify duplication candidates is as following.
    # We first compute the squared euclidean distance to a random point
    # in interval arithmetic. We then compute all overlapping intervals and return this
    # list of tuples.
    # The complexity of this is O(n log(n)) (we need to perform an array sort)

    # 1 Compute squared euclidean distance to a random point
    original_indices = Int[]

    n = length(solution_candidate(first(certificates)))
    pt = randn(ComplexF64, n)

    sq_eucl_distances = Interval{Float64}[]

    a, b = Arblib.Arf(prec = 53), Arblib.Arf(prec = 53)
    for (i, cert) in enumerate(certificates)
        if is_certified(cert)
            push!(original_indices, i)

            d = zero(Interval{Float64})
            for i = 1:n
                yᵢ = IComplexF64(cert.I′[i], a, b)
                d +=
                    IntervalArithmetic.sqr(real(yᵢ) - real(pt[i])) +
                    IntervalArithmetic.sqr(imag(yᵢ) - imag(pt[i]))
            end
            push!(sq_eucl_distances, d)
        end
    end

    # 2) Find all overlapping intervals

    # strategy:
    # Look at close and open of an interval separately and merge all these in one large vector
    # (keeping track if open or close of an interval),
    # sort this (first by value, then open, then close).
    # Now iterate through sorted array and keep track of open intervals

    # value, index, and false if opening and true if closing
    interval_bounds = Tuple{Float64,Int,Bool}[]
    for (i, d) in enumerate(sq_eucl_distances)
        push!(interval_bounds, (d.lo, i, false), (d.hi, i, true))
    end

    sort!(interval_bounds; lt = (a, b) -> a[1] != b[1] ? a[1] < b[1] : !a[3])

    current_open = Int[]
    # dict to help group duplicates:
    # Assume (1,3) and (3,5) overlap then we want to have entries 1=>1, 3=>1 and 5=>1
    duplicate_grouping = Dict{Int,Int}()
    # final duplicates. From our example: 1 => [1,3,5]
    duplicates_dict = Dict{Int,Vector{Int}}()

    for (_v, index, is_closing) in interval_bounds
        if is_closing
            if index == current_open[1]
                popfirst!(current_open)
            else
                for (k, idx) in enumerate(current_open)
                    if idx == index
                        deleteat!(current_open, k)
                        break
                    end
                end
            end
        else
            if !isempty(current_open)
                for idx in current_open
                    # 3) Interval (idx, index) is overlapping. Check for duplicate
                    check_duplicate_candidate!(
                        duplicate_grouping,
                        duplicates_dict,
                        original_indices[idx],
                        original_indices[index],
                        certificates,
                    )
                end
            end
            push!(current_open, index)
        end
    end

    isempty(duplicates_dict) ? Vector{Int}[] : collect(values(duplicates_dict))
end

function check_duplicate_candidate!(duplicate_grouping, duplicates_dict, i, j, certificates)
    i_is_dupl = haskey(duplicate_grouping, i)
    j_is_dupl = haskey(duplicate_grouping, j)
    #If i and j are already in a duplicate cluster then there is no need to
    # check again overlapping by transitivity of our clustering
    i_is_dupl && j_is_dupl && return true
    if Bool(Arblib.overlaps(certificates[i].I′, certificates[j].I′))
        if i_is_dupl
            duplicate_grouping[j] = duplicate_grouping[i]
            push!(duplicates_dict[duplicate_grouping[i]], j)
        elseif j_is_dupl
            duplicate_grouping[i] = duplicate_grouping[j]
            push!(duplicates_dict[duplicate_grouping[j]], i)
        else
            duplicate_grouping[i] = i
            duplicate_grouping[j] = i
            duplicates_dict[i] = [i, j]
        end
        true
    else
        false
    end
end

function certify(
    F::AbstractVector{Expression},
    X,
    p = nothing;
    parameters = Variable[],
    variables = setdiff(variables(F), parameters),
    variable_ordering = variables,
    variable_groups = nothing,
    kwargs...,
)
    sys = System(
        F,
        variables = variable_ordering,
        parameters = parameters,
        variable_groups = variable_groups,
    )
    certify(sys, X, p; kwargs...)
end
function certify(
    F::AbstractVector{<:MP.AbstractPolynomial},
    X,
    p = nothing;
    parameters = similar(MP.variables(F), 0),
    variables = setdiff(MP.variables(F), parameters),
    variable_ordering = variables,
    variable_groups = nothing,
    target_parameters = nothing,
    kwargs...,
)
    if !isnothing(p)
        target_parameters = p
    end
    # as in solver_startsolutions():
    # handle special case that we have no parameters
    # to shift the coefficients of the polynomials to the parameters
    # this was the behaviour of HC.jl v1
    if isnothing(target_parameters) && isempty(parameters)
        sys, target_parameters = ModelKit.system_with_coefficents_as_params(
            F,
            variables = variable_ordering,
            variable_groups = variable_groups,
        )
    else
        sys = System(
            F,
            variables = variable_ordering,
            parameters = parameters,
            variable_groups = variable_groups,
        )
    end
    certify(sys, X, p; target_parameters = target_parameters, kwargs...)
end

function certify(F::System, args...; compile::Union{Bool,Symbol} = false, kwargs...)
    certify(fixed(F; compile = compile), args...; kwargs...)
end

function certify(
    F::AbstractSystem,
    X::MonodromyResult,
    cache::CertificationCache = CertificationCache(F);
    max_precision::Int = 256,
    kwargs...,
)
    _certify(
        F,
        solutions(X),
        certification_parameters(parameters(X)),
        cache;
        max_precision = max_precision,
        kwargs...,
    )
end

function certify(
    F::AbstractSystem,
    r::PathResult,
    p::Union{Nothing,AbstractArray} = nothing,
    cache::CertificationCache = CertificationCache(F);
    target_parameters = nothing,
    max_precision::Int = 256,
    kwargs...,
)
    cert_params =
        certification_parameters(isnothing(p) ? target_parameters : p; prec = max_precision)
    _certify(F, [solution(r)], cert_params, cache; max_precision = max_precision, kwargs...)
end

function certify(
    F::AbstractSystem,
    r::AbstractVector{PathResult},
    p::Union{Nothing,AbstractArray} = nothing,
    cache::CertificationCache = CertificationCache(F);
    target_parameters = nothing,
    max_precision::Int = 256,
    kwargs...,
)
    cert_params =
        certification_parameters(isnothing(p) ? target_parameters : p; prec = max_precision)
    _certify(F, solution.(r), cert_params, cache; max_precision = max_precision, kwargs...)
end

function certify(
    F::AbstractSystem,
    x::AbstractVector{<:Number},
    p::Union{Nothing,AbstractArray} = nothing,
    cache::CertificationCache = CertificationCache(F);
    target_parameters = nothing,
    max_precision::Int = 256,
    kwargs...,
)
    cert_params =
        certification_parameters(isnothing(p) ? target_parameters : p; prec = max_precision)
    _certify(F, [x], cert_params, cache; max_precision = max_precision, kwargs...)
end

function certify(
    F::AbstractSystem,
    X::Result,
    p::Union{Nothing,AbstractArray} = nothing,
    cache::CertificationCache = CertificationCache(F);
    target_parameters = nothing,
    max_precision::Int = 256,
    kwargs...,
)
    params = isnothing(p) ? target_parameters : p
    cert_params = certification_parameters(params; prec = max_precision)
    _certify(
        F,
        solutions(X; only_nonsingular = true),
        cert_params,
        cache;
        max_precision = max_precision,
        kwargs...,
    )
end

function certify(
    F::AbstractSystem,
    X,
    p::Union{Nothing,AbstractArray} = nothing,
    cache::CertificationCache = CertificationCache(F);
    target_parameters = nothing,
    max_precision::Int = 256,
    kwargs...,
)
    cert_params =
        certification_parameters(isnothing(p) ? target_parameters : p; prec = max_precision)
    _certify(F, X, cert_params, cache; max_precision = max_precision, kwargs...)
end


text:
abstract type OrgComponent end

include("objects.jl")
include("elements.jl")

include("cache.jl")

mutable struct Org
    settings::Dict
    contents::Vector{Union{Heading, Section}}
    cache::OrgCache
    function Org(settings::Dict, contents::Vector{<:Union{Heading, Section}})
        o = new(settings, contents, OrgCache())
        setfield!(o.cache, :doc, o)
        o
    end
end

Org(contents::Vector) = Org(Dict(), contents)
Org() = Org(Union{Heading, Section}[])

ext: jl


text:
x = 2
function f() global x=5; @assert x == 5 end
f()
@assert x == 5

max_stars_repo_name: idrougge/ple

max_issues_repo_name: idrougge/ple

max_forks_repo_name: idrougge/ple


text:
using ToeplitzMatrices
import LinearAlgebra: Symmetric, cholesky, *, mul!

import Base: transpose, adjoint, copy
import ToeplitzMatrices: Toeplitz

function copy(T::Toeplitz)
    return Toeplitz(copy(T.vc), copy(T.vr), copy(T.vcvr_dft), copy(T.tmp), T.dft)
end
function copy(T::SymmetricToeplitz)
    return SymmetricToeplitz(copy(T.vc), copy(T.vcvr_dft), copy(T.tmp), T.dft)
end

"""
    cholesky(T::SymmetricToeplitz)

Implementation adapted from "On the stability of the Bareiss and related
Toeplitz factorization algorithms", Bojanczyk et al, 1993.
"""
function cholesky!(L::AbstractMatrix, T::SymmetricToeplitz)

    # Initialize.
    L[:, 1] .= T.vc ./ sqrt(T.vc[1])
    v = copy(L[:, 1])
    N = size(T, 1)

    # Iterate.
    @inbounds for n in 1:N-1
        sinθn = v[n + 1] / L[n, n]
        cosθn = sqrt(1 - sinθn^2)

        for n′ in n+1:N
            v[n′] = (v[n′] - sinθn * L[n′ - 1, n]) / cosθn
            L[n′, n + 1] = -sinθn * v[n′] + cosθn * L[n′ - 1, n]
        end
    end
    return Cholesky(L, 'L', 0)
end
function cholesky(T::SymmetricToeplitz)
    return cholesky!(Matrix{eltype(T)}(undef, size(T, 1), size(T, 1)), T)
end

function +(T::SymmetricToeplitz, u::UniformScaling)
    Tu = copy(T)
    Tu.vc[1] += u.λ
    Tu.vcvr_dft .+= u.λ
    return Tu
end
+(u::UniformScaling, T::SymmetricToeplitz) = T + u

transpose(T::Toeplitz) = Toeplitz(T.vr, T.vc)
adjoint(T::Toeplitz) = Toeplitz(conj.(T.vr), conj.(T.vc))
transpose(T::SymmetricToeplitz) = T
adjoint(T::SymmetricToeplitz) = T

@inline LinearAlgebra.Symmetric(T::SymmetricToeplitz) = T

Toeplitz(vc::AbstractVector, vr::AbstractVector) = Toeplitz(Vector(vc), Vector(vr))

# """
#     mul!(C::Matrix, A::AbstractToeplitz, B::AbstractToeplitz)

# `O(prod(size(C)))` matrix multiplication for Toeplitz matrices. Follows from a skeleton of
# an algorithm on stackoverflow:
# https://stackoverflow.com/questions/15889521/product-of-two-toeplitz-matrices
# """
# function mul!(C::Matrix, A::AbstractToeplitz, B::AbstractToeplitz)
#     for q in 1:size(C, 2)
#         for p in 1:size(C, 1)
#             C[p, q]
#         end
#     end
#     return C
# end

# # function *(A::AbstractToeplitz, B::AbstractToeplitz)

# # end


text:
using LinearAlgebra
using DifferentialEquations
using Distributions
using Random
using Printf
using Plots
using RCall
using CPUTime
using DynamicalSystems
using StatsPlots
# using StatsPlots
include("utils.jl")
include("lorenz_ode.jl")
include("abc.jl")
include("abcMCMC.jl")
include("BayesianCalibration.jl")
include("abcSMC.jl")

# sigma, rho, beta, x0
# sigma_prior = Uniform(5,15)
# rho_prior = Uniform(20, 30)
# beta_prior = Uniform(8/3 - 1.5, 8/3 +1.5)
# x0_prior = Uniform(7, 13)

# sigma_prior = truncated(Normal(10., 5), 0, 20)
# rho_prior = truncated(Normal(26, 10), 0, 40)
# beta_prior = truncated(Normal(8/3., 2), 1, 4)

sigma_prior = Uniform(5,15)
rho_prior = Uniform(20, 30)
beta_prior = Uniform(8/3 - 1.5, 8/3 +1.5)

sd = 0.01
sigma_true = truncated(Normal(10., sd), 0, Inf)
rho_true = truncated(Normal(26., sd), 0, Inf)
beta_true = truncated(Normal(8/3., sd), 0, Inf)

function data_generator_lorenz(p)

    u0 = [10., 10., 10.]
    time_window = (0., 20.)

    y = solve_lorenz(u0,time_window,p)
end

function data_generator_lorenz10(p)

    u0 = [10., 10., 10.]
    time_window = (0., 10.)

    y = solve_lorenz(u0,time_window,p)
end

true_prior=(sigma_true, rho_true, beta_true)
true_p = map(x -> rand(x),true_prior)
y = data_generator_lorenz(true_p)

y1 = data_generator_lorenz((10., 26., 8/3))
y2 = data_generator_lorenz((10, 26, 8/3 - 1e-3))

sol1=hcat(y1[:]...)'
sol2=hcat(y2[:]...)'
l2_error = sqrt(sum(sum((sol1.-sol2).^2,dims=2)*0.01,dims=1))[1]

#
#
#
# ds_1 = Systems.lorenz([10., 10., 10.];σ= 10., ρ= 26., β= 8/3)
# ds_2 = Systems.lorenz([10., 10., 10.];σ= 10. , ρ= 26., β= 8/3 + 1)
# lyap_1 = ChaosTools.lyapunov(ds_1, 1000.0, dt=1.0, Ttr=10.0)
# lyap_2 = ChaosTools.lyapunov(ds_2, 1000.0, dt=1.0, Ttr=10.0)
# dist_pert = abs(lyap_1 - lyap_2)

#abc
algo_parameter_ABC = (prior = ( sigma_prior, rho_prior, beta_prior),epsilon = 45.,
eta= identity_mapping, d= compute_l2_norm_decay)

#  ABC MCMC
# algo_parameter_mcmc = (prior = ( sigma_prior, rho_prior, beta_prior) ,epsilon = 60., eta = identity_mapping,
#                         d= compute_l2_norm_decay, proposal = proposal_Normal,
#                         proposalRatio = proposalRatio_Normal,sd = (0.5,0.5,0.5),
#                         thinning = 100, burn_in = 100,verbose=true)

algo_parameter_mcmc = (prior = ( sigma_prior, rho_prior, beta_prior) ,epsilon = 60., eta = identity_mapping,
                        d= compute_l2_norm_decay, proposal = proposal_Normal,
                        proposalRatio = proposalRatio_Normal,sd = (0.5,0.5,0.5),
                        thinning = 20, burn_in = 100,verbose=true)

#smc
algo_parameter_smc = (prior = ( sigma_prior, rho_prior, beta_prior), time_final=8, eps_list = [80,70, 65,60,55,50,44,45],
eta= identity_mapping, d= compute_l2_norm_decay, kernel=proposal_Normal,
kernel_density=proposal_Normal_density,
sd=(0.3,0.3, 0.3), resample_method=systematic_resample, verbose=true)

output_abc = ABC(y, data_generator_lorenz, algo_parameter_ABC, 250)
output_smc = ABC_SMC(y, data_generator_lorenz, algo_parameter_smc, 250)
output_mcmc = ABC_MCMC(y, data_generator_lorenz, algo_parameter_mcmc, 250)

n=100

sigma_abc=output_abc[1][:,1]; rho_abc=output_abc[1][:,2]; beta_abc=output_abc[1][:,3]

ds_observed = Systems.lorenz([10., 10., 10.];σ= true_p[1], ρ= true_p[2], β= true_p[3])
lyap_observed = ChaosTools.lyapunov(ds_observed, 1000.0, dt=10.0, Ttr=100.0)
dist_lyapunov_abc=zeros(n)
lya_abc=zeros(n)
for j=1:n
    u0 = [10., 10., 10.]
    ds_fitted = Systems.lorenz(u0;σ= sigma_abc[j], ρ= rho_abc[j], β= beta_abc[j])
    lyap_fit = ChaosTools.lyapunov(ds_fitted, 1000.0, dt=10.0, Ttr=100.0)
    dist_lyapunov_abc[j] = abs(lyap_fit-lyap_observed)/lyap_observed
    lya_abc[j] = lyap_fit
end

sigma_mcmc=output_mcmc[1][:,1]; rho_mcmc=output_mcmc[1][:,2]; beta_mcmc=output_mcmc[1][:,3]
dist_lyapunov_mcmc=zeros(n)
lya_mcmc=zeros(n)
for j=1:n
    u0 = [10., 10., 10.]
    ds_fitted = Systems.lorenz(u0;σ= sigma_mcmc[j], ρ= rho_mcmc[j], β= beta_mcmc[j])
    lyap_fit = ChaosTools.lyapunov(ds_fitted, 1000.0, dt=1.0, Ttr=10.0)
    dist_lyapunov_mcmc[j] = abs(lyap_observed-lyap_fit)/lyap_observed
    lya_mcmc[j]=lyap_fit
end

sigma_smc=output_smc[1][:,1]; rho_smc=output_smc[1][:,2]; beta_smc=output_smc[1][:,3]
dist_lyapunov_smc=zeros(n)
lyapunov_smc=zeros(n)
for j=1:n
    u0 = [10., 10., 10.]
    ds_fitted = Systems.lorenz(u0;σ= sigma_smc[j], ρ= rho_smc[j], β= beta_smc[j])
    lyap_fit = ChaosTools.lyapunov(ds_fitted, 1000.0, dt=10.0, Ttr=10.0)
    dist_lyapunov_smc[j] = abs(lyap_fit-lyap_observed)/lyap_observed
    lyapunov_smc[j] = lyap_fit
end
using LaTeXStrings
lyap_dist_plot = plot(plot(1:n,dist_lyapunov_abc, seriestype=:scatter, markersize=3, alpha=0.5, label="ABC"),
plot(1:n,dist_lyapunov_mcmc,seriestype=:scatter, markersize=3, alpha=0.5, label="MCMC"),
plot(1:n,dist_lyapunov_smc,seriestype=:scatter,  markersize=3, alpha=0.5,label="SMC"),
ylim=(0,1.1), xlabel=L"n", ylabel=L"$|\lambda^* - \lambda^n|$",
legend=:outertopright,
 dpi=350, layout=(3,1),
 foreground_color_legend=nothing)
plot(lya_abc,zeros(n),seriestype=:scatter)
plot!([lyap_observed],[0], seriestype=:scatter, color="red")

lyap_density_plot = density(lya_abc); vline!([lyap_observed])

n= size(output_abc[1])[1]
p_abc=plot(y, vars=(1,3),w=1,legend=false,title=L"\textrm{ABC}", dpi=300)
abc_best_ind = sortperm(dist_lyapunov_abc)
abc_best_ind = randperm(n)
# abc_best_ind = 1:1:n
output_abc_best = output_abc[1][abc_best_ind,:]
for k in 1:20
    p_pred = (output_abc_best[k,1], output_abc_best[k,2], output_abc_best[k,3])
    y_fitted = data_generator_lorenz(p_pred)
    plot!(p_abc,y_fitted, vars=(1,3),alpha=0.3, color="#BBBBBB", legend=false, xlim=(-20, 20),xlabel=L"x", ylabel=L"z")
end
p_abc

# mcmc_best_ind = sortperm(dist_lyapunov_mcmc)
mcmc_best_ind =randperm(n)
output_mcmc_best = output_mcmc[1][mcmc_best_ind,:]
p_mcmc=plot(y, vars=(1,3),w=1,legend=false,title=L"\textrm{ABC-MCMC}", dpi=300)
for k in 1:20
    p_pred = (output_mcmc_best[k,1], output_mcmc_best[k,2], output_mcmc_best[k,3])
    y_fitted = data_generator_lorenz(p_pred)
    plot!(p_mcmc,y_fitted, vars=(1,3),alpha=0.3, color="#BBBBBB", legend=false, xlim=(-20, 20),xlabel=L"x", ylabel=L"z")
end
p_mcmc

smc_best_ind = sortperm(dist_lyapunov_smc)
mcmc_best_ind =randperm(n)
# smc_best_ind = 1:1:n
output_smc_best = output_smc[1][smc_best_ind,:]
p_smc=plot(y, vars=(1,3),w=1,legend=false,title=L"\textrm{ABC-SMC}",dpi=300)
for k in 1:20
    p_pred = (output_smc_best[k,1], output_smc_best[k,2], output_smc_best[k,3])
    y_fitted = data_generator_lorenz(p_pred)
    plot!(p_smc,y_fitted, vars=(1,3),alpha=0.3, color="#BBBBBB", legend=false, xlim=(-20, 20),xlabel=L"x", ylabel=L"z")
end
p_smc

sigma_density=density(output_abc[1][:,1], xlabel="σ", label="ABC")
density!(output_mcmc[1][:,1], label="MCMC")
density!(output_smc[1][:,1], label="SMC")

rho_density=density(output_abc[1][:,2], xlabel="ρ", label="ABC")
density!(output_mcmc[1][:,2], label="MCMC")
density!(output_smc[1][:,2], label="SMC")

b_density=density(output_abc[1][:,3], xlabel="b", label="ABC")
density!(output_mcmc[1][:,3], label="MCMC")
density!(output_smc[1][:,3], label="SMC")


text:
@testset "SemialgebraicSets" begin
    solver = SemialgebraicSetsHCSolver(; compile = false)
    @polyvar x y
    V = SemialgebraicSets.@set x^2 == 1 && y^2 == 2 solver
    S = collect(V)
    S = sort!(map(s -> round.(s, digits = 2), S))
    @test S == [[-1.0, -1.41], [-1.0, 1.41], [1.0, -1.41], [1.0, 1.41]]

    # Inspired from https://jump.dev/SumOfSquares.jl/v0.4.6/generated/Polynomial%20Optimization/
    function f(ε, tol)
        @testset "ε=$ε tol=$tol" begin
            solver = SemialgebraicSetsHCSolver(;
                excess_residual_tol = tol,
                real_tol = tol,
                compile = false,
            )
            @test sprint(show, solver) ==
                  "SemialgebraicSetsHCSolver(; excess_residual_tol = $tol, real_tol = $tol, compile = false)"
            o = 1 + ε
            V = SemialgebraicSets.algebraicset(
                [-x - y + o, -o * x * y + o * y^2 - y, -o * x^2 + y^2 - 2y + o],
                solver,
            )
            S = collect(V)
            S = sort!(map(s -> round.(s, digits = 2), S))
            @test length(S) == 2
            @test S[1] ≈ [0.0, 1.0] atol = tol
            @test S[2] ≈ [1.0, 0.0] atol = tol
        end
    end
    f(1e-5, 1e-4)
    f(1e-4, 1e-3)
    f(1e-3, 1e-2)
    f(1e-2, 1e-1)
end


text:
########## Parameters for scripts ##########

###### bspline-MLE ######

#### cross-timescale ####

cscl = Int(Npd*sfct)  # central scale
wscl = 10  # half scale window width, e.g. <=10
minscl, maxscl = max(1, cscl-wscl), cscl+wscl
# maxscl = Npd - 0*Npd÷4
# # minscl = maxscl ÷ 2
# minscl = maxscl - Npd÷2

wfct = 40  # factor of moving window

parms = (
    estimator = "bspline_MLE_timescale",
    method = :optim,  # method of optimizer

    # Parameters for B-Spline wavelet transform
    scalefactor = sfct,
    sclrng = 2*(minscl:maxscl),  # index range of inertial scales
    # sclrng = 2*[minscl,maxscl],  # index of working scale
    vm = 1,  # vanishing moment of the wavelet
    fmode = :causal,  # mode of filtration
    # ρ = 0,  # threshold for truncation of small eigen values in the covariance matrix
    partial = true,  # true for whittle estimator

    dfct = cscl,  # down-sampling factor
    ssize = wfct,  # size of sub window, set ssize=wsize/dfct to disable sub window
    dlen = 1,  # length of decorrelation of sub window

    windowfactor = wfct,
    wsize = cscl*wfct,  # size of moving window
    pov = Npd÷1,  # step of moving window
    boundary = :hard,  # boundary condition
)

base_estimator = x -> FracFin.bspline_MLE_estim(x, parms.sclrng, parms.vm, parms.ssize, parms.dfct, parms.dlen; partial=parms.partial)

max_stars_repo_name: yanncalec/FracFin

max_issues_repo_name: yanncalec/FracFin

max_forks_repo_name: yanncalec/FracFin


text:
#!/usr/bin/env julia

#=
 Licensing Information: You are free to use or extend these projects for
 education or reserach purposes provided that (1) you retain this notice
 and (2) you provide clear attribution to UC Berkeley, including a link
 to http://barc-project.com

 Attibution Information: The barc project ROS code-base was developed
 at UC Berkeley in the Model Predictive Control (MPC) lab by Jon Gonzales
 (jon.gonzales@berkeley.edu). The cloud services integation with ROS was developed
 by Kiet Lam  (kiet.lam@berkeley.edu). The web-server app Dator was
 based on an open source project by Bruce Wootton
=#

using RobotOS
@rosimport barc.msg: ECU, pos_info, Encoder, Ultrasound, Z_KinBkMdl#, Logging
@rosimport geometry_msgs.msg: Vector3
rostypegen()
using barc.msg
using geometry_msgs.msg
using JuMP
using Ipopt

# define model parameters
L_a     = 0.18         # distance from CoG to front axel
L_b     = 0.14         # distance from CoG to rear axel
dt      = 0.1           # time step of system
coeffCurvature   = [0,0,0,0]

# preview horizon
N       = 15

# define targets [generic values]
v_ref   = 1.0

# define objective function values
c_ey = 50
c_ev = 10
c_epsi = 50
c_df = 1
c_a = 0.36
c_ey_f = 0
c_ev_f = 0
c_epsi_f = 0

# define decision defVars
# states: position (x,y), yaw angle, and velocity
# inputs: acceleration, steering angle
println("Creating kinematic bicycle model ....")
# This turns out to be very restrictive for the Arduino at some time steps
# -> think about increasing the computation time and reducing the frequency
mdl     = Model(solver = IpoptSolver(print_level=0,max_cpu_time=0.25))

@variable( mdl, s[1:(N+1)] )
@variable( mdl, ey[1:(N+1)] )
@variable( mdl, epsi[1:(N+1)] )
@variable( mdl, 0.0 <= v[1:(N+1)] <= 3.0 )
@variable( mdl, 2.0 >= a[1:N] >= -1.0 )
@variable( mdl, -0.2 <= d_f[1:N] <= 0.2 )

# define objective function
@NLobjective(mdl, Min, sum(c_ey*ey[i]^2+c_ev*(v[i]-v_ref)^2+c_epsi*epsi[i]^2+c_df*d_f[i]^2+(a[i])^2 for i=1:N) + c_ey_f*ey[N+1]^2 + c_ev_f*(v[N+1]-v_ref)^2 + c_epsi_f*epsi[N+1]^2)

# define constraints
# define system dynamics
# Reference: R.Rajamani, Vehicle Dynamics and Control, set. Mechanical Engineering Series,
#               Spring, 2011, page 26

@NLparameter(mdl, s0      == 0); @NLconstraint(mdl, s[1]      == s0);
@NLparameter(mdl, ey0     == 0); @NLconstraint(mdl, ey[1]     == ey0);
@NLparameter(mdl, epsi0   == 0); @NLconstraint(mdl, epsi[1]   == epsi0 );
@NLparameter(mdl, v0      == 0); @NLconstraint(mdl, v[1]      == v0);
@NLparameter(mdl, coeff[i=1:length(coeffCurvature)]==coeffCurvature[i]);
@NLexpression(mdl, c[i = 1:N],    coeff[1]*s[i]^3+coeff[2]*s[i]^2+coeff[3]*s[i]+coeff[4])
@NLexpression(mdl, bta[i = 1:N],    atan( L_b / (L_a + L_b) * tan( d_f[i]) ) )

@NLexpression(mdl, dsdt[i = 1:N], v[i]*cos(epsi[i]+bta[i])/(1-ey[i]*c[i]))
for i in 1:N
    @NLconstraint(mdl, s[i+1]     == s[i]       + dt*dsdt[i]  )
    @NLconstraint(mdl, ey[i+1]    == ey[i]      + dt*v[i]*sin(epsi[i]+bta[i])  )
    @NLconstraint(mdl, epsi[i+1]  == epsi[i]    + dt*(v[i]/L_b*sin(bta[i])-dsdt[i]*c[i])  )
    @NLconstraint(mdl, v[i+1]     == v[i]       + dt*(a[i]  - 0.63 *abs(v[i])*v[i])  )
end
# restrict change of steering angle between two time steps
for i in 1:N-1
    @constraint(mdl, -0.06 <= d_f[i+1]-d_f[i] <= 0.06)
end

# status update
println("initial solve ...")
solve(mdl)
println("finished initial solve!")

# write info in logfile
# currenttime = now()
# f = open("/home/odroid/rosbag/MPClog_$currenttime.txt","w")
# byteswritten=write(f,"MPC LOG, created $currenttime\n=========================\nv_ref = $v_ref\nN = $N\nc_ey = $c_ey\nc_ev = $c_ev\nc_epsi = $c_epsi\nc_df = $c_df\nc_a = $c_a\nc_ey_f = $c_ey_f\nc_ev_f = $c_ev_f\nc_epsi_f = $c_epsi_f")
# close(f)

# if byteswritten>0
#     println("Successfully written to Logfile MPClog_$currenttime.txt")
# else
#     println("Problem with Logfile")
# end

function SE_callback(msg::pos_info)
    # update mpc initial condition
    # println("initial condition updated")
    # println(msg.coeffCurvature)
    setvalue(s0,     msg.s)
    setvalue(ey0,    msg.ey)
    setvalue(epsi0,  msg.epsi)
    setvalue(v0,     msg.v)
    setvalue(coeff,  msg.coeffCurvature)
end

function main()
    # initiate node, set up publisher / subscriber topics
    init_node("controller_mpc_traj")
    pub = Publisher{ECU}("ecu", queue_size=10)      # publishes acceleration at this time step
    pub3 = Publisher{ECU}("ecu2", queue_size=10)    # publishes velocity at next time step at rear wheels
    pub0 = Publisher{ECU}("ecu0", queue_size=10)    # publishes velocity at current time step
    # pub2 = Publisher("logging", Logging, queue_size=10)
    s1  = Subscriber{pos_info}("pos_info", SE_callback, queue_size=1)
    loop_rate = Rate(10)    # with the steering angle change bounds, we might need to adapt this (solving became harder)
    cmdcount = 0
    failcount = 0
    while ! is_shutdown()
        # run mpc, publish command
        tic()
        status = solve(mdl)
        solvetime = toc()
        if status == Symbol("Optimal")
            # get optimal solutions
            a_opt   = getvalue(a[1])
            d_f_opt = getvalue(d_f[1])
            v_opt   = getvalue(v[2])*cos(getvalue(bta[2]))
            cmd = ECU(a_opt, d_f_opt)
            cmd2 = ECU(v_opt, d_f_opt)
            cmd0 = ECU(getvalue(v[1]), d_f_opt)
            # println(d_f_opt)
            # publish commands
            if cmdcount>10      # ignore first 10 commands since MPC often stagnates during the first seconds (why?)
                publish(pub0, cmd0)
                publish(pub, cmd)
                publish(pub3, cmd2)
            end
            failcount = 0
        else
            failcount = failcount + 1
            if failcount >= 5       # if at least 5 unsolved problems in a row
                cmd = ECU(0,0)      # stop car
                publish(pub0, cmd)
                publish(pub, cmd)
                publish(pub3, cmd)
                println("Car stopped because failcount hit.")
            end
        end
        println("Solve Status: ", string(status), "\nSolve Time: ", solvetime,"\n\n")
        # loginfo = Logging(getobjectivevalue(mdl),string(status),solvetime,getvalue(s),getvalue(ey),getvalue(epsi),getvalue(v),getvalue(a),getvalue(d_f))
        # publish(pub2, loginfo)
        cmdcount = cmdcount + 1
        rossleep(loop_rate)
    end
end

if ! isinteractive()
    main()
end

max_stars_repo_name: RoyaFiroozi/barc

max_issues_repo_name: RoyaFiroozi/barc

max_forks_repo_name: RoyaFiroozi/barc


text:
using VLOptionsModeling

# create an emprt contract set for now -
contract_set = Set{VLAbstractAsset}()
assetSymbol = "XYZ"
strikePrice = 50.0
expirationDate = nothing
premimumValue = 0.0
numberOfContracts = 1
sense = :buy
contractMultiplier = 100.0
put_contract = VLPutOptionContract(assetSymbol, strikePrice, premimumValue; sense=sense, contractMultiplier=1.0)
push!(contract_set, put_contract)

# setup the lattice (Example 19.1/Fig. 19.3 Hull)
𝝙t = 0.0833 # years
μ = log(1.0084) * (1 / 𝝙t)
σ = (1 / √𝝙t) * log(1.1224)
numberOfLevels = 6
underlyingPrice = 50.0

# Create the lattice model -
latticeModel = VLBinomialLattice(μ, σ, 𝝙t, numberOfLevels)

# compute the price -
price_tree = binomial_price(contract_set, latticeModel, underlyingPrice) |> check


text:
using DoubleFloats
using LinearAlgebra

using SemiLagrangian:
    Advection,
    standardsplit,
    strangsplit,
    triplejumpsplit,
    sizeall,
    AdvectionData,
    getdata,
    advection!,
    UniformMesh,
    gettranslationvar,
    AbstractInterpolation,
    interpolate!,
    Lagrange,
    B_SplineLU,
    B_SplineFFT,
    Hermite
"""

   exact(tf, mesh1, mesh2)

   Julia function to compute exact solution "

```math
\frac{d f}{dt} +  (y \frac{df}{delta1} - x \frac{df}{delta2}) = 0
```

"""
function exact!(f, decall::NTuple{Nsum,T}, t::T) where {T,Nsum}
    for ind in CartesianIndices(f)
        f[ind] = exp(
            -sum([
                sin(2T(pi) * (ind.I[i] + decall[i] * t) / size(f, i) + T(big"0.56") * i^2) for i = 1:Nsum
            ]),
        )
    end
end

function test_translation(
    sz::NTuple{2,Int},
    interp_sp::AbstractInterpolation{T},
    interp_v::AbstractInterpolation{T},
    dt::T;
    tc::Vector{T}=strangsplit(dt),
) where {T}
@time begin
    spmin, spmax, nsp = T(-5), T(5), sz[1]
    vmin, vmax, nv = -T(5.55), T(5), sz[2]

    mesh_sp = UniformMesh(spmin, spmax, nsp)
    mesh_v = UniformMesh(vmin, vmax, nv)
end
println("trace1")
@time begin
    nbdt = Int(round(T(1)/dt))
    adv = Advection(
        (mesh_sp, mesh_v),
        [interp_sp, interp_v],
        dt,
        [([1, 2], 1, 1, true), ([2, 1], 1, 2, true)],
        tab_coef= tc
    )

    tabref = zeros(T, sz)

    v1 = T(big"0.83545655467782872872782870029282982828737872878776717190927267611111")
    v2 = T(big"0.945678101929276765616767176761671771766717828781828998101092981877817176")
end
println("trace2")
@time begin

    exact!(tabref, (v1, v2), T(0))

    pvar = gettranslationvar((v1, v2))

    advdata = AdvectionData(adv, tabref, pvar)

    diffmax = 0
    data = getdata(advdata)
end
println("trace3")
    for ind = 1:nbdt
        data = getdata(advdata)
        while advection!(advdata)
        end
        #        GC.gc()
        exact!(tabref, (v1, v2), T(ind)*dt)
        data = getdata(advdata)
        diff = norm(data .- tabref, Inf)
        @show v1, v2, ind, diff
        diffmax = max(diffmax, diff)
    end
    println(
        "length(tab_coef)=$(length(tc)) test_translation sz=$sz interp=$interp_sp, $interp_v nbdt=$nbdt diffmax=$diffmax",
    )
    return diffmax

end
function test_translation(
    sz::NTuple{N,Int},
    interp_t::Vector{I},
    nbdt::Int,
) where {T,N,I<:AbstractInterpolation{T}}

    dt = T(1)

    vall = ntuple(x -> T(10rand(BigFloat) - 5), N)

    vall = T(1) / 10 .* vall

    @show vall

    tabref = zeros(T, sz)
    exact!(tabref, vall, T(0))

    data = copy(tabref)

    buf = zeros(T, sz)

    # dec1 = fill(v1, sz)
    # dec2 = fill(v2, sz)

    diffmax = 0
    for ind = 1:nbdt
        interpolate!(buf, data, (i, _) -> vall[i], interp_t)
        #        interpolate!(buf, data, (dec1, dec2), interp)
        copyto!(data, buf)
        exact!(tabref, vall, T(ind))
        diff = norm(data .- tabref, Inf)
        diffmax = max(diffmax, diff)
        @show ind, diff
    end
    println("test_translation sz=$sz interp=$interp_t nbdt=$nbdt diffmax=$diffmax")
    return diffmax

end

@testset "test translation" begin
    T = Float64
    dt = T(1)/11
    @time @test test_translation(
        (200, 300),
        Lagrange(5, T),
        Lagrange(5, T),
        dt,
    ) < 1e-3
    @time @test test_translation(
        (200, 300),
        Lagrange(5, T),
        Lagrange(5, T),
        dt, tc=triplejumpsplit(dt)
    ) < 1e-3
    @time @test test_translation(
        (200, 300),
        Lagrange(5, T),
        Lagrange(5, T),
        dt, tc=standardsplit(dt)
    ) < 1e-3
    @time @test test_translation(
        (200, 300),
        Hermite(5, T),
        Hermite(5, T),
        dt,
    ) < 1e-3
    @time @test test_translation(
        (200, 300),
        Hermite(5, T),
        Hermite(5, T),
        dt, tc=triplejumpsplit(dt)
    ) < 1e-3
    @time @test test_translation(
        (200, 300),
        Hermite(5, T),
        Hermite(5, T),
        dt,tc=standardsplit(dt)
    ) < 1e-3
    @time @test test_translation(
        (128, 64),
        B_SplineLU(5, 128, T),
        B_SplineLU(5, 64, T),
        dt,
    ) < 1e-6
    @time @test test_translation(
        (128, 64),
        B_SplineLU(5, 128, T),
        B_SplineLU(5, 64, T),
        dt, tc=triplejumpsplit(dt)
    ) < 1e-6
    @time @test test_translation(
        (128, 64),
        B_SplineFFT(5, 128, T),
        B_SplineFFT(5, 64, T),
        dt,
    ) < 1e-6
    T = Double64
    dt = T(1)/11
    @time @test test_translation((100, 120), Lagrange(15, T), Lagrange(15, T), dt) < 1e-14
    @time @test test_translation(
        (128, 64),
        B_SplineLU(15, 128, T),
        B_SplineLU(15, 64, T),
        dt,
    ) < 1e-18
    @time @test test_translation(
        (128, 64),
        B_SplineFFT(15, 128, T),
        B_SplineFFT(15, 64, T),
        dt,
    ) < 1e-17

    T = BigFloat
    dt = T(1)/11
@time @test test_translation((100, 90), Lagrange(25, T), Lagrange(25, T), dt) < 1e-20
@time @test test_translation(
        (128, 64),
        B_SplineLU(25, 128, T),
        B_SplineLU(25, 64, T),
        dt,
    ) < 1e-22
   @time @test test_translation(
        (128, 64),
        B_SplineFFT(25, 128, T),
        B_SplineFFT(25, 64, T),
        dt,
    ) < 1e-22

end


text:
@safetestset "Control flow rendering" begin
  @safetestset "IF conditional rendering" begin
    @safetestset "IF true" begin
      using Genie
      using Genie.Renderer.Html

      view = raw"""

                  <% if true; [ %>

Hello

Welcome

                  <% ]end %>

"""

      @test String(html(view).body) == raw"

Hello

Welcome

"
    end;

    @safetestset "IF false" begin
      using Genie
      using Genie.Renderer.Html

      view = raw"""

                  <% if false; [ %>

Hello

Welcome

                  <% ]end %>

"""

      @test String(html(view).body) == raw"""

"""
    end;
  end;
end;


text:
using Test
import PortfolioMngmt

@test cmp(PortfolioMngmt.test_me(), "PortfolioMngmt") == 0


text:
###
#   An implementation of Transitional Markov Chain Monte Carlo in Julia
#
#            Institute for Risk and Uncertainty, Uni of Liverpool
#
#                       Authors: Ander Gray, Adolphus Lye
#
#                       Email: Ander.Gray@liverpool.ac.uk,
#                              Adolphus.Lye@liverpool.ac.uk
#
#
#   This Transitional MCMC algorithm is inspirted by OpenCOSSAN:
#                   https://github.com/cossan-working-group/OpenCossan
#
#
#   Algorithm originally proposed by:
#           J. Ching, and Y. Chen (2007). Transitional Markov Chain Monte Carlo method
#           for Bayesian model updating, Model class selection, and Model averaging.
#           Journal of Engineering Mechanics, 133(7), 816-832.
#           doi:10.1061/(asce)0733-9399(2007)133:7(816)
#
###

using TimerOutputs
const to = TimerOutput()

function tmcmc(
    log_fD_T::Function,
    log_fT::Function,
    sample_fT::Function,
    Nsamples::Integer,
    burnin::Integer=20,
    thin::Integer=3,
    beta2::Float64=0.01)

    j1 = 0;                     # Iteration number
    βj = 0;                     # Tempering parameter
    θ_j = sample_fT(Nsamples);  # Samples of prior
    Lp_j = zeros(Nsamples, 1);   # Log liklihood of first iteration

    Log_ev = 0                  # Log Evidence

    Ndims = size(θ_j, 2)         # Number of dimensions (input)

    covariance_method = LinearShrinkage(DiagonalUnitVariance(), :lw) # * DiagonalUnitVariance is always SPD!

    @timeit_debug to "Main while loop" while βj < 1

        j1 = j1 + 1

        @debug "Beginning iteration $j1"

        ###
        # Parallel evaluation of the likelihood
        ###

        @debug "Computing likelihood with $(nworkers()) workers..."

        Lp_j = pmap(log_fD_T, eachrow(θ_j))
        Lp_j = reduce(vcat, Lp_j)
        Lp_adjust = maximum(Lp_j)

        ###
        #   Computing new βj
        ###
        @debug "Computing β_j and weights..."

        @timeit_debug to "Compute β and weights " βj1, w_j = _beta_and_weights(βj, Lp_j .- Lp_adjust)

        @info "β_$(j1) = $(βj1)"

        @timeit_debug to "Log evidence" Log_ev = log(mean(w_j)) + (βj1 - βj) * Lp_adjust + Log_ev   # Log evidence in current iteration

        prop = mu -> proprnd(mu, Σ_j, log_fT) # Anonymous function for proposal
        target = x -> log_fD_T(x) .* βj1 .+ log_fT(x) # Anonymous function for transitional distribution

        # (Normalized) Weighted resampling of θj (indices with replacement)
        @timeit_debug to "Normalised weights" wn_j = w_j ./ sum(w_j); # normalize weights
        @timeit_debug to "Compute indices" indices = sample(1:Nsamples, Weights(wn_j), Nsamples, replace=true)
        @timeit_debug to "Update θ_j1" θ_j1 = θ_j[indices, :]

        # Estimate covariance matrix
        @timeit_debug to "Compute covariance" Σ_j = beta2 * cov(covariance_method, θ_j1)

        @debug "Markov chains with $(nworkers()) workers..."

        # pmap is a parallel map
        @timeit_debug to "Run markov chains" θ_j1 = pmap(x -> metropolis_hastings_simple(target, prop, x, 1, burnin, thin)[1], eachrow(θ_j1))
        θ_j1 = reduce(vcat, θ_j1)

        βj = βj1
        θ_j = θ_j1
    end
    return θ_j, Log_ev
end

"""
    _beta_and_weights(β, likelihood)

Compute the next value for `β` and the nominal weights `w` using bisection.
"""
function _beta_and_weights(β::Real, adjusted_likelihood::AbstractVector{<:Real})
    low = β
    high = 2

    local x, w # Declare variables so they are visible outside the loop

    while (high - low) / ((high + low) / 2) > 1e-6
        x = (high + low) / 2
        w = exp.((x .- β) .* adjusted_likelihood)

        if std(w) / mean(w) > 1
            high = x
        else
            low = x
        end
    end

    return min(1, x), w
end

function proprnd(mu::AbstractVector, covMat::AbstractMatrix, prior::Function)
    samp = rand(MvNormal(mu, covMat), 1)
    while isinf(prior(samp))
        samp = rand(MvNormal(mu, covMat), 1)
    end
    return samp[:]
end


text:
module Special
using  ...PatternStructure.Trees
using  ...PatternStructure.Checks
using  ...PatternStructure.SlurpTypes
import Base.Meta.isexpr
export special_heads, special_shortcuts, is_slurp, is_special,
       slurptype, patterntype, special_name, is_special_expr,
       is_binding_name

const special_heads =
  [:literal, :autobinding, :binding, :type, :consistent, :predicate, :equals, :iterable, :raw]

const special_shortcuts =
  Dict(:(:A) => :(:autobinding),
       :(:L) => :(:literal),
       :(:B) => :(:binding),
       :(:T) => :(:type),
       :(:C) => :(:consistent),
       :(:P) => :(:predicate),
       :(:I) => :(:iterable),
       :(:E) => :(:equals),
       :(:EQ)=> :(:equals),
       :(:R) => :(:raw)
       )

is_slurp(node::PatternNode) = isa(node.head, SlurpHead)
is_slurp(x) = false

is_special(node::PatternNode) = node.head in special_heads || is_slurp(node)

slurptype(node::PatternNode) =
  isa(node.head, LazySlurp)   ? :? :
  isa(node.head, GreedySlurp) ? :* :
  throw(ArgumentError(
    "Couldn't determine the slurp type of the pattern node $(node.head)"))

patterntype(node::PatternNode) =
  is_slurp(node) ? slurptype(node) : node.head

extract_name(x::QuoteNode) = x.value
extract_name(x) = x

special_name(x)       = extract_name(get(special_shortcuts, x, x))
special_name(x::Expr) = special_name(QuoteNode(x.args[1]))

function is_binding_name(x::Symbol)
  firstchar = string(x)[1]

  isalpha(firstchar) || firstchar in ['#', '@', '_']
end
is_binding_name(x) = false

function is_special_expr(ex::Expr)
  if ex.head == :curly && ex.args[1] == :?
    warn("?{...} is deprecated, use :?{...}")
    ex.args[1] = :(:?)
  end

  ex.head == :curly &&
    (isa(ex.args[1], QuoteNode) ||
    (isexpr(ex.args[1], :quote) ||
    (ex.args[1] in [:(:?), :(*)])))
end

end


text:
export simple_eval, getfield_vec, abserror, relerror

abserror(sol, refsol) = abs(sol-refsol)
relerror(sol, refsol) = abs((sol-refsol)./refsol)
abserror(tr::TestResults) = abserror(tr.solution[tr.testrun.tc.refsolinds],
                                     tr.testrun.tc.refsol)
relerror(tr::TestResults) = relerror(tr.solution[tr.testrun.tc.refsolinds],
                                     tr.testrun.tc.refsol)

@doc """Calculates the scd error estimate from Mazzia & Magherini p.II-ii.
        It is the minimum significant correct digits in the solution.
        (Note this is does not work if any component of refsol==0)""" ->
calc_error_scd(sol, refsol) = -log10(norm(relerror(sol,refsol), Inf))

@doc """Same as two argument call but use absolute error instead
        of relative error at absinds""" ->
function calc_error_scd(sol, refsol, absinds)
    out = relerror(sol,refsol)
    out[absinds] = abserror(sol[absinds],refsol[absinds])
    return -log10(norm(out, Inf))
end
calc_error_scd(sol, tr::TestRun) = calc_error_scd(sol[tr.tc.refsolinds],
                                                 tr.tc.refsol,
                                                 tr.tc.scd_absinds)
calc_error_scd(res::TestResults) = calc_error_scd(res.yend[res.testrun.tc.refsolinds],
                                                 res.testrun.tc.refsol,
                                                 res.testrun.tc.scd_absinds)

@doc """Calculates the mescd error estimate from Mazzia & Magherini p.II-ii.
        Mixed significant digit of the approximate solution.  (Very similar to
        the scd error estimate, except when refsol is close to zero?)
        """ ->
function calc_error_mescd(sol, refsol, abstol, reltol)
    -log10( norm( abserror(sol,refsol)./(abstol./reltol + refsol) ,Inf))
end
calc_error_mescd(sol, tr::TestRun) = calc_error_mescd(sol[tr.tc.refsolinds],
                                                      tr.tc.refsol,
                                                      tr.abstol, tr.reltol)
calc_error_mescd(res::TestResults) = calc_error_mescd(res.yend[res.testrun.tc.refsolinds],
                                                     res.testrun.tc.refsol,
                                                     res.testrun.abstol, res.testrun.reltol)

function simple_eval{N,T}(tend::T, yend::Vector{T}, stats, walltime, mem, gc_time, tr::TestRun{N,T})
    tc = tr.tc
    if tend!=tc.tspan[2]
        st = "Integration did not run to specified end time.  tend=$tend but should be $(tc.tspan[2])."
        warn(st)
        e = ErrorException(st)
        return simple_eval(e, tr)
    end
    scd = calc_error_scd(yend, tr)
    #scd = scd<0 ? NaN : scd # return NaN if significant digits less than 0
    if isadaptive(tr)
        mescd = calc_error_mescd(yend, tr)
    else
        mescd = calc_error_mescd(yend[tc.refsolinds], tc.refsol, 1, 1) #
    end
    #mescd = mescd<0 ? NaN : scd # return NaN if significant digits less than 0

    return TestResults(tr, tend, yend, stats..., walltime, mem, gc_time, scd, mescd, nothing)
end
function simple_eval{N,T}(e::Exception, tr::TestRun{N,T})
    return TestResults(tr, tr.tc.tspan[1], T[], -1, -1, -1, -1, -1,     NaN,  -1,    NaN, NaN, NaN, e)
end

@doc """For getting metrics (or any field out of) out of an Vector filled with
        one type having fields""" ->
function getfield_vec(trs::Vector, field::Symbol)
    tr = trs[1]
    out = zeros(typeof(tr.(field)), length(trs))
    for (i,tr) in enumerate(trs)
        out[i] = tr.(field)
    end
    out
end


text:
"""
    HighOrderCBF

High Order Control barrier function for a control affine system.
# Fields
- `CBF`: CBF used to derive HOCBF
- `ψ`: highest order derivative of CBF
- `∇ψ`: gradient of ψ
"""
struct HighOrderCBF <: BarrierFunction
    ψ
    ∇ψ
    α
end

"""
	HighOrderCBF(CBF::ControlBarrierFunction, Σ::ControlAffineSystem, degree::Int)

Construct a High Order Control Barrier Function of a specified degree.
"""
function HighOrderCBF(CBF::ControlBarrierFunction, Σ::ControlAffineSystem, degree::Int)
    if degree == 2
        return second_order_hocbf(CBF, Σ)
    end
end

"""
	(HOCBF::HighOrderCBF)(x)

Evaluate HOCBF at state x.
"""
function (HOCBF::HighOrderCBF)(x)
    return HOCBF.ψ(x)
end

"""
	second_order_hocbf(CBF::ControlBarrierFunction, Σ::ControlAffineSystem)

Construct a HOCBF of relative degree 2.
"""
function second_order_hocbf(CBF::ControlBarrierFunction, Σ::ControlAffineSystem)
    # Need to do these autodiffs separately for now...maybe we can automate this later.
    dfdx(x) = jacobian(Σ.f, x)[1]
    dh2dx2(x) = hessian(CBF.h, x)
    dα(s) = gradient(CBF.α, s)[1]
    ψ1(x) = CBF.∇h(x) * Σ.f(x) + CBF.α(CBF(x))
    ∇ψ1(x) = CBF.∇h(x) * dfdx(x) + (dh2dx2(x) * Σ.f(x))' + dα(CBF(x)) * CBF.∇h(x)

    return HighOrderCBF(ψ1, ∇ψ1, CBF.α)
end


text:
# Hacks in a new Gadfly GuideElement that behaves like an Annotation but
# sits under (behind) the main plot rather than over it.

immutable Underguide <: Gadfly.GuideElement
    ctx::Context
end

const underguide = Underguide

import Gadfly.render

function render(guide::Underguide, theme::Gadfly.Theme,
                aes::Gadfly.Aesthetics)
    ctx = compose(context(), svgclass("geometry"), guide.ctx)
    return [Gadfly.Guide.PositionedGuide([ctx], 0, Gadfly.Guide.under_guide_position)]
end


text:
module SimpleBufferStream

include("BufferStream.jl")

end # module SimpleBufferStream

max_stars_repo_name: staticfloat/SimpleBufferStream

max_issues_repo_name: staticfloat/SimpleBufferStream

max_forks_repo_name: staticfloat/SimpleBufferStream


text:
#
println("Ec) Three-step cascade computations and simulations for the decay of the neon 1s^-1 3p hole states: AverageSCA model.")
using JLD
#
setDefaults("method: continuum, asymptotic Coulomb")    ## setDefaults("method: continuum, Galerkin")
setDefaults("method: normalization, pure sine")         ## setDefaults("method: normalization, pure Coulomb")    setDefaults("method: normalization, pure sine")
#
if false

elseif  false    ## Stepwise decay cascade
setDefaults("print summary: open", "zzz-old-Cascade-computation-following-decay.sum")

name = "Cascade after neon 1s --> 3p excitation"
grid = Radial.Grid(false)
wa   = Cascade.Computation(Cascade.Computation(); name=name, nuclearModel=Nuclear.Model(10.), grid=grid, approach=Cascade.AverageSCA(),
                           scheme=Cascade.StepwiseDecayScheme([Auger(), Radiative()], 1, Dict{Int64,Float64}(), 0, Shell[], Shell[]),
                           initialConfigs=[Configuration("1s^1 2s^2 2p^6 3p")] )
println(wa)
wb = perform(wa; output=true)
setDefaults("print summary: close", "")

else
setDefaults("print summary: open", "zzz-Cascade-simulation.sum")

data = [JLD.load("zzz-cascade-decay-computations-2020-07-13T08.jld")]
## data = [wb]
name = "Simulation of the neon 1s^-1 3p decay"

wc   = Cascade.Simulation(Cascade.Simulation(), name=name, property=Cascade.IonDistribution(),
                          settings=Cascade.SimulationSettings(0., 0., 0., 0., 0., [(1, 2.0), (2, 1.0), (3, 0.5)]), computationData=data )
println(wc)
wd = perform(wc; output=true)
setDefaults("print summary: close", "")

end


text:
######
# This file is part of the MomentArithmetic.jl package, for performing arithmetic operations
# between moments of uncertain numbers, when the moments and the uncertain numbers are only
# partially known.
#
#   Definition of Moment variable class
#
#           University of Liverpool,
#           Institute for Risk and Unertainty
#
#                                           Authors: Ander Gray, Scott Ferson
#                                           Email:  ander.gray@liverpool.ac.uk
#
#
#           (look in downs)
#
#               -> Doesn't exits: bound and a variance
#
#               -> From moments, we can't get the range... unless we know something about the shapes
#               -> However possible option to truncate pba style. Give option. (Consitent throughout)
#
#               -> Arithmetic when dependency is known (?) Perf, op, correlation. Interval correlation coeff.
#               -> Dependeny. Calulate the covarience between Z and X. when z = x + y. Langewisch paper has some
#
#
######

#global correlationDict = Dict()
#global unitDict = Dict()

#interval(x :: AbstractInterval, y :: AbstractInterval) = interval(x.lo, y.hi);

mutable struct Moments <: AbstractMoment

    mean ::  Union{Real, Interval}                   # Can also make it a pbox as an input
    var  ::  Union{Real, Interval}
    range :: Union{Missing, Interval}

    #CorreltedVars :: Union{Missing, Array{UInt64,1}} # For defining which variables this one is correlted too, for dependency tracking
    #Unit ::  Symbol   # Will add this later

    function Moments( mean = missing,  var = missing, range = missing)

        mean, var, range = make_consistent(mean, var, range)

        return new(mean, var, range)

    end
end

Moments(;mean = missing, var =missing, range =missing) = Moments(mean, var , range)

intervalM =  Moments

function make_consistent(mean = missing,  var = missing, range = missing)

    if ismissing(mean);     mean    = interval(-Inf,Inf);   end
    if ismissing(var);      var     = interval(0,Inf);      end
    if ismissing(range);    range   = interval(-Inf,Inf);   end

    # mean constraint
    if (mean ∩ range == ∅); throw(ArgumentError("Provided information not valid. Mean ∩ Range = ∅.\n       $mean ∩ $range = ∅")); end
    if !( mean ⊆ range);
        ml = max(left(mean), left(range));
        mh = min(right(mean), right(range));
        mean = interval(ml,mh);
    end

    # Var constraint
    maxVar = interval(0,Inf)
    if isfinite(range.hi) && isfinite(range.lo)
        maxRange = range.hi; minRange = range.lo;
        meanMax = right(mean); meanMin = left(mean);

        v1 = (maxRange-minRange)*(maxRange-meanMin)-(maxRange-meanMin)*(maxRange-meanMin);
        v2 = (maxRange-minRange)*(maxRange-meanMax)-(maxRange-meanMax)*(maxRange-meanMax);
        v3 = 0; mid = (range.hi + range.lo)/2;

        if (mid ∈ mean); v3 =  (maxRange-minRange)*(maxRange-mid)-(maxRange-mid)*(maxRange-mid); end

        #println("v1: $v1");println("v2: $v2");println("v3: $v3")

        vh = max(v1, v2, v3); vl = 0;
        maxVar = interval(vl,vh);
    end

    if (var ∩ maxVar == ∅); throw(ArgumentError("Provided information not valid. Variance ∩ VarBounds = ∅.\n       $var ∩ $maxVar = ∅")); end
    if !(var ⊆ maxVar)
        vl = max(left(var), left(maxVar));
        vh = min(right(var), right(maxVar));
        var = interval(vl, vh)
    end

    return mean, var, range

end

function makepbox(x :: Moments)

    if isvacuous(x.range); return meanVar(x.mean, x.var); end

    return MinMaxMeanVar(x.range.lo, x.range.hi, x.mean, x.var)

end

plot(x::Moments) = plot(makepbox(x))
plot(x::Moments, c...) = plot(makepbox(x), c...)

function Base.show(io::IO, z::Moments)

    if (typeof(z.mean) <: AbstractInterval); statement1 = "mean = [$(z.mean.lo),$(z.mean.hi)]"; else statement1= "mean  = $(z.mean)";end
    if (typeof(z.var) <: AbstractInterval); statement2 = "var = [$(z.var.lo),$(z.var.hi)]"; else statement2 = "var = $(z.var)";end
    if ( !ismissing(z.range) ); statement3 = ", range = [$(z.range.lo),$(z.range.hi)]"; else; statement3 = "";end

    print(io, "moment: \t  ~ ( $statement1, $statement2 $statement3 )");

end


text:
function edit{T <: Union(AbstractFixedVector, Real)}(style::Style{:Default}, numbertex::Texture{T, 1, 2}, customization::Dict{Symbol,Any})

  backgroundcolor = customization[:backgroundcolor]
  screen          = customization[:screen]
  camera          = screen.orthographiccam
  color           = customization[:color]
  maxdigits       = customization[:maxdigits]
  maxlength       = customization[:maxlength]
  model           = customization[:model]
  font            = customization[:font]

  numbers         = data(numbertex) # get data from texture/video memory
  text            = Array(GLGlyph{Uint16}, int(size(numbers,1)*maxdigits), size(numbers, 2))

  fill!(text, GLGlyph()) # Fill text array with blanks, as we don't need all of them

  # handle real values
  Base.stride(x::Real, i)       = 1
  # remove f0
  makestring(x::Integer)        = string(int(x))
  makestring(x::FloatingPoint)  = string(float64(x))
  makestring{T}(x::Vector1{T})  = string(float64(x[1]))
  makestring{T}(x::Vector1{T}, maxlen) = makestring(float64(x[1]), maxlen)

  makestring(x::Integer, maxlen) = begin
    tmp = string(int(x))
    len = length(tmp)
    if len > maxlen
        tmp = tmp[1:maxlen]
    elseif len < maxlen
      tmp = rpad(tmp, maxlen, " ")
    end
    tmp
  end
  makestring(x::FloatingPoint, maxlen) = begin
    tmp = string(float64(x))
    len = length(tmp)
    if len > maxlen
        tmp = tmp[1:maxlen]
    elseif len < maxlen
      tmp = rpad(tmp, maxlen, "0")
    end
    tmp
  end
  textgpu     = Texture(text, keepinram=true)
  customization[:style_group] = Texture([color])
  customization[:textlength]  = length(textgpu)
  obj         = visualize(style, textgpu, customization)

  startposition   = Uint16[0,0]
  positionrunner  = startposition

  maxlength = 5
  for i=1:size(numbers,1)
    for j=1:size(numbers, 2)
      number = numbers[i,j]
      i3 = ((i-1)*maxlength) + 1
      textgpu[i3:i3+maxlength-1, j:j] = GLGlyph{Uint16}[GLGlyph(c, positionrunner[1], positionrunner[2]+k-1, 0) for (k,c) in enumerate(makestring(number, maxlength))]
      positionrunner += [0, maxlength+1]
    end
    positionrunner = startposition + [i,0]
  end
  selectiondata = lift(first, SELECTION[:mouse_hover])
  # We allocated more space on the gpu then needed (length(numbers)*maxdigits)
  # So we need to update the render method, to render only length(numbers) * maxlength
  obj[:postrender, renderinstanced] = (obj.vertexarray, length(textgpu))
  # We allocated more space on the gpu then needed (length(numbers)*maxdigits)
  # So we need to update the render method, to render only length(numbers) * maxlength
  # ([numbers], zero(eltype(numbers)), -1, -1, -1, Vector2(0.0))
  number_lift = foldl(([numbers;], zero(eltype(numbers)), -1, -1, -1, Vector2(0.0)), screen.inputs[:mouseposition], screen.inputs[:mousebuttonspressed], selectiondata) do v0, mposition, mbuttons, selection
    numbers0, value0, inumbers0, igpu0, mbutton0, mposition0 = v0
    # if over a number           && nothing selected &&         only           left mousebutton clicked
    if selection[1] == obj.id && inumbers0 == -1 && length(mbuttons) == 1 && in(0, mbuttons)
      iorigin   = selection[2]
      inumbers  = div(iorigin, maxlength) + 1
      igpu      = int((iorigin - (iorigin%maxlength)) + 1) # get the first index of the number
      return (numbers0, numbers0[inumbers], inumbers, igpu, 0, mposition)
    end
    # if a number is selected && previous click was left && still only left button ist clicked
    if inumbers0 > 0 && mbutton0 == 0 && length(mbuttons) == 1 && in(0, mbuttons)
      xdiff                    = mposition[1] - mposition0[1]
      numbers0[inumbers0]      = value0 + (float32(xdiff)/ 50.0f0)
      numbertex[inumbers0]     = numbers0[inumbers0]
      a = mod(inumbers0-1, size(numbers0, 1))
      b = div(inumbers0-1, size(numbers0, 1))
      for (k,c) in enumerate(makestring(numbers0[inumbers0], maxlength))
        textgpu[igpu0+k-1] = GLGlyph(c, a, (b-1)+b*maxlength+k, 0)
      end
      return (numbers0, value0, inumbers0, igpu0, 0, mposition0)
    end
    return (numbers0, zero(eltype(numbers0)), -1, -1, -1, Vector2(0.0))
  end

  obj, numbertex
end

# remove f0
makestring(x::Integer)        = string(int(x))
makestring(x::FloatingPoint)  = string(float64(x))
makestring{T}(x::Vector1{T})  = string(float64(x[1]))
makestring{T}(x::Vector1{T}, maxlen) = makestring(float64(x[1]), maxlen)
makestring(x::Integer, maxlen) = begin
  tmp = string(int(x))
  len = length(tmp)
  if len > maxlen
      tmp = tmp[1:maxlen]
  elseif len < maxlen
    tmp = rpad(tmp, maxlen, " ")
  end
  tmp
end
makestring(x::FloatingPoint, maxlen) = begin
  tmp = string(float64(x))
  len = length(tmp)
  if len > maxlen
      tmp = tmp[1:maxlen]
  elseif len < maxlen
    tmp = rpad(tmp, maxlen, "0")
  end
  tmp
end

function edit{T <: Union(AbstractVector, Real)}(style::Style, numbers::T, customization::Dict{Symbol,Any})

  backgroundcolor = customization[:backgroundcolor]
  screen          = customization[:screen]
  camera          = screen.orthographiccam
  color           = customization[:color]
  maxdigits       = customization[:maxdigits]
  maxlength       = customization[:maxlength]
  model           = customization[:model]
  font            = customization[:font]

  if ndims(numbers) == 2
    xn, yn = size(numbers)
  elseif ndims(numbers) == 1 || ndims(numbers) == 0
    xn, yn = length(numbers), 1
  else
    error("Dimension Missmatch. Wanted: 0D,1D, 2D Vector. Got: ", typeof(numbers), " ndims: ", ndims(numbers))
  end

  text        = Array(GLGlyph{Uint16}, int(xn*maxdigits), yn)
  fill!(text, GLGlyph()) # Fill text array with blanks, as we don't need all of them

  textgpu     = Texture(text)
  customization[:style_group] = Texture([color])
  customization[:textlength]  = length(textgpu)
  obj         = visualize(style, textgpu, customization)

  startposition   = Uint16[0,0]
  positionrunner  = startposition

  maxlength = 5
  for i=1:xn
    for j=1:yn
      number = numbers[i,j]
      i3 = ((i-1)*maxlength) + 1
      textgpu[i3:i3+maxlength-1, j:j] = GLGlyph{Uint16}[GLGlyph(c, positionrunner[1], positionrunner[2]+k-1, 0) for (k,c) in enumerate(makestring(number, maxlength))]
      positionrunner += [0, maxlength+1]
    end
    positionrunner = startposition + [i,0]
  end
  selectiondata = lift(first, SELECTION[:mouse_hover])
  # We allocated more space on the gpu then needed (length(numbers)*maxdigits)
  # So we need to update the render method, to render only length(numbers) * maxlength
  obj[:postrender, renderinstanced] = (obj.vertexarray, length(textgpu))
  # We allocated more space on the gpu then needed (length(numbers)*maxdigits)
  # So we need to update the render method, to render only length(numbers) * maxlength
  # ([numbers], zero(eltype(numbers)), -1, -1, -1, Vector2(0.0))
  number_lift = foldl(([numbers;], zero(eltype(numbers)), -1, -1, -1, Vector2(0.0)), screen.inputs[:mouseposition], screen.inputs[:mousebuttonspressed], selectiondata) do v0, mposition, mbuttons, selection
    numbers0, value0, inumbers0::Int, igpu0, mbutton0, mposition0 = v0
    # if over a number           && nothing selected &&         only           left mousebutton clicked
    if selection[1] == obj.id && inumbers0::Int == -1 && length(mbuttons) == 1 && in(0, mbuttons)
      iorigin   = selection[2]
      inumbers  = div(iorigin, maxlength) + 1
      igpu      = int((iorigin - (iorigin%maxlength)) + 1) # get the first index of the number
      return (numbers0, numbers0[inumbers], inumbers, igpu, 0, mposition)
    end
    # if a number is selected && previous click was left && still only left button ist clicked
    if inumbers0::Int > 0 && mbutton0 == 0 && length(mbuttons) == 1 && in(0, mbuttons)
      xdiff                    = mposition[1] - mposition0[1]
      numbers0[inumbers0::Int]      = value0 + (float32(xdiff)/ 50.0f0)
      #numbertex[inumbers0]     = numbers0[inumbers0]
      a = mod(inumbers0::Int-1, size(numbers0::Int, 1))
      b = div(inumbers0::Int-1, size(numbers0::Int, 1))
      for (k,c) in enumerate(makestring(numbers0[inumbers0::Int], maxlength))
        textgpu[igpu0+k-1] = GLGlyph(c, a, (b-1)+b*maxlength+k, 0)
      end
      return (numbers0, value0, inumbers0::Int, igpu0, 0, mposition0)
    end
    return (numbers0, zero(eltype(numbers0)), -1, -1, -1, Vector2(0.0))
  end
  convert2{XX}(::Type{XX}, x::Vector{XX}) = first(x)
  convert2(t::Type, x) = convert(t, x)
  obj, lift(x->convert2(T,first(x)), number_lift)
end


text:
@testset "Constant input" begin
    c = ConstantInput(1.0)

    @test eltype(c) == Float64
    @test Base.IteratorSize(typeof(c)) == Base.IsInfinite()
    @test Base.IteratorEltype(typeof(c)) == Base.HasEltype()

    i = 0
    for ci in c
        @test ci == 1.0
        i += 1
        i <= 10 || break
    end
    @test i == 11

    @test  length(nextinput(c, 2)) == 2
    i = 1
    for val in nextinput(c, 2)
        @test i <= 2 && val == 1.0
        i += 1
    end

    c2 = map(x->x^2, c)
    @test c2.U == (c.U)^2
 end

@testset "Varying input" begin
    v = VaryingInput([1.0, 2.0, 3.0])

    @test eltype(v) == Float64
    @test Base.IteratorSize(typeof(v)) == Base.HasLength()
    @test Base.IteratorEltype(typeof(v)) == Base.HasEltype()

    @test length(v) == 3
    @test collect(v) == [1.0, 2.0, 3.0]

    @test  length(nextinput(v, 2)) == 2
    i = 1
    for val in nextinput(v, 2)
        @test i <= 2 && val == v.U[i]
        i += 1
    end

    v2 = map(x->x.^2, v)
    @test all([v2.U[i] == vi^2 for (i, vi) in enumerate(v)])

    v = VaryingInput([1:3., 4:6.])
    v3 = map(x->norm(x, Inf), v)
    @test all([v3.U[i] == maximum(vi) for (i, vi) in enumerate(v)])
 end


text:
struct EnumerableDefaultIfEmpty{T,S} <: Enumerable
    source::S
    default_value::T
end

Base.eltype(iter::Type{EnumerableDefaultIfEmpty{T,S}}) where {T,S} = T

_default_value_expr(::Type{T}) where {T} = :( DataValues.DataValue{$T}() )

_default_value_expr(::Type{T}) where {T<:DataValues.DataValue} = :( $T() )

function _default_value_expr(::Type{T}) where {T<:NamedTuple}
    return :( NamedTuple{$(fieldnames(T))}( ($( (_default_value_expr(fieldtype(T,i)) for i in 1:length(fieldnames(T)))...   ),)) )
end

@generated function default_if_empty(source::S) where {S}
    T_source = eltype(source)

    default_value_expr = _default_value_expr(T_source)

    q = quote
        default_value = $default_value_expr

        T = typeof(default_value)

        return EnumerableDefaultIfEmpty{T,$S}(source, default_value)
    end

    return q
end

function default_if_empty(source::S, default_value::TD) where {S,TD}
    T = eltype(source)
    if T!=TD
        error("The default value must have the same type as the elements from the source.")
    end
    return EnumerableDefaultIfEmpty{T,S}(source, default_value)
end

function Base.iterate(iter::EnumerableDefaultIfEmpty{T,S}) where {T,S}
    s = iterate(iter.source)

    if s===nothing
        return iter.default_value, nothing
    else
        return convert(T,s[1]), s[2]
    end
end

function Base.iterate(iter::EnumerableDefaultIfEmpty{T,S}, state) where {T,S}
    state===nothing && return nothing

    s = iterate(iter.source, state)
    if s===nothing
        return nothing
    else
        return convert(T, s[1]), s[2]
    end
end


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

mutable struct IMU <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function IMU(; kwargs...)
        obj = new(meta(IMU), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct IMU
const __meta_IMU = Ref{ProtoMeta}()
function meta(::Type{IMU})
    ProtoBuf.metalock() do
        if !isassigned(__meta_IMU)
            __meta_IMU[] = target = ProtoMeta(IMU)
            fnum = Int[1,2,3,5,6,7,8]
            allflds = Pair{Symbol,Union{Type,String}}[:acc_x => Float64, :acc_y => Float64, :acc_z => Float64, :gyr_x => Float64, :gyr_y => Float64, :gyr_z => Float64, :time => Float64]
            meta(target, IMU, allflds, ProtoBuf.DEF_REQ, fnum, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_IMU[]
    end
end
function Base.getproperty(obj::IMU, name::Symbol)
    if name === :acc_x
        return (obj.__protobuf_jl_internal_values[name])::Float64
    elseif name === :acc_y
        return (obj.__protobuf_jl_internal_values[name])::Float64
    elseif name === :acc_z
        return (obj.__protobuf_jl_internal_values[name])::Float64
    elseif name === :gyr_x
        return (obj.__protobuf_jl_internal_values[name])::Float64
    elseif name === :gyr_y
        return (obj.__protobuf_jl_internal_values[name])::Float64
    elseif name === :gyr_z
        return (obj.__protobuf_jl_internal_values[name])::Float64
    elseif name === :time
        return (obj.__protobuf_jl_internal_values[name])::Float64
    else
        getfield(obj, name)
    end
end

mutable struct VICON <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function VICON(; kwargs...)
        obj = new(meta(VICON), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct VICON
const __meta_VICON = Ref{ProtoMeta}()
function meta(::Type{VICON})
    ProtoBuf.metalock() do
        if !isassigned(__meta_VICON)
            __meta_VICON[] = target = ProtoMeta(VICON)
            allflds = Pair{Symbol,Union{Type,String}}[:pos_x => Float32, :pos_y => Float32, :pos_z => Float32, :quat_w => Float32, :quat_x => Float32, :quat_y => Float32, :quat_z => Float32, :time => Float64]
            meta(target, VICON, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_VICON[]
    end
end
function Base.getproperty(obj::VICON, name::Symbol)
    if name === :pos_x
        return (obj.__protobuf_jl_internal_values[name])::Float32
    elseif name === :pos_y
        return (obj.__protobuf_jl_internal_values[name])::Float32
    elseif name === :pos_z
        return (obj.__protobuf_jl_internal_values[name])::Float32
    elseif name === :quat_w
        return (obj.__protobuf_jl_internal_values[name])::Float32
    elseif name === :quat_x
        return (obj.__protobuf_jl_internal_values[name])::Float32
    elseif name === :quat_y
        return (obj.__protobuf_jl_internal_values[name])::Float32
    elseif name === :quat_z
        return (obj.__protobuf_jl_internal_values[name])::Float32
    elseif name === :time
        return (obj.__protobuf_jl_internal_values[name])::Float64
    else
        getfield(obj, name)
    end
end

mutable struct IMU_VICON <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function IMU_VICON(; kwargs...)
        obj = new(meta(IMU_VICON), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct IMU_VICON
const __meta_IMU_VICON = Ref{ProtoMeta}()
function meta(::Type{IMU_VICON})
    ProtoBuf.metalock() do
        if !isassigned(__meta_IMU_VICON)
            __meta_IMU_VICON[] = target = ProtoMeta(IMU_VICON)
            allflds = Pair{Symbol,Union{Type,String}}[:imu => IMU, :vicon => VICON]
            meta(target, IMU_VICON, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_IMU_VICON[]
    end
end
function Base.getproperty(obj::IMU_VICON, name::Symbol)
    if name === :imu
        return (obj.__protobuf_jl_internal_values[name])::IMU
    elseif name === :vicon
        return (obj.__protobuf_jl_internal_values[name])::VICON
    else
        getfield(obj, name)
    end
end

export IMU, VICON, IMU_VICON


text:
# This script checks sparse arrays can be used as caches for Γ₅, Γ₆, Λ_sss, and Σ_sss
using BenchmarkTools, RiskAdjustedLinearizations, Test
include(joinpath(dirname(@__FILE__), "..", "..", "examples", "rbc_cc", "rbc_cc.jl"))
include(joinpath(dirname(@__FILE__), "..", "..", "examples", "crw", "crw.jl"))

@testset "Sparse caching for Γ₅, Γ₆, Λ, and Σ (RBC-CC)" begin
    m_rbc_cc = RBCCampbellCochraneHabits()
    m = rbc_cc(m_rbc_cc, 0; sparse_arrays = true)

    @test issparse(m[:Γ₅])
    @test issparse(m[:Γ₆])
    @test issparse(m[:Λ_sss])
    @test issparse(m[:Σ_sss])

    # Solve!
    solve!(m; algorithm = :relaxation, verbose = :none)

    sssout = JLD2.jldopen(joinpath(dirname(@__FILE__), "..", "reference", "rbccc_sss_iterative_output.jld2"), "r")

    @test isapprox(sssout["z_rss"], m.z, atol=1e-4)
    @test isapprox(sssout["y_rss"], m.y, atol=1e-4)
    @test isapprox(sssout["Psi_rss"], m.Ψ, atol=1e-4)
end

@testset "Sparse caching for Γ₅, Γ₆, Λ, and Σ with jump-dependent shock matrices (CRW)" begin
m_crw = CoeurdacierReyWinant()
m = crw(m_crw; sparse_arrays = true)

solve!(m; algorithm = :homotopy, step = .5, verbose = :none)
sssout = JLD2.jldopen(joinpath(dirname(@__FILE__), "..", "reference", "crw_sss.jld2"), "r")

    @test isapprox(sssout["z_rss"], m.z, atol=1e-4)
    @test isapprox(sssout["y_rss"], m.y, atol=1e-4)
    @test isapprox(sssout["Psi_rss"], m.Ψ, atol=1e-4)
end