text:
using Revise
using JuMP, EAGO, Ipopt

println(" ")
for i=1:6
    println("----- CASE NOTE: NEW NOT TEST AFFINE SPEC OFF recent -----")
end

println(" ")

m = Model(optimizer_with_attributes(EAGO.Optimizer, "verbosity" => 1,
                                                    "output_iterations" => 1,
                                                    "iteration_limit" => 200,

                                                    "cp_depth" => -1,
                                                    "cp_repetitions" => 2,
                                                    "cp_forward_reverse_limit" => 2,

                                                    "cut_min_iterations" => 0,
                                                    "cut_max_iterations" => 0,
                                                    "objective_cut_on" => false,
                                                    "subgrad_tighten" => false,

                                                    "obbt_depth" => -1,
                                                    "obbt_repetitions" => 4,
                                                    "obbt_aggressive_on" => true,

                                                    "upper_bounding_depth" => 12,

                                                    "fbbt_lp_depth" => -1,
                                                    "fbbt_lp_repetitions" => 3,
                                                    "relax_tag" => NS()))

#m = Model(Ipopt.Optimizer)

# ----- Variables ----- #
x_Idx = Any[2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]
@variable(m, x[x_Idx])

#JuMP.set_lower_bound(x[1], -30.0)
#JuMP.set_upper_bound(x[1], 19.0)

JuMP.set_lower_bound(x[2], 0.0)
JuMP.set_upper_bound(x[2], 2.0)

JuMP.set_lower_bound(x[3], 0.0)
JuMP.set_upper_bound(x[3], 1.6)

JuMP.set_lower_bound(x[4], 0.0)
JuMP.set_upper_bound(x[4], 1.2)

JuMP.set_lower_bound(x[5], 0.0)
JuMP.set_upper_bound(x[5], 5.0)

JuMP.set_lower_bound(x[6], 0.0)
JuMP.set_upper_bound(x[6], 2.0)

JuMP.set_lower_bound(x[7], 0.85)
JuMP.set_upper_bound(x[7], 0.93)

JuMP.set_lower_bound(x[8], 0.9)
JuMP.set_upper_bound(x[8], 0.95)

JuMP.set_lower_bound(x[9], 3.0)
JuMP.set_upper_bound(x[9], 12.0)

JuMP.set_lower_bound(x[10], 1.2)
JuMP.set_upper_bound(x[10], 4.0)

JuMP.set_lower_bound(x[11], 1.45)
JuMP.set_upper_bound(x[11], 1.62)

JuMP.set_lower_bound(x[12], 0.99)
JuMP.set_upper_bound(x[12], 1.01010101010101)

JuMP.set_lower_bound(x[13], 0.99)
JuMP.set_upper_bound(x[13], 1.01010101010101)

JuMP.set_lower_bound(x[14], 0.9)
JuMP.set_upper_bound(x[14], 1.11111111111111)

JuMP.set_lower_bound(x[15], 0.99)
JuMP.set_upper_bound(x[15], 1.01010101010101)

# ----- Constraints ----- #
#@constraint(m, e2, -0.819672131147541*x[2]+x[5]-0.819672131147541*x[6] == 0.0)

#@NLconstraint(m, e3, 0.98*x[4]-(0.01*(x[5]*x[10])+x[4])*x[7] == 0.0) # first

#@constraint(m, e4, -x[2]*x[9] + 10*x[3] + x[6] == 0.0)

#@NLconstraint(m, e7, x[10]*x[14] + 22.2*x[11] == 35.82)

#@NLconstraint(m, e8, x[11]*x[15]-3*x[8] == -1.33)

#@NLconstraint(m, e5, x[5]*x[12]-(1.12+0.13167*x[9]-0.0067*x[9]^2)*x[2] == 0.0)

#@NLconstraint(m, e6, x[8]*x[13]-0.01*(1.098*x[9]-0.038*x[9]^2)-0.325*x[7] == 0.57425)

#=
@NLconstraint(m, e3a, 0.98*x[4]-x[7]*(0.01*x[5]*x[10]+x[4]) <= 0.0)
@NLconstraint(m, e4a, -x[2]*x[9]+10*x[3]+x[6] <= 0.0)
@NLconstraint(m, e7a, x[10]*x[14]+22.2*x[11] <= 35.82)
@NLconstraint(m, e8a, x[11]*x[15]-3*x[8] <= -1.33)
@NLconstraint(m, e5a, x[5]*x[12]-x[2]*(1.12+0.13167*x[9]-0.0067*x[9]^2) <= 0.0)
@NLconstraint(m, e6a, x[8]*x[13]-0.01*(1.098*x[9]-0.038*x[9]^2)-0.325*x[7] <= 0.57425)
=#
#=
@NLconstraint(m, e3b, -(0.98*x[4]-x[7]*(0.01*x[5]*x[10]+x[4])) <= 0.0)
@NLconstraint(m, e4b, -(-x[2]*x[9]+10*x[3]+x[6]) <= 0.0)
@NLconstraint(m, e7b, -(x[10]*x[14]+22.2*x[11]) <= -35.82)
@NLconstraint(m, e8b, -(x[11]*x[15]-3*x[8]) <= 1.33)
@NLconstraint(m, e5b, -(x[5]*x[12]-x[2]*(1.12+0.13167*x[9]-0.0067*x[9]^2)) <= 0.0)
@NLconstraint(m, e6b, -(x[8]*x[13]-0.01*(1.098*x[9]-0.038*x[9]^2)-0.325*x[7]) <= -0.57425)
=#
# EQ TO LEQ/LEQ fixes issue....
# OPTION 1: CC BOUND NOT ADDED CORRECTLY...
# OPTION 2: CC BOUND INCORRECT FOR SOME TERMS...

# ALL EQ -> (0.0, 0.0) BOUNDS
# ALL LEQ -> (-Inf, 0.0) BOUNDs

# NOT A SQR VS MULT ISSUE
#@NLconstraint(m, e5, x[5]*x[12]-x[2]*(1.12+0.13167*x[9]-0.0067*x[9]*x[9]) == 0.0)
#@NLconstraint(m, e6, x[8]*x[13]-0.01*(1.098*x[9]-0.038*x[9]*x[9])-0.325*x[7] == 0.57425)

# NOT A BOUND ISSUE
#@NLconstraint(m, e6, x[8]*x[13]-0.01*(1.098*x[9]-0.038*x[9]*x[9])-0.325*x[7] - 0.57425 == 0.0)
#@NLconstraint(m, e7, x[10]*x[14]+22.2*x[11] - 35.82 == 0.0)
#@NLconstraint(m, e8, x[11]*x[15]-3*x[8] + 1.33 == 0.0)

#=
-1.7650
=#
# ----- Objective ----- #
#@NLconstraint(m, (5.04*x[2] + 0.35*x[3] + x[4] + 3.36*x[6] - 6.3*x[5]*x[8]) - x[1] == 0.0)
#@objective(m, Min, x[2])

@objective(m, Min, 5.04*x[2] + 0.35*x[3] + x[4] + 3.36*x[6] - 6.3*x[5]*x[8] + x[4]*(x[6] -1.0)  + x[8]*x[4] - (x[4]+2.0)^2)

#@NLobjective(m, Min, 5.04*x[2] + 0.35*x[3] + x[4] + 3.36*x[6] - 6.3*x[5]*x[8] + x[7] +  x[10]*x[14] + 22.2*x[11] + x[11]*x[15]-3*x[8]
#                       + x[5]*x[12]-(1.12+0.13167*x[9]-0.0067*x[9]^2)*x[2] + x[13])
JuMP.optimize!(m)

println("primal status = $(primal_status(m))")
println("termination status = $(termination_status(m))")

x2 = Interval(0.0, 2.0)
x3 = Interval(0.0, 1.6)
x4 = Interval(0.0, 1.2)

x5 = Interval(0.0, 5.0)
x6 = Interval(0.0, 2.0)
x8 = Interval(0.9, 0.95)
nbnds = 5.04*x2 + 0.35*x3 + x4 + 3.36*x6 - 6.3*x5*x8
println("nbnds = $nbnds")


text:
function isuniqueroot(f,
                      a::arb,
                      b::arb;
                      evaltype = :taylor)
    if evaltype == :ball
        @warn ("Cannot determine if a root is unique with only "
               *"evaluation of f, it requires an evaluation of the Taylor "
               *"expansion")

        return false, true
    elseif evaltype == :taylor
        x = setinterval(a, b)
        PP = ArbPolyRing(parent(x), :x)

        # Evaluate f' on the interval and check that it's non-zero
        dy = f(arb_series(PP([x, parent(x)(1)])))

        if contains_zero(dy[1])
            return false, true
        end

        # Evaluate f at the endpoints and check that they have different signs
        ya = f(a)
        yb = f(b)

        if ispositive(ya) && isnegative(yb) || isnegative(ya) && ispositive(yb)
            return true, true
        end

        if ispositive(ya) && ispositive(yb) || isnegative(ya) && isnegative(yb)
            return false, false
        end

        return false, true
    end

    return false, true
end

function refine_root(f,
                     a::arb,
                     b::arb)
    root = setinterval(a, b)
    PP = ArbPolyRing(parent(root), :x)

    iterations = 1
    max_iterations = 5div(precision(parent(root)), 64)
    previous_precision = -typemin(Int)
    current_precision = rel_accuracy_bits(root)
    # Perform a maximum of max_iterations iterations or until the
    # precision no longer increases
    while iterations < 20 && previous_precision < current_precision
        iterations += 1
        y = f(midpoint(root))
        dy = f(arb_series(PP([root, parent(root)(1)])))[1]

        root = setintersection(midpoint(root) - y/dy, root)
        if isnan(root)
            return (a, b)
        end

        previous_precision = current_precision
        current_precision = rel_accuracy_bits(root)
    end

    return getinterval(root)
end

function refine_root(poly::arb_poly,
                     (a, b)::Tuple{arb, arb})
    root = setinterval(a, b)
    # Perform a couple of simple newton iterations to gain at least 20
    # bits of relative accuracy
    iterations = 1
    while iterations < 5 && rel_accuracy_bits(root) < 20
        iterations += 1
        y = evaluate(poly, midpoint(root))
        dy = evaluate(derivative(poly), root)

        root = setintersection(midpoint(root) - y/dy, root)

        if isnan(root)
            return (a, b)
        end
    end

    # Refine the root using method in arb
    convergence_interval = parent(root)(root)
    convergence_factor_ball = parent(root)(0)

    GC.@preserve convergence_factor_ball begin
        convergence_factor = ccall((:arb_mid_ptr, Nemo.libarb), Ptr{Nemo.arf_struct},
                                   (Ref{arb}, ), convergence_factor_ball)

        ccall(("_arb_poly_newton_convergence_factor", Nemo.libarb), Cvoid,
              (Ptr{Nemo.arf_struct}, Ptr{arb}, Int, Ref{arb}, Int),
              convergence_factor, poly.coeffs, length(poly), convergence_interval, precision(parent(root)))

        # FIXME: The extra bits used in the computation should be
        # looked over. I just chose it randomly.
        ccall(("_arb_poly_newton_refine_root", Nemo.libarb), Cvoid,
              (Ref{arb}, Ptr{arb}, Int, Ref{arb}, Ref{arb}, Ref{Nemo.arf_struct}, Int, Int),
              root, poly.coeffs, length(poly), root, convergence_interval,
              convergence_factor, div(precision(parent(root)), 4), precision(parent(root)))
    end

    if a < root
        return getinterval(root)
    else
        return (a, b)
    end
end

function isolateroots(f,
                      a::arb,
                      b::arb;
                      evaltype = :ball,
                      refine = false,
                      atol = sqrt(eps()),
                      rtol = sqrt(eps()),
                      maxevals = 10^3,
                      maxfound = typemax(Int),
                      store_trace = false,
                      show_trace = false,
                      extended_trace = false)
    @assert isfinite(a) && isfinite(b) && a < b

    intervals = [(a, b)]
    numfound = 0
    numevals = 0
    iteration = 0

    found = Vector{eltype(intervals)}()
    flags = Vector{Bool}()

    while !isempty(intervals) && numfound < maxfound && numevals < maxevals
        iteration += 1
        nextintervals = Vector{eltype(intervals)}()

        for (lower, upper) in intervals
            numevals += 1

            mayberoot = contains_zero(f(setinterval(lower, upper)))
            unique = false

            if mayberoot
                if evaltype == :taylor
                    # Try to prove uniqueness
                    unique, mayberoot = isuniqueroot(f, lower, upper, evaltype = evaltype)
                end

                if unique && refine
                    # Try to refine the root
                    (lower, upper) = refine_root(f, lower, upper)
                end
            end

            if mayberoot
                x = setinterval(lower, upper)

                if (radius(x) < atol) || (radius(x)/x < rtol)
                    push!(found, (lower, upper))
                    push!(flags, unique)
                else
                    midpoint = 0.5*(lower + upper)
                    push!(nextintervals, (lower, midpoint))
                    push!(nextintervals, (midpoint, upper))
                end
            end
        end

        if show_trace
            println("$iteration      $(length(nextintervals))")
        end

        intervals = nextintervals
    end

    found = [found; intervals]
    flags = [flags; zeros(Bool, length(intervals))]

    p = sortperm(found, by=x -> BigFloat(getindex(x, 1)))

    found[p], flags[p]
end

function isolateroots(poly::arb_poly,
                      a::arb,
                      b::arb;
                      atol = sqrt(eps()),
                      rtol = sqrt(eps()),
                      maxevals = 10^3,
                      maxfound = typemax(Int),
                      store_trace = false,
                      show_trace = false,
                      extended_trace = false)
    @assert isfinite(a) && isfinite(b) && a < b

    intervals = [(a, b)]
    numfound = 0
    numevals = 0
    iteration = 0

    found = Vector{eltype(intervals)}()
    flags = Vector{Bool}()

    while !isempty(intervals) && numfound < maxfound && numevals < maxevals
        iteration += 1
        nextintervals = Vector{eltype(intervals)}()

        for (lower, upper) in intervals
            numevals += 1

            x = setinterval(lower, upper)
            y = evaluate(poly, x)

            if contains_zero(y)
                dy = evaluate(derivative(poly), x)

                if contains_zero(dy)
                    ylower = evaluate(poly, lower)
                    yupper = evaluate(poly, upper)

                    if ylower > 0 && yupper < 0 || ylower < 0 && yupper > 0
                        (lower, upper) = refine_root(poly, (lower, upper))
                        push!(found, (lower, upper))
                        push!(flags, true)
                    elseif !(ylower > 0 && yupper > 0 || ylower < 0 && yupper < 0) &&
                        if (radius(x) < atol) || (radius(x)/x < rtol)
                            push!(found, (lower, upper))
                            push!(flags, false)
                        else
                            midpoint = 0.5*(lower + upper)
                            push!(nextintervals, (lower, midpoint))
                            push!(nextintervals, (midpoint, upper))
                        end
                    end
                elseif (radius(x) < atol) || (radius(x)/x < rtol)
                    push!(found, (lower, upper))
                    push!(flags, false)
                else
                    midpoint = 0.5*(lower + upper)
                    push!(nextintervals, (lower, midpoint))
                    push!(nextintervals, (midpoint, upper))
                end
            end
        end

        if show_trace
            println("$iteration      $(length(nextintervals))")
        end

        intervals = nextintervals
    end

    found = [found; intervals]
    flags = [flags; zeros(Bool, length(intervals))]

    p = sortperm(found, by=x -> BigFloat(getindex(x, 1)))

    found[p], flags[p]
end

ext: jl


text:
# Utility function which returns bitstring mask for one output of the
# packed representation.
# Examples:  output_mask(2) returns  0b1111 = 0XF
#        and output_mask(3) returns  0b11111111 = 0XFF
function output_mask(num_inputs)
    mask = convert(BitString, 0x1)
    for i in 0:(num_inputs-1)
        mask |= (mask << 2^i)
    end
    return mask
end

# Returns an array of BitStrings which are the standard inputs to a chromosome
#    when it is executed.
# num_inputs  is the number of inputs of the chromosome.
function std_input_context(num_inputs)
   if 2^num_inputs > 8*sizeof(BitString)
      error("number of inputs is too large for size of BitString in std_input_context")
   end
   std_in = Array(BitString,num_inputs)
   std_in[1] = 0b10
   for i in 2:num_inputs
      for j in i:-1:2
         std_in[j] = std_in[j-1] $ std_in[j-1] << 2^(i-1)
      end
      std_in[1] = output_mask(i-1) << 2^(i-1)
   end
   return std_in
end


text:
################################################################################
# Copyright (C) Marcus Greiff 2020
#
# @file example_FSF_attitude_utils.jl
# @author Marcus Greiff
# @date July 2020
#
# Utility functions for initialization and plotting in the FSF example.
#
################################################################################

function get_errors_SO3(qvec, qrvec, wvec, wrvec)
    NT = size(qvec, 2);
    eRvec = zeros(3, NT);
    ewvec = zeros(3, NT);
    zvec  = zeros(2, NT);
    znorm = zeros(1, NT);
    for t = 1:NT
        R  = quat_2_SO3(reshape(qvec[:,t],  (4,1)))
        Rr = quat_2_SO3(reshape(qrvec[:,t], (4,1)))
        w  = wvec[:,t];
        wr = wrvec[:,t];

        eRvec[:,t]  = SO3_vee(Rr'*R - R'*Rr) / 2.0;
        ewvec[:,t]  = w -  R'*Rr * wr;
        zvec[1,t]   = norm(eRvec[:,t])
        zvec[2,t]   = norm(ewvec[:,t])
        znorm[1,t]  = zvec[:,t]'*zvec[:,t]

    end
    return eRvec, ewvec, zvec, znorm
end

function get_errors_SU2(qvec, qrvec, wvec, wrvec)
    NT = size(qvec, 2);
    eXvec = zeros(3, NT);
    ewvec = zeros(3, NT);
    zvec  = zeros(2, NT);
    znorm = zeros(1, NT);
    for t = 1:NT
        X  = quat_2_SU2(reshape(qvec[:,t],  (4,1)))
        Xr = quat_2_SU2(reshape(qrvec[:,t], (4,1)))
        w  = wvec[:,t];
        wr = wrvec[:,t];

        Xe = Xr'*X;

        eXvec[:,t]  = SU2_vee(Xe - tr(Xe)*Matrix{Float64}(I, 2, 2)) / 2.0;
        ewvec[:,t]  = w -  SU2_vee(Xe'*SU2_hat(wr)*Xe);
        zvec[1,t]   = norm(eXvec[:,t])
        zvec[2,t]   = norm(ewvec[:,t])
        znorm[1,t]  = zvec[:,t]'*zvec[:,t]

    end
    return eXvec, ewvec, zvec, znorm
end

function get_bounds(M1, M2, W, zvec)
    NT = size(zvec, 2);
    Vmin = zeros(1, NT);
    Vmax = zeros(1, NT);
    dV   = zeros(1, NT);
    for t = 1:NT
        z = zvec[:,t];
        Vmin[1,t] = +z'*M1*z;
        Vmax[1,t] = +z'*M2*z;
        dV[1,t]   = -z'*W*z;
    end
    return Vmin, Vmax, dV
end

function maximal_feasible_eps_SU2(
    C::con_state_qw_fsf_t,
    X0::Array{Complex{T},2},
    Xr0::Array{Complex{T},2},
    w0::Array{T,2},
    wr0::Array{T,2};
    isContinuous=true
) where {T<:Number}

    # Extract parameters
    kX, kc, kw  = C.gain_kR, C.gain_kc, C.gain_kw;
    J           = ntuple_2_mat(C.inertia, 3, 3);

    if isContinuous
        # Compute maximum permissible phi
        Gamma     = abs(0.5 * tr(Matrix{Float64}(I, 2, 2) - Xr0'*X0));
        Xe = Xr0'*X0;
        eX = SU2_vee(Xe - tr(Xe)*Matrix{Float64}(I, 2, 2)) / 2.0
        ew = w0 -  SU2_vee(Xe'*SU2_hat(wr0)*Xe)
        attpart   = kX * Gamma;
        crosspart = kc * ew' * eX;
        ratepart  = ew'*J*ew / 2.0;
        Lyap0     = attpart + crosspart[1,1] + ratepart[1,1];
        phi       = Lyap0/kX;
    else
        # phi is fixed at 1.0 in the discontinuous case
        phi = 1.0
    end

    M1, M2, W   = get_M1M2W_SU2(C, phi)

    M1max, M1min = maximum(eigvals(M1)), minimum(eigvals(M1));
    M2max, M2min = maximum(eigvals(M2)), minimum(eigvals(M2));
    Wmax,  Wmin  = maximum(eigvals(W)),  minimum(eigvals(W));

    mindecayrate = Wmin/M2max;
    maxeps       = ((M1min*Wmin) / M2max)*phi*(2.0 - phi);
    uniboundgain = (M2max / (M1min*Wmin));

    return maxeps, mindecayrate, uniboundgain, phi, M1, M2, W
end

function get_info(
    C::con_state_qw_fsf_t,
    controllerType::Integer,
    q0::Array{T,2},
    qr0::Array{T,2},
    w0::Array{T,2},
    wr0::Array{T,2}
) where {T<:Number}

    controllerStrings = [
        "continuous FSF controller on SO(3) (called from C)",
        "robust FSF controller on SO(3) (called from C)",
        "continuous FSF controller on SU(2) (called from C)",
        "discontinuous FSF controller on SU(2) (called from C)",
        "robust FSF controller on SU(2) (called from C)"
    ]

    nstr = length(controllerStrings[controllerType]);
    println("\n$('~'^(34+nstr))")
    println("~~~ Simulation done with the $(controllerStrings[controllerType]) ~~~")
    println("$('~'^(34+nstr))")

    if controllerType in [1,2]
        R0  = quat_2_SO3(q0);
        Rr0 = quat_2_SO3(qr0);

        phi0 = tr(I-Rr0'*R0)/2.0;
        println("* The initial attitude error is: Psi(Rr(t0), R(t0)) = $(phi0)")
        if (phi0 > 1.5)
            println("* Warning: This is relatively large, requiring very small initial attitude rate errors")
        end
        maxeps, mindecayrate, uniboundgain, phi, M1, M2, W = maximal_feasible_eps_SO3(C, R0, Rr0, w0, wr0)
    end
    if controllerType in [3,4,5]
        X0  = quat_2_SU2(q0);
        Xr0 = quat_2_SU2(qr0);

        phi0 = abs(0.5 * tr(Matrix{Float64}(I, 2, 2) - Xr0'*X0));
        println("* The initial attitude error is: Gamma(Xr(t0), X(t0)) = $(phi0)")
        if (phi0 > 1.5)
            println("* Warning: This is relatively large, requiring very small initial attitude rate errors")
        end
        if controllerType in [3,5]
            maxeps, mindecayrate, uniboundgain, phi, M1, M2, W = maximal_feasible_eps_SU2(C, X0, Xr0, w0, wr0; isContinuous=true)
        end
        if controllerType in [4]
            maxeps, mindecayrate, uniboundgain, phi, M1, M2, W = maximal_feasible_eps_SU2(C, X0, Xr0, w0, wr0; isContinuous=false)
        end
    end

    println("* Any V(t0)/kR = $(phi) < phi < 2 can be used in the stability proof. We let phi = $(phi).")
    if (phi > 2);
        println("* Waring. The system may move outside of the domain of exponential attraction given phi")
        println("* Lyap0/kR = $(phi), please retune the controller, or lower the initial rate error")
    end

    println("* Worst case decay rate of the Lyapunov function : $(mindecayrate)")
    if controllerType in [2,5]
        println("* Upper bound on the allowed epsilon given phi   : $(maxeps)")
        println("* Uniform ultimate bound K*eps, where K is       : $(uniboundgain)")
    end
    return M1, M2, W, uniboundgain
end

function initialize_FSF_attitude_example(controllerType::Integer)

    # Define the controller parameters
    J  = rand_PSD(3);
    J  = J./opnorm(J);
    if controllerType in [1,2]
        kR = 5.0;
    else
        kR = 10.0;
    end
    kw = 2.0;

    # Compute a feasible cross gain kc
    if controllerType in [1,2]
        # In case of an SO(3) configured controller (continuous or robust)
        maxkc = maximal_feasible_kc_SO3(kR, kw, J);
    else
        maxkc = maximal_feasible_kc_SU2(kR, kw, J); # here we let kR=kX for simplicity
    end
    kc = 0.1 * maxkc;

    # Define the controller state object
    C  = con_state_qw_fsf_t(J, kR, kc, kw)

    # Initial states of the reference system
    qr0 = rand(4,1).-0.5;
    qr0 = qr0./norm(qr0);
    wr0 = zeros(3,1);

    # Initial states of the physical system
    q0 = rand(4,1).-0.5;
    q0 = q0./norm(q0);
    w0 = rand(3,1).-0.5;

    # Select a feasible epsilon gain, with the uptimate bound characterize by alpha
    alpha = 1e-2;
    if controllerType == 2
        # In case of an SO(3) configured controller (only applicable in the robust case)
        maxeps, _, _, _, _, _ = maximal_feasible_eps_SO3(C, quat_2_SO3(q0), quat_2_SO3(qr0), w0, wr0);
        C.gain_eps = alpha*maxeps;
    end
    if controllerType == 5
        # In case of an SU(2) configured controller (only applicable in the robust case)
        maxeps, _, _, _, _, _ = maximal_feasible_eps_SU2(C, quat_2_SU2(q0), quat_2_SU2(qr0), w0, wr0);
        C.gain_eps = alpha*maxeps;
    end

    # Initial states for the debugging signals
    output = zeros(14,1);
    output[13,1] = 1.0;
    output[14,1] = 1.0;

    x0    = vcat(q0, w0, qr0, wr0, output)
    tspan = (0.0,20.0)
    return x0, tspan, C
end

function solution_analysis(
    sol::ODESolution,
    C::con_state_qw_fsf_t,
    controllerType::Integer;
    showPlot=[1,2,3],
    savePlot=false,
    namePlot="unnamed",
    showDisturbance=false
)

    # Posterior solution analysis
    t     = sol.t;
    data  = zeros(length(t),length(sol.u[1]));
    for i = 1:length(t); data[i,:] = sol.u[i]; end

    q     = data[:,1:4]';
    w     = data[:,5:7]';
    qr    = data[:,8:11]';
    wr    = data[:,12:14]';
    tau   = data[:,15:17]';
    taur  = data[:,18:20]';
    psi   = data[:,21]';
    gamma = data[:,22]';
    dist  = data[:,24:26]';
    lyap  = data[:,23]';
    lyapd = numerical_differentiation(lyap, t);

    # The first element of certain logged variables (computed in the C-code)
    # are not set at t0, hence we let tis value be the one at the next time step
    # to tidy up and simplify the plotting.
    psi[:,1], gamma[:,1] = psi[:,2], gamma[:,2];
    tau[:,1] = tau[:,2];
    taur[:,1] = taur[:,2];
    dist[:,1] = dist[:,2];
    lyap[:,1] = lyap[:,2];

    q0  = reshape(q[:,1],  (4,1));
    qr0 = reshape(qr[:,1], (4,1));
    w0  = reshape(w[:,1],  (3,1));
    wr0 = reshape(wr[:,1], (3,1));

    M1, M2, W, K = get_info(C, controllerType, q0, qr0, w0, wr0);

    # In case we are dealine with an SO(3) controller (continuous, or robust)
    if controllerType in [1,2]
        eAtt, eW, z, znorm = get_errors_SO3(q, qr, w, wr);
    end
    # In case we are dealine with an SU(2) controller
    if controllerType in [3,4,5]
        eAtt, eW, z, znorm = get_errors_SU2(q, qr, w, wr)
    end

    Vmin, Vmax, dV   = get_bounds(M1, M2, W, z);
    zbound           = ones(length(t),1).*K.*C.gain_eps;

    if 1 in showPlot
        pqqr = plot( t, q',    color=:black, linewidth=2, xaxis="Time (t)",yaxis="Attitude [.]",          label=[L"q(t)" nothing nothing nothing])
               plot!(t, qr',   color=:blue,  linewidth=2, xaxis="Time (t)",yaxis="Attitude [.]",          label=[L"q_r(t)" nothing nothing nothing])
        pwwr = plot( t, w',    color=:black, linewidth=2, xaxis="Time (t)",yaxis="Attitude rate [rad/s]", label=[L"\omega(t)" nothing nothing])
               plot!(t, wr',   color=:blue,  linewidth=2, xaxis="Time (t)",yaxis="Attitude rate [rad/s]", label=[L"\omega_{r}(t)" nothing nothing])
        pttr = plot( t, tau',  color=:black, linewidth=2, xaxis="Time (t)",yaxis="Torque [Nm]",           label=[L"\tau(t)" nothing nothing])
               plot!(t, taur', color=:blue,  linewidth=2, xaxis="Time (t)",yaxis="Torque [Nm]",           label=[L"\tau_{r}(t)" nothing nothing])
        pdis = plot( t, dist', color=:black, linewidth=2, xaxis="Time (t)",yaxis="Load disturbance [Nm]", label=[L"d(t)" nothing nothing])

        if showDisturbance
            plot(pqqr, pwwr, pttr, pdis, layout=(4,1), size=(1000,1000))
        else
            plot(pqqr, pwwr, pttr, layout=(3,1), size=(1000,750))
        end

        gui()
        if savePlot
            savefig(joinpath(AerialVehicleControl.CONT_LIB_DOCS, "$(namePlot)_states.png"))
        end
    end

    if 2 in showPlot

        attitudeErrors = [L"e_{R}(t)", L"e_{R}(t)", L"e_{X}(t)", L"e_{X}^{\pm}(t)", L"e_{X}(t)"]
        plotattd    = plot( t, psi',        color=:black,             linewidth=2,                            label=L"\Psi(R_r(t), R(t))")
                      plot!(t, gamma',      color=:blue,               linewidth=2,                            label=L"\Gamma(X_r(t), X(t))")
                      plot!(t, 2 .- gamma', color=:blue,  style=:dash, linewidth=2, yaxis="Attitude distance", label=L"\bar{\Gamma}(X_r(t), X(t))")
        plotatterr  = plot( t, eAtt', color=:black, linewidth=2,                   yaxis="Attitude error",    label=[attitudeErrors[controllerType] nothing nothing])
        plotrateerr = plot( t, eW',   color=:black, linewidth=2, xaxis="Time (t)", yaxis="Rate error",        label=[L"e_{\omega}(t)" nothing nothing])
        plot(plotattd, plotatterr,  plotrateerr, layout=(3,1), size=(1000,750))
        gui()
        if savePlot
            savefig(joinpath(AerialVehicleControl.CONT_LIB_DOCS, "$(namePlot)_errors.png"))
        end
    end

    if 3 in showPlot
        ubstyle    = :dash;
        lbstyle    = :dot;
        sigcolor   = :black;
        boundcolor = :blue;
        boundalpha = 0.1;
        fillcolor  = :black;

        bounds     = (Vmin[:].*0, Vmax[:]-Vmin[:]);
        plotlyap   = plot(t,  Vmin[:], ribbon=bounds, color=boundcolor, linewidth=2, style=lbstyle, fillcolor=fillcolor, fillalpha=boundalpha, label=L"z^{\top}M_1z")
                     plot!(t, Vmax[:],                color=boundcolor, linewidth=2, style=ubstyle,                                            label=L"z^{\top}M_2z")
                     plot!(t, lyap[:],                color=sigcolor,   linewidth=2, xaxis="Time (t)",yaxis=L"\mathcal{V}(t)",                 label=L"$\mathcal{V}(t)$")

        logVmin = log10.(abs.(Vmin[:]))
        logVmax = log10.(abs.(Vmax[:]))
        logLyap = log10.(abs.(lyap[:]))
        bounds  = (logVmin.*0, logVmax - logVmin)
        plotlyaplog = plot(t,  logVmin, ribbon=bounds, color=boundcolor, linewidth=2, style=lbstyle, fillcolor=fillcolor, fillalpha=boundalpha, label=L"$log_{10}(z^{\top}M_1z)$")
                      plot!(t, logVmax,                color=boundcolor, linewidth=2, style=ubstyle,                                            label=L"$log_{10}(z^{\top}M_2z)$")
                      plot!(t, logLyap,                color=sigcolor,   linewidth=2, xaxis="Time (t)",yaxis=L"log_{10}(\mathcal{V}(t))",       label=L"$log_{10}(\mathcal{V}(t))$")

        ind           = 3;
        tshort        = t[ind:end];
        dVshort       = dV[1,ind:end];
        lyapshort     = lyapd[1,ind:end];
        lbval         = 1.1*minimum([dVshort; lyapshort])
        bounds        = (dVshort .- lbval, 0 .* dVshort)
        plotlyapderiv = plot(tshort,    dVshort, ribbon=bounds, color=boundcolor, linewidth=2, style=ubstyle, fillcolor=fillcolor, fillalpha=boundalpha,   label=L"-z^{\top}Wz")
                        plot!(tshort, lyapshort,                color=sigcolor,   linewidth=2, xaxis="Time (t)",yaxis=L"\frac{d}{dt}\mathcal{V}(t)",       label=L"$({d}/{dt}){\mathcal{V}}(t)$")

        log10zbound   = zbound[:]; #log10.(abs.(zbound[:]))
        log10znorm    = znorm[:];  #log10.(abs.(znorm[:]))
        lbval         = minimum([log10zbound; log10znorm]) - 1.0;
        bounds        = (log10zbound .- lbval, 0 .*log10zbound);
        plotznorm     = plot(t,  log10zbound, ribbon=bounds, color=boundcolor, linewidth=2, style=ubstyle, fillcolor=fillcolor, fillalpha=boundalpha, label=L"$\lambda_M(M_2)(\lambda_m(M_1)\lambda_m(W))^{-1}\epsilon$")
                        plot!(t, log10znorm,                 color=sigcolor,   linewidth=2, xaxis="Time (t)",yaxis=L"$\|\|z(t)\|\|_2^2$", label=L"$\|\|z(t)\|\|_2^2$")

        if showDisturbance
            plot(plotlyap, plotlyaplog, plotlyapderiv, plotznorm, layout=(4,1), size=(1000,1000))
        else
            plot(plotlyap, plotlyaplog, plotlyapderiv, layout=(3,1), size=(1000,750))
        end
        gui()
        if savePlot
            savefig(joinpath(AerialVehicleControl.CONT_LIB_DOCS, "$(namePlot)_analysis.png"))
        end
    end
end


text:
export ForwardSensRStepAdapter, AbstractForwardSensSpaceAdapterSparse

abstract type AbstractForwardSensSpaceAdapterSparse <: AbstractForwardSensSpaceAdapter end

"""
Simple adapter based on reachability. Whenever the current FSP solution error is found to exceed the acceptable tolerance, the adapter will expand the state space by exploring all states that could be reached from the current state space within a set number of reaction events.

This adapter only works with `StateSpaceSparse`.
"""
struct ForwardSensRStepAdapter <: AbstractForwardSensSpaceAdapterSparse
    initial_step_count::Int
    max_step_count::Int
    dropstates::Bool
end

"""
`init!(statespace::StateSpaceSparse, adapter::RStepAdapter, p::Vector{RealT}, S::Vector{Vector{RealT}}, t::RealT, fsptol::RealT) where {RealT <: AbstractFloat}`

Initialize the state space, solution vector and sensitivity indices for FSP integration.
"""
function init!(statespace::StateSpaceSparse, adapter::ForwardSensRStepAdapter, p::Vector{RealT}, S::Vector{Vector{RealT}}, t::AbstractFloat, fsptol::AbstractFloat) where {RealT <: AbstractFloat}
    nold = get_state_count(statespace)
    expand!(statespace, adapter.initial_step_count)
    append!(p, zeros(RealT, get_state_count(statespace) - nold))
    for svec in S
        append!(svec, zeros(RealT, get_state_count(statespace) - nold))
    end
    nothing
end

"""
`adapt!(statespace::StateSpaceSparse, adapter::RStepAdapter, p::Vector{RealT}, sinks::Vector{RealT}, t::RealT, fsptol::RealT; integrator::Union{DEIntegrator, Nothing}=nothing) where {RealT <: AbstractFloat}`

Adapt the state space, probability vector and sensitivity vectors based on current error recorded in the `sinks` vector.
"""
function adapt!(statespace::StateSpaceSparse, adapter::ForwardSensRStepAdapter, p::Vector{RealT}, S::Vector{Vector{RealT}},
    sinks::Vector{RealT},
    dsinks::Vector{Vector{RealT}},
    t::RealT, tend::RealT,
    fsptol::RealT; integrator::Union{DEIntegrator,Nothing} = nothing) where {RealT<:AbstractFloat}

    if adapter.dropstates
        pids = sortperm(p)
        dropcount = sum((sum(p) .- cumsum(p[pids])) .>= (1.0 - t * fsptol / tend))
        dropids = sort(pids[1:dropcount])
        deleteat!(statespace, dropids)
        deleteat!(p, dropids)
        for svec in S
            deleteat!(svec, dropids)
        end
    end

    nold = get_state_count(statespace)
    expand!(statespace, adapter.max_step_count)
    append!(p, zeros(RealT, get_state_count(statespace) - nold))
    for svec in S
        append!(svec, zeros(RealT, get_state_count(statespace) - nold))
    end
    nothing
end


text:
"""
    LogitNormal

* Parametrization
parametrized in \\mu and \\sigma^2

* Score

* Fisher Information

* `time_varying_params` map.

* Default link
"""
LogitNormal

logit(x) = log(x) - log(1-x)

function score!(score_til::Matrix{T}, y::T, ::Type{LogitNormal}, param::Matrix{T}, t::Int) where T
    score_til[t, 1] = (logit(y) - param[t, 1])/param[t, 2]
    score_til[t, 2] = -(0.5/param[t, 2]) * (1 - ((logit(y) - param[t, 1])^2)/param[t, 2])
    return
end

function fisher_information!(aux::AuxiliaryLinAlg{T}, ::Type{LogitNormal}, param::Matrix{T}, t::Int) where T
    aux.fisher[1, 1] = 1/(param[t, 2])
    aux.fisher[2, 2] = 1/(2*(param[t, 2]^2))
    aux.fisher[2, 1] = 0
    aux.fisher[1, 2] = 0
    return
end

function log_likelihood(::Type{LogitNormal}, y::Vector{T}, param::Matrix{T}, n::Int) where T
    loglik = -0.5*log(2*pi)*n
    for t in 1:n
        loglik -= log((y[t] * (1 - y[t])) * sqrt(param[t, 2])) + 0.5*(logit(y[t]) - param[t, 1])^2/param[t, 2]
    end
    return -loglik
end

# Links
function link!(param_tilde::Matrix{T}, ::Type{LogitNormal}, param::Matrix{T}, t::Int) where T
    param_tilde[t, 1] = link(IdentityLink, param[t, 1])
    param_tilde[t, 2] = link(LogLink, param[t, 2], zero(T))
    return
end
function unlink!(param::Matrix{T}, ::Type{LogitNormal}, param_tilde::Matrix{T}, t::Int) where T
    param[t, 1] = unlink(IdentityLink, param_tilde[t, 1])
    param[t, 2] = unlink(LogLink, param_tilde[t, 2], zero(T))
    return
end
function jacobian_link!(aux::AuxiliaryLinAlg{T}, ::Type{LogitNormal}, param::Matrix{T}, t::Int) where T
    aux.jac[1] = jacobian_link(IdentityLink, param[t, 1])
    aux.jac[2] = jacobian_link(LogLink, param[t, 2], zero(T))
    return
end

# utils
function update_dist(::Type{LogitNormal}, param::Matrix{T}, t::Int) where T
    # LogitNormal here is parametrized as sigma^2
    return LogitNormal(param[t, 1], sqrt(param[t, 2]))
end

function params_sdm(d::LogitNormal)
    pars = Distributions.params(d)
    return (pars[1], pars[2]^2)
end

function num_params(::Type{LogitNormal})
    return 2
end


text:
# Hebrew char set
map(Char, 1488:1514)
ncodeunits('ש')
isascii('ש')
 codepoint('ש')

 snake2camel_case(s::AbstractString) = join(map(uppercasefirst, split(s, '_')))

 function snake2camel_case()
    s =snake2camel_case(clipboard())
    clipboard(s)
    s
 end

function camel2snake_case(s::AbstractString)
    words = String[]
    i = 1
    for j in findall(isuppercase, s)
        push!(words, lowercase(s[i:j-1]))
        i = j
    end
    push!(words, lowercase(s[i:end]))
    join(words, "_")
 end

max_stars_repo_name: FourMInfo/JuliaforBeginners

max_issues_repo_name: FourMInfo/JuliaforBeginners

max_forks_repo_name: FourMInfo/JuliaforBeginners


text:
export elem_in_algebra

parent_type(::Type{AlgAssAbsOrdElem{S, T}}) where {S, T} = AlgAssAbsOrd{S, T}

parent_type(::AlgAssAbsOrdElem{S, T}) where {S, T} = AlgAssAbsOrd{S, T}

@doc Markdown.doc"""
    parent(x::AlgAssAbsOrdElem) -> AlgAssAbsOrd

Returns the order containing $x$.
"""
@inline parent(x::AlgAssAbsOrdElem) = x.parent

################################################################################
#
#  Parent check
#
################################################################################

function check_parent(x::AlgAssAbsOrdElem{S, T}, y::AlgAssAbsOrdElem{S, T}) where {S, T}
  return parent(x) === parent(y)
end

################################################################################
#
#  Parent object overloading
#
################################################################################

(O::AlgAssAbsOrd{S, T})(a::T, check::Bool = true) where {S, T} = begin
  if check
    (x, y) = _check_elem_in_order(a, O)
    !x && error("Algebra element not in the order")
    return AlgAssAbsOrdElem{S, T}(O, deepcopy(a), y)
  else
    return AlgAssAbsOrdElem{S, T}(O, deepcopy(a))
  end
end

(O::AlgAssAbsOrd{S, T})(a::T, arr::Vector{fmpz}, check::Bool = false) where {S, T} = begin
  if check
    (x, y) = _check_elem_in_order(a, O)
    (!x || arr != y) && error("Algebra element not in the order")
    return AlgAssAbsOrdElem{S, T}(O, deepcopy(a), y)
  else
    return AlgAssAbsOrdElem{S, T}(O, deepcopy(a), deepcopy(arr))
  end
end

(O::AlgAssAbsOrd{S, T})(arr::Vector{fmpz}) where {S, T} = begin
  M = basis_matrix(O, copy = false)
  N = matrix(FlintZZ, 1, degree(O), arr)
  NM = N*M
  x = elem_from_mat_row(algebra(O), NM.num, 1, NM.den)
  return AlgAssAbsOrdElem{S, T}(O, x, deepcopy(arr))
end

(O::AlgAssAbsOrd{S, T})(a::AlgAssAbsOrdElem{S, T}, check::Bool = true) where {S, T} = begin
  b = elem_in_algebra(a) # already a copy
  if check
    (x, y) = _check_elem_in_order(b, O)
    !x && error("Algebra element not in the order")
    return AlgAssAbsOrdElem{S, T}(O, b, y)
  else
    return AlgAssAbsOrdElem{S, T}(O, b)
  end
end

(O::AlgAssAbsOrd)(a::T, check::Bool = true) where T = O(algebra(O)(a), check)

################################################################################
#
#  Deepcopy
#
################################################################################

function Base.deepcopy_internal(a::AlgAssAbsOrdElem, dict::IdDict)
  b = parent(a)()
  b.elem_in_algebra = Base.deepcopy_internal(a.elem_in_algebra, dict)
  if a.has_coord
    b.has_coord = true
    b.coordinates = Base.deepcopy_internal(a.coordinates, dict)
  end
  return b
end

################################################################################
#
#  Special elements
#
################################################################################

(O::AlgAssAbsOrd{S, T})() where {S, T} = AlgAssAbsOrdElem{S, T}(O)

one(O::AlgAssAbsOrd) = O(one(algebra(O)))

zero(O::AlgAssAbsOrd) = O(algebra(O)())

################################################################################
#
#  Element in algebra
#
################################################################################

@doc Markdown.doc"""
    elem_in_algebra(x::AlgAssAbsOrdElem; copy::Bool = true) -> AbsAlgAssElem
    elem_in_algebra(x::AlgAssRelOrdElem; copy::Bool = true) -> AbsAlgAssElem

Returns $x$ as an element of the algebra containing it.
"""
function elem_in_algebra(x::AlgAssRelOrdElem{S, T, U}; copy::Bool = true) where {S, T, U}
  if copy
    return deepcopy(x.elem_in_algebra)::elem_type(U)
  else
    return x.elem_in_algebra::elem_type(U)
  end
end

function elem_in_algebra(x::AlgAssAbsOrdElem{S, T}; copy::Bool = true) where {S, T}
  if copy
    return deepcopy(x.elem_in_algebra)
  else
    return x.elem_in_algebra
  end
end

_elem_in_algebra(x::Union{ AlgAssAbsOrdElem, AlgAssRelOrdElem }; copy::Bool = true) = elem_in_algebra(x, copy = copy)

################################################################################
#
#  "Assure" functions for fields
#
################################################################################

function assure_has_coord(x::Union{ AlgAssAbsOrdElem, AlgAssRelOrdElem })
  if x.has_coord
    return nothing
  end

  a, b = _check_elem_in_order(elem_in_algebra(x, copy = false), parent(x))
  !a && error("Not a valid order element")
  x.coordinates = b
  x.has_coord = true
  return nothing
end

################################################################################
#
#  Coordinates
#
################################################################################

@doc Markdown.doc"""
    coordinates(x::AlgAssAbsOrdElem; copy::Bool = true) -> Vector{fmpz}
    coordinates(x::AlgAssRelOrdElem; copy::Bool = true) -> Vector{NumFieldElem}

Returns the coordinates of $x$ in the basis of `parent(x)`.
"""
function coordinates(x::Union{ AlgAssAbsOrdElem, AlgAssRelOrdElem }; copy::Bool = true)
  assure_has_coord(x)
  if copy
    return deepcopy(x.coordinates)
  else
    return x.coordinates
  end
end

################################################################################
#
#  Unary operations
#
################################################################################

@doc Markdown.doc"""
    -(x::AlgAssAbsOrdElem) -> AlgAssAbsOrdElem
    -(x::AlgAssRelOrdElem) -> AlgAssRelOrdElem

Returns $-x$.
"""
function -(x::Union{ AlgAssAbsOrdElem, AlgAssRelOrdElem })
  return parent(x)(-elem_in_algebra(x, copy = false))
end

###############################################################################
#
#  Binary operations
#
###############################################################################

@doc Markdown.doc"""
    *(x::AlgAssAbsOrdElem, y::AlgAssAbsOrdElem) -> AlgAssAbsOrdElem
    *(x::AlgAssRelOrdElem, y::AlgAssRelOrdElem) -> AlgAssRelOrdElem

Return $x \cdot y$.
"""
function *(x::T, y::T) where { T <: Union{ AlgAssAbsOrdElem, AlgAssRelOrdElem } }
  !check_parent(x, y) && error("Wrong parents")
  return parent(x)(elem_in_algebra(x, copy = false)*elem_in_algebra(y, copy = false))
end

@doc Markdown.doc"""
    +(x::AlgAssAbsOrdElem, y::AlgAssAbsOrdElem) -> AlgAssAbsOrdElem
    +(x::AlgAssRelOrdElem, y::AlgAssRelOrdElem) -> AlgAssRelOrdElem

Return $x + y$.
"""
function +(x::T, y::T) where { T <: Union{ AlgAssAbsOrdElem, AlgAssRelOrdElem } }
  !check_parent(x, y) && error("Wrong parents")
  z = parent(x)(elem_in_algebra(x, copy = false) + elem_in_algebra(y, copy = false))
  if x.has_coord && y.has_coord
    z.coordinates = [ x.coordinates[i] + y.coordinates[i] for i = 1:degree(parent(x)) ]
    z.has_coord = true
  end
  return z
end

@doc Markdown.doc"""
    -(x::AlgAssAbsOrdElem, y::AlgAssAbsOrdElem) -> AlgAssAbsOrdElem
    -(x::AlgAssRelOrdElem, y::AlgAssRelOrdElem) -> AlgAssRelOrdElem

Return $x - y$.
"""
function -(x::T, y::T) where { T <: Union{ AlgAssAbsOrdElem, AlgAssRelOrdElem } }
  !check_parent(x, y) && error("Wrong parents")
  z = parent(x)(elem_in_algebra(x, copy = false) - elem_in_algebra(y, copy = false))
  if x.has_coord && y.has_coord
    z.coordinates = [ x.coordinates[i] - y.coordinates[i] for i = 1:degree(parent(x)) ]
    z.has_coord = true
  end
  return z
end

function *(n::Union{Integer, fmpz}, x::AlgAssAbsOrdElem)
  O=x.parent
  y=Array{fmpz,1}(undef, O.dim)
  z=coordinates(x, copy = false)
  for i=1:O.dim
    y[i] = z[i] * n
  end
  return O(y)
end

*(x::AlgAssAbsOrdElem, n::Union{Integer, fmpz}) = n*x

# Computes a/b if action is :right and b\a if action is :left (and if this is possible)
function divexact(a::T, b::T, action::Symbol, check::Bool = true) where { T <: Union{ AlgAssAbsOrdElem, AlgAssRelOrdElem } }
  !check_parent(a, b) && error("Wrong parents")
  O = parent(a)
  c = divexact(elem_in_algebra(a, copy = false), elem_in_algebra(b, copy = false), action)
  if check
    (x, y) = _check_elem_in_order(c, O)
    !x && error("Quotient not an element of the order")
    return typeof(a)(O, c, y) # Avoid unneccessary copies
  end
  return typeof(a)(O, c)
end

@doc Markdown.doc"""
    divexact_right(a::AlgAssAbsOrdElem, b::AlgAssAbsOrdElem, check::Bool = true)
    divexact_right(a::AlgAssRelOrdElem, b::AlgAssRelOrdElem, check::Bool = true)
      -> AlgAssRelOrdElem

Returns an element $c \in O$ such that $a = c \cdot b$ where $O$ is the order
containing $a$.
If `check` is `false`, it is not checked whether $c$ is an element of $O$.
"""
divexact_right(a::T, b::T, check::Bool = true) where { T <: Union{ AlgAssAbsOrdElem, AlgAssRelOrdElem } } = divexact(a, b, :right, check)

@doc Markdown.doc"""
    divexact_left(a::AlgAssAbsOrdElem, b::AlgAssAbsOrdElem, check::Bool = true)
    divexact_left(a::AlgAssRelOrdElem, b::AlgAssRelOrdElem, check::Bool = true)
      -> AlgAssRelOrdElem

Returns an element $c \in O$ such that $a = b \cdot c$ where $O$ is the order
containing $a$.
If `check` is `false`, it is not checked whether $c$ is an element of $O$.
"""
divexact_left(a::T, b::T, check::Bool = true) where { T <: Union{ AlgAssAbsOrdElem, AlgAssRelOrdElem } } = divexact(a, b, :left, check)

################################################################################
#
#  Conversion from matrix
#
################################################################################

function elem_from_mat_row(O::AlgAssAbsOrd, M::fmpz_mat, i::Int)
  return O(fmpz[ M[i, j] for j = 1:degree(O) ])
end

function elem_to_mat_row!(M::fmpz_mat, i::Int, a::AlgAssAbsOrdElem)
  for c = 1:ncols(M)
    M[i, c] = deepcopy(coordinates(a; copy = false))[c]
  end
  return nothing
end

################################################################################
#
#  Exponentiation
#
################################################################################

@doc Markdown.doc"""
    ^(x::AlgAssAbsOrdElem, y::Union{ Int, fmpz }) -> AlgAssAbsOrdElem
    ^(x::AlgAssRelOrdElem, y::Union{ Int, fmpz }) -> AlgAssRelOrdElem

Returns $x^y$.
"""
function ^(x::Union{ AlgAssAbsOrdElem, AlgAssRelOrdElem }, y::fmpz)
  z = parent(x)()
  z.elem_in_algebra = elem_in_algebra(x, copy = false)^y
  return z
end

function ^(x::Union{ AlgAssAbsOrdElem, AlgAssRelOrdElem }, y::Integer)
  z = parent(x)()
  z.elem_in_algebra = elem_in_algebra(x, copy = false)^y
  return z
end

################################################################################
#
#  Equality
#
################################################################################

@doc Markdown.doc"""
    ==(a::AlgAssAbsOrdElem, b::AlgAssAbsOrdElem) -> Bool
    ==(a::AlgAssRelOrdElem, b::AlgAssRelOrdElem) -> Bool

Returns `true` if $a$ and $b$ are equal and `false` otherwise.
"""
function ==(a::T, b::T) where { T <: Union{ AlgAssAbsOrdElem, AlgAssRelOrdElem } }
  if parent(a) !== parent(b)
    return false
  end
  return elem_in_algebra(a, copy = false) == elem_in_algebra(b, copy = false)
end

################################################################################
#
#  Unsafe operations
#
################################################################################

function add!(z::T, x::T, y::T) where { T <: Union{ AlgAssAbsOrdElem, AlgAssRelOrdElem } }
  z.elem_in_algebra = add!(elem_in_algebra(z, copy = false), elem_in_algebra(x, copy = false), elem_in_algebra(y, copy = false))
  z.has_coord = false
  return z
end

function mul!(z::T, x::T, y::T) where { T <: Union{ AlgAssAbsOrdElem, AlgAssRelOrdElem } }
  z.elem_in_algebra = mul!(elem_in_algebra(z, copy = false), elem_in_algebra(x, copy = false), elem_in_algebra(y, copy = false))
  z.has_coord = false
  return z
end

function mul!(z::AlgAssAbsOrdElem, x::Union{ Int, fmpz }, y::AlgAssAbsOrdElem)
  z.elem_in_algebra = mul!(elem_in_algebra(z, copy = false), x, elem_in_algebra(y, copy = false))
  if z.has_coord && y.has_coord
    x = fmpz(x)
    for i = 1:degree(parent(y))
      z.coordinates[i] = mul!(z.coordinates[i], x, coordinates(y, copy = false)[i])
    end
  end
  return z
end

mul!(z::AlgAssAbsOrdElem, y::AlgAssAbsOrdElem, x::Union{ Int, fmpz }) = mul!(z, x, y)

################################################################################
#
#  String I/O
#
################################################################################

function show(io::IO, a::AlgAssAbsOrdElem)
  print(io, a.elem_in_algebra)
end

################################################################################
#
#  Representation matrices
#
################################################################################

@doc Markdown.doc"""
    representation_matrix(x::AlgAssAbsOrdElem, action::Symbol = :left) -> fmpz_mat

Returns a matrix representing multiplication with $x$ with respect to the basis
of `order(x)`.
The multiplication is from the left if `action == :left` and from the right if
`action == :right`.
"""
function representation_matrix(x::AlgAssAbsOrdElem, action::Symbol = :left)

  O = parent(x)
  M = basis_matrix(O, copy = false)
  M1 = basis_mat_inv(O, copy = false)

  B = FakeFmpqMat(representation_matrix(elem_in_algebra(x, copy = false), action))
  B = mul!(B, M, B)
  B = mul!(B, B, M1)

  @assert B.den == 1
  return B.num
end

function representation_matrix_mod(x::AlgAssAbsOrdElem, d::fmpz, action::Symbol = :left)
  O = parent(x)
  M = basis_matrix(O, copy = false)
  M1 = basis_mat_inv(O, copy = false)

  A = FakeFmpqMat(representation_matrix(elem_in_algebra(x, copy = false), action))
  d2 = M.den * M1.den*A.den
  d2c, d2nc = ppio(d2, d)
  d1 = d * d2c
  A1 = A.num
  mod!(A.num, d1)
  S1 = mod(M.num, d1)
  mul!(A1, S1, A1)
  S2 = mod(M1.num, d1)
  mul!(A1, A1, S2)
  mod!(A1, d1)
  divexact!(A1, A1, d2c)
  inver = invmod(d2nc, d1)
  mul!(A1, A1, inver)
  mod!(A1, d)
  return A1
end

################################################################################
#
#  Trace
#
################################################################################

@doc Markdown.doc"""
    tr(x::AlgAssAbsOrdElem) -> fmpz

Returns the trace of $x$.
"""
function tr(x::AlgAssAbsOrdElem)
  return FlintZZ(tr(x.elem_in_algebra))
end

@doc Markdown.doc"""
    trred(x::AlgAssAbsOrdElem) -> fmpz

Returns the reduced trace of $x$.
"""
function trred(x::AlgAssAbsOrdElem)
  return FlintZZ(trred(x.elem_in_algebra))
end

################################################################################
#
#  Modular exponentiation and division
#
################################################################################

@doc Markdown.doc"""
    powermod(a::AlgAssAbsOrdElem, i::Union{ fmpz, Int }, m::AlgAssAbsOrdIdl)
      -> AlgAssAbsOrdElem

Returns an element $b$ of `order(a)` such that $a^i \equiv b \mod m$.
"""
function powermod(a::AlgAssAbsOrdElem, i::Union{fmpz, Int}, m::AlgAssAbsOrdIdl)
  if i < 0
    b, a = isdivisible_mod_ideal(one(parent(a)), a, m)
    @assert b "Element is not invertible modulo the ideal"
    return powermod(a, -i, m)
  end

  if i == 0
    return one(parent(a))
  end

  if i == 1
    b = mod(a, m)
    return b
  end

  if mod(i, 2) == 0
    j = div(i, 2)
    b = powermod(a, j, m)
    b = b^2
    b = mod(b, m)
    return b
  end

  b = mod(a*powermod(a, i - 1, m), m)
  return b
end

# This is mostly isdivisible in NfOrd/ResidueRing.jl
function isdivisible_mod_ideal(x::AlgAssAbsOrdElem, y::AlgAssAbsOrdElem, a::AlgAssAbsOrdIdl)

  iszero(y) && error("Dividing by zero")

  if iszero(x)
    return true, zero(parent(x))
  end

  O = parent(x)
  d = degree(O)
  V = zero_matrix(FlintZZ, 2*d + 1, 2*d + 1)
  V[1, 1] = fmpz(1)

  for i = 1:d
    V[1, 1 + i] = coordinates(x, copy = false)[i]
  end

  A = representation_matrix(y)
  B = integral_basis_matrix_wrt(a, O, copy = false)

  _copy_matrix_into_matrix(V, 2, 2, A)
  _copy_matrix_into_matrix(V, 2 + d, 2, B)

  for i = 1:d
    V[1 + i, d + 1 + i] = 1
  end

  V = hnf(V)

  for i = 2:(d + 1)
    if !iszero(V[1, i])
      return false, O()
    end
  end

  z = -O([ V[1, i] for i = (d + 2):(2*d + 1) ])
  return true, z
end

################################################################################
#
#  isone/iszero
#
################################################################################

iszero(a::AlgAssAbsOrdElem) = iszero(elem_in_algebra(a, copy = false))

isone(a::AlgAssAbsOrdElem) = isone(elem_in_algebra(a, copy = false))


text:
struct SRA1ConstantCache <: StochasticDiffEqConstantCache end
alg_cache(alg::SRA1,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}}) = SRA1ConstantCache()

@cache struct SRA1Cache{randType,rateType,uType,rateNoiseType} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType
  chi2::randType
  tmp1::uType
  E₁::rateType
  E₂::rateType
  gt::rateNoiseType
  k₁::rateType
  k₂::rateType
  gpdt::rateNoiseType
  tmp::uType
end

function alg_cache(alg::SRA1,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  if typeof(ΔW) <: Union{SArray,Number}
    chi2 = copy(ΔW)
  else
    chi2 = zero(ΔW)
  end
  tmp1 = zero(u)
  E₁ = zero(rate_prototype); gt = zero(noise_rate_prototype); gpdt = zero(noise_rate_prototype)
  E₂ = zero(rate_prototype); k₁ = zero(rate_prototype); k₂ = zero(rate_prototype)
  tmp = zero(u)
  SRA1Cache(u,uprev,chi2,tmp1,E₁,E₂,gt,k₁,k₂,gpdt,tmp)
end

struct SRA2ConstantCache{T,T2} <: StochasticDiffEqConstantCache
  a21::T
  b21::T
  c02::T2
  c11::T2
  c12::T2
  α1::T
  α2::T
  beta12::T
  beta21::T
  beta22::T
end

function SRA2ConstantCache(T::Type, T2::Type)
  a21 = convert(T, 3//4)
  b21 = convert(T, 3//2)
  c02 = convert(T2, 3//4)
  c11 = convert(T2, 1//3)
  c12 = convert(T2, 1)
  α1 = convert(T, 1//3)
  α2 = convert(T, 2//3)
  beta12 = convert(T, 1)
  beta21 = convert(T, 3//2)
  beta22 = convert(T, -3//2)
  SRA2ConstantCache(a21,b21,c02,c11,c12,α1,α2,beta12,beta21,beta22)
end

function alg_cache(alg::SRA2,prob,u,ΔW,ΔZ,p,rate_prototype,
                   noise_rate_prototype,uEltypeNoUnits,
                   uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  SRA2ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
end

@cache struct SRA2Cache{uType,randType,tabType,rateNoiseType,T} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType
  chi2::randType
  tab::tabType
  g1::rateNoiseType
  g2::rateNoiseType
  k1::T
  k2::T
  E₁::T
  E₂::T
  tmp::T
end

function alg_cache(alg::SRA2,prob,u,ΔW,ΔZ,p,rate_prototype,
                   noise_rate_prototype,uEltypeNoUnits,
                   uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  if typeof(ΔW) <: Union{SArray,Number}
    chi2 = copy(ΔW)
  else
    chi2 = zero(ΔW)
  end
  tab = SRA2ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
  g1 = zero(noise_rate_prototype); g2 = zero(noise_rate_prototype)
  k1 = zero(rate_prototype); k2 = zero(rate_prototype)
  E₁ = zero(rate_prototype); E₂ = zero(rate_prototype)
  tmp = k1
  SRA2Cache(u,uprev,chi2,tab,g1,g2,k1,k2,E₁,E₂,tmp)
end

struct ThreeStageSRAConstantCache{T,T2} <: StochasticDiffEqConstantCache
  a21::T
  a31::T
  a32::T
  b21::T
  b31::T
  b32::T
  c02::T2
  c03::T2
  c11::T2
  c12::T2
  c13::T2
  α1::T
  α2::T
  α3::T
  beta11::T
  beta12::T
  beta13::T
  beta21::T
  beta22::T
  beta23::T
end

function SRA3ConstantCache(T::Type, T2::Type)
  a21 = convert(T, 1)
  a31 = convert(T, 1//4)
  a32 = convert(T, 1//4)
  b21 = convert(T, 0)
  b31 = convert(T, 1)
  b32 = convert(T, 1//2)
  c02 = convert(T2, 1)
  c03 = convert(T2, 1//2)
  c11 = convert(T2, 1)
  c12 = convert(T2, 0)
  c13 = convert(T2, 0)
  α1 = convert(T, 1//6)
  α2 = convert(T, 1//6)
  α3 = convert(T, 2//3)
  beta11 = convert(T, 1)
  beta12 = convert(T, 0)
  beta13 = convert(T, 0)
  beta21 = convert(T, -1)
  beta22 = convert(T, 1)
  beta23 = convert(T, 0)
  ThreeStageSRAConstantCache(a21,a31,a32,b21,b31,b32,c02,c03,c11,c12,c13,
                             α1,α2,α3,beta11,beta12,beta13,beta21,beta22,beta23)
end

function SOSRAConstantCache(T::Type, T2::Type)
  α1 = convert(T, 0.2889874966892885)
  α2 = convert(T, 0.6859880440839937)
  α3 = convert(T, 0.025024459226717772)
  c02 = convert(T2, 0.6923962376159507)
  c03 = convert(T2, 1)
  c11 = convert(T2, 0)
  c12 = convert(T2, 0.041248171110700504)
  c13 = convert(T2, 1)
  beta11 = convert(T, -16.792534242221663)
  beta12 = convert(T, 17.514995785380226)
  beta13 = convert(T, 0.27753845684143835)
  beta21 = convert(T, 0.4237535769069274)
  beta22 = convert(T, 0.6010381474428539)
  beta23 = convert(T, -1.0247917243497813)
  a21 = convert(T, 0.6923962376159507)
  a31 = convert(T, -3.1609142252828395)
  a32 = convert(T, 4.1609142252828395)
  b21 = convert(T, 1.3371632704399763)
  b31 = convert(T, 1.442371048468624)
  b32 = convert(T, 1.8632741501139225)
  ThreeStageSRAConstantCache(a21,a31,a32,b21,b31,b32,c02,c03,c11,c12,c13,
                             α1,α2,α3,beta11,beta12,beta13,beta21,beta22,beta23)
end

function SOSRA2ConstantCache(T::Type, T2::Type)
  α1 = convert(T, 0.4999999999999998)
  α2 = convert(T, -0.9683897375354181)
  α3 = convert(T, 1.4683897375354185)
  c02 = convert(T2, 1)
  c03 = convert(T2, 1)
  c11 = convert(T2, 0)
  c12 = convert(T2, 1)
  c13 = convert(T2, 1)
  beta11 = convert(T, 0)
  beta12 = convert(T, 0.92438032145683)
  beta13 = convert(T, 0.07561967854316998)
  beta21 = convert(T, 1)
  beta22 = convert(T, -0.8169981105823436)
  beta23 = convert(T, -0.18300188941765633)
  a21 = convert(T, 1)
  a31 = convert(T, 0.9511849235504364)
  a32 = convert(T, 0.04881507644956362)
  b21 = convert(T, 0.7686101171003622)
  b31 = convert(T, 0.43886792994934987)
  b32 = convert(T, 0.7490415909204886)
  ThreeStageSRAConstantCache(a21,a31,a32,b21,b31,b32,c02,c03,c11,c12,c13,
                             α1,α2,α3,beta11,beta12,beta13,beta21,beta22,beta23)
end

function alg_cache(alg::SRA3,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,
                   uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,
                   f,t,dt,::Type{Val{false}})
  SRA3ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
end

function alg_cache(alg::SOSRA,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,
                   uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,
                   f,t,dt,::Type{Val{false}})
  SOSRAConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
end

function alg_cache(alg::SOSRA2,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,
                   uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,
                   f,t,dt,::Type{Val{false}})
  SOSRA2ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
end

@cache struct ThreeStageSRACache{uType,randType,tabType,rateNoiseType,rateType,GT} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType
  chi2::randType
  tab::tabType
  g1::rateNoiseType
  g2::rateNoiseType
  g3::rateNoiseType
  k1::rateType
  k2::rateType
  k3::rateType
  E₁::rateType
  E₂::rateType
  tmp::rateType
  gtmp::GT
end

function alg_cache(alg::SRA3,prob,u,ΔW,ΔZ,p,rate_prototype,
                   noise_rate_prototype,uEltypeNoUnits,
                   uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  if typeof(ΔW) <: Union{SArray,Number}
    chi2 = copy(ΔW)
  else
    chi2 = zero(ΔW)
  end
  tab = SRA3ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
  g1 = zero(noise_rate_prototype); g2 = zero(noise_rate_prototype)
  g3 = zero(noise_rate_prototype)
  k1 = zero(rate_prototype); k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  E₁ = zero(rate_prototype); E₂ = zero(rate_prototype)
  tmp = k1
  if typeof(noise_rate_prototype) == typeof(rate_prototype)
    gtmp = nothing
  else
    gtmp = zero(noise_rate_prototype)
  end
  ThreeStageSRACache(u,uprev,chi2,tab,g1,g2,g3,k1,k2,k3,E₁,E₂,tmp,gtmp)
end

function alg_cache(alg::SOSRA,prob,u,ΔW,ΔZ,p,rate_prototype,
                   noise_rate_prototype,uEltypeNoUnits,
                   uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  if typeof(ΔW) <: Union{SArray,Number}
    chi2 = copy(ΔW)
  else
    chi2 = zero(ΔW)
  end
  tab = SOSRAConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
  g1 = zero(noise_rate_prototype); g2 = zero(noise_rate_prototype)
  g3 = zero(noise_rate_prototype)
  k1 = zero(rate_prototype); k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  E₁ = zero(rate_prototype); E₂ = zero(rate_prototype)
  tmp = k1
  if typeof(noise_rate_prototype) == typeof(rate_prototype)
    gtmp = nothing
  else
    gtmp = zero(noise_rate_prototype)
  end
  ThreeStageSRACache(u,uprev,chi2,tab,g1,g2,g3,k1,k2,k3,E₁,E₂,tmp,gtmp)
end

function alg_cache(alg::SOSRA2,prob,u,ΔW,ΔZ,p,rate_prototype,
                   noise_rate_prototype,uEltypeNoUnits,
                   uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  if typeof(ΔW) <: Union{SArray,Number}
    chi2 = copy(ΔW)
  else
    chi2 = zero(ΔW)
  end
  tab = SOSRA2ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
  g1 = zero(noise_rate_prototype); g2 = zero(noise_rate_prototype)
  g3 = zero(noise_rate_prototype)
  k1 = zero(rate_prototype); k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  E₁ = zero(rate_prototype); E₂ = zero(rate_prototype)
  tmp = k1
  if typeof(noise_rate_prototype) == typeof(rate_prototype)
    gtmp = nothing
  else
    gtmp = zero(noise_rate_prototype)
  end
  ThreeStageSRACache(u,uprev,chi2,tab,g1,g2,g3,k1,k2,k3,E₁,E₂,tmp,gtmp)
end

struct SRAConstantCache{VType1,VType2,MType,uType} <: StochasticDiffEqConstantCache
  c₀::VType1
  c₁::VType1
  A₀::MType
  B₀::MType
  α::VType2
  β₁::VType2
  β₂::VType2
  stages::Int
  H0::Vector{uType}
end

function SRAConstantCache(tableau,rate_prototype)
  @unpack c₀,c₁,A₀,B₀,α,β₁,β₂ = tableau
  stages = length(α)
  H0 = Vector{typeof(rate_prototype)}(undef,stages)
  SRAConstantCache(c₀,c₁,A₀',B₀',α,β₁,β₂,stages,H0)
end

function alg_cache(alg::SRA,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  SRAConstantCache(alg.tableau,rate_prototype)
end

@cache struct SRACache{uType,rateType,tabType,randType,rateNoiseType} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType
  H0::Vector{uType}
  A0temp::rateType
  B0temp::rateType
  ftmp::rateType
  gtmp::rateNoiseType
  chi2::randType
  atemp::rateType
  btemp::rateType
  E₁::rateType
  E₁temp::rateType
  E₂::rateType
  tmp::uType
  tab::tabType
end

u_cache(c::SRACache) = ()
du_cache(c::SRACache) = (c.A0temp,c.B0temp,c.ftmp,c.gtmp,c.chi2,c.chi2,c.atemp,
                          c.btemp,c.E₁,c.E₁temp,c.E₂)
user_cache(c::SRACache) = (c.u,c.uprev,c.tmp,c.H0...)

function alg_cache(alg::SRA,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  H0 = Vector{typeof(u)}()
  tab = SRAConstantCache(alg.tableau,rate_prototype)
  for i = 1:tab.stages
    push!(H0,zero(u))
  end
  A0temp = zero(rate_prototype); B0temp = zero(rate_prototype)
  ftmp = zero(rate_prototype); gtmp = zero(noise_rate_prototype);
  if typeof(ΔW) <: Union{SArray,Number}
    chi2 = copy(ΔW)
  else
    chi2 = zero(ΔW)
  end
  atemp = zero(rate_prototype); btemp = zero(rate_prototype); E₂ = zero(rate_prototype); E₁temp = zero(rate_prototype)
  E₁ = zero(rate_prototype)
  tmp = zero(u)
  SRACache(u,uprev,H0,A0temp,B0temp,ftmp,gtmp,chi2,atemp,btemp,E₁,E₁temp,E₂,tmp,tab)
end


text:
# precompile is now the default
VERSION < v"0.7.0-rc2" && __precompile__()

module LibSndFile

import SampledSignals
using SampledSignals: SampleSource, SampleSink, SampleBuf
using SampledSignals: PCM16Sample, PCM32Sample
using SampledSignals: nframes, nchannels, samplerate
using FileIO: File, Stream, filename, stream
using FileIO: add_format, add_loader, add_saver, @format_str
if VERSION >= v"0.7-"
    using Printf: @printf
    using LinearAlgebra: transpose!
else
    using Compat: Cvoid, @cfunction, @warn
end

const supported_formats = (format"WAV", format"FLAC", format"OGG")

include("libsndfile_h.jl")
include("lengthIO.jl")
include("sourcesink.jl")
include("loadsave.jl")
include("readwrite.jl")

function __init__()
    # ogg currently not in the registry
    add_format(format"OGG", "OggS", [".ogg", ".oga"], [:LibSndFile])
    for fmt in supported_formats
        add_loader(fmt, :LibSndFile)
        add_saver(fmt, :LibSndFile)
    end
end

depsjl = joinpath(@__DIR__, "..", "deps", "deps.jl")
if isfile(depsjl)
    include(depsjl)
else
    error("LibSndFile not properly installed. Please run Pkg.build(\"LibSndFile\")")
end

end # module LibSndFile