text:
function FieldVector3D(arg0::Field, arg1::Vector3D)
    return FieldVector3D((Field, Vector3D), arg0, arg1)
end

function FieldVector3D(arg0::RealFieldElement, arg1::FieldVector3D)
    return FieldVector3D((RealFieldElement, FieldVector3D), arg0, arg1)
end

function FieldVector3D(arg0::RealFieldElement, arg1::FieldVector3D, arg2::RealFieldElement, arg3::FieldVector3D)
    return FieldVector3D((RealFieldElement, FieldVector3D, RealFieldElement, FieldVector3D), arg0, arg1, arg2, arg3)
end

function FieldVector3D(arg0::RealFieldElement, arg1::FieldVector3D, arg2::RealFieldElement, arg3::FieldVector3D, arg4::RealFieldElement, arg5::FieldVector3D)
    return FieldVector3D((RealFieldElement, FieldVector3D, RealFieldElement, FieldVector3D, RealFieldElement, FieldVector3D), arg0, arg1, arg2, arg3, arg4, arg5)
end

function FieldVector3D(arg0::RealFieldElement, arg1::FieldVector3D, arg2::RealFieldElement, arg3::FieldVector3D, arg4::RealFieldElement, arg5::FieldVector3D, arg6::RealFieldElement, arg7::FieldVector3D)
    return FieldVector3D((RealFieldElement, FieldVector3D, RealFieldElement, FieldVector3D, RealFieldElement, FieldVector3D, RealFieldElement, FieldVector3D), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7)
end

function FieldVector3D(arg0::RealFieldElement, arg1::RealFieldElement)
    return FieldVector3D((RealFieldElement, RealFieldElement), arg0, arg1)
end

function FieldVector3D(arg0::RealFieldElement, arg1::RealFieldElement, arg2::RealFieldElement)
    return FieldVector3D((RealFieldElement, RealFieldElement, RealFieldElement), arg0, arg1, arg2)
end

function FieldVector3D(arg0::RealFieldElement, arg1::Vector3D)
    return FieldVector3D((RealFieldElement, Vector3D), arg0, arg1)
end

function FieldVector3D(arg0::RealFieldElement, arg1::Vector3D, arg2::RealFieldElement, arg3::Vector3D)
    return FieldVector3D((RealFieldElement, Vector3D, RealFieldElement, Vector3D), arg0, arg1, arg2, arg3)
end

function FieldVector3D(arg0::RealFieldElement, arg1::Vector3D, arg2::RealFieldElement, arg3::Vector3D, arg4::RealFieldElement, arg5::Vector3D)
    return FieldVector3D((RealFieldElement, Vector3D, RealFieldElement, Vector3D, RealFieldElement, Vector3D), arg0, arg1, arg2, arg3, arg4, arg5)
end

function FieldVector3D(arg0::RealFieldElement, arg1::Vector3D, arg2::RealFieldElement, arg3::Vector3D, arg4::RealFieldElement, arg5::Vector3D, arg6::RealFieldElement, arg7::Vector3D)
    return FieldVector3D((RealFieldElement, Vector3D, RealFieldElement, Vector3D, RealFieldElement, Vector3D, RealFieldElement, Vector3D), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7)
end

function FieldVector3D(arg0::Vector{RealFieldElement})
    return FieldVector3D((Vector{RealFieldElement},), arg0)
end

function FieldVector3D(arg0::jdouble, arg1::FieldVector3D)
    return FieldVector3D((jdouble, FieldVector3D), arg0, arg1)
end

function FieldVector3D(arg0::jdouble, arg1::FieldVector3D, arg2::jdouble, arg3::FieldVector3D)
    return FieldVector3D((jdouble, FieldVector3D, jdouble, FieldVector3D), arg0, arg1, arg2, arg3)
end

function FieldVector3D(arg0::jdouble, arg1::FieldVector3D, arg2::jdouble, arg3::FieldVector3D, arg4::jdouble, arg5::FieldVector3D)
    return FieldVector3D((jdouble, FieldVector3D, jdouble, FieldVector3D, jdouble, FieldVector3D), arg0, arg1, arg2, arg3, arg4, arg5)
end

function FieldVector3D(arg0::jdouble, arg1::FieldVector3D, arg2::jdouble, arg3::FieldVector3D, arg4::jdouble, arg5::FieldVector3D, arg6::jdouble, arg7::FieldVector3D)
    return FieldVector3D((jdouble, FieldVector3D, jdouble, FieldVector3D, jdouble, FieldVector3D, jdouble, FieldVector3D), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7)
end

function add(obj::FieldVector3D, arg0::FieldVector3D)
    return jcall(obj, "add", FieldVector3D, (FieldVector3D,), arg0)
end

function add(obj::FieldVector3D, arg0::RealFieldElement, arg1::FieldVector3D)
    return jcall(obj, "add", FieldVector3D, (RealFieldElement, FieldVector3D), arg0, arg1)
end

function add(obj::FieldVector3D, arg0::RealFieldElement, arg1::Vector3D)
    return jcall(obj, "add", FieldVector3D, (RealFieldElement, Vector3D), arg0, arg1)
end

function add(obj::FieldVector3D, arg0::Vector3D)
    return jcall(obj, "add", FieldVector3D, (Vector3D,), arg0)
end

function add(obj::FieldVector3D, arg0::jdouble, arg1::FieldVector3D)
    return jcall(obj, "add", FieldVector3D, (jdouble, FieldVector3D), arg0, arg1)
end

function add(obj::FieldVector3D, arg0::jdouble, arg1::Vector3D)
    return jcall(obj, "add", FieldVector3D, (jdouble, Vector3D), arg0, arg1)
end

function angle(::Type{FieldVector3D}, arg0::FieldVector3D, arg1::FieldVector3D)
    return jcall(FieldVector3D, "angle", RealFieldElement, (FieldVector3D, FieldVector3D), arg0, arg1)
end

function angle(::Type{FieldVector3D}, arg0::FieldVector3D, arg1::Vector3D)
    return jcall(FieldVector3D, "angle", RealFieldElement, (FieldVector3D, Vector3D), arg0, arg1)
end

function angle(::Type{FieldVector3D}, arg0::Vector3D, arg1::FieldVector3D)
    return jcall(FieldVector3D, "angle", RealFieldElement, (Vector3D, FieldVector3D), arg0, arg1)
end

function cross_product(::Type{FieldVector3D}, arg0::FieldVector3D, arg1::FieldVector3D)
    return jcall(FieldVector3D, "crossProduct", FieldVector3D, (FieldVector3D, FieldVector3D), arg0, arg1)
end

function cross_product(::Type{FieldVector3D}, arg0::FieldVector3D, arg1::Vector3D)
    return jcall(FieldVector3D, "crossProduct", FieldVector3D, (FieldVector3D, Vector3D), arg0, arg1)
end

function cross_product(::Type{FieldVector3D}, arg0::Vector3D, arg1::FieldVector3D)
    return jcall(FieldVector3D, "crossProduct", FieldVector3D, (Vector3D, FieldVector3D), arg0, arg1)
end

function cross_product(obj::FieldVector3D, arg0::FieldVector3D)
    return jcall(obj, "crossProduct", FieldVector3D, (FieldVector3D,), arg0)
end

function cross_product(obj::FieldVector3D, arg0::Vector3D)
    return jcall(obj, "crossProduct", FieldVector3D, (Vector3D,), arg0)
end

function distance(::Type{FieldVector3D}, arg0::FieldVector3D, arg1::FieldVector3D)
    return jcall(FieldVector3D, "distance", RealFieldElement, (FieldVector3D, FieldVector3D), arg0, arg1)
end

function distance(::Type{FieldVector3D}, arg0::FieldVector3D, arg1::Vector3D)
    return jcall(FieldVector3D, "distance", RealFieldElement, (FieldVector3D, Vector3D), arg0, arg1)
end

function distance(::Type{FieldVector3D}, arg0::Vector3D, arg1::FieldVector3D)
    return jcall(FieldVector3D, "distance", RealFieldElement, (Vector3D, FieldVector3D), arg0, arg1)
end

function distance(obj::FieldVector3D, arg0::FieldVector3D)
    return jcall(obj, "distance", RealFieldElement, (FieldVector3D,), arg0)
end

function distance(obj::FieldVector3D, arg0::Vector3D)
    return jcall(obj, "distance", RealFieldElement, (Vector3D,), arg0)
end

function distance1(::Type{FieldVector3D}, arg0::FieldVector3D, arg1::FieldVector3D)
    return jcall(FieldVector3D, "distance1", RealFieldElement, (FieldVector3D, FieldVector3D), arg0, arg1)
end

function distance1(::Type{FieldVector3D}, arg0::FieldVector3D, arg1::Vector3D)
    return jcall(FieldVector3D, "distance1", RealFieldElement, (FieldVector3D, Vector3D), arg0, arg1)
end

function distance1(::Type{FieldVector3D}, arg0::Vector3D, arg1::FieldVector3D)
    return jcall(FieldVector3D, "distance1", RealFieldElement, (Vector3D, FieldVector3D), arg0, arg1)
end

function distance1(obj::FieldVector3D, arg0::FieldVector3D)
    return jcall(obj, "distance1", RealFieldElement, (FieldVector3D,), arg0)
end

function distance1(obj::FieldVector3D, arg0::Vector3D)
    return jcall(obj, "distance1", RealFieldElement, (Vector3D,), arg0)
end

function distance_inf(::Type{FieldVector3D}, arg0::FieldVector3D, arg1::FieldVector3D)
    return jcall(FieldVector3D, "distanceInf", RealFieldElement, (FieldVector3D, FieldVector3D), arg0, arg1)
end

function distance_inf(::Type{FieldVector3D}, arg0::FieldVector3D, arg1::Vector3D)
    return jcall(FieldVector3D, "distanceInf", RealFieldElement, (FieldVector3D, Vector3D), arg0, arg1)
end

function distance_inf(::Type{FieldVector3D}, arg0::Vector3D, arg1::FieldVector3D)
    return jcall(FieldVector3D, "distanceInf", RealFieldElement, (Vector3D, FieldVector3D), arg0, arg1)
end

function distance_inf(obj::FieldVector3D, arg0::FieldVector3D)
    return jcall(obj, "distanceInf", RealFieldElement, (FieldVector3D,), arg0)
end

function distance_inf(obj::FieldVector3D, arg0::Vector3D)
    return jcall(obj, "distanceInf", RealFieldElement, (Vector3D,), arg0)
end

function distance_sq(::Type{FieldVector3D}, arg0::FieldVector3D, arg1::FieldVector3D)
    return jcall(FieldVector3D, "distanceSq", RealFieldElement, (FieldVector3D, FieldVector3D), arg0, arg1)
end

function distance_sq(::Type{FieldVector3D}, arg0::FieldVector3D, arg1::Vector3D)
    return jcall(FieldVector3D, "distanceSq", RealFieldElement, (FieldVector3D, Vector3D), arg0, arg1)
end

function distance_sq(::Type{FieldVector3D}, arg0::Vector3D, arg1::FieldVector3D)
    return jcall(FieldVector3D, "distanceSq", RealFieldElement, (Vector3D, FieldVector3D), arg0, arg1)
end

function distance_sq(obj::FieldVector3D, arg0::FieldVector3D)
    return jcall(obj, "distanceSq", RealFieldElement, (FieldVector3D,), arg0)
end

function distance_sq(obj::FieldVector3D, arg0::Vector3D)
    return jcall(obj, "distanceSq", RealFieldElement, (Vector3D,), arg0)
end

function dot_product(::Type{FieldVector3D}, arg0::FieldVector3D, arg1::FieldVector3D)
    return jcall(FieldVector3D, "dotProduct", RealFieldElement, (FieldVector3D, FieldVector3D), arg0, arg1)
end

function dot_product(::Type{FieldVector3D}, arg0::FieldVector3D, arg1::Vector3D)
    return jcall(FieldVector3D, "dotProduct", RealFieldElement, (FieldVector3D, Vector3D), arg0, arg1)
end

function dot_product(::Type{FieldVector3D}, arg0::Vector3D, arg1::FieldVector3D)
    return jcall(FieldVector3D, "dotProduct", RealFieldElement, (Vector3D, FieldVector3D), arg0, arg1)
end

function dot_product(obj::FieldVector3D, arg0::FieldVector3D)
    return jcall(obj, "dotProduct", RealFieldElement, (FieldVector3D,), arg0)
end

function dot_product(obj::FieldVector3D, arg0::Vector3D)
    return jcall(obj, "dotProduct", RealFieldElement, (Vector3D,), arg0)
end

function equals(obj::FieldVector3D, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_alpha(obj::FieldVector3D)
    return jcall(obj, "getAlpha", RealFieldElement, ())
end

function get_delta(obj::FieldVector3D)
    return jcall(obj, "getDelta", RealFieldElement, ())
end

function get_minus_i(::Type{FieldVector3D}, arg0::Field)
    return jcall(FieldVector3D, "getMinusI", FieldVector3D, (Field,), arg0)
end

function get_minus_j(::Type{FieldVector3D}, arg0::Field)
    return jcall(FieldVector3D, "getMinusJ", FieldVector3D, (Field,), arg0)
end

function get_minus_k(::Type{FieldVector3D}, arg0::Field)
    return jcall(FieldVector3D, "getMinusK", FieldVector3D, (Field,), arg0)
end

function get_na_n(::Type{FieldVector3D}, arg0::Field)
    return jcall(FieldVector3D, "getNaN", FieldVector3D, (Field,), arg0)
end

function get_negative_infinity(::Type{FieldVector3D}, arg0::Field)
    return jcall(FieldVector3D, "getNegativeInfinity", FieldVector3D, (Field,), arg0)
end

function get_norm(obj::FieldVector3D)
    return jcall(obj, "getNorm", RealFieldElement, ())
end

function get_norm1(obj::FieldVector3D)
    return jcall(obj, "getNorm1", RealFieldElement, ())
end

function get_norm_inf(obj::FieldVector3D)
    return jcall(obj, "getNormInf", RealFieldElement, ())
end

function get_norm_sq(obj::FieldVector3D)
    return jcall(obj, "getNormSq", RealFieldElement, ())
end

function get_plus_i(::Type{FieldVector3D}, arg0::Field)
    return jcall(FieldVector3D, "getPlusI", FieldVector3D, (Field,), arg0)
end

function get_plus_j(::Type{FieldVector3D}, arg0::Field)
    return jcall(FieldVector3D, "getPlusJ", FieldVector3D, (Field,), arg0)
end

function get_plus_k(::Type{FieldVector3D}, arg0::Field)
    return jcall(FieldVector3D, "getPlusK", FieldVector3D, (Field,), arg0)
end

function get_positive_infinity(::Type{FieldVector3D}, arg0::Field)
    return jcall(FieldVector3D, "getPositiveInfinity", FieldVector3D, (Field,), arg0)
end

function get_x(obj::FieldVector3D)
    return jcall(obj, "getX", RealFieldElement, ())
end

function get_y(obj::FieldVector3D)
    return jcall(obj, "getY", RealFieldElement, ())
end

function get_z(obj::FieldVector3D)
    return jcall(obj, "getZ", RealFieldElement, ())
end

function get_zero(::Type{FieldVector3D}, arg0::Field)
    return jcall(FieldVector3D, "getZero", FieldVector3D, (Field,), arg0)
end

function hash_code(obj::FieldVector3D)
    return jcall(obj, "hashCode", jint, ())
end

function is_infinite(obj::FieldVector3D)
    return jcall(obj, "isInfinite", jboolean, ())
end

function is_na_n(obj::FieldVector3D)
    return jcall(obj, "isNaN", jboolean, ())
end

function negate(obj::FieldVector3D)
    return jcall(obj, "negate", FieldVector3D, ())
end

function normalize(obj::FieldVector3D)
    return jcall(obj, "normalize", FieldVector3D, ())
end

function orthogonal(obj::FieldVector3D)
    return jcall(obj, "orthogonal", FieldVector3D, ())
end

function scalar_multiply(obj::FieldVector3D, arg0::RealFieldElement)
    return jcall(obj, "scalarMultiply", FieldVector3D, (RealFieldElement,), arg0)
end

function scalar_multiply(obj::FieldVector3D, arg0::jdouble)
    return jcall(obj, "scalarMultiply", FieldVector3D, (jdouble,), arg0)
end

function subtract(obj::FieldVector3D, arg0::FieldVector3D)
    return jcall(obj, "subtract", FieldVector3D, (FieldVector3D,), arg0)
end

function subtract(obj::FieldVector3D, arg0::RealFieldElement, arg1::FieldVector3D)
    return jcall(obj, "subtract", FieldVector3D, (RealFieldElement, FieldVector3D), arg0, arg1)
end

function subtract(obj::FieldVector3D, arg0::RealFieldElement, arg1::Vector3D)
    return jcall(obj, "subtract", FieldVector3D, (RealFieldElement, Vector3D), arg0, arg1)
end

function subtract(obj::FieldVector3D, arg0::Vector3D)
    return jcall(obj, "subtract", FieldVector3D, (Vector3D,), arg0)
end

function subtract(obj::FieldVector3D, arg0::jdouble, arg1::FieldVector3D)
    return jcall(obj, "subtract", FieldVector3D, (jdouble, FieldVector3D), arg0, arg1)
end

function subtract(obj::FieldVector3D, arg0::jdouble, arg1::Vector3D)
    return jcall(obj, "subtract", FieldVector3D, (jdouble, Vector3D), arg0, arg1)
end

function to_array(obj::FieldVector3D)
    return jcall(obj, "toArray", Vector{RealFieldElement}, ())
end

function to_string(obj::FieldVector3D)
    return jcall(obj, "toString", JString, ())
end

function to_string(obj::FieldVector3D, arg0::NumberFormat)
    return jcall(obj, "toString", JString, (NumberFormat,), arg0)
end

function to_vector3d(obj::FieldVector3D)
    return jcall(obj, "toVector3D", Vector3D, ())
end


text:
@doc raw"""
    conjugate_gradient_descent(M, F, ∇F, x)

perform a conjugate gradient based descent
````math
x_{k+1} = \operatorname{retr}_{x_k} \bigl( s_k\delta_k \bigr),
````
where $\operatorname{retr}$ denotes a retraction on the `Manifold` `M`
and one can employ different rules to update the descent direction $\delta_k$ based on
the last direction $\delta_{k-1}$ and both gradients $\nabla f(x_k)$,$\nabla f(x_{k-1})$.
The [`Stepsize`](@ref) $s_k$ may be determined by a [`Linesearch`](@ref).

Available update rules are [`SteepestDirectionUpdateRule`](@ref), which yields a [`steepest_descent`](@ref),
[`ConjugateDescentCoefficient`](@ref), [`DaiYuanCoefficient`](@ref), [`FletcherReevesCoefficient`](@ref),
[`HagerZhangCoefficient`](@ref), [`HeestenesStiefelCoefficient`](@ref),
[`LiuStoreyCoefficient`](@ref), and [`PolakRibiereCoefficient`](@ref).

They all compute $\beta_k$ such that this algorithm updates the search direction as
````math
\delta_k=\nabla f(x_k) + \beta_k \delta_{k-1}
````

# Input
* `M` : a manifold $\mathcal M$
* `F` : a cost function $F\colon\mathcal M\to\mathbb R$ to minimize
* `∇F`: the gradient $∇ F\colon\mathcal M\to T\mathcal M$ of F
* `x` : an initial value $x\in\mathcal M$

# Optional
* `coefficient` : ([`SteepestDirectionUpdateRule`](@ref) `<:` [`DirectionUpdateRule`](@ref)
  rule to compute the descent direction update coefficient $\beta_k$,
  as a functor, i.e. the resulting function maps `(p,o,i) -> β`, where
  `p` is the current [`GradientProblem`](@ref), `o` are the
  [`ConjugateGradientDescentOptions`](@ref) `o` and `i` is the current iterate.
* `retraction_method` - (`ExponentialRetraction`) a retraction method to use, by default the exponntial map
* `return_options` – (`false`) – if actiavated, the extended result, i.e. the
    complete [`Options`](@ref) re returned. This can be used to access recorded values.
    If set to false (default) just the optimal value `xOpt` if returned
* `stepsize` - (`Constant(1.)`) A [`Stepsize`](@ref) function applied to the
  search direction. The default is a constant step size 1.
* `stopping_criterion` : (`stopWhenAny( stopAtIteration(200), stopGradientNormLess(10.0^-8))`)
  a function indicating when to stop.
* `vector_transport_method` – (`ParallelTransport()`) vector transport method to transport
  the old descent direction when computing the new descent direction.

# Output
* `xOpt` – the resulting (approximately critical) point of gradientDescent
OR
* `options` - the options returned by the solver (see `return_options`)
"""
function conjugate_gradient_descent(
    M::Manifold,
    F::Function,
    ∇F::Function,
    x;
    coefficient::DirectionUpdateRule = SteepestDirectionUpdateRule(),
    stepsize::Stepsize = ConstantStepsize(1.),
    retraction_method::AbstractRetractionMethod = ExponentialRetraction(),
    stopping_criterion::StoppingCriterion = StopWhenAny(
        StopAfterIteration(500), StopWhenGradientNormLess(10^(-8))
    ),
    transport_method = ParallelTransport(),
    return_options=false,
    kwargs...
)
    p = GradientProblem(M,F,∇F)
    o = ConjugateGradientDescentOptions(
        x,
        stopping_criterion,stepsize,
        coefficient,
        retraction_method,
        transport_method,
        )
    o = decorate_options(o; kwargs...)
    resultO = solve(p,o)
    if return_options
        return resultO
    end
    return get_solver_result(resultO)
end
function initialize_solver!(p::GradientProblem,o::ConjugateGradientDescentOptions)
    o.∇ = get_gradient(p,o.x)
    o.δ = -o.∇
    o.β = 0.
end
function step_solver!(p::GradientProblem, o::ConjugateGradientDescentOptions, i)
    xOld = o.x
    o.x = retract(p.M, o.x, get_stepsize(p, o, i, o.δ) * o.δ, o.retraction_method)
    o.∇ = get_gradient(p, o.x)
    o.β = o.coefficient(p,o,i)
    o.δ = -o.∇ + o.β * vector_transport_to(p.M, xOld, o.δ, o.x, o.vector_transport_method)
end
get_solver_result(o::ConjugateGradientDescentOptions) = o.x

ext: jl


text:
#=
Table of quadrature points and weights for the reference tet.

values adapted from:

Lee Shunn, Frank Ham,
Symmetric quadrature rules for tetrahedra based on a cubic close-packed lattice arrangement,
Journal of Computational and Applied Mathematics,
Volume 236, Issue 17,
2012,
Pages 4348-4364,
ISSN 0377-0427,
https://doi.org/10.1016/j.cam.2012.03.032.

=#

function tetrahedron_quadrature_nodes_weights(order)

    if order > 6 || order < 1
        printerr("unsupported order for tet quadrature nodes. Choose 1 to 5.")
        return zeros(0,0);
    end

    cubTet = Array{Array{Float64,2},1}(undef,6);

    cubTet[1] = [
        -0.5 -0.5 -0.5 2.0
    ];

    cubTet[2] = [
        -0.7236067977499786 -0.7236067977499785 -0.7236067977499785 0.5;
        0.17082039324993548 -0.7236067977499785 -0.7236067977499785 0.5;
        -0.7236067977499785 0.17082039324993545 -0.7236067977499785 0.5;
        -0.7236067977499784 -0.7236067977499785 0.1708203932499354 0.5
    ];

    cubTet[3] = [
        -0.8523301965475533 -0.8523301965475534 -0.8523301965475534 0.0952662696864178;
        0.5569905896426599 -0.8523301965475534 -0.8523301965475534 0.0952662696864178;
        -0.8523301965475534 0.5569905896426598 -0.8523301965475533 0.0952662696864178;
        -0.8523301965475534 -0.8523301965475533 0.5569905896426598 0.0952662696864178;
        -0.1875113122318981 -0.8124886877681019 -0.8124886877681019 0.269822486875722;
        -0.8124886877681019 -0.18751131223189815 -0.8124886877681019 0.269822486875722;
        -0.8124886877681017 -0.812488687768102 -0.1875113122318981 0.269822486875722;
        -0.18751131223189807 -0.18751131223189815 -0.8124886877681019 0.269822486875722;
        -0.18751131223189804 -0.812488687768102 -0.1875113122318981 0.269822486875722;
        -0.8124886877681019 -0.18751131223189815 -0.1875113122318981 0.269822486875722
    ];

    cubTet[4] = [
        -0.9352948105455121 -0.935294810545512 -0.9352948105455121 0.014134149588939;
        0.8058844316365361 -0.935294810545512 -0.9352948105455121 0.014134149588939;
        -0.935294810545512 0.8058844316365361 -0.9352948105455121 0.014134149588939;
        -0.935294810545512 -0.9352948105455121 0.8058844316365362 0.014134149588939;
        0.23319306612387392 -0.879279116949716 -0.8792791169497158 0.0939973379437754;
        -0.4746348322244419 -0.879279116949716 -0.8792791169497158 0.0939973379437754;
        -0.8792791169497159 0.23319306612387394 -0.8792791169497158 0.0939973379437754;
        -0.8792791169497156 -0.47463483222444225 -0.8792791169497159 0.0939973379437754;
        -0.8792791169497158 -0.8792791169497159 0.23319306612387386 0.0939973379437754;
        -0.8792791169497156 -0.8792791169497159 -0.47463483222444214 0.0939973379437754;
        -0.47463483222444214 0.23319306612387394 -0.8792791169497158 0.0939973379437754;
        0.23319306612387397 -0.47463483222444225 -0.8792791169497159 0.0939973379437754;
        -0.4746348322244421 -0.8792791169497159 0.23319306612387386 0.0939973379437754;
        0.23319306612387392 -0.8792791169497159 -0.47463483222444214 0.0939973379437754;
        -0.8792791169497158 -0.4746348322244422 0.23319306612387386 0.0939973379437754;
        -0.8792791169497158 0.2331930661238738 -0.4746348322244421 0.0939973379437754;
        -0.38046139145427543 -0.3804613914542755 -0.8586158256371735 0.203873836579736;
        -0.3804613914542754 -0.8586158256371736 -0.3804613914542755 0.203873836579736;
        -0.8586158256371735 -0.38046139145427543 -0.3804613914542755 0.203873836579736;
        -0.3804613914542755 -0.38046139145427543 -0.3804613914542755 0.203873836579736
    ];

    cubTet[5] = [
        -0.946526448891253 -0.9465264488912534 -0.9465264488912533 0.0043800927930776;
        0.8395793466737598 -0.9465264488912534 -0.9465264488912533 0.0043800927930776;
        -0.9465264488912533 0.8395793466737599 -0.9465264488912533 0.0043800927930776;
        -0.9465264488912531 -0.9465264488912533 0.8395793466737596 0.0043800927930776;
        0.49551977696361765 -0.9217955187287027 -0.9217955187287028 0.028679134035533;
        -0.6519287395062123 -0.9217955187287027 -0.9217955187287028 0.028679134035533;
        -0.9217955187287027 0.49551977696361776 -0.9217955187287028 0.028679134035533;
        -0.9217955187287027 -0.6519287395062123 -0.9217955187287028 0.028679134035533;
        -0.9217955187287028 -0.9217955187287027 0.49551977696361765 0.028679134035533;
        -0.9217955187287028 -0.9217955187287027 -0.6519287395062123 0.028679134035533;
        -0.6519287395062122 0.49551977696361776 -0.9217955187287028 0.028679134035533;
        0.4955197769636177 -0.6519287395062123 -0.9217955187287028 0.028679134035533;
        -0.6519287395062123 -0.9217955187287027 0.49551977696361765 0.028679134035533;
        0.49551977696361765 -0.9217955187287027 -0.6519287395062123 0.028679134035533;
        -0.9217955187287027 -0.6519287395062123 0.49551977696361765 0.028679134035533;
        -0.9217955187287027 0.49551977696361776 -0.6519287395062123 0.028679134035533;
        -0.09049080003103432 -0.9095091999689656 -0.9095091999689658 0.0500610791373492;
        -0.9095091999689657 -0.09049080003103421 -0.9095091999689658 0.0500610791373492;
        -0.9095091999689657 -0.9095091999689658 -0.09049080003103428 0.0500610791373492;
        -0.09049080003103424 -0.09049080003103421 -0.9095091999689658 0.0500610791373492;
        -0.09049080003103421 -0.9095091999689658 -0.09049080003103428 0.0500610791373492;
        -0.9095091999689657 -0.09049080003103421 -0.09049080003103428 0.0500610791373492;
        -0.5535979240753701 -0.55359792407537 -0.899041441878456 0.0959678666115108;
        0.006237290029196096 -0.55359792407537 -0.899041441878456 0.0959678666115108;
        -0.55359792407537 0.006237290029196047 -0.899041441878456 0.0959678666115108;
        -0.5535979240753701 -0.8990414418784561 -0.55359792407537 0.0959678666115108;
        0.006237290029196152 -0.8990414418784561 -0.55359792407537 0.0959678666115108;
        -0.5535979240753701 -0.8990414418784561 0.006237290029196064 0.0959678666115108;
        -0.899041441878456 -0.55359792407537 -0.55359792407537 0.0959678666115108;
        -0.8990414418784561 0.006237290029196047 -0.55359792407537 0.0959678666115108;
        -0.8990414418784562 -0.55359792407537 0.006237290029196064 0.0959678666115108;
        0.006237290029196041 -0.55359792407537 -0.55359792407537 0.0959678666115108;
        -0.55359792407537 0.006237290029196047 -0.55359792407537 0.0959678666115108;
        -0.5535979240753699 -0.55359792407537 0.006237290029196064 0.0959678666115108;
        -0.5 -0.5 -0.5 0.186349146239068
    ];

    cubTet[6] = [
        -0.9700958696938813 -0.9700958696938815 -0.9700958696938814 0.002074622467228;
        0.9102876090816443 -0.9700958696938815 -0.9700958696938814 0.002074622467228;
        -0.9700958696938815 0.9102876090816445 -0.9700958696938814 0.002074622467228;
        -0.9700958696938815 -0.9700958696938816 0.9102876090816444 0.002074622467228;
        -0.6963361016681259 -0.9318079576074771 -0.931807957607477 0.019203329079896;
        0.55995201688308 -0.9318079576074771 -0.931807957607477 0.019203329079896;
        -0.9318079576074769 -0.696336101668126 -0.931807957607477 0.019203329079896;
        -0.931807957607477 0.55995201688308 -0.931807957607477 0.019203329079896;
        -0.9318079576074768 -0.9318079576074771 -0.696336101668126 0.019203329079896;
        -0.9318079576074769 -0.931807957607477 0.55995201688308 0.019203329079896;
        0.55995201688308 -0.696336101668126 -0.931807957607477 0.019203329079896;
        -0.696336101668126 0.55995201688308 -0.931807957607477 0.019203329079896;
        0.55995201688308 -0.9318079576074771 -0.696336101668126 0.019203329079896;
        -0.696336101668126 -0.931807957607477 0.55995201688308 0.019203329079896;
        -0.931807957607477 0.55995201688308 -0.696336101668126 0.019203329079896;
        -0.9318079576074769 -0.696336101668126 0.55995201688308 0.019203329079896;
        0.10531128621203428 -0.9075896991699963 -0.9075896991699962 0.0328987953596464;
        -0.2901318878720418 -0.9075896991699963 -0.9075896991699962 0.0328987953596464;
        -0.9075896991699962 0.10531128621203427 -0.9075896991699962 0.0328987953596464;
        -0.9075896991699961 -0.2901318878720418 -0.9075896991699962 0.0328987953596464;
        -0.9075896991699963 -0.9075896991699964 0.10531128621203434 0.0328987953596464;
        -0.9075896991699961 -0.9075896991699963 -0.29013188787204175 0.0328987953596464;
        -0.29013188787204175 0.10531128621203427 -0.9075896991699962 0.0328987953596464;
        0.10531128621203428 -0.2901318878720418 -0.9075896991699962 0.0328987953596464;
        -0.29013188787204175 -0.9075896991699964 0.10531128621203434 0.0328987953596464;
        0.10531128621203434 -0.9075896991699963 -0.29013188787204175 0.0328987953596464;
        -0.9075896991699963 -0.29013188787204175 0.10531128621203434 0.0328987953596464;
        -0.9075896991699963 0.10531128621203434 -0.29013188787204175 0.0328987953596464;
        -0.5436190778624788 -0.5436190778624789 -0.9889704900510458 0.030749553302662;
        0.07620864577600328 -0.5436190778624789 -0.9889704900510458 0.030749553302662;
        -0.5436190778624788 0.07620864577600327 -0.9889704900510458 0.030749553302662;
        -0.5436190778624788 -0.9889704900510458 -0.5436190778624788 0.030749553302662;
        0.07620864577600323 -0.9889704900510458 -0.5436190778624788 0.030749553302662;
        -0.5436190778624787 -0.9889704900510459 0.07620864577600328 0.030749553302662;
        -0.9889704900510456 -0.5436190778624788 -0.5436190778624788 0.030749553302662;
        -0.9889704900510458 0.07620864577600334 -0.5436190778624788 0.030749553302662;
        -0.9889704900510456 -0.5436190778624788 0.07620864577600328 0.030749553302662;
        0.07620864577600328 -0.5436190778624788 -0.5436190778624788 0.030749553302662;
        -0.5436190778624789 0.07620864577600334 -0.5436190778624788 0.030749553302662;
        -0.5436190778624788 -0.5436190778624788 0.07620864577600328 0.030749553302662;
        -0.29538947982401254 -0.29538947982401254 -0.8015885595010944 0.058704023675046;
        -0.6076324808508804 -0.29538947982401254 -0.8015885595010944 0.058704023675046;
        -0.2953894798240124 -0.6076324808508806 -0.8015885595010944 0.058704023675046;
        -0.2953894798240124 -0.8015885595010945 -0.29538947982401254 0.058704023675046;
        -0.6076324808508804 -0.8015885595010945 -0.29538947982401254 0.058704023675046;
        -0.2953894798240124 -0.8015885595010945 -0.6076324808508805 0.058704023675046;
        -0.8015885595010943 -0.29538947982401254 -0.29538947982401254 0.058704023675046;
        -0.8015885595010943 -0.6076324808508805 -0.29538947982401254 0.058704023675046;
        -0.8015885595010943 -0.2953894798240126 -0.6076324808508805 0.058704023675046;
        -0.6076324808508804 -0.29538947982401254 -0.29538947982401254 0.058704023675046;
        -0.2953894798240125 -0.6076324808508805 -0.29538947982401254 0.058704023675046;
        -0.29538947982401254 -0.2953894798240126 -0.6076324808508805 0.058704023675046;
        -0.7310433304140116 -0.7310433304140114 -0.7310433304140115 0.0732582732810216;
        0.1931299912420344 -0.7310433304140114 -0.7310433304140115 0.0732582732810216;
        -0.7310433304140115 0.1931299912420344 -0.7310433304140115 0.0732582732810216;
        -0.7310433304140116 -0.7310433304140115 0.19312999124203437 0.0732582732810216
    ];

    return cubTet[order];
end


text:
using Plots
using Random
Random.seed!(1)

include_ddp()

# Model
include_model("double_integrator")

function f(model::DoubleIntegratorContinuous, x, u, w)
    [x[2]; w[1] * u[1]]
end

function fd(model::DoubleIntegratorContinuous{Midpoint, FixedTime}, x, u, w, h, t)
	if t == 1
		w = u[2:2]
    	return [x + h * f(model, x + 0.5 * h * f(model, x, view(u, 1:1), w), view(u, 1:1), w); w]
	else
		w = x[3:3]
		return [view(x, 1:2) + h * f(model, view(x, 1:2) + 0.5 * h * f(model, view(x, 1:2), u, w), u, w); w]
	end
end

# Time
T = 11
h = 0.1

model = DoubleIntegratorContinuous{Midpoint, FixedTime}(2, 1, 1)
n = [t == 1 ? model.n : model.n + 1 for t = 1:T]
m = [t == 1 ? model.m + 1 : model.m for t = 1:T]

# Initial conditions, controls, disturbances
x1 = [1.0; 0.0]
xT = [0.0; 0.0]
ū = [t == 1 ? [1.0 * randn(1); 1.0] : 1.0 * randn(1) for t = 1:T-1]
w = [zeros(model.d) for t = 1:T-1]

# Rollout
x̄ = rollout(model, x1, ū, w, h, T)

# Objective
Q = [(t < T ? h : 1.0) * (t > 1 ? Diagonal([1.0; 1.0; 0.0]) : Diagonal([1.0; 1.0])) for t = 1:T]
q = [-2.0 * Q[t] * (t > 1 ? [xT; 1.0] : xT) for t = 1:T]
R = h * [Diagonal(1.0 * (t == 1 ? ones(2) : ones(1))) for t = 1:T-1]
r = [zeros(m[t]) for t = 1:T-1]
obj = StageCosts([QuadraticCost(Q[t], q[t],
	t < T ? R[t] : nothing, t < T ? r[t] : nothing) for t = 1:T], T)

function g(obj::StageCosts, x, u, t)
	T = obj.T
    if t < T
		Q = obj.cost[t].Q
		q = obj.cost[t].q
	    R = obj.cost[t].R
		r = obj.cost[t].r
        return x' * Q * x + q' * x + u' * R * u + r' * u
    elseif t == T
		Q = obj.cost[T].Q
		q = obj.cost[T].q
        return x' * Q * x + q' * x
    else
        return 0.0
    end
end

# Constraints
p = [t == 1 ? 2 + 2 * model.m : (t == T ? 2 : 2 * model.m) for t = 1:T]
ul = [-5.0]
uu = [5.0]
info_1 = Dict(:pl => [0.0], :pu => [5.0],
 	:ul => ul, :uu => uu, :inequality => (1:(2 + 2 * model.m)))
info_t = Dict(:ul => ul, :uu => uu, :inequality => (1:(2 * model.m)))
info_T = Dict(:xT => xT)
con_set = [StageConstraint(p[t], t < T ? (t == 1 ? info_1 : info_t) : info_T) for t = 1:T]

function c!(c, cons::StageConstraints, x, u, t)
	T = cons.T
	if t == 1
		c[1:1] = view(u, 1:1) - cons.con[t].info[:uu]
		c[2:2] = cons.con[t].info[:ul] - view(u, 1:1)
		c[3:3] = view(u, 2:2) - cons.con[t].info[:pu]
		c[4:4] = cons.con[t].info[:pl] - view(u, 2:2)
	elseif t < T
		c[1:1] = view(u, 1:1) - cons.con[t].info[:uu]
		c[2:2] = cons.con[t].info[:ul] - view(u, 1:1)
	end
	if t == T
		c[1:2] .= view(x, 1:2) - cons.con[t].info[:xT]
	end
end

prob = problem_data(model, obj, con_set, copy(x̄), copy(ū), w, h, T, n = n, m = m)

# Solve
@time constrained_ddp_solve!(prob,
    max_iter = 1000, max_al_iter = 6,
	ρ_init = 1.0, ρ_scale = 10.0)

x, u = current_trajectory(prob)
x̄, ū = nominal_trajectory(prob)

plt = plot(hcat([xt[1:2] for xt in x]...)',
	width = 2.0,
	color = [:cyan :orange],
	label = ["x" "ẋ"],
	xlabel = "time step",
	ylabel = "state")

savefig(plt,
	joinpath("/home/taylor/Research/parameter_optimization_manuscript/figures/di_codesign_state.png"))

plt = plot(hcat([ut[1] for ut in u]..., u[end][1])',
	width = 2.0,
	color = :magenta,
	linetype = :steppost,
	xlabel = "time step",
	ylabel = "control",
	label = "")

savefig(plt,
	joinpath("/home/taylor/Research/parameter_optimization_manuscript/figures/di_codesign_control.png"))

@show u[1][2]

max_issues_repo_name: thowell/DirectMotionPlanning


text:
# This file includes testset for SDESystem

import DifferentialEquations: LambaEM

@testset "SDESystemTestSet" begin
    @info "Running SDESystemTestSet ..."

    # SDESystem construction
    f(dx, x, u, t) = (dx[1] = -x[1])
    h(dx, x, u, t) = (dx[1] = -x[1])
    g(x, u, t) = x
    ds = SDESystem(drift=f, diffusion=h, readout=g, state=[1.], input=nothing, output=Outport(1), alg=LambaEM{true}())
    ds = SDESystem(drift=f, diffusion=h, readout=g, state=[1.], input=nothing, output=Outport(1))
    @test typeof(ds.trigger) == Inpin{Float64}
    @test typeof(ds.handshake) == Outpin{Bool}
    @test ds.input === nothing
    @test isa(ds.output, Outport)
    @test length(ds.output) == 1
    @test ds.integrator.sol.prob.p === nothing

    # Driving SDESystem
    iport = Inport(1)
    trg = Outpin()
    hnd = Inpin{Bool}()
    connect!(ds.output, iport)
    connect!(trg, ds.trigger)
    connect!(ds.handshake, hnd)
    tsk = launch(ds)
    tsk2 = @async while true
        all(take!(iport) .=== NaN) && break
    end
    for t in 1. : 10.
        put!(trg, t)
        take!(hnd)
        @test ds.t == t
        @test [read(pin.link.buffer) for pin in iport] == ds.state
    end
    put!(trg, NaN)
    sleep(0.1)
    @test istaskdone(tsk)
    put!(ds.output, NaN * ones(length(ds.output)))
    sleep(0.1)
    @test istaskdone(tsk2)

    # SDESystem with input
    function f2(dx, x, u, t)
        dx[1] = -x[1] + sin(u[2](t))
        dx[2] = -x[2] + cos(u[1](t))
        dx[3] = -x[3] + cos(u[1](t))
    end
    function h2(dx, x, u, t)
        dx[1] = -x[1]
        dx[2] = -x[2] + cos(u[1](t))
        dx[3] = -x[3] + cos(u[2](t))
    end
    g2(x, u, t) = x
    ds = SDESystem(drift=f2, diffusion=h2, readout=g2, state=ones(3), input=Inport(2), output=Outport(3))
    oport = Outport(2)
    iport = Inport(3)
    trg = Outpin()
    hnd = Inpin{Bool}()
    connect!(oport, ds.input)
    connect!(ds.output, iport)
    connect!(trg, ds.trigger)
    connect!(ds.handshake, hnd)
    tsk = launch(ds)
    tsk2 = @async while true
        all(take!(iport) .=== NaN) && break
    end
    for t in 1 : 10
        put!(trg, t)
        put!(oport, [t, 2t])
        take!(hnd)
        @test ds.t == t
        @test [read(pin.link.buffer) for pin in iport] == ds.state
    end
    put!(trg, NaN)
    sleep(0.1)
    @test istaskdone(tsk)
    put!(ds.output, NaN * ones(length(ds.output)))
    sleep(0.1)
    @test istaskdone(tsk2)

    # Test defining new SDESystem types
    # The type must be mutable
    @test_throws Exception @eval @def_sde_system struct SDESystem{DR, DF, RO, ST, IP, OP}
        drift::DR
        diffusion::DF
        readout::RO
        state::ST
        input::IP
        output::OP
    end

    # The type must be subtype of AbstractSDESystem
    @test_throws Exception @eval @def_sde_system mutable struct SDESystem{DR, DF, RO, ST, IP, OP}
        drift::DR
        diffusion::DF
        readout::RO
        state::ST
        input::IP
        output::OP
    end

    # The type must be subtype of AbstractSDESystem
    @test_throws Exception @eval @def_sde_system mutable struct SDESystem{DR, DF, RO, ST, IP, OP} <: MyDummyAbstractSDESystem
        drift::DR
        diffusion::DF
        readout::RO
        state::ST
        input::IP
        output::OP
    end

    @info "Running SDESystemTestSet ..."
end  # testset


text:
# Utility function to count the properly colored edges in a graph coloring (i.e. endpoints have different colors)
function properly_colored_edges(graph::Graph, coloring::Vector{Int})
    length(coloring) == graph.n || throw(ArgumentError("Length of coloring must equal the number of vertices."))
    return count(coloring[a] != coloring[b] for (a, b) ∈ graph.edges)
end

# Utility function to compute the probabilities of the outcomes represented by a wavefunction ψ, sorted by descending probability
function wavefunction_distribution(ψ::Vector{ComplexF64}; as_bitstrings::Bool = true,
        include_zero = false)::Union{Vector{Tuple{Int, Float64}}, Vector{Tuple{Vector{Int}, Float64}}}
    distribution = [(i-1, abs(amplitude)^2) for (i, amplitude) ∈ enumerate(ψ) if abs(amplitude) > 0 || include_zero]
    if as_bitstrings
        N = Int(log2(length(ψ)))
        distribution = [(digits(i, base=2, pad=N) |> reverse, p) for (i, p) ∈ distribution]
    end

    return sort(distribution, by=(t -> -t[2]))
end

# Utility function to compute the probabilities of the outcome wavefunction of a circuit applied to ψ, sorted by descending probability
function output_distribution(circ::Circuit, ψ::Vector{ComplexF64}; as_bitstrings::Bool = true,
        include_zero = false)::Union{Vector{Tuple{Int, Float64}}, Vector{Tuple{Vector{Int}, Float64}}}
    ψ_out = apply(ψ, circ.moments)
    return wavefunction_distribution(ψ_out, as_bitstrings=as_bitstrings, include_zero=include_zero)
end

# Utility function to compute the probabilities of the output colorings of a circuit applied to ψ, sorted by descending probability
function output_colorings_distribution(circ::Circuit{N}, n::Int;
        include_zero = false)::Vector{Tuple{Dict{Int, Vector{Int}}, Float64}} where {N}
    κ = Int(N / n)
    ψ = ψ_initial(n, κ)
    ψ_out = apply(ψ, circ.moments)
    distribution = wavefunction_distribution(ψ_out, as_bitstrings=false, include_zero=include_zero)

    return [(decode_basis_state(state, n, κ), p) for (state, p) ∈ distribution]
end

# Utility function to compute the probabilities of the output colorings of a circuit applied to ψ, sorted by descending probability
function output_colorings_distribution_scored(circ::Circuit{N}, graph::Graph;
        include_zero = false)::Vector{Tuple{Vector{Int}, Int, Float64}} where {N}
    dist = output_colorings_distribution(circ, graph.n, include_zero = include_zero)

    coloring_dict_to_list(dict) = [length(dict[i]) == 1 ? dict[i][1] :
        throw(ArgumentError("At least one vertex has not one unique color.")) for i ∈ 1:graph.n]

    return [(coloring_dict_to_list(coloring), properly_colored_edges(graph, coloring_dict_to_list(coloring)), p) for (coloring, p) ∈ dist]
end

# Decodes the coloring represented by a single computational basis state, represented as integer
function decode_basis_state(basis_state::Int, n::Int, κ::Int)::Dict{Int, Vector{Int}}
    (n > 0 && κ > 0) || throw(DomainError("Parameters n and κ must be positive integers."))

    N = n * κ
    bits = digits(basis_state, base=2, pad=N) |> reverse
    vertex_bits = vertex -> bits[((vertex - 1) * κ + 1):(vertex * κ)]
    colors_by_vertex = Dict([(v, findall(!iszero, vertex_bits(v))) for v ∈ 1:n])

    all(!isnothing, colors_by_vertex) || throw(ArgumentError("The state `basis_state` has at least one vertex without a color."))

    return colors_by_vertex
end

# Prototypical log function. Implementations should have the @Zygote.nograd attribute and use
# a data structure as first argument which can save log data. params is ::Zygote.Params.
function log_qaoa(logger::T, round::Int, ψ_out::Vector{ComplexF64}, objective::Float64, params) where {T}
    throw(ArgumentError("The type `$(T)` does not implement the `log_qaoa`
        function and is not a valid logger."))
end

function plot_circuit_output_distribution(circ::Circuit, ψ)
    data = output_distribution(circ, ψ, as_bitstrings=true, include_zero=false)
    sort!(data, by=d -> d[1])
    bar(data .|> (d -> "|$(join(d[1]))⟩"),
                            data .|> (d -> d[2]),
                        ylims=(0, 1), xrotation = 60, color="#ff8080", x_ticks=:all,
                        bottom_margin=5Plots.mm,
                        xtickfontsize=5, legend=:topleft,
                        label = "outcome probabilities"
                    )
end


text:
module RegisteredFixturesTests

using Fixtures
using FactCheck

facts("Registered Fixtures tests") do
  context("test add_fixture adds to fixtures") do
    patch(Fixtures, :fixtures, Dict()) do
      function somefunction(fn::Function)
        fn()
      end

      @fact Fixtures.fixtures --> Dict()
      add_fixture(:somecontext, somefunction)
      @fact length(Fixtures.fixtures)-->1
      @fact Fixtures.fixtures[:somecontext][1].name --> nothing
      @fact Fixtures.fixtures[:somecontext][1].fn --> anything_of_type(Function)
    end
  end

  context("test add_fixture with single fixture function") do
    patch(Fixtures, :fixtures, Dict()) do
      global x = 0

      function fixture_func(fn::Function)
        global x = 100
        try
          fn()
        finally
          global x = 200
        end
      end

      @fact Fixtures.fixtures --> Dict()
      add_fixture(:somecontext, fixture_func)
      @fact x-->0
      apply_fixtures(:somecontext) do
        @fact x-->100
      end
      @fact x-->200
    end
  end

  context("test add_simple_fixture with separate setup and teardown functions") do
    patch(Fixtures, :fixtures, Dict()) do
      global x = 0

      function setup_func()
        global x = 100
      end

      function teardown_func()
        global x = 200
      end

      @fact Fixtures.fixtures --> Dict()
      add_simple_fixture(:somecontext, setup_func, teardown_func)
      @fact x-->0
      apply_fixtures(:somecontext) do
        @fact x-->100
      end
      @fact x-->200
    end
  end

  context("Test that fixtures are scoped") do
    patch(Fixtures, :fixtures, Dict()) do
      global x = 0

      function fixture_func(fn::Function)
        global x = 100
        try
          fn()
        finally
          global x = 200
        end
      end

      apply_fixtures(:toplevel) do
        @fact Fixtures.fixtures --> Dict()
        add_fixture(:somecontext, fixture_func)

        @fact Fixtures.fixtures --> not(Dict())

        @fact x-->0
        apply_fixtures(:somecontext) do
          @fact x-->100
        end
        @fact x-->200
      end

      @fact Fixtures.fixtures --> Dict()
    end
  end

  context("A 1-part-function fixture", using_fixtures) do
    global x = 0

    function changex(fn::Function)
      global x = 100
      fn()
    end

    add_fixture(:examplescope, changex)

    @fact x --> 0
    apply_fixtures(:examplescope) do
      @fact x --> 100
    end
    @fact x --> 100
  end

  context("A 2-part-function fixture", using_fixtures) do
    global x = 0

    function changex(fn::Function)
      global x = 100
      try
        fn()
      finally
        global x = 200
      end
    end

    add_fixture(:examplescope, changex)

    @fact x --> 0
    apply_fixtures(:examplescope) do
      @fact x --> 100
    end
    @fact x --> 200
  end

  context("combine two 2-part-function fixtures", using_fixtures) do
    global x = 0
    global y = 0

    function changex(fn::Function)
      global x = 100
      try
        fn()
      finally
        global x = 200
      end
    end

    function changey(fn::Function)
      global y = 100
      try
        fn()
      finally
        global y = 200
      end
    end

    add_fixture(:examplescope, changex)
    add_fixture(:examplescope, changey)

    @fact x --> 0
    @fact y --> 0
    apply_fixtures(:examplescope) do
      @fact x --> 100
      @fact y --> 100
    end
    @fact x --> 200
    @fact y --> 200
  end

  context("combine two 2-part-function fixtures check order", using_fixtures) do
    global x = ""

    function appenda(fn::Function)
      global x = x*"a"
      try
        fn()
      finally
        global x = x*"A"
      end
    end

    function appendb(fn::Function)
      global x = x*"b"
      try
        fn()
      finally
        global x = x*"B"
      end
    end

    add_fixture(:examplescope, appenda)
    add_fixture(:examplescope, appendb)

    @fact x --> ""
    apply_fixtures(:examplescope) do
      @fact x --> "ab"
    end
    @fact x --> "abBA"
  end

  context("combine 1-part-function and 2-part-function fixtures", using_fixtures) do
    global x = 0
    global y = 0

    function changex(fn::Function)
      global x = 100
      fn()
    end

    function changey(fn::Function)
      global y = 100
      try
        fn()
      finally
        global y = 200
      end
    end

    add_fixture(:examplescope, changex)
    add_fixture(:examplescope, changey)

    @fact x --> 0
    @fact y --> 0
    apply_fixtures(:examplescope) do
      @fact x --> 100
      @fact y --> 100
    end
    @fact x --> 100
    @fact y --> 200
  end

  context("A fixture with no values", using_fixtures) do
    function nadda(fn::Function)
      fn()
    end

    add_fixture(:examplescope, :nowt, nadda)

    apply_fixtures(:examplescope, fixture_values=true) do values
      @fact collect(keys(values)) --> [:nowt]
      @fact values[:nowt] --> nothing
    end
  end

  context("A fixture with one value", using_fixtures) do
    function produce_n(fn::Function, n::Integer)
      fn(n)
    end

    add_fixture(:examplescope, :one_hundred, produce_n, 100)

    apply_fixtures(:examplescope, fixture_values=true) do values
      @fact collect(keys(values)) --> [:one_hundred]
      @fact values[:one_hundred] --> 100
    end
  end

  context("A fixture with one 0-tuple value", using_fixtures) do
    function produce_empty(fn::Function)
      fn(())
    end

    add_fixture(:examplescope, :empty, produce_empty)

    apply_fixtures(:examplescope, fixture_values=true) do values
      @fact collect(keys(values)) --> [:empty]
      @fact values[:empty] --> ()
    end
  end

  context("A fixture with one 1-tuple value", using_fixtures) do
    function produce_n(fn::Function, n::Integer)
      fn((n,))
    end

    add_fixture(:examplescope, :one_hundred, produce_n, 100)

    apply_fixtures(:examplescope, fixture_values=true) do values
      @fact collect(keys(values)) --> [:one_hundred]
      @fact values[:one_hundred] --> (100,)
    end
  end

  context("A fixture with one 2-tuple value", using_fixtures) do
    function produce_n(fn::Function, n::Any)
      fn(n)
    end

    add_fixture(:examplescope, :some_nums, produce_n, (100, 200))

    apply_fixtures(:examplescope, fixture_values=true) do values
      @fact collect(keys(values)) --> [:some_nums]
      @fact values[:some_nums] --> (100, 200)
    end
  end

  context("A fixture with two values", using_fixtures) do
    function produce_n_and_m(fn::Function, n::Integer, m::Integer)
      fn(n, m)
    end

    add_fixture(:examplescope, :some_nums, produce_n_and_m, 100, 200)

    apply_fixtures(:examplescope, fixture_values=true) do values
      @fact collect(keys(values)) --> [:some_nums]
      @fact values[:some_nums] --> (100, 200)
    end
  end

  context("A fixture with three values", using_fixtures) do
    function produce_n_m_and_o(fn::Function, n::Integer, m::Integer, o::Integer)
      fn(n, m, o)
    end

    add_fixture(:examplescope, :some_nums, produce_n_m_and_o, 100, 200, 300)

    apply_fixtures(:examplescope, fixture_values=true) do values
      @fact collect(keys(values)) --> [:some_nums]
      @fact values[:some_nums] --> (100, 200, 300)
    end
  end

  context("An anonymous fixture with one value", using_fixtures) do
    function produce_n(fn::Function, n::Integer)
      fn(n)
    end

    add_fixture(:examplescope, produce_n, 100)

    apply_fixtures(:examplescope, fixture_values=true) do values
      @fact collect(keys(values)) --> []
    end
  end

  context("Two fixture with one value each", using_fixtures) do
    function produce_n(fn::Function, n::Integer)
      fn(n)
    end

    add_fixture(:examplescope, :peter, produce_n, 100)
    add_fixture(:examplescope, :paul, produce_n, 200)

    apply_fixtures(:examplescope, fixture_values=true) do values
      @fact sort(collect(keys(values))) --> sort([:peter, :paul])
      @fact values[:peter] --> 100
      @fact values[:paul] --> 200
    end
  end

  context("Two fixture with the same name", using_fixtures) do
    function produce_n(fn::Function, n::Integer)
      fn(n)
    end

    name = :same_name

    add_fixture(:examplescope, name, produce_n, 100)
    @fact_throws add_fixture(:examplescope, name, produce_n, 100)
  end
end

end


text:
"""
    valid_reduce2(A, fun)

Reduce a 2d matrix across its columns using a given function, but ignoring
NaN, Inf, and -Inf.

# Arguments

- A = 2d array of floats
- fun = function with which to reduce across the columns of A

# Value

2d array of floats

"""
function valid_reduce2(A::Array{Float64,3}, fun::Function=mean)

    # Initialize array for storing output
    out = Array{Float64, 2}(undef, size(A, 1), size(A, 2))
    # Iterate through rows of A and reduce across columns
    for i in 1:size(A, 1), j in 1:size(A, 2)
        # Note that the following line also drops NaNs, which are not numbers
        out[i, j] = fun(A[i, j, (A[i, j, :] .< Inf) .& (A[i, j, :] .> -Inf)])
    end
    return out
end

"""
    calc_avg_mse(MLMNet_cv)

Calculates average test MSE across folds.

# Arguments

- MLMNet_cv = MLMNet_cv object

# Value

2d array of floats

"""
function calc_avg_mse(MLMNet_cv::Mlmnet_cv)

    # return valid_reduce2(MLMNet_cv.mse, mean)[:, :, 1]
    return valid_reduce2(MLMNet_cv.mse, mean)
end

"""
    calc_avg_prop_zero(MLMNet_cv)

Calculates average proportion of zero interaction coefficients across folds.

# Arguments

- MLMNet_cv = MLMNet_cv object

# Value

1d array of floats

"""
function calc_avg_prop_zero(MLMNet_cv::Mlmnet_cv)

    # return valid_reduce2(MLMNet_cv.propZero, mean)[:,:,1]
    return valid_reduce2(MLMNet_cv.propZero, mean)
end

"""
    mlmnet_cv_summary(MLMNet_cv)

Summarizes results of cross-validation by returning a table with:

- Lambdas used
- Average MSE across folds for each lambda
- Average proportion of zero interaction coefficeints across folds for each
  lambda

# Arguments

- MLMNet_cv = MLMNet_cv object

# Value

DataFrame summarizing average MSE and proportion of zero interactions across
folds for each lambda.

"""
function mlmnet_cv_summary(MLMNet_cv::Mlmnet_cv)

    # Calculate summary information across folds

    lenVecLambdas = length(MLMNet_cv.lambdas)
    lenVecAlphas = length(MLMNet_cv.alphas)

    λαIndex =  map((x,y) -> (x,y),
            repeat(collect(1:lenVecLambdas), lenVecAlphas),
            vec(repeat(permutedims(collect(1:lenVecAlphas)), lenVecLambdas)) )

    # Second column
    # a 2-d matrix of average MSEs w.r.p.t each (l1, l2)
    avg_mse = vec(calc_avg_mse(MLMNet_cv))

    # Third column
    mseStd = vec(valid_reduce2(MLMNet_cv.mse, std))

    # Fourth column
    # a 2-d matrix of average zero proportions w.r.p.t each (l1, l2)
    avg_prop_zero = vec(calc_avg_prop_zero(MLMNet_cv))

    out_df = DataFrame(𝜆_𝛼_Index = λαIndex,
                       Lambda = repeat(MLMNet_cv.lambdas, lenVecAlphas),
                       Alpha =  vec(repeat(permutedims(MLMNet_cv.alphas), lenVecLambdas)),
                       AvgMSE = avg_mse,
                       StdMSE = mseStd,
                       AvgPercentZero =  avg_prop_zero);

    return out_df
end

"""
    lambda_min(MLMNet_cv)

Returns summary information for lambdas corresponding to the minimum average
test MSE across folds and the MSE one that is standard error greater.

# Arguments

- MLMNet_cv = MLMNet_cv object

# Value

DataFrame from mlmnet_cv_summary restricted to the lambdas and alphas that correspond to
the minimum average test MSE across folds and the MSE that is one standard
error greater.

"""
function lambda_min(MLMNet_cv::Mlmnet_cv)
    # Calculate average proportion of zeros
    prop_zeroMean = calc_avg_prop_zero(MLMNet_cv)

    # Calculate average test MSE across folds.
    mseMean = calc_avg_mse(MLMNet_cv)
    mseStd = valid_reduce2(MLMNet_cv.mse, std)

    # Find index of minimum average test MSE
    minIdx = argmin(mseMean)[1]
    minIdy = argmin(mseMean)[2]

    # Compute standard error across folds for the minimum MSE
    mse1StdErr = mseMean[minIdx, minIdy] + mseStd[minIdx, minIdy]
    # Find the index of the lambda that is closest to being 1 SE greater than
    # the lowest lambda, in the direction of the bigger lambdas
    min1StdErrIdx = argmin(abs.(mseMean[1:minIdx[1], 1:minIdy[1]].-mse1StdErr))[1]
    min1StdErrIdy = argmin(abs.(mseMean[1:minIdx[1], 1:minIdy[1]].-mse1StdErr))[2]

    # Pull out summary information for these two lambdas
    out = hcat(MLMNet_cv.lambdas[minIdx], MLMNet_cv.alphas[minIdy],
               mseMean[minIdx, minIdy], prop_zeroMean[minIdx, minIdy])
    out2 = hcat(MLMNet_cv.lambdas[min1StdErrIdx], MLMNet_cv.alphas[min1StdErrIdy],
                mseMean[min1StdErrIdx, min1StdErrIdy],
                prop_zeroMean[min1StdErrIdx, min1StdErrIdy])
    out = DataFrame(vcat(out, out2), :auto)

    colnames = ["Lambda", "Alpha", "AvgMSE", "AvgPercentZero"];
    rename!(out, Symbol.(colnames))

    # Add names first column
    insertcols!(out, 1, :Name => ["(𝜆, 𝛼)_min", "(𝜆, 𝛼)_min1se"])
    # Add indices as second column
    insertcols!(out, 2, :Index => [Tuple([minIdx, minIdy]), Tuple([min1StdErrIdx, min1StdErrIdy])])

    return out
end


text:
# ---
# title: 1025. Divisor Game
# id: problem1025
# author: Indigo
# date: 2022-04-14
# difficulty: Easy
# categories: Math, Dynamic Programming
# link: <https://leetcode.com/problems/divisor-game/description/>
# hidden: true
# ---
#
# Alice and Bob take turns playing a game, with Alice starting first.
#
# Initially, there is a number `N` on the chalkboard.  On each player's turn,
# that player makes a _move_  consisting of:
#
#   * Choosing any `x` with `0 < x < N` and `N % x == 0`.
#   * Replacing the number `N` on the chalkboard with `N - x`.
#
# Also, if a player cannot make a move, they lose the game.
#
# Return `True` if and only if Alice wins the game, assuming both players play
# optimally.
#
#
#
# **Example 1:**
#
#
#
#     Input: 2
#     Output: true
#     Explanation: Alice chooses 1, and Bob has no more moves.
#
#
# **Example 2:**
#
#
#
#     Input: 3
#     Output: false
#     Explanation: Alice chooses 1, Bob chooses 1, and Alice has no more moves.
#
#
#
#
# **Note:**
#
#   1. `1 <= N <= 1000`
#
#
## @lc code=start
using LeetCode

divisor_game(n::Int) = iseven(n)
## @lc code=end


text:
## MATLAB interfcae to find_zero
## Main functions are
## * fzero(f, ...) to find _a_ zero of f, a univariate function
## * fzeros(f, ...) to attempt to find all zeros of f, a univariate function
## unlike `find_zero` these do not specialize on f, so
## will be faster the first use, and slower for subsequent uses (for the same f)
struct FnWrapper
    f
end
(F::FnWrapper)(x::Number) = first(F.f(x))

"""
    fzero(f, x0; order=0; kwargs...)
    fzero(f, x0, M; kwargs...)
    fzero(f, x0, M, N; kwargs...)
    fzero(f, x0; kwargs...)
    fzero(f, a::Number, b::Number; kwargs...)
    fzero(f, a::Number, b::Number; order=?, kwargs...)
    fzero(f, fp, a::Number; kwargs...)

Find zero of a function using one of several iterative algorithms.

* `f`: a scalar function or callable object

* `x0`: an initial guess, a scalar value or tuple of two values

* `order`: An integer, symbol, or string indicating the algorithm to
   use for `find_zero`. The `Order0` default may be specified directly
   by `order=0`, `order=:0`, or `order="0"`; `Order1()` by `order=1`,
   `order=:1`, `order="1"`, or `order=:secant`; `Order1B()` by
   `order="1B"`, etc.

* `M`: a specific method, as would be passed to `find_zero`, bypassing
  the use of the `order` keyword

* `N`: a specific bracketing method. When given, if a bracket is
  identified, method `N` will be used to finish instead of method `M`.

* `a`, `b`: When two values are passed along, if no `order` value is
  specified, `Bisection` will be used over the bracketing interval
  `(a,b)`. If an `order` value is specified, the value of `x0` will be set to
  `(a,b)` and the specified method will be used.

* `fp`: when `fp` is specified (assumed to compute the derivative of `f`),
  Newton's method will be used

* `kwargs...`: See `find_zero` for the specification of tolerances and other keyword arguments

Examples:
```
fzero(sin, 3)                  # use Order0() method, the default
fzero(sin, 3, order=:secant)   # use secant method (also just `order=1`)
fzero(sin, 3, Roots.Order1B()) # use secant method variant for multiple roots.
fzero(sin, 3, 4)               # use bisection method over (3,4)
fzero(sin, 3, 4, xatol=1e-6)   # use bisection method until |x_n - x_{n-1}| <= 1e-6
fzero(sin, 3, 3.1, order=1)    # use secant method with x_0=3.0, x_1 = 3.1
fzero(sin, (3, 3.1), order=2)  # use Steffensen's method with x_0=3.0, x_1 = 3.1
fzero(sin, cos, 3)             # use Newton's method
```

!!! note
    Unlike `find_zero`, `fzero` does not specialize on the type of the function argument.
    This has the advantage of making the first use of the function `f` faster, but subsequent uses slower.

"""
function fzero(f, x0::Number; kwargs...)
    x = float(x0)
    isinf(x) && throw(ConvergenceFailed("An initial value must be finite"))
    derivative_free(f, x; kwargs...)
end

function fzero(f, x0, M::AbstractUnivariateZeroMethod; kwargs...)
    find_zero(FnWrapper(f), x0, M; kwargs...)
end

function fzero(f, x0, M::AbstractUnivariateZeroMethod, N::AbstractBracketing; kwargs...)
    find_zero(FnWrapper(f), x0, M, N; kwargs...)
end

function fzero(f, bracket::Tuple{T,S}; kwargs...)  where {T <: Number, S<:Number}
    d = Dict(kwargs)
    if haskey(d, :order)
        find_zero(FnWrapper(f), bracket, _method_lookup[d[:order]]; kwargs...)
    else
        find_zero(FnWrapper(f), bracket, Bisection();kwargs...)
    end
end

fzero(f, a::Number, b::Number,args...; kwargs...) = fzero(f, (a,b), args...; kwargs...)
fzero(f, bracket::Vector{T}, args...; kwargs...)  where {T <: Number} = fzero(f,(bracket[1],bracket[2]),args...; kwargs...)

fzero(f::Function, fp::Function, x0::Real; kwargs...) = find_zero((f,fp), x0, Newton(); kwargs...)

# match fzero up with find_zero
_method_lookup = Dict(0   => Order0(),
                     :0  => Order0(),
                     "0" => Order0(),

                      1    => Order1(),
                      :1   => Order1(),
                      "1"  => Order1(),
                      :secant => Order1(),
                      :Secant => Order1(),
                      "1B" => Order1B(),
                      :king => Order1B(),
                      :King => Order1B(),

                     2    => Order2(),
                     :2   => Order2(),
                     :steffensen   => Order2(),
                     :Steffensen   => Order2(),
                     "2"  => Order2(),
                     "2B" => Order2B(),
                     :esser   => Order2B(),
                     :Esser   => Order2B(),

                      5  => Order5(),
                     :5  => Order5(),
                     "5" => Order5(),

                     8   => Order8(),
                     :8  => Order8(),
                     "8" => Order8(),

                     16   => Order16(),
                     :16  => Order16(),
                     "16" => Order16(),
)

@noinline function derivative_free(f, x0; order=0, kwargs...)

    if haskey(_method_lookup, order)
        M = _method_lookup[order]
    else
        warn("Invalid order specified. See ?fzero.")
        throw(ArgumentError())
    end

    # d = (kv[1] == :ftol ? :atol=>kv[2] :
    #      kv[1] == :ftolrel ? :rtol=>kv[2] :
    #      kv[1] == :xtol ? :xatol=>kv[2] :
    #      kv[1] == :xtolrel ? xrtol=>kv[2] :
    #      kv[1] => kv[1] for kv in kwargs)

    d = Dict(kwargs)
     for (o, n) in ((:ftol, :atol), (:ftolrel, :rtol),
                    (:xtol, :xatol), (:xtolrel, :xrtol))
         if haskey(d, o)
             d[n] = d[o]
         end
     end

    find_zero(FnWrapper(f), x0, M; d...)
end

## fzeros
"""

`fzeros(f, a, b; kwargs...)`

Searches for all zeros of `f` within an interval `(a,b)`. Assume neither `a` or `b` is a zero.

Dispatches to `find_zeros(f, a, b; kwargs...)`.
"""
function fzeros(f, a::Number, b::Number; kwargs...)
    find_zeros(FnWrapper(f), float(a), float(b); kwargs...)
end
fzeros(f, bracket::Vector{T}; kwargs...) where {T <: Number} = fzeros(f, bracket[1], bracket[2]; kwargs...)
fzeros(f, bracket::Tuple{T,S}; kwargs...) where {T <: Number, S<:Number} = fzeros(f, bracket[1], bracket[2]; kwargs...)


text:
export simplify_constants

function simplify_constants(O::Operation, shorten_tree)
    while true
        O′ = _simplify_constants(O, shorten_tree)
        if is_operation(O′)
            O′ = Operation(O′.op, simplify_constants.(O′.args, shorten_tree))
        end
        isequal(O, O′) && return O
        O = O′
    end
end
simplify_constants(x, shorten_tree) = x
simplify_constants(x) = simplify_constants(x, true)

const AC_OPERATORS = (*, +)

function _simplify_constants(O::Operation, shorten_tree)
    # Tree shrinking
    if shorten_tree && O.op ∈ AC_OPERATORS
        # Flatten tree
        idxs = findall(x -> is_operation(x) && x.op === O.op, O.args)
        if !isempty(idxs)
            keep_idxs = eachindex(O.args) .∉ (idxs,)
            args = Vector{Expression}[O.args[i].args for i in idxs]
            push!(args, O.args[keep_idxs])
            return Operation(O.op, vcat(args...))
        end

        # Collapse constants
        idxs = findall(is_constant, O.args)
        if length(idxs) > 1
            other_idxs = eachindex(O.args) .∉ (idxs,)
            new_const = Constant(mapreduce(get, O.op, O.args[idxs]))
            args = push!(O.args[other_idxs], new_const)

            length(args) == 1 && return first(args)
            return Operation(O.op, args)
        end
    end

    if O.op === (*)
        # If any variable is `Constant(0)`, zero the whole thing
        any(iszero, O.args) && return Constant(0)

        # If any variable is `Constant(1)`, remove that `Constant(1)` unless
        # they are all `Constant(1)`, in which case simplify to a single variable
        if any(isone, O.args)
            args = filter(!isone, O.args)

            isempty(args)     && return Constant(1)
            length(args) == 1 && return first(args)
            return Operation(O.op, args)
        end

        return O
    end

    if O.op === (+) && any(iszero, O.args)
        # If there are Constant(0)s in a big `+` expression, get rid of them
        args = filter(!iszero, O.args)

        isempty(args)     && return Constant(0)
        length(args) == 1 && return first(args)
        return Operation(O.op, args)
    end

    (O.op, length(O.args)) === (identity, 1) && return O.args[1]

    (O.op, length(O.args)) === (-, 1) && return Operation(*, Expression[-1, O.args[1]])

    return O
end
_simplify_constants(x, shorten_tree) = x
_simplify_constants(x) = _simplify_constants(x, true)


text:
module econ628ps7

using ForwaDiff, LinearAlgebra

function newton(f, f_prime, x_0; tolerance = 1E-7, maxiter = 1000)
    #Setup the algorithm
    x_old = x_0
    error = Inf
    iter = 1
    while error > tolerance && iter <= maxiter
        x_new = x_old - (f(x_old)/f_prime(x_old))
        error = sqrt((x_new - x_old)^2)
        x_old = x_new
        iter = iter + 1
    end
    if iter == maxiter + 1
        return nothing
    else
        return x_old
end
end

f(x) = (x-1)^3
f_prime(x) = 3*(x-1)^2

function newton(f, x_0; tolerance = 1E-7, maxiter = 1000)
f_prime = x -> ForwardDiff.derivative(f, x)
return newton(f,f_prime,x_0;tolerance = tolerance, maxiter = maxiter)
end

export newton, f, f_prime
end


text:
using FinancialDerivatives
using Test

eu_put = EuropeanOption(100.0, 90.0, 0.05, 0.3, 180/365, -1)
eu_call = EuropeanOption(100.0, 90.0, 0.05, 0.3, 180/365, 1)

am_put = AmericanOption(100.0, 90.0, 0.05, 0.3, 180/365, -1)
am_call = AmericanOption(100.0, 90.0, 0.05, 0.3, 180/365, 1)

@testset "Black-Scholes" begin
    @test isapprox(evaluate(eu_put, BlackScholes()), 3.22, atol=0.1)
    @test isapprox(evaluate(eu_call, BlackScholes()), 15.42, atol=0.1)
end

@testset "Cox-Ross-Rubinstein" begin
    @test isapprox(evaluate(am_put, CoxRossRubinstein()), 3.22, atol=0.1)
    @test isapprox(evaluate(am_call, CoxRossRubinstein()), 15.42, atol=0.1)
end

@testset "Jarrow-Rudd" begin
    @test isapprox(evaluate(am_put, JarrowRudd()), 3.22, atol=0.25)
    @test isapprox(evaluate(am_call, JarrowRudd()), 15.42, atol=0.25)
    @test isapprox(evaluate(am_put, JarrowRudd(), false), 3.22, atol=0.25)
    @test isapprox(evaluate(am_call, JarrowRudd(), false), 15.42, atol=0.25)
end

@testset "Longstaff-Schwartz" begin
    @test isapprox(evaluate(am_put, LongstaffSchwartz()), 3.22, atol=1)
    @test isapprox(evaluate(am_call, LongstaffSchwartz()), 15.42, atol=1)
end

@testset "Tian" begin
    @test isapprox(evaluate(am_put, Tian()), 5.15, atol=0.25)
    @test isapprox(evaluate(am_call, Tian()), 11.87, atol=0.25)
end

@testset "Garman–Kohlhagen" begin
    fxp = FXOption(100.0, 90.0, 0.05, 0.025, 0.3, 180/365, -1)
    fxc = FXOption(100.0, 90.0, 0.05, 0.025, 0.3, 180/365, 1)
    @test isapprox(evaluate(fxp, GarmanKohlhagen()), 3.51, atol=0.25)
    @test isapprox(evaluate(fxc, GarmanKohlhagen()), 14.48, atol=0.25)
end

ird = InterestRateDerivative(0.01875, 0.20, 0.01, 0.012, 180/365)

@testset "Black-Karasinski" begin
    @test isapprox(evaluate(ird, BlackKarasinski(), 2), [0.2, 0.2, 0.2, 0.2], atol=0.25)
end

@testset "Brennan-Schwartz" begin
    @test isapprox(evaluate(ird, BrennanSchwartz(), 2), [0.2, 0.2, 0.2, 0.2], atol=0.25)
end

@testset "Rendleman-Bartter" begin
    @test isapprox(evaluate(am_put, RendlemanBartter()), 3.22, atol=0.25)
    @test isapprox(evaluate(am_call, RendlemanBartter()), 15.42, atol=0.25)
    @test isapprox(evaluate(ird, RendlemanBartter(), 2), [0.2, 0.2, 0.2, 0.2], atol=0.025)
end

@testset "Vasicek" begin
    @test isapprox(evaluate(ird, Vasicek(), 2), [0.2, 0.2, 0.2, 0.2], atol=0.025)
end


text:
using CircuitComponentRounding
using Documenter

DocMeta.setdocmeta!(CircuitComponentRounding, :DocTestSetup, :(using CircuitComponentRounding); recursive=true)

makedocs(;
    modules=[CircuitComponentRounding],
    authors="KronosTheLate",
    repo="https://github.com/KronosTheLate/CircuitComponentRounding.jl/blob/{commit}{path}#{line}",
    sitename="CircuitComponentRounding.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://KronosTheLate.github.io/CircuitComponentRounding.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/KronosTheLate/CircuitComponentRounding.jl",
)


text:
# Autogenerated wrapper script for Help2man_jll for any
export help2man

JLLWrappers.@generate_wrapper_header("Help2man")
JLLWrappers.@declare_file_product(help2man)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_file_product(
        help2man,
        "bin/help2man",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
export ClassicalAMGsetup;
include("coloring.jl")
include("interpolation.jl")

function ClassicalAMGsetup(AT::SparseMatrixCSC{VAL,IND},param::MGparam{VAL,IND},
							symm::Bool = true,nrhs::Int64 = 1,verbose::Bool=false) where {VAL,IND}

Ps = Array{SparseMatrixCSC{real(VAL),IND}}(undef, param.levels-1);
Rs = Array{SparseMatrixCSC{real(VAL),IND}}(undef, param.levels-1);
As = Array{SparseMatrixCSC{VAL,IND}}(undef, param.levels);

if symm == false
	error("not supported yet...");
end

N = size(AT,2);
As[1] = AT;
relaxPrecs = Array{Vector{VAL}}(undef, param.levels-1);
Cop = nnz(AT);
time = 0;
for l = 1:(param.levels-1)
	if verbose
		time = time_ns();
	end
	AT = As[l];
	if param.relaxType=="Jac" || param.relaxType=="Jac-GMRES" || param.relaxType=="SPAI"
		relaxPrecs[l] = getRelaxPrec(AT,param.relaxType,param.relaxParam,[],false);
	else
		error("Unknown relaxation type !!!!");
	end

	S = getStrengthMatrix(AT,param.strongConnParam, N);
	#println("Time of first coloring pass:")
	coloring = getColoringFirst(S, N)
	checkColoring(coloring,N)
	#println("Time of second coloring pass:")
	coloring = getColoringSecond(S, coloring, N)
	checkColoring(coloring,N)
	#println("Time of interpolation:")
	(P, PT) = getInterpolation(AT, S, coloring, N)
	Nc = size(P,2);
	## P should also be row-wise, so we hold P^T.
	if (size(P,1)==size(P,2))
		if verbose; println(string("Stopped Coarsening at level ",l)); end
		param.levels = l;
		As = As[1:l];
		Ps = Ps[1:l-1];
		Rs = Rs[1:l-1];
		relaxPrecs = relaxPrecs[1:l-1];
		break;
	else
		Rs[l] = P; # this is becasue we hold the transpose of the matrices and P = R' anyway here....
		Ps[l] = PT;

		Act = Ps[l]*AT*Rs[l];
		As[l+1] = Act;
		Cop = Cop + nnz(Act);

		if verbose; println("MG setup: ",N," took:", (time_ns()-time)/1e+9); end;
		N = Nc;
	end
end
if verbose
	time = time_ns();
	println("MG Setup: Operator complexity = ",Cop/nnz(As[1]));
end
As[end] = As[end] + 1e-8*norm(As[end],1)*sparse(1.0I,size(As[end],2),size(As[end],2));
defineCoarsestAinv(param,As[end]);

if verbose
	println("MG setup coarsest: ",N,", done LU in ",(time_ns()-time)/1e+9);
end
param.As = As;
param.Ps = Ps;
param.Rs = Rs;
param.relaxPrecs = relaxPrecs;
param = adjustMemoryForNumRHS(param,nrhs,verbose);
return;
end

function getStrengthMatrix(AT::SparseMatrixCSC,strengthConnParam::Float64, n::Int64)
S = -AT;
mm = 1e-16*maximum(S.nzval);
for j = 1:n
	maxVal_j = mm;
	for gidx = S.colptr[j] : S.colptr[j+1]-1 ### here we go over the j-th row
		if S.nzval[gidx] > maxVal_j
			maxVal_j = S.nzval[gidx];
		end
	end
	scal_k = 1.0./maxVal_j;
	for gidx = S.colptr[j] : S.colptr[j+1]-1
		S.nzval[gidx]*=scal_k;
	end
	for gidx = S.colptr[j] : S.colptr[j+1]-1
		if S.nzval[gidx] <  strengthConnParam
			S.nzval[gidx] = 0.0;
		end
	end
	for gidx = S.colptr[j] : S.colptr[j+1]-1
		if S.rowval[gidx] == j
			# We have a diagonal element.
			S.nzval[gidx] = 1.0;
		end
	end
end
S = (S + S')/2;
dropzeros!(S)
end

function checkColoring(coloring::Array, n::Int64)
	coarse = 0
	fine = 0
	for i=1:n
		if coloring[i] == 1
			coarse+=1
		else
			fine+=1
		end
	end
	#println("fine : $fine , coarse : $coarse")
end


text:
#=
Filename: aiyagari_compute_equilibrium.jl
Author: Victoria Gregory
Date: 8/30/2016

Draws the aggregate supply and demand curves for
the Aiyagari model.
=#

using QuantEcon
include("aiyagari_household.jl")
using Plots
pyplot()

# Firms' parameters
A = 1
N = 1
alpha = 0.33
beta = 0.96
delta = 0.05

"""
Compute wage rate given an interest rate, r
"""
function r_to_w(r::Float64)
    return A * (1 - alpha) * (A * alpha / (r + delta)) ^ (alpha / (1 - alpha))
end

"""
Inverse demand curve for capital. The interest rate
associated with a given demand for capital K.
"""
function rd(K::Float64)
    return A * alpha * (N / K) ^ (1 - alpha) - delta
end

"""
Map prices to the induced level of capital stock.

##### Arguments
- `am::Household` : Household instance for problem we want to solve
- `r::Float64` : interest rate

##### Returns
- The implied level of aggregate capital
"""
function prices_to_capital_stock(am::Household, r::Float64)

    # Set up problem
    w = r_to_w(r)
    setup_R!(am, r, w)
    aiyagari_ddp = DiscreteDP(am.R, am.Q, am.beta)

    # Compute the optimal policy
    results = solve(aiyagari_ddp, PFI)

    # Compute the stationary distribution
    stationary_probs = stationary_distributions(results.mc)[:, 1][1]

    # Return K
    return sum(am.s_vals[:, 1] .* stationary_probs)
end

# Create an instance of Household
z_chain = MarkovChain([0.67 0.33; 0.33 0.67], [0.5, 1.5])
am = Household(z_chain=z_chain, beta=beta, a_max=20.0)

# Create a grid of r values at which to compute demand and supply of capital
num_points = 20
r_vals = linspace(0.02, 1/beta - 1, num_points)

# Compute supply of capital
k_vals = Array(Float64, num_points, 1)
for i in 1:num_points
    k_vals[i] = prices_to_capital_stock(am, r_vals[i])
end

# Plot against demand for capital by firms
demand = [rd(k) for k in k_vals]
labels =  ["demand for capital"; "supply of capital"]
plot(k_vals, [demand r_vals], label=labels', lw=2, alpha=0.6)
plot!(xlabel="capital", ylabel="interest rate")

max_stars_repo_name: mwhchen/quantecon

max_issues_repo_name: mwhchen/quantecon

max_forks_repo_name: mwhchen/quantecon


text:
using ApproxFunOrthogonalPolynomials, ApproxFunBase, LazyArrays, FillArrays, LinearAlgebra, SpecialFunctions, Test
    import ApproxFunBase: interlace, Multiplication, ConstantSpace, PointSpace,
                        ArraySpace, testblockbandedoperator, blocklengths, ∞,
                        testraggedbelowoperator, Vec

@testset "Vector" begin
    @testset "Construction" begin
        f = Fun(x->[1.,0.])
        @test f(0.) ≈ [1.,0.]
        @test f[end] ≈ 0
        @test f[firstindex(f)] ≈ 1
    end

    @testset "Broadcast" begin
        x = Fun()

        @test x .+ [1,2] ≈ [x+1,x+2]
        @test [1,2] .+ x ≈ [x+1,x+2]

        @test x.*[1,2] ≈ [x,2x]
        @test [1,2].*x ≈ [x,2x]
        @test (x*[1,2])(0.1) ≈ [0.1,0.2]
        @test ([1,2]*x)(0.1) ≈ [0.1,0.2]

        @test (x*[1,2])(0.1) ≈ [0.1,0.2]
        @test ([1,2]*x)(0.1) ≈ [0.1,0.2]

        @test (x*Fun([1,2]))(0.1) ≈ [0.1,0.2]
        @test (Fun([1,2])*x)(0.1) ≈ [0.1,0.2]

        @test (Fun(x->[1., 2.]) + [2, 2])(0.) ≈ [3., 4.]
    end

    @testset "Vector*Vector{Fun}" begin
        f = Fun(x->[exp(x),cos(x)])
        @test f(0.1) ≈ [exp(0.1),cos(0.1)]
        @test ([1 2]*f)(0.1) ≈ [1 2]*f(0.1)
        @test (f*3)(0.1) ≈ f(0.1)*3
        @test (3*f)(0.1) ≈ f(0.1)*3

        @test_broken transpose(f)*[1,2] ≈ transpose(f(0.1))*[1,2]

        @test norm(f) ≈ sqrt(sinh(2)+1+cos(1)sin(1))
        @test norm(f,2) ≈ sqrt(sinh(2)+1+cos(1)sin(1))
        @test norm(f,1) ≈ sqrt(sinh(2))+sqrt(1+cos(1)sin(1))
    end

    @testset "Chebyshev Vector" begin
        a = [1 2; 3 4]
        f = Fun(x->[exp(x),cos(x)])
        @test (a*f)(0.1) ≈ [exp(0.1)+2cos(0.1); 3exp(0.1)+4cos(0.1)]
        @test (a*f)(0.1) ≈ a*f(0.1)
        @test Fun(a)*f ≈ a*f
        @test Fun(a*Array(f)) ≈ a*f
    end

    @testset "Chebyshev Matrix" begin
        a = [1 2; 3 4]
        m = Fun(x->[exp(x) cos(x); sin(x) airyai(x)])
        @test (a*m)(0.1) ≈ [exp(0.1)+2sin(0.1) cos(0.1)+2airyai(0.1);
                            3exp(0.1)+4sin(0.1) 3cos(0.1)+4airyai(0.1)]
        @test (a*m)(0.1) ≈ a*m(0.1)
        @test (m*a)(0.1) ≈ m(0.1)*a
        @test Fun(a)*m   ≈ a*m
        @test m*Fun(a)   ≈ m*a
        @test Fun(a*Array(m))   ≈ a*m
        @test Fun(Array(m)*a)   ≈ m*a

        @test (a+m)(0.1) ≈ a+m(0.1)
        @test (m+a)(0.1) ≈ m(0.1)+a

        @test (m+I)(0.1) ≈ m(0.1)+I
    end

    @testset "Matrix{Fun}*Matrix{Fun}" begin
    # note that 2x2 and 3x3 mult are special cases
        x=Fun()
        A = [x x; x x]

        @test A isa Fun

        @test norm(map(norm,A*A-[2x^2 2x^2; 2x^2 2x^2])) <eps()
        @test norm(map(norm,A*Array(A)-[2x^2 2x^2; 2x^2 2x^2])) < eps()
        @test norm(map(norm,Array(A)*A-[2x^2 2x^2; 2x^2 2x^2])) < eps()
        @test norm(map(norm,Array(A)*Array(A)-Array([2x^2 2x^2; 2x^2 2x^2]))) < eps()

        @test [x;x] == A[:,1]
        @test [[x;2x] [3x;4x]] == [x 3x; 2x 4x]
        @test [x x] == A[1:1,:]

        A = fill(x,3,3)
        @test norm(map(norm,A^3-fill(9x^3,3,3))) <eps()
        @test norm((A^2*Fun(A)-fill(9x^3,3,3)).coefficients) < eps()
        @test norm((A*Fun(A)^2-fill(9x^3,3,3)).coefficients) < eps()
        @test norm((Fun(A)^3-fill(9x^3,3,3)).coefficients) < eps()

        A = fill(x,4,4)
        @test norm(map(norm,A^2-fill(4x^2,4,4))) <eps()
        @test norm((A*Fun(A)-fill(4x^2,4,4)).coefficients) < eps()
        @test norm((Fun(A)*A-fill(4x^2,4,4)).coefficients) < eps()
        @test norm((Fun(A)^2-fill(4x^2,4,4)).coefficients) < eps()
    end

    @testset "Vector ODE" begin
        d=ChebyshevInterval()
        D=Derivative(d);
        B=ldirichlet();
        Bn=lneumann();
        A=[B 0;
           0 B;
           D-I 2I;
           0I D+I];

        f=Fun(x->[exp(x),cos(x)],d)

        b=Any[0.,0.,f...]
        f1,f2=vec(f)
        @time u=A\b
        u1=vec(u)[1];u2=vec(u)[2];

        @test norm(u1'-u1+2u2-f1)<10eps()
        @test norm(u2'+u2-f2)<10eps()

        A=[B 0;
           Bn 0;
           0 B;
           D^2-I 2I;
           0 D+I];

        b=Any[0.,0.,0.,f...]

        @time u=A\b
        u1=vec(u)[1];u2=vec(u)[2];

        @test norm(differentiate(u1,2)-u1+2u2-f1) < 4eps()
        @test norm(u2'+u2-f2) < 4eps()
    end

    @testset "Matrix exponential" begin
        n=4
        d=fill(Interval(0.,1.),n)
        B=Evaluation(d,0.)
        D=Derivative(d)
        A=rand(n,n)
        L=[B;D-A]
        @time u=L\Matrix(I,2n,n)
        @test norm(u(1.)-exp(A))<eps(1000.)

        n=4
        d=fill(Interval(0.,1.),n)
        B=Evaluation(d,0.)
        D=Derivative(d)
        A=rand(n,n)
        L=[B;D-A]
        @time u=L\Matrix(I,2n,2)
        @test norm(u(1.)-exp(A)[:,1:2])<eps(1000.)
    end

    @testset "Multiplication" begin
        d = ChebyshevInterval()
        t=Fun(identity,d)
        f = Fun([t^2, sin(t)])
        @test norm(((Derivative(space(f))*f)-Fun(t->[2t,cos(t)])).coefficients)<100eps()
        @test norm((([1 2;3 4]*f)-Fun(t->[t^2+2sin(t),3t^2+4sin(t)])).coefficients)<100eps()
    end

    @testset "Interlace" begin
        S1 = Chebyshev()^2
        S2 = Chebyshev()
        TS = ArraySpace([ConstantSpace(),S1,ConstantSpace(),S2,PointSpace([1.,2.])])
        f = Fun(TS,collect(1:13))
        @test f[1] == Fun(TS[1],[1.])
        @test f[2] == Fun(TS[2],[2.,7.,8.,10.,11.,12.])
        @test f[3] == Fun(TS[3],[3.])
        @test f[4] == Fun(TS[4],[4.,9.,13.])
        @test f[5] == Fun(TS[5],[5.,6.])
    end

    @testset "Operator * Matrix" begin
        f = [Fun(exp) Fun(cos)]
        @test f(0.1) ≈ [exp(0.1) cos(0.1)]

        D=Derivative()
        u=D*f
        @test u(0.1) ≈ [exp(0.1) -sin(0.1)]
        u=D*Array(f)
        @test u(0.1) ≈ [exp(0.1) -sin(0.1)]
    end

    @testset "Floquet" begin
        # TODO: fix when IntervalSets.jl is fixed
        T = Float64(π, RoundUp); a=0.15
        t = Fun(identity,0..T)
        d=domain(t)
        D=Derivative(d)
        B=ldirichlet(d)

        B_row = [D             -I  0I            0I]

        f=Fun(exp,d)
        @test norm((B_row*[f;f;f;f])[1]) ≤ 1000eps()
        @test B_row isa ApproxFunBase.MatrixInterlaceOperator

        @test size([B_row;B_row].ops) == (2,4)

        @test norm(([B_row;B_row]*[f;f;f;f])[1]) ≤ 1000eps()

        n=4
        Dg = Operator(diagm(0 => fill(ldirichlet(d),n)))

        @test Dg isa ApproxFunBase.MatrixInterlaceOperator
        @test size([Dg; B_row].ops) == (5,4)

        L = [Dg; B_row]
        testraggedbelowoperator(L)
        @test ([Dg; B_row]*[f;f;f;f])(0.1) ≈ [fill(1.0,4);0]

        @test hcat(Dg).ops == Dg.ops

        B_row2 = [(2+a*cos(2t))   D  -I            0I]
        @test ([B_row;B_row2]*[f;f;f;f])(0.1) ≈ [0.,(2+a*cos(2*0.1))*f(0.1) + f'(0.1) - f(0.1)]

        A=[ Operator(diagm(0 => fill(ldirichlet(d),n)));
            D             -I  0I            0I;
           (2+a*cos(2t))   D  -I            0I;
           0I             0I   D            -I;
           -I             0I  (2+a*cos(2t))  D]

        Φ = A\Matrix(I,2n,n);

        @test Φ(π) ≈ [-0.170879 -0.148885 -0.836059 0.265569;
                 0.732284 -0.170879 -0.612945 -0.836059;
                 -0.836059 0.265569 -0.170879 -0.148885;
                 -0.612945 -0.836059 0.732284 -0.170879] atol=1E-3
    end

    @testset "diagonal" begin
        x = Fun(0.0..1)
        @test isambiguous(domain(zero(typeof(x))))
        @test diagm(0 => [x,x]) * [1,2] == [x,2x]
    end
end


text:
include("../server/simulation/simulation.jl")

using Printf

using .Simulation

const SAMPLE_DEST_PATH = "/media/RAMDisk/"

synapses = 10
time_delta = 10 # microseconds/(step|pass)
duration = 20 * time_delta # = 2000ms * 100us/step
firing_rate = 0.02
chunk_count = 3

function test_write_poisson_samples()
    println("Test: generate poisson samples")

    for chunk in 1:chunk_count
        file = SAMPLE_DEST_PATH * "poi_samples_c" * string(chunk) * ".data"

        samples = zeros(UInt8, synapses, duration)
        streams = Array{Simulation.AbstractBitStream,1}()

        # Create and collect a stream for each synapse
        for id in 1:synapses
            seed = Int64(round(rand(1)[1] * 10000.0))

            stream = Simulation.PoissonStream(UInt64(seed), firing_rate)

            push!(streams, stream)
        end

        # Exercise each stream and capture outputs
        for t in 1:duration
            for id in 1:synapses
                stream = streams[id]
                # Exercise stream
                Simulation.step!(stream)
                # Capture output
                samples[id, t] = Simulation.output(stream)
            end
        end

        # Now write all each stream/synpase-input
        # Format:
        # id 1010001011010...::
        open(file, "w") do f
            for id in 1:synapses
                # Write stream id
                print(f, @sprintf("%03d ", id))
                # write all spikes and non-spikes
                for t in 1:duration
                    print(f, samples[id, t])
                end
                # Terminate synapse stream with "::" marker
                println(f, "::")
            end
        end
    end

    println("Passed")
end

function test_read_poisson_samples()
    samples = zeros(UInt8, synapses, duration)
    file = SAMPLE_DEST_PATH * "poi_samples_c1.data"

    idx = 1
    open(file, "r") do f
        syn_samples = readlines(f)
        for syn_sample in syn_samples
            # Format:
            # id 1010001011010...::

            # Parse out "id" field
            idx_range = findfirst(" ", syn_sample)
            id = parse(Int64, SubString(syn_sample, 1, idx_range[1]))
            # println(id)

            # EOL marker
            bits_end = findfirst("::", syn_sample)

            # Extract just spike data
            bits = SubString(syn_sample, idx_range[1] + 1, bits_end[1] - 1)

            # load spikes into samples
            t = 1
            for bit in bits
                samples[idx, t] = parse(UInt8, bit)
                t += 1
            end

            idx += 1
        end

        # println("samples")
        # println(samples[8, 1:100])
    end
end

test_write_poisson_samples()
# test_read_poisson_samples()


text:
# implement julia_commander
# for the moment, only help is implemented

function filterranges(u::Vector{UnitRange{Int}})
    v = Vector{Int}()
    while length(u) > 1
        if all(m -> (u[1] ⊈ m), u[2:end])
            push!(v, u[1][end])
        end
        u = u[findall(m -> m ⊈ u[1], u[2:end]).+1]
    end
    if length(u) == 1
        push!(v,u[1][end])
    end
    return v
end

const STYLES = [
    r"```.+?```"s, r"`.+?`", r"~~.+?~~", r"_.+?_", r"__.+?__",
    r"\*.+?\*", r"\*\*.+?\*\*", r"\*\*\*.+?\*\*\*",
]

function split_message_fixed(text::AbstractString, chunk_limit::Int=2000; extraregex::Vector{Regex}=Vector{Regex}())
    length(text) <= chunk_limit && return String[text]
    chunks = String[]

    while !isempty(text)
        mranges = vcat(findall.(union(STYLES,extraregex),Ref(text))...)

        stop = maximum(filter(i -> length(text[1:i]) ≤ chunk_limit, filterranges(mranges)))

        # give up if first chunk cannot be broken down
        if stop == 0
            push!(chunks, strip(text))
            @warn "message could not be broken down into chunks smaller than the desired length $chunk_limit"
            return chunks
        end

        push!(chunks, strip(text[1:stop]))
        text = strip(text[stop+1:end])
    end

    return chunks
end

function parse_doc(doc::AbstractString)
    doc = replace(doc, r"\n\n\n+" => "\n\n")
    for m in eachmatch(r"(^|\n)(#+ |!!! )(.*)\n",doc)
        if m.captures[2] == "# "
            doc = replace(doc, m.match => m.captures[1]*"**"*m.captures[3]*"**\n"*"≡"^(length(m.captures[3])), count = 1)
        elseif m.captures[2] == "!!! "
            doc = replace(doc, m.match => m.captures[1]*"__"*m.captures[3]*"__\n", count = 1)
        else
            doc = replace(doc, m.match => m.captures[1]*"*"*m.captures[3]*"*\n"*"-"^(length(m.captures[3])), count = 1)
        end
    end
    doc = replace(doc, r"(```.+)\n" => "```julia\n")
    doc = replace(doc, "```\n" => "```")
    for m in eachmatch(r"\[([^ ]*)\]\(@ref\)",doc)
        doc = replace(doc, m.match => "`"*m.captures[1]*"`", count = 1)
    end
    return doc
end

function commander(c::Client, m::Message, ::Val{:julia_doc})
    # @info "julia_commander called"
    # @info "Message content" m.content m.author.username m.author.discriminator
    startswith(m.content, COMMAND_PREFIX * "j") || return
    regex = Regex(COMMAND_PREFIX * raw"j(\?| help| doc| packages| stats)? *(.*)$")
    matches = match(regex, m.content)
    if matches === nothing || matches.captures[1] in (" help", nothing)
        help_commander(c, m, Val(:julia_doc))
    elseif matches.captures[1] in ("?", " doc")
        handle_julia_help_commander(c, m, matches.captures[2])
    elseif matches.captures[1] == " packages"
        handle_julia_package_list(c, m)
    elseif matches.captures[1] == " stats"
        handle_doc_stats(c, m, matches.captures[2])
    else
        reply(c, m, "Sorry, I don't understand the request; use `j help` for help")
    end
    return nothing
end

function help_commander(c::Client, m::Message, ::Val{:julia_doc})
    # @info "Sending help for message" m.id m.author
    reply(c, m, """
        The `j` commands shows the docstring of names in the Base and other selected packages,
        as well as statistics of its use.

        How to use the `j` command:
        ```
        j help
        j? <name>
        j doc <name>
        j packages
        j stats <name>
        j stats <place> <number>
        ```
        `j help` returns this help
        `j packages` shows which packages are available for showing their docstrings (incomplete)
        `j? <name>` and `j doc <name>` return the documentation for `<name>`
        `j stats <name>` return how many times the docstring for `name` has been queried.
        `j stats <place> <number>` return the top (if `place` is equal to either "top" or "head") or the bottom (if `place` is equal to either "bottom" or "tail") `number` names that have been queried.
        """)
    return nothing
end

function handle_julia_help_commander(c::Client, m::Message, name)
    # @info "julia_help_commander called"
    try
        doc = string(eval(:(Base.Docs.@doc $(Symbol(name)))))
        doc = parse_doc(doc)
        user = @discord retrieve(c, User, m.author.id)
        channel = @discord get_channel(c, m.channel_id)
        if !occursin("No documentation found", doc)
            count = update_names_count(
                "namescount", name,
                m.channel_id, channel.name)
            doc *= "\n*(Count number for `$name`: $count)*"
        end
        docs = split_message_fixed(doc, extraregex = [r"\n≡.+\n", r"n-.+\n"])
        for doc_chunck in docs
            reply(c, m, doc_chunck)
        end
    catch ex
        @show ex
        reply(c, m, "Sorry, it didn't work.")
    end
    return nothing
end

function handle_doc_stats(c::Client, m::Message, captured::AbstractString)
    # @info "handle_doc_stats called"
    try
        r = match(r"(top|head|high|bottom|tail|low) +(\d*) *$", captured)
        if captured == ""
            statsmgs = stats_namescount("namescount")
        elseif r === nothing
            statsmgs = stats_namescount("namescount", name=strip(captured))
        else
            place = r.captures[1] in ("top", "head", "high") ? "top" : "bottom"
            statsmgs = stats_namescount("namescount", place=place, number=parse(Int,r.captures[2]))
        end
        reply(c, m, statsmgs)
    catch ex
        @show ex
        reply(c, m, "Sorry, it didn't work.")
    end
    return nothing
end

function handle_julia_package_list(c::Client, m::Message)
    reply(c, m, "Packages available with docstrings:\n≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡\n `Base`")
end


text:
"""
    scaleInflate!(L::AbstractMatrix, A::AbstractMatrix, Λ::AbstractTerm)

Overwrite a diagonal block of `L` with the corresponding block of `Λ'AΛ + I` except when Λ
is a [`MatrixTerm`]{@ref}, in which case this becomes `copyto!(L, A)`.
"""
function scaleInflate! end

function scaleInflate!(Ljj::Matrix{T}, Ajj::Matrix{T}, Λj::MatrixTerm{T}) where {T}
    @argcheck(size(Ljj) == size(Ajj), DimensionMismatch)
    copyto!(Ljj, Ajj)
end

function scaleInflate!(Ljj::Diagonal{T}, Ajj::Diagonal{T},
        Λj::ScalarFactorReTerm{T}) where {T}
    broadcast!((x,k) -> k * x + one(T), Ljj.diag, Ajj.diag, abs2(Λj.Λ))
    Ljj
end

function scaleInflate!(Ljj::Matrix{T}, Ajj::Diagonal{T},
        Λj::ScalarFactorReTerm{T}) where {T}
    n = LinearAlgebra.checksquare(Ljj)
    Ad = Ajj.diag
    @argcheck(length(Ad) == n, DimensionMismatch)
    lambsq = abs2(Λj.Λ)
    fill!(Ljj, zero(T))
    for (j, dj) in enumerate(Ad)
        Ljj[j,j] = lambsq * dj + one(T)
    end
    Ljj
end

function scaleInflate!(Ljj::UniformBlockDiagonal{T}, Ajj::UniformBlockDiagonal{T},
        Λj::VectorFactorReTerm{T}) where {T}
    Ljjdd = Ljj.data
    Ajjdd = Ajj.data
    if ((m, n, l) = size(Ljjdd)) != size(Ajjdd)
        throw(DimensionMismatch("size(Ljj.data) = $(size(Ljjdd)) != $(size(Ajjdd)) = size(Ajj.data)"))
    end
    copyto!(Ljjdd, Ajjdd)
    m, n, l = size(Ljjdd)
    λ = Λj.Λ
    Lfv = Ljj.facevec
    dind = diagind(Lfv[1])
    @inbounds for Lf in Lfv
        lmul!(adjoint(λ), rmul!(Lf, λ))
    end
    @inbounds for k in 1:l, i in 1:m
        Ljjdd[i, i, k] += one(T)
    end
    Ljj
end

function scaleInflate!(Ljj::Matrix{T}, Ajj::UniformBlockDiagonal{T},
        Λj::VectorFactorReTerm{T}) where {T}
    if size(Ljj) != size(Ajj)
        throw(DimensionMismatch("size(Ljj) = $(size(Ljj)) != $(size(Ajj)) = size(Ajj)"))
    end
    λ = Λj.Λ
    Afv = Ajj.facevec
    m, n, l = size(Ajj.data)
    m == n || throw(ArgumentError("Diagonal blocks of Ajj must be square"))
    fill!(Ljj, zero(T))
    tmp = Array{T}(undef, m, m)
    @inbounds for (k, Af) in enumerate(Afv)
        lmul!(adjoint(λ), rmul!(copyto!(tmp, Af), λ))
        offset = (k - 1)*m
        for j in 1:m, i in 1:m
            Ljj[offset + i, offset + j] = tmp[i, j] + (i == j)
        end
    end
    Ljj
end


text:
module FixedEffectModels

##############################################################################
##
## Dependencies
##
##############################################################################
import Base: size, copyto!, getindex, length, fill!, eltype, length, view, adjoint
import LinearAlgebra: mul!, rmul!, norm, Matrix, Diagonal, cholesky, cholesky!, Symmetric, Hermitian, rank, dot, eigen, axpy!, svd, I, Adjoint, diag, qr
import LinearAlgebra.BLAS: gemm!
import Statistics: mean, quantile
import Printf: @sprintf
import Distributions: ccdf, TDist, FDist, Chisq
import DataFrames: DataFrame, AbstractDataFrame, completecases, names!, ismissing, groupby, groupindices
import Combinatorics: combinations
using CategoricalArrays
using FillArrays
import StatsBase: coef, nobs, coeftable, vcov, predict, residuals, var, RegressionModel, model_response, stderror, confint, fit, CoefTable, dof_residual, r2, adjr2, deviance, mss, rss, islinear, response
import StatsModels: @formula,  Formula, ModelFrame, ModelMatrix, Terms, coefnames, evalcontrasts, check_non_redundancy!
using Reexport
@reexport using StatsBase
@reexport using StatsModels
@reexport using FixedEffects

##############################################################################
##
## Exported methods and types
##
##############################################################################

export reg,
partial_out,
allvars,
fes,
WeightFormula,

AbstractRegressionResult,
title,
top,
RegressionResult,
RegressionResultIV,
RegressionResultFE,
RegressionResultFEIV,

AbstractVcovFormula,
VcovSimpleFormula,
VcovRobustFormula,
VcovClusterFormula,
VcovFormula,

AbstractVcovMethod,
VcovMethod,
VcovSimpleMethod,
VcovWhiteMethod,
VcovClusterMethod,

vcov!,
shat!,
VcovData,

Model,
@model

##############################################################################
##
## Load files
##
##############################################################################
include("utils/model.jl")
include("utils/weights.jl")
include("utils/fixedeffects.jl")
include("utils/basecol.jl")
include("utils/tss.jl")

include("formula/formula_iv.jl")
include("formula/formula_fe.jl")

include("RegressionResult.jl")

include("vcov/types.jl")
include("vcov/vcovsimple.jl")
include("vcov/vcovrobust.jl")
include("vcov/vcovcluster.jl")
include("vcov/utils.jl")

include("reg.jl")
include("partial_out.jl")

end  # module FixedEffectModels