text:
"""
    AbstractVectorTransportMethod

Abstract type for methods for transporting vectors. Such vector transports are not
necessarily linear.

# See also

[`AbstractLinearVectorTransportMethod`](@ref)
"""
abstract type AbstractVectorTransportMethod end

"""
    AbstractLinearVectorTransportMethod <: AbstractVectorTransportMethod

Abstract type for linear methods for transporting vectors, that is transport of a linear
combination of vectors is a linear combination of transported vectors.
"""
abstract type AbstractLinearVectorTransportMethod <: AbstractVectorTransportMethod end

@doc raw"""
    DifferentiatedRetractionVectorTransport{R<:AbstractRetractionMethod} <:
        AbstractVectorTransportMethod

A type to specify a vector transport that is given by differentiating a retraction.
This can be introduced in two ways. Let ``\mathcal M`` be a Riemannian manifold,
``p∈\mathcal M`` a point, and ``X,Y∈ T_p\mathcal M`` denote two tangent vectors at ``p``.

Given a retraction (cf. [`AbstractRetractionMethod`](@ref)) ``\operatorname{retr}``,
the vector transport of `X` in direction `Y` (cf. [`vector_transport_direction`](@ref))
by differentiation this retraction, is given by

```math
\mathcal T^{\operatorname{retr}}_{p,Y}X
= D_Y\operatorname{retr}_p(Y)[X]
= \frac{\mathrm{d}}{\mathrm{d}t}\operatorname{retr}_p(Y+tX)\Bigr|_{t=0}.
```
see [^AbsilMahonySepulchre2008], Section 8.1.2 for more details.

This can be phrased similarly as a [`vector_transport_to`](@ref) by introducing
``q=\operatorname{retr}_pX`` and defining

```math
\mathcal T^{\operatorname{retr}}_{q \gets p}X = \mathcal T^{\operatorname{retr}}_{p,Y}X
```

which in practice usually requires the [`inverse_retract`](@ref) to exists in order to
compute ``Y = \operatorname{retr}_p^{-1}q``.

# Constructor

    DifferentiatedRetractionVectorTransport(m::AbstractRetractionMethod)
"""
struct DifferentiatedRetractionVectorTransport{R<:AbstractRetractionMethod} <:
       AbstractLinearVectorTransportMethod
    retraction::R
end

@doc raw"""
    ParallelTransport <: AbstractVectorTransportMethod

Compute the vector transport by parallel transport, see
[`parallel_transport_to`](@ref)
"""
struct ParallelTransport <: AbstractLinearVectorTransportMethod end

"""
    ProjectionTransport <: AbstractVectorTransportMethod

Specify to use projection onto tangent space as vector transport method within
[`vector_transport_to`](@ref), [`vector_transport_direction`](@ref), or
[`vector_transport_along`](@ref). See [`project`](@ref) for details.
"""
struct ProjectionTransport <: AbstractLinearVectorTransportMethod end

@doc raw"""
    PoleLadderTransport <: AbstractVectorTransportMethod

Specify to use [`pole_ladder`](@ref) as vector transport method within
[`vector_transport_to`](@ref), [`vector_transport_direction`](@ref), or
[`vector_transport_along`](@ref), i.e.

Let $X∈ T_p\mathcal M$ be a tangent vector at $p∈\mathcal M$ and $q∈\mathcal M$ the
point to transport to. Then $x = \exp_pX$ is used to call
`y = `[`pole_ladder`](@ref)`(M, p, x, q)` and the resulting vector is obtained by computing
$Y = -\log_qy$.

The [`PoleLadderTransport`](@ref) posesses two advantages compared to
[`SchildsLadderTransport`](@ref):
* it is cheaper to evaluate, if you want to transport several vectors, since the
  mid point $c$ then stays unchanged.
* while both methods are exact if the curvature is zero, pole ladder is even exact in
  symmetric Riemannian manifolds[^Pennec2018]

The pole ladder was was proposed in [^LorenziPennec2014]. Its name stems from the fact that
it resembles a pole ladder when applied to a sequence of points usccessively.

# Constructor
````julia
PoleLadderTransport(
    retraction = ExponentialRetraction(),
    inverse_retraction = LogarithmicInverseRetraction(),
)
````
Construct the classical pole ladder that employs exp and log, i.e. as proposed
in[^LorenziPennec2014]. For an even cheaper transport the inner operations can be
changed to an [`AbstractRetractionMethod`](@ref) `retraction` and an
[`AbstractInverseRetractionMethod`](@ref) `inverse_retraction`, respectively.

[^LorenziPennec2014]:
    > Lorenzi, M. and Pennec, X: Efficient parallel transport of deformations in time
    > series of images: From Schild’s to pole ladder.
    > Journal of Mathematical Imaging and Vision (2014), 50(1), pp. 5–17
    > doi [10.1007/s10851-013-0470-3](https://doi.org/10.1007/s10851-013-0470-3),
    > hal: [hal-00870489](https://hal.inria.fr/hal-00870489)
[^Pennec2018]:
    > Pennec, X: Parallel Transport with Pole Ladder: a Third Order Scheme in Affine
    > Connection Spaces which is Exact in Affine Symmetric Spaces.
    > arXiv: [1805.11436](https://arxiv.org/abs/1805.11436)

"""
struct PoleLadderTransport{
    RT<:AbstractRetractionMethod,
    IRT<:AbstractInverseRetractionMethod,
} <: AbstractLinearVectorTransportMethod
    retraction::RT
    inverse_retraction::IRT
    function PoleLadderTransport(
        retraction = ExponentialRetraction(),
        inverse_retraction = LogarithmicInverseRetraction(),
    )
        return new{typeof(retraction),typeof(inverse_retraction)}(
            retraction,
            inverse_retraction,
        )
    end
end

@doc raw"""
    ScaledVectorTransport{T} <: AbstractVectorTransportMethod

Introduce a scaled variant of any [`AbstractVectorTransportMethod`](@ref) `T`,
as introduced in [^SatoIwai2013] for some ``X∈ T_p\mathcal M`` as

```math
    \mathcal T^{\mathrm{S}}(X) = \frac{\lVert X\rVert_p}{\lVert \mathcal T(X)\rVert_q}\mathcal T(X).
```

Note that the resulting point `q` has to be known, i.e. for [`vector_transport_direction`](@ref)
the curve or more precisely its end point has to be known (via an exponential map or a
retraction). Therefore a default implementation is only provided for the [`vector_transport_to`](@ref)

# Constructor

    ScaledVectorTransport(m::AbstractVectorTransportMethod)

[^SatoIwai2013]:
    > Sato, H., Iwai, T.: _A new, globally convergent Riemannian conjugate gradient method_,
    > Optimization, 2013, Volume 64(4), pp. 1011–1031.
    > doi: [10.1080/02331934.2013.836650](https://doi.org/10.1080/02331934.2013.836650),
    > arXiv: [1302.0125](https://arxiv.org/abs/1302.0125).
"""
struct ScaledVectorTransport{T<:AbstractVectorTransportMethod} <:
       AbstractVectorTransportMethod
    method::T
end

@doc raw"""
    SchildsLadderTransport <: AbstractVectorTransportMethod

Specify to use [`schilds_ladder`](@ref) as vector transport method within
[`vector_transport_to`](@ref), [`vector_transport_direction`](@ref), or
[`vector_transport_along`](@ref), i.e.

Let $X∈ T_p\mathcal M$ be a tangent vector at $p∈\mathcal M$ and $q∈\mathcal M$ the
point to transport to. Then

````math
P^{\mathrm{S}}_{q\gets p}(X) =
    \log_q\bigl( \operatorname{retr}_p ( 2\operatorname{retr}_p^{-1}c ) \bigr),
````

where $c$ is the mid point between $q$ and $d=\exp_pX$.

This method employs the internal function [`schilds_ladder`](@ref)`(M, p, d, q)` that avoids
leaving the manifold.

The name stems from the image of this paralleltogram in a repeated application yielding the
image of a ladder. The approximation was proposed in [^EhlersPiraniSchild1972].

# Constructor
````julia
SchildsLadderTransport(
    retraction = ExponentialRetraction(),
    inverse_retraction = LogarithmicInverseRetraction(),
)
````
Construct the classical Schilds ladder that employs exp and log, i.e. as proposed
in[^EhlersPiraniSchild1972]. For an even cheaper transport these inner operations can be
changed to an [`AbstractRetractionMethod`](@ref) `retraction` and an
[`AbstractInverseRetractionMethod`](@ref) `inverse_retraction`, respectively.

[^EhlersPiraniSchild1972]:
    > Ehlers, J., Pirani, F.A.E., Schild, A.: The geometry of free fall and light
    > propagation. In: O’Raifeartaigh, L. (ed.) General Relativity: Papers in Honour of
    > J. L. Synge, pp. 63–84. Clarendon Press, Oxford (1972).
    > reprint doi: [10.1007/s10714-012-1353-4](https://doi.org/10.1007/s10714-012-1353-4)

"""
struct SchildsLadderTransport{
    RT<:AbstractRetractionMethod,
    IRT<:AbstractInverseRetractionMethod,
} <: AbstractLinearVectorTransportMethod
    retraction::RT
    inverse_retraction::IRT
    function SchildsLadderTransport(
        retraction = ExponentialRetraction(),
        inverse_retraction = LogarithmicInverseRetraction(),
    )
        return new{typeof(retraction),typeof(inverse_retraction)}(
            retraction,
            inverse_retraction,
        )
    end
end

@doc raw"""
    VectorTransportDirection{VM<:AbstractVectorTransportMethod,RM<:AbstractRetractionMethod}
        <: AbstractVectorTransportMethod

Specify a [`vector_transport_direction`](@ref) using a [`AbstractVectorTransportMethod`](@ref)
with explicitly using the [`AbstractRetractionMethod`](@ref) to determine the point in
the specified direction where to transsport to.
Note that you only need this for the non-default (non-implicit) second retraction method
associated to a vector transport, i.e. when a first implementation assumed
an implicit associated retraction.
"""
struct VectorTransportDirection{
    VM<:AbstractVectorTransportMethod,
    RM<:AbstractRetractionMethod,
} <: AbstractVectorTransportMethod
    retraction::RM
    vector_transport::VM
    function VectorTransportDirection(
        vector_transport = ParallelTransport(),
        retraction = ExponentialRetraction(),
    )
        return new{typeof(vector_transport),typeof(retraction)}(
            retraction,
            vector_transport,
        )
    end
end

@doc raw"""
    VectorTransportTo{VM<:AbstractVectorTransportMethod,RM<:AbstractRetractionMethod}
        <: AbstractVectorTransportMethod

Specify a [`vector_transport_to`](@ref) using a [`AbstractVectorTransportMethod`](@ref)
with explicitly using the [`AbstractInverseRetractionMethod`](@ref) to determine the direction
that transports from  in `p`to `q`.
Note that you only need this for the non-default (non-implicit) second retraction method
associated to a vector transport, i.e. when a first implementation assumed
an implicit associated retraction.
"""
struct VectorTransportTo{
    VM<:AbstractVectorTransportMethod,
    IM<:AbstractInverseRetractionMethod,
} <: AbstractVectorTransportMethod
    inverse_retraction::IM
    vector_transport::VM
    function VectorTransportTo(
        vector_transport = ParallelTransport(),
        inverse_retraction = LogarithmicInverseRetraction(),
    )
        return new{typeof(vector_transport),typeof(inverse_retraction)}(
            inverse_retraction,
            vector_transport,
        )
    end
end

"""
    default_vector_transport_method(M::AbstractManifold)

The [`AbstractVectorTransportMethod`](@ref) that is used when calling
[`vector_transport_along`](@ref), [`vector_transport_to`](@ref), or
[`vector_transport_direction`](@ref) without specifying the vector transport method.
By default, this is [`ParallelTransport`](@ref).
"""
function default_vector_transport_method(::AbstractManifold)
    return ParallelTransport()
end

@doc raw"""
    pole_ladder(
        M,
        p,
        d,
        q,
        c = mid_point(M, p, q);
        retraction=default_retraction_method(M),
        inverse_retraction=default_inverse_retraction_method(M)
    )

Compute an inner step of the pole ladder, that can be used as a [`vector_transport_to`](@ref).
Let $c = \gamma_{p,q}(\frac{1}{2})$ mid point between `p` and `q`, then the pole ladder is
given by

````math
    \operatorname{Pl}(p,d,q) = \operatorname{retr}_d (2\operatorname{retr}_d^{-1}c)
````

Where the classical pole ladder employs $\operatorname{retr}_d=\exp_d$
and $\operatorname{retr}_d^{-1}=\log_d$ but for an even cheaper transport these can be set
to different [`AbstractRetractionMethod`](@ref) and [`AbstractInverseRetractionMethod`](@ref).

When you have $X=log_pd$ and $Y = -\log_q \operatorname{Pl}(p,d,q)$,
you will obtain the [`PoleLadderTransport`](@ref). When performing multiple steps, this
method avoids the switching to the tangent space. Keep in mind that after $n$ successive
steps the tangent vector reads $Y_n = (-1)^n\log_q \operatorname{Pl}(p_{n-1},d_{n-1},p_n)$.

It is cheaper to evaluate than [`schilds_ladder`](@ref), sinc if you want to form multiple
ladder steps between `p` and `q`, but with different `d`, there is just one evaluation of a geodesic
each., since the center `c` can be reused.
"""
function pole_ladder(
    M,
    p,
    d,
    q,
    c = mid_point(M, p, q);
    retraction = default_retraction_method(M),
    inverse_retraction = default_inverse_retraction_method(M),
)
    return retract(M, d, 2 * inverse_retract(M, d, c, inverse_retraction), retraction)
end
@doc raw"""
    pole_ladder(
        M,
        pl,
        p,
        d,
        q,
        c = mid_point(M, p, q),
        X = allocate_result_type(M, log, d, c);
        retraction = default_retraction_method(M),
        inverse_retraction = default_inverse_retraction_method(M),
    )

Compute the [`pole_ladder`](@ref), i.e. the result is saved in `pl`.
`X` is used for storing intermediate inverse retraction.
"""
function pole_ladder!(
    M,
    pl,
    p,
    d,
    q,
    c = mid_point(M, p, q),
    X = allocate_result(M, log, d, c);
    retraction = default_retraction_method(M),
    inverse_retraction = default_inverse_retraction_method(M),
)
    inverse_retract!(M, X, d, c, inverse_retraction)
    X *= 2
    return retract!(M, pl, d, X, retraction)
end

@doc raw"""
    schilds_ladder(
        M,
        p,
        d,
        q,
        c = mid_point(M, q, d);
        retraction = default_retraction_method(M),
        inverse_retraction = default_inverse_retraction_method(M),
    )

Perform an inner step of schilds ladder, which can be used as a
[`vector_transport_to`](@ref), see [`SchildsLadderTransport`](@ref).
Let $c = \gamma_{q,d}(\frac{1}{2})$ denote the mid point
on the shortest geodesic connecting $q$ and the point $d$. Then Schild's ladder reads as

````math
\operatorname{Sl}(p,d,q) = \operatorname{retr}_x( 2\operatorname{retr}_p^{-1} c)
````

Where the classical Schilds ladder employs $\operatorname{retr}_d=\exp_d$
and $\operatorname{retr}_d^{-1}=\log_d$ but for an even cheaper transport these can be set
to different [`AbstractRetractionMethod`](@ref) and [`AbstractInverseRetractionMethod`](@ref).

In consistency with [`pole_ladder`](@ref) you can change the way the mid point is computed
using the optional parameter `c`, but note that here it's the mid point between `q` and `d`.

When you have $X=log_pd$ and $Y = \log_q \operatorname{Sl}(p,d,q)$,
you will obtain the [`PoleLadderTransport`](@ref).
Then the approximation to the transported vector is given by $\log_q\operatorname{Sl}(p,d,q)$.

When performing multiple steps, this method avoidsd the switching to the tangent space.
Hence after $n$ successive steps the tangent vector reads
$Y_n = \log_q \operatorname{Pl}(p_{n-1},d_{n-1},p_n)$.
"""
function schilds_ladder(
    M,
    p,
    d,
    q,
    c = mid_point(M, q, d);
    retraction = default_retraction_method(M),
    inverse_retraction = default_inverse_retraction_method(M),
)
    return retract(M, p, 2 * inverse_retract(M, p, c, inverse_retraction), retraction)
end
@doc raw"""
    schilds_ladder!(
        M,
        sl
        p,
        d,
        q,
        c = mid_point(M, q, d),
        X = allocate_result_type(M, log, d, c);
        retraction = default_retraction_method(M),
        inverse_retraction = default_inverse_retraction_method(M),
    )

Compute [`schilds_ladder`](@ref) and return the value in the parameter `sl`.
If the required mid point `c` was computed before, it can be passed using `c`,
and the allocation of new memory can be avoided providing a tangent vector `X`
for the interims result.
"""
function schilds_ladder!(
    M,
    sl,
    p,
    d,
    q,
    c = mid_point(M, q, d),
    X = allocate_result(M, log, d, c);
    retraction = default_retraction_method(M),
    inverse_retraction = default_inverse_retraction_method(M),
)
    inverse_retract!(M, X, p, c, inverse_retraction)
    X *= 2
    return retract!(M, sl, p, X, retraction)
end

function show(io::IO, ::ParallelTransport)
    return print(io, "ParallelTransport()")
end
function show(io::IO, m::ScaledVectorTransport)
    return print(io, "ScaledVectorTransport($(m.method))")
end

"""
    vector_transport_along(M::AbstractManifold, p, X, c)
    vector_transport_along(M::AbstractManifold, p, X, c, m::AbstractVectorTransportMethod)

Transport a vector `X` from the tangent space at a point `p` on the [`AbstractManifold`](@ref) `M`
along the curve represented by `c` using the `method`, which defaults to
[`default_vector_transport_method`](@ref)`(M)`.
"""
function vector_transport_along(
    M::AbstractManifold,
    p,
    X,
    c,
    m::AbstractVectorTransportMethod = default_vector_transport_method(M),
)
    return _vector_transport_along(M, p, X, c, m)
end
function _vector_transport_along(M::AbstractManifold, p, X, c, ::ParallelTransport)
    return parallel_transport_along(M, p, X, c)
end
function _vector_transport_along(
    M::AbstractManifold,
    p,
    X,
    c,
    m::DifferentiatedRetractionVectorTransport,
)
    return vector_transport_along_diff(M, p, X, c, m.retraction)
end
@doc raw"""
    vector_transport_along_diff(M::AbstractManifold, p, X, c, m::AbstractRetractionMethod)

Compute the vector transport of `X` from ``T_p\mathcal M`` along the curve `c`
using the differential of the [`AbstractRetractionMethod`](@ref) `m`.
"""
function vector_transport_along_diff(
    M::AbstractManifold,
    p,
    X,
    c,
    m::AbstractRetractionMethod,
)
    Y = allocate_result(M, vector_transport_along, X, p)
    return vector_transport_along_diff!(M, Y, p, X, c, m)
end
@doc raw"""
    vector_transport_along_diff!(M::AbstractManifold, Y, p, X, c, m::AbstractRetractionMethod)

Compute the vector transport of `X` from ``T_p\mathcal M`` along the curve `c`
using the differential of the [`AbstractRetractionMethod`](@ref) `m` in place of `Y`.
"""
vector_transport_along_diff!(M::AbstractManifold, Y, p, X, c, m)

function vector_transport_along_diff! end

function _vector_transport_along(M::AbstractManifold, p, X, c, ::ProjectionTransport)
    return vector_transport_along_project(M, p, X, c)
end
@doc raw"""
    vector_transport_along_project(M::AbstractManifold, p, X, c::AbstractVector)

Compute the vector transport of `X` from ``T_p\mathcal M`` along the curve `c`
using a projection.
"""
function vector_transport_along_project(M::AbstractManifold, p, X, c::AbstractVector)
    Y = allocate_result(M, vector_transport_along, X, p)
    return vector_transport_along_project!(M, Y, p, X, c)
end
@doc raw"""
    vector_transport_along_project!(M::AbstractManifold, Y, p, X, c::AbstractVector)

Compute the vector transport of `X` from ``T_p\mathcal M`` along the curve `c`
using a projection. The result is computed in place of `Y`.
"""
vector_transport_along_project!(M::AbstractManifold, Y, p, X, c::AbstractVector)

function vector_transport_along_project! end

function _vector_transport_along(
    M::AbstractManifold,
    p,
    X,
    c::AbstractVector,
    m::PoleLadderTransport,
)
    Y = allocate_result(M, vector_transport_along, X, p)
    return _vector_transport_along!(M, Y, p, X, c, m)
end
function _vector_transport_along(
    M::AbstractManifold,
    p,
    X,
    c::AbstractVector,
    m::SchildsLadderTransport,
)
    Y = allocate_result(M, vector_transport_along, X, p)
    return _vector_transport_along!(M, Y, p, X, c, m)
end

"""
    vector_transport_along!(M::AbstractManifold, Y, p, X, c::AbstractVector)
    vector_transport_along!(M::AbstractManifold, Y, p, X, c::AbstractVector, m::AbstractVectorTransportMethod)

Transport a vector `X` from the tangent space at a point `p` on the [`AbstractManifold`](@ref) `M`
along the curve represented by `c` using the `method`, which defaults to
[`default_vector_transport_method`](@ref)`(M)`. The result is saved to `Y`.
"""
function vector_transport_along!(
    M::AbstractManifold,
    Y,
    p,
    X,
    c::AbstractVector,
    m::AbstractVectorTransportMethod = default_vector_transport_method(M),
)
    return _vector_transport_along!(M, Y, p, X, c, m)
end

function parallel_transport_along!(M::AbstractManifold, Y, p, X, c::AbstractVector)
    n = length(c)
    if n == 0
        copyto!(Y, X)
    else
        parallel_transport_to!(M, Y, p, X, c[1])
        for i in 1:(length(c) - 1)
            parallel_transport_to!(M, Y, c[i], Y, c[i + 1])
        end
    end
    return Y
end

function _vector_transport_along!(
    M::AbstractManifold,
    Y,
    p,
    X,
    c::AbstractVector,
    ::ParallelTransport,
)
    return parallel_transport_along!(M, Y, p, X, c)
end
function _vector_transport_along!(
    M::AbstractManifold,
    Y,
    p,
    X,
    c::AbstractVector,
    m::DifferentiatedRetractionVectorTransport,
)
    return vector_transport_along_diff!(M, Y, p, X, c, m.retraction)
end
function _vector_transport_along!(
    M::AbstractManifold,
    Y,
    p,
    X,
    c::AbstractVector,
    ::ProjectionTransport,
)
    return vector_transport_along_project!(M, Y, p, X, c)
end

@doc raw"""
    function vector_transport_along(
        M::AbstractManifold,
        p,
        X,
        c::AbstractVector,
        m::PoleLadderTransport
    )

Compute the vector transport along a discretized curve using
[`PoleLadderTransport`](@ref) succesively along the sampled curve.
This method is avoiding additional allocations as well as inner exp/log by performing all
ladder steps on the manifold and only computing one tangent vector in the end.
"""
vector_transport_along(
    M::AbstractManifold,
    Y,
    p,
    X,
    c::AbstractVector,
    m::PoleLadderTransport,
)

function _vector_transport_along!(
    M::AbstractManifold,
    Y,
    p,
    X,
    c::AbstractVector,
    m::PoleLadderTransport,
)
    clen = length(c)
    if clen == 0
        copyto!(Y, X)
    else
        d = retract(M, p, X, m.retraction)
        mp = mid_point(M, p, c[1])
        pole_ladder!(
            M,
            d,
            p,
            d,
            c[1],
            mp,
            Y;
            retraction = m.retraction,
            inverse_retraction = m.inverse_retraction,
        )
        for i in 1:(clen - 1)
            # precompute mid point inplace
            ci = c[i]
            cip1 = c[i + 1]
            mid_point!(M, mp, ci, cip1)
            # compute new ladder point
            pole_ladder!(
                M,
                d,
                ci,
                d,
                cip1,
                mp,
                Y;
                retraction = m.retraction,
                inverse_retraction = m.inverse_retraction,
            )
        end
        inverse_retract!(M, Y, c[clen], d, m.inverse_retraction)
        Y *= (-1)^clen
    end
    return Y
end

@doc raw"""
    vector_transport_along(
        M::AbstractManifold,
        p,
        X,
        c::AbstractVector,
        m::SchildsLadderTransport
    )

Compute the vector transport along a discretized curve using
[`SchildsLadderTransport`](@ref) succesively along the sampled curve.
This method is avoiding additional allocations as well as inner exp/log by performing all
ladder steps on the manifold and only computing one tangent vector in the end.
"""
vector_transport_along(
    M::AbstractManifold,
    p,
    X,
    c::AbstractVector,
    m::SchildsLadderTransport,
)
function _vector_transport_along!(
    M::AbstractManifold,
    Y,
    p,
    X,
    c::AbstractVector,
    m::SchildsLadderTransport,
)
    clen = length(c)
    if clen == 0
        copyto!(Y, X)
    else
        d = retract(M, p, X, m.retraction)
        mp = mid_point(M, c[1], d)
        schilds_ladder!(
            M,
            d,
            p,
            d,
            c[1],
            mp,
            Y;
            retraction = m.retraction,
            inverse_retraction = m.inverse_retraction,
        )
        for i in 1:(clen - 1)
            ci = c[i]
            cip1 = c[i + 1]
            # precompute mid point inplace
            mid_point!(M, mp, cip1, d)
            # compute new ladder point
            schilds_ladder!(
                M,
                d,
                ci,
                d,
                cip1,
                mp,
                Y;
                retraction = m.retraction,
                inverse_retraction = m.inverse_retraction,
            )
        end
        inverse_retract!(M, Y, c[clen], d, m.inverse_retraction)
    end
    return Y
end

@doc raw"""
    vector_transport_direction(M::AbstractManifold, p, X, d)
    vector_transport_direction(M::AbstractManifold, p, X, d, m::AbstractVectorTransportMethod)

Given an [`AbstractManifold`](@ref) ``\mathcal M`` the vector transport is a generalization of the
[`parallel_transport_direction`](@ref) that identifies vectors from different tangent spaces.

More precisely using [^AbsilMahonySepulchre2008], Def. 8.1.1, a vector transport
``T_{p,d}: T_p\mathcal M \to T_q\mathcal M``, ``p∈ \mathcal M``, ``Y∈ T_p\mathcal M`` is a smooth mapping
associated to a retraction ``\operatorname{retr}_p(Y) = q`` such that

1. (associated retraction) ``\mathcal T_{p,d}X ∈ T_q\mathcal M`` if and only if ``q = \operatorname{retr}_p(d)``.
2. (consistency) ``\mathcal T_{p,0_p}X = X`` for all ``X∈T_p\mathcal M``
3. (linearity) ``\mathcal T_{p,d}(αX+βY) = α\mathcal T_{p,d}X + β\mathcal T_{p,d}Y``

For the [`AbstractVectorTransportMethod`](@ref) we might even omit the third point.
The [`AbstractLinearVectorTransportMethod`](@ref)s are linear.

# Input Parameters
* `M` a manifold
* `p` indicating the tangent space of
* `X` the tangent vector to be transported
* `d` indicating a transport direction (and distance through its length)
* `m` an [`AbstractVectorTransportMethod`](@ref), by default [`default_vector_transport_method`](@ref), so usually [`ParallelTransport`](@ref)

Usually this method requires a [`AbstractRetractionMethod`](@ref) as well.
By default this is assumed to be the [`default_retraction_method`](@ref) or
implicitly given (and documented) for a vector transport.
To explicitly distinguish different retractions for a vector transport,
see [`VectorTransportDirection`](@ref).

Instead of spcifying a start direction `d` one can equivalently also specify a target tanget space
``T_q\mathcal M``, see [`vector_transport_to`](@ref).
By default [`vector_transport_direction`](@ref) falls back to using [`vector_transport_to`](@ref),
using the [`default_retraction_method`](@ref) on `M`.

[^AbsilMahonySepulchre2008]:
    > Absil, P.-A., Mahony, R. and Sepulchre R.,
    > _Optimization Algorithms on Matrix Manifolds_
    > Princeton University Press, 2008,
    > doi: [10.1515/9781400830244](https://doi.org/10.1515/9781400830244)
    > [open access](http://press.princeton.edu/chapters/absil/)
"""
function vector_transport_direction(
    M::AbstractManifold,
    p,
    X,
    d,
    m::AbstractVectorTransportMethod = default_vector_transport_method(M),
)
    return _vector_transport_direction(M, p, X, d, m)
end
function _vector_transport_direction(
    M::AbstractManifold,
    p,
    X,
    d,
    m::AbstractVectorTransportMethod = default_vector_transport_method(M),
)
    r = default_retraction_method(M)
    return vector_transport_to(M, p, X, retract(M, p, d, r), m)
end
function _vector_transport_direction(
    M::AbstractManifold,
    p,
    X,
    d,
    m::VectorTransportDirection,
)
    return vector_transport_to(M, p, X, retract(M, p, d, m.retraction), m.vector_transport)
end
function _vector_transport_direction(
    M::AbstractManifold,
    p,
    X,
    d,
    m::DifferentiatedRetractionVectorTransport{R},
) where {R<:AbstractRetractionMethod}
    return vector_transport_direction_diff(M, p, X, d, m.retraction)
end
@doc raw"""
    vector_transport_direction_diff(M::AbstractManifold, p, X, d, m::AbstractRetractionMethod)

Compute the vector transport of `X` from ``T_p\mathcal M`` into the direction `d`
using the differential of the [`AbstractRetractionMethod`](@ref) `m`.
"""
function vector_transport_direction_diff(
    M::AbstractManifold,
    p,
    X,
    d,
    r::AbstractRetractionMethod,
)
    Y = allocate_result(M, vector_transport_direction, p, X, d)
    return vector_transport_direction_diff!(M, Y, p, X, d, r)
end
function _vector_transport_direction(M::AbstractManifold, p, X, d, ::ParallelTransport)
    return parallel_transport_direction(M, p, X, d)
end

"""
    vector_transport_direction!(M::AbstractManifold, Y, p, X, d)
    vector_transport_direction!(M::AbstractManifold, Y, p, X, d, m::AbstractVectorTransportMethod)

Transport a vector `X` from the tangent space at a point `p` on the [`AbstractManifold`](@ref) `M`
in the direction indicated by the tangent vector `d` at `p`. By default, [`retract`](@ref) and
[`vector_transport_to!`](@ref) are used with the `m` and `r`, which default
to [`default_vector_transport_method`](@ref)`(M)` and [`default_retraction_method`](@ref)`(M)`, respectively.
The result is saved to `Y`.

See [`vector_transport_direction`](@ref) for more details.
"""
function vector_transport_direction!(
    M::AbstractManifold,
    Y,
    p,
    X,
    d,
    m::AbstractVectorTransportMethod = default_vector_transport_method(M),
)
    return _vector_transport_direction!(M, Y, p, X, d, m)
end
function _vector_transport_direction!(
    M::AbstractManifold,
    Y,
    p,
    X,
    d,
    m::AbstractVectorTransportMethod = default_vector_transport_method(M),
)
    r = default_retraction_method(M)
    return vector_transport_to!(M, Y, p, X, retract(M, p, d, r), m)
end
function _vector_transport_direction!(
    M::AbstractManifold,
    Y,
    p,
    X,
    d,
    m::VectorTransportDirection,
)
    return vector_transport_to!(
        M,
        Y,
        p,
        X,
        retract(M, p, d, m.retraction),
        m.vector_transport,
    )
end
function _vector_transport_direction!(
    M::AbstractManifold,
    Y,
    p,
    X,
    d,
    m::DifferentiatedRetractionVectorTransport,
)
    return vector_transport_direction_diff!(M, Y, p, X, d, m.retraction)
end
@doc raw"""
    vector_transport_direction_diff!(M::AbstractManifold, Y, p, X, d, m::AbstractRetractionMethod)

Compute the vector transport of `X` from ``T_p\mathcal M`` into the direction `d`
using the differential of the [`AbstractRetractionMethod`](@ref) `m` in place of `Y`.
"""
vector_transport_direction_diff!(M, Y, p, X, d, m)

function vector_transport_direction_diff! end

function _vector_transport_direction!(M::AbstractManifold, Y, p, X, d, ::ParallelTransport)
    return parallel_transport_direction!(M, Y, p, X, d)
end

@doc raw"""
    vector_transport_to(M::AbstractManifold, p, X, q)
    vector_transport_to(M::AbstractManifold, p, X, q, m::AbstractVectorTransportMethod)
    vector_transport_to(M::AbstractManifold, p, X, q, m::AbstractVectorTransportMethod)

Transport a vector `X` from the tangent space at a point `p` on the [`AbstractManifold`](@ref) `M`
along a curve implicitly given by an [`AbstractRetractionMethod`](@ref) associated to `m`.
By default `m` is the [`default_vector_transport_method`](@ref)`(M)`.
To explicitly specify a (different) retraction to the implicitly assumeed retraction, see [`VectorTransportTo`](@ref).
Note that some vector transport methods might also carry their own retraction they are associated to,
like the  [`DifferentiatedRetractionVectorTransport`](@ref) and some are even independent of the retraction, for example the [`ProjectionTransport`](@ref).

This method is equivalent to using ``d = \operatorname{retr}^{-1}_p(q)`` in [`vector_transport_direction`](@ref)`(M, p, X, q, m, r)`,
where you can find the formal definition. This is the fallback for [`VectorTransportTo`](@ref).
"""
function vector_transport_to(
    M::AbstractManifold,
    p,
    X,
    q,
    m::AbstractVectorTransportMethod = default_vector_transport_method(M),
)
    return _vector_transport_to(M, p, X, q, m)
end
function _vector_transport_to(M::AbstractManifold, p, X, q, m::VectorTransportTo)
    d = inverse_retract(M, p, q, m.inverse_retraction)
    return vector_transport_direction(M, p, X, d, m.vector_transport)
end
function _vector_transport_to(M::AbstractManifold, p, X, q, ::ParallelTransport)
    return parallel_transport_to(M, p, X, q)
end
function _vector_transport_to(
    M::AbstractManifold,
    p,
    X,
    q,
    m::DifferentiatedRetractionVectorTransport,
)
    return vector_transport_to_diff(M, p, X, q, m.retraction)
end
@doc raw"""
    vector_transport_to_diff(M::AbstractManifold, p, X, q, r)

Compute a vector transport by using a [`DifferentiatedRetractionVectorTransport`](@ref) `r`.
"""
function vector_transport_to_diff(M::AbstractManifold, p, X, q, r)
    Y = allocate_result(M, vector_transport_to, X, p)
    return vector_transport_to_diff!(M, Y, p, X, q, r)
end
@doc raw"""
    vector_transport_to_diff(M::AbstractManifold, p, X, q, r)

Compute a vector transport by using a [`DifferentiatedRetractionVectorTransport`](@ref) `r` in place of `Y`.
"""
vector_transport_to_diff!(M::AbstractManifold, Y, p, X, q, r)

function vector_transport_to_diff! end

function _vector_transport_to(M::AbstractManifold, p, X, q, ::ProjectionTransport)
    return vector_transport_to_project(M, p, X, q)
end
@doc raw"""
    vector_transport_to_project(M::AbstractManifold, p, X, q)

Compute a vector transport by projecting ``X\in T_p\mathcal M`` onto the tangent
space ``T_q\mathcal M`` at ``q``.
"""
function vector_transport_to_project(M::AbstractManifold, p, X, q)
    Y = allocate_result(M, vector_transport_to, X, p)
    return vector_transport_to_project!(M, Y, p, X, q)
end
@doc raw"""
    vector_transport_to_project!(M::AbstractManifold, Y, p, X, q)

Compute a vector transport by projecting ``X\in T_p\mathcal M`` onto the tangent
space ``T_q\mathcal M`` at ``q`` in place of `Y`.
"""
function vector_transport_to_project!(M::AbstractManifold, Y, p, X, q)
    # Note that we have to use embed (not embed!) since we do not have memory to store this embedded value in
    return project!(M, Y, q, embed(M, p, X))
end

"""
    vector_transport_to!(M::AbstractManifold, Y, p, X, q)
    vector_transport_to!(M::AbstractManifold, Y, p, X, q, m::AbstractVectorTransportMethod)

Transport a vector `X` from the tangent space at a point `p` on the [`AbstractManifold`](@ref) `M`
to `q` using the [`AbstractVectorTransportMethod`](@ref) `m` and the [`AbstractRetractionMethod`](@ref) `r`.

The result is computed in `Y`.
See [`vector_transport_to`](@ref) for more details.
"""
function vector_transport_to!(
    M::AbstractManifold,
    Y,
    p,
    X,
    q,
    m::AbstractVectorTransportMethod = default_vector_transport_method(M),
)
    return _vector_transport_to!(M, Y, p, X, q, m)
end
function _vector_transport_to!(M::AbstractManifold, Y, p, X, q, m::VectorTransportTo)
    d = inverse_retract(M, p, q, m.inverse_retraction)
    return vector_transport_direction!(M, Y, p, X, d, m.vector_transport)
end
function _vector_transport_to!(M::AbstractManifold, Y, p, X, q, ::ParallelTransport)
    return parallel_transport_to!(M, Y, p, X, q)
end
function _vector_transport_to!(
    M::AbstractManifold,
    Y,
    p,
    X,
    q,
    m::DifferentiatedRetractionVectorTransport,
)
    return vector_transport_to_diff!(M, Y, p, X, q, m.retraction)
end

function _vector_transport_to!(M::AbstractManifold, Y, p, X, q, ::ProjectionTransport)
    return vector_transport_to_project!(M, Y, p, X, q)
end

function _vector_transport_to(M::AbstractManifold, p, X, q, m::PoleLadderTransport)
    Y = allocate_result(M, vector_transport_to, X, p)
    return _vector_transport_to!(M, Y, p, X, q, m)
end
function _vector_transport_to!(M::AbstractManifold, Y, p, X, q, m::PoleLadderTransport)
    inverse_retract!(
        M,
        Y,
        q,
        pole_ladder(
            M,
            p,
            retract(M, p, X, m.retraction),
            q;
            retraction = m.retraction,
            inverse_retraction = m.inverse_retraction,
        ),
        m.inverse_retraction,
    )
    copyto!(Y, -Y)
    return Y
end

function _vector_transport_to(M::AbstractManifold, p, X, c, m::ScaledVectorTransport)
    Y = allocate_result(M, vector_transport_to, X, p)
    return _vector_transport_to!(M, Y, p, X, c, m)
end
function _vector_transport_to!(M::AbstractManifold, Y, p, X, q, m::ScaledVectorTransport)
    vector_transport_to!(M, Y, p, X, q, m.method)
    Y .*= norm(M, p, X) / norm(M, q, Y)
    return Y
end
function _vector_transport_to(M::AbstractManifold, p, X, c, m::SchildsLadderTransport)
    Y = allocate_result(M, vector_transport_to, X, p)
    return _vector_transport_to!(M, Y, p, X, c, m)
end
function _vector_transport_to!(M::AbstractManifold, Y, p, X, q, m::SchildsLadderTransport)
    return inverse_retract!(
        M,
        Y,
        q,
        schilds_ladder(
            M,
            p,
            retract(M, p, X, m.retraction),
            q;
            retraction = m.retraction,
            inverse_retraction = m.inverse_retraction,
        ),
        m.inverse_retraction,
    )
end


text:
# Install PyCall
import Pkg
Pkg.add("PyCall")
Pkg.add("IJulia")

ext: jl

max_stars_repo_name: APIOxy/APIOxy

max_issues_repo_name: APIOxy/APIOxy

max_forks_repo_name: APIOxy/APIOxy


text:
using PackageCompiler, Test

@testset "compilage_package" begin
    @testset "FixedPointNumbers" begin
        sysimage = PackageCompiler.compile_package("FixedPointNumbers", verbose = true)
        test_code = """
        using FixedPointNumbers; N0f8(0.5); println("no segfaults, yay")
        """
        cmd = PackageCompiler.julia_code_cmd(test_code, J = sysimage)
        @test read(cmd, String) == "no segfaults, yay\n"
    end
    @testset "FixedPointNumbers ColorTypes" begin
        sysimage = PackageCompiler.compile_package("FixedPointNumbers", "ColorTypes", verbose = true)
        test_code = """
        using FixedPointNumbers, ColorTypes; N0f8(0.5); RGB(0.0, 0.0, 0.0); println("no segfaults, yay")
        """
        cmd = PackageCompiler.julia_code_cmd(test_code, J = sysimage)
        @test read(cmd, String) == "no segfaults, yay\n"
    end
end

@testset "compile_incremental" begin
    @testset "unregistered package with runtests.jl" begin
        path = joinpath(@__DIR__, "TestPackage")
        push!(Base.LOAD_PATH, path)
        syso, syso_old = PackageCompiler.compile_incremental(:TestPackage)
        test_code = """
        push!(Base.LOAD_PATH, $(repr(path)))
        using TestPackage; TestPackage.greet()
        """
        cmd = PackageCompiler.julia_code_cmd(test_code, J = syso)
        @test read(cmd, String) == "Hello World!"
        pop!(Base.LOAD_PATH)
    end
    @testset "unregistered package with snoopfile.jl" begin
        path = joinpath(@__DIR__, "TestPackage2")
        push!(Base.LOAD_PATH, path)
        syso, syso_old = PackageCompiler.compile_incremental(:TestPackage2)
        test_code = """
        push!(Base.LOAD_PATH, $(repr(path)))
        using TestPackage2; TestPackage2.greet()
        """
        cmd = PackageCompiler.julia_code_cmd(test_code, J = syso)
        @test read(cmd, String) == "Hello World!"
        pop!(Base.LOAD_PATH)
    end
    @testset "FixedPointNumbers" begin
        # This is the new compile_package
        syso, syso_old = PackageCompiler.compile_incremental(:FixedPointNumbers)
        test_code = """
        using FixedPointNumbers; N0f8(0.5); println("no segfaults, yay")
        """
        cmd = PackageCompiler.julia_code_cmd(test_code, J = syso)
        @test read(cmd, String) == "no segfaults, yay\n"
    end
    @testset "FixedPointNumbers ColorTypes" begin
        syso, syso_old = PackageCompiler.compile_incremental(:FixedPointNumbers, :ColorTypes)
        test_code = """
        using FixedPointNumbers, ColorTypes; N0f8(0.5); RGB(0.0, 0.0, 0.0); println("no segfaults, yay")
        """
        cmd = PackageCompiler.julia_code_cmd(test_code, J = syso)
        @test read(cmd, String) == "no segfaults, yay\n"
    end
    @testset "JSON with Distributed blacklisted" begin
        # This is the new compile_package
        syso, syso_old = PackageCompiler.compile_incremental(:JSON, blacklist=[:Distributed])
        test_code = """
        using JSON
        s = \"{\\"a_number\\" : 5.0, \\"an_array\\" : [\\"string\\", 9]}\"
        j = JSON.parse(s)
        println(\"no segfaults, yay\")
        """
        cmd = PackageCompiler.julia_code_cmd(test_code, J = syso)
        @test read(cmd, String) == "no segfaults, yay\n"
    end
end

julia = Base.julia_cmd().exec[1]

@testset "build_executable" begin
    jlfile = joinpath(@__DIR__, "..", "examples", "hello.jl")
    basedir = mktempdir()
    relativebuilddir = "build"
    cd(basedir) do
        mkdir(relativebuilddir)
        snoopfile = "snoop.jl"
        open(snoopfile, "w") do io
            write(io, open(read, jlfile))
            println(io)
            println(io, "using .Hello; Hello.julia_main(String[])")
        end
        build_executable(
            jlfile, snoopfile = snoopfile, builddir = relativebuilddir, verbose = true
        )
    end
    builddir = joinpath(basedir, relativebuilddir)
    @test isfile(joinpath(builddir, "hello.$(PackageCompiler.Libdl.dlext)"))
    @test isfile(joinpath(builddir, "hello$executable_ext"))
    @test startswith(read(`$(joinpath(builddir, "hello$executable_ext"))`, String), "hello, world")
    for i = 1:100
        # Windows seems to have problems with removing files - it can error
        # making this test fail.
        try rm(basedir, recursive = true) catch end
        sleep(1/100)
    end
end

@testset "program.c" begin
    @testset "args" begin
        basedir = mktempdir();
        argsjlfile = mktemp(basedir)[1];
        write(argsjlfile, raw"""
            Base.@ccallable function julia_main(argv::Vector{String})::Cint
                println("@__FILE__: $(@__FILE__)")
                println("PROGRAM_FILE: $(PROGRAM_FILE)")
                println("argv: $(argv)")
                # Sometimes code accesses ARGS directly, as a global
                println("ARGS: $ARGS")
                println("Base.ARGS: $(Base.ARGS)")
                println("Core.ARGS: $(Core.ARGS)")
                return 0
            end
        """)
        builddir = joinpath(basedir, "builddir")
        outname = "args"
        build_executable(
            argsjlfile, outname; builddir = builddir
        )
        # Check that the output from the program is as expected:
        exe = joinpath(builddir, outname*executable_ext)
        output = read(`$exe a b c`, String)
        println(output)
        @test output ==
            "@__FILE__: $argsjlfile\n" *
            "PROGRAM_FILE: $exe\n" *
            "argv: [\"a\", \"b\", \"c\"]\n" *
            "ARGS: [\"a\", \"b\", \"c\"]\n" *
            "Base.ARGS: [\"a\", \"b\", \"c\"]\n" *
            "Core.ARGS: Any[\"$(Sys.iswindows() ? replace(exe, "\\" => "\\\\") : exe)\", \"a\", \"b\", \"c\"]\n"
    end
end

@testset "juliac" begin
    mktempdir() do builddir
        juliac = joinpath(@__DIR__, "..", "juliac.jl")
        jlfile = joinpath(@__DIR__, "..", "examples", "hello.jl")
        cfile = joinpath(@__DIR__, "..", "examples", "program.c")
        @test success(`$julia $juliac -vaej $jlfile $cfile --builddir $builddir`)
        @test isfile(joinpath(builddir, "hello.$(PackageCompiler.Libdl.dlext)"))
        @test isfile(joinpath(builddir, "hello$executable_ext"))
        @test success(`$(joinpath(builddir, "hello$executable_ext"))`)
        @testset "--cc-flag" begin
            # Try passing `--help` to $cc. This should work for any system compiler.
            # Then grep the output for "-g", which should be present on any system.
            @test occursin("-g", read(`$julia $juliac -se --cc-flag="--help" $jlfile $cfile --builddir $builddir`, String))
            # Just as a control, make sure that without passing '--help', we don't see "-g"
            @test !occursin("-g", read(`$julia $juliac -se $jlfile $cfile --builddir $builddir`, String))
        end
    end
end


text:
# System module.
module System

LIB_SKL_AVAILABLE = false
LIB_CRT_AVAILABLE = false

export LIB_SKL_AVAILABLE,
       LIB_CRT_AVAILABLE

#using RCall
#using Conda
#
#import PyCall: pyimport, pycall

#function check_py_dep(package::AbstractString)
#  is_available = true
#  try
#    pyimport(package)
#  catch
#    try
#      Conda.add(package)
#      is_available = true
#    catch
#      is_available = false
#    end
#  end
#  return is_available
#end
#
#function check_r_dep()
#  is_available = true
#  try
#    R"library(caret)"
#  catch
#    try
#      R"install.packages('caret',repos='https://cloud.r-project.org',type='binary')"
#      R"install.packages('earth',repos='https://cloud.r-project.org',type='binary')"
#      R"install.packages('mda',repos='https://cloud.r-project.org',type='binary')"
#      R"install.packages('e1071',repos='https://cloud.r-project.org',type='binary')"
#      R"install.packages('gam',repos='https://cloud.r-project.org',type='binary')"
#      R"install.packages('randomForest',repos='https://cloud.r-project.org',type='binary')"
#      R"install.packages('nnet',repos='https://cloud.r-project.org',type='binary')"
#      R"install.packages('kernlab',repos='https://cloud.r-project.org',type='binary')"
#      R"install.packages('grid',repos='https://cloud.r-project.org',type='binary')"
#      R"install.packages('MASS',repos='https://cloud.r-project.org',type='binary')"
#      R"install.packages('pls',repos='https://cloud.r-project.org',type='binary')"
#      R"install.packages('xgboost',repos='https://cloud.r-project.org',type='binary')"
#      is_available = true
#    catch
#      is_available = false
#    end
#  end
#  return is_available
#end

# disable support for caret/scikitlearn for faster building/development
# use julia native libraries only

### Check system for python dependencies.
#if "LOAD_SK_CARET" in keys(ENV) && ENV["LOAD_SK_CARET"] == "true"
#  LIB_SKL_AVAILABLE = check_py_dep("scikit-learn")
#  LIB_CRT_AVAILABLE = check_r_dep()
#elseif "LOAD_SK_CARET" in keys(ENV) && ENV["LOAD_SK_CARET"] == "false"
#  LIB_SKL_AVAILABLE = false
#  LIB_CRT_AVAILABLE = false
#else
#  LIB_SKL_AVAILABLE = check_py_dep("scikit-learn")
#  LIB_CRT_AVAILABLE = check_r_dep()
#end

end # module


text:
## change working directory

cd(joinpath(ENV["HOME"],"Documents","debugging","racecarimgs","5000","easy"))

## Load packages

using FreeTypeAbstraction
using AprilTags
using Caesar
using Images, ImageView, ImageDraw
using FileIO
using LinearAlgebra
using Cairo, Fontconfig, RoMEPlotting
Gadfly.set_default_plot_size(25cm, 20cm)

## make multiprocess

using Distributed

addprocs(4)
using FreeTypeAbstraction, AprilTags, Images, ImageView, FileIO, Caesar
@everywhere using FreeTypeAbstraction, AprilTags, Images, ImageView, FileIO, Caesar

## which images

imgFiles = [
  "1_4cb1b54b-e297-41ac-8a81-f0ddab02dac2.png",
  "2_ddaf9899-4c6a-4182-a7ad-4cc9be6f127d.png",
  "3_077227d4-e6cd-4884-befd-a0d73a07c172.png",
  "4_7545b43e-0731-42ee-8b57-df03755abbd6.png",
  "5_50ce6f66-1546-4f90-a7e5-336e595f14a2.png",
  "6_7676be16-8f57-4a2e-ac2b-1528747213d1.png",
  "7_c7c38417-eecd-4e37-91fc-4628e065ac66.png",
  "8_e42c4cc8-25b1-4df9-8790-61f8212e9ca2.png"
]

## load images into memory

imgs = load.(imgFiles)

rows, cols = size(imgs[1])
c_width = cols/2
c_height = rows/2
# f_width = 500.0
# f_height = f_width
taglength = 0.172

# f_width = 496.70305287384673
# c_height = 311.6691866287029

f_width = 504.91911888349216
# c_height = 324.2630820675616

f_height = f_width

K = Caesar._defaultCameraCalib(f_width=f_width, f_height=f_height, c_width=c_width, c_height=c_height)[1:2,:]

## test to see

imshow(vcat(hcat(imgs[1:3]...),hcat(imgs[4:6]...)))

## Run AprilTag detector (app)

detector = AprilTagDetector()
allTags = deepcopy(detector.(imgs))

## init the factor graph

fg = generateCanonicalFG_ZeroPose2()

getSolverParams(fg).useMsgLikelihoods = true

#

posecount = 0
prevSym = :x0

for tags in allTags[1:5]
  # new variable in factor graph
  posecount += 1
  poseSym = Symbol("x$posecount")
  addVariable!(fg, poseSym, Pose2, tags=[:POSE;])

  # add basic odo relative
  addFactor!(fg, [prevSym; poseSym], Pose2Pose2(MvNormal(zeros(3),diagm([1.0;0.5;0.5]))), tags=[:ODOMETRY;])
  for tag in tags
    @show tag.id
    tagSym = Symbol("tag$(tag.id)")
    if !exists(fg, tagSym)
      addVariable!(fg, tagSym, Pose2, tags=[:APRILTAG;])
    end

    apt4 = Pose2AprilTag4Corners(corners=tag.p, homography=tag.H, f_width=f_width, c_width=c_width, c_height=c_height, taglength=taglength, id=tag.id)
    addFactor!(fg, [poseSym; tagSym], apt4, tags=[:SIGHTING;])
  end
  prevSym = poseSym
end

## solve the graph

for i in 1:1
  solveTree!(fg, storeOld=true);
end

## solving and looking at the subgraph

# listSolveKeys(fg)

plp = plotSLAM2DPoses(fg, contour=false, drawEllipse=true)

# development on plot landmarks

pll = plotSLAM2DLandmarks(fg, regexLandmark=r"tag", contour=false, drawEllipse=true)

#

union!(plp.layers, pll.layers)
# co = Coord.Cartesian(xmax=4)
# plp.coord = co
plp

##

pl = plotSLAM2D(fg, regexLandmark=r"tag", contour=false, drawEllipse=true)

##

freeDetector!(detector)

## dev

imshow( drawTags(imgs[1], K) )

## look at one sighting factor

pred, meas = manikde!.(approxDeconv(fg, :x1tag16f1), Pose2)

plotPose(Pose2(), [meas;pred], c=["red";"green"])

## draw back for validation about what is going on

imgs_ = deepcopy(imgs);

for i in 1:5
  @show fSym = intersect(lsf(fg, Pose2AprilTag4Corners), ls(fg, Symbol("x$i")))
  for fct in getFactorType.(fg, fSym)
    drawBearingLinesAprilTags!( imgs_[i], fct,
                                f_width=f_width, c_width=c_width, taglength=taglength);
    #
  end
end

imshow( vcat(hcat( imgs_[1:3]... ), hcat( imgs_[4:6]... )) )

##

# img = deepcopy(imgs[1])

# drawBearingLinesAprilTags!(img, f_width=f_width, c_height=c_height, taglength=taglength)

# drawBearingLinesAprilTags!(deepcopy(imgs[1]), allTags[1], f_width=f_width, c_width=c_width, taglength=taglength) |> imshow
# drawBearingLinesAprilTags!(deepcopy(imgs[1]), f_width=f_width, c_width=c_width, taglength=taglength) |> imshow

# drawBearingLinesAprilTags!(deepcopy(imgs[2]), f_width=f_width, c_width=c_width, taglength=taglength) |> imshow
# drawBearingLinesAprilTags!(deepcopy(imgs[3]), f_width=f_width, c_width=c_width, taglength=taglength) |> imshow
# drawBearingLinesAprilTags!(deepcopy(imgs[4]), f_width=f_width, c_width=c_width, taglength=taglength) |> imshow

##

@show getFactorType(fg, :x1tag10f1).Zij.z.μ

allTags[1][2].id
Pose2AprilTag4Corners( corners=allTags[1][2].p, homography=allTags[1][2].H, f_width=f_width, c_width=c_width, c_height=c_height, taglength=taglength).Zij.z.μ

##

ls(fg, :x4)

##

plotLocalProduct(fg, :x4, dims=[1;2])

plotLocalProduct(fg, :x5, dims=[1;2])

plotLocalProduct(fg, :tag10, dims=[1;2])
plotLocalProduct(fg, :tag10, dims=[3])

## Some development tests

# accept default cost function argument mapping `args = (x)->[x[1]; x[1]; c_width; x[2]; taglength]`
#   such that `cost([f_width; c_height])`
cost = generateCostAprilTagsPreimageCalib(fg, f_width=f_width, f_height=f_height, c_width=c_width, c_height=c_height)

cost([f_width; c_height])

##

using Optim

##

result = optimize(cost, [f_width; c_height], BFGS())

result.minimizer

f_width, c_width

##

PL = []
for sk in [Symbol("default_$i") for i in 0:5]
  push!(PL, plotSLAM2D(fg,xmin=-2,xmax=7,ymin=-3,ymax=3, solveKey=sk))
end

plimgs = convert.(Matrix{RGB}, PL);

writevideo("pl100.avi", plimgs, fps=5)

##

## build subgraph [:x0; .. :x3;]

varfct = setdiff( union( ls.(fg, [:x0, :x1, :x2, :x3])...),  [:x4, :x3x4f1] )
vars = union(ls.(fg, varfct)...)
union!(varfct, vars)

sfg = buildSubgraph(fg,varfct)

##

for i in 1:5
  solveTree!(sfg, storeOld=true);
end

## subgraph 2

sfg = buildSubgraph(fg, [:x4; :x5; :tag10; :tag16], 1)

ls(sfg, :tag16)

## freeze the tags

# setMarginalize!()

##


text:
# Autogenerated wrapper script for libcgal_julia_jll for i686-w64-mingw32-cxx11
export libcgal_julia_exact, libcgal_julia_inexact

using CGAL_jll
using libcxxwrap_julia_jll
JLLWrappers.@generate_wrapper_header("libcgal_julia")
JLLWrappers.@declare_library_product(libcgal_julia_exact, "libcgal_julia_exact.dll")
JLLWrappers.@declare_library_product(libcgal_julia_inexact, "libcgal_julia_inexact.dll")
function __init__()
    JLLWrappers.@generate_init_header(CGAL_jll, libcxxwrap_julia_jll)
    JLLWrappers.@init_library_product(
        libcgal_julia_exact,
        "bin\\libcgal_julia_exact.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcgal_julia_inexact,
        "bin\\libcgal_julia_inexact.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# this is a better implementatoin of Pkg.clone
import Base.Pkg
import Base.Pkg.Reqs
using Base.Pkg.Types
using Base.Pkg.Git

"""
  Adds a package.  This clones the package, installs its dependencies, and
  executes the package's build script.

  This function uses Julia Pkg.resolve() to install dependencies, which is
  known to be unreliable, but its the only option at the moment.

  **Inputs**

   * repo: repository URL

  **KeyWord Argument**
   * branch_ish: see [`clone`](@ref).  Note that if branch_ish is not specified,
     it is possible that installing dependencies may change the version of
     the package in the repository `repo`
"""
function add(repo::AbstractString; branch_ish="")

  repo_name = getRepoName(repo)

  if PkgFix.installed(repo_name) != nothing
    println("Package $repo_name already installed, not installing...")
    return nothing
  end

  clone(repo, branch_ish)
  if branch_ish != ""
    pin(repo_name)  # make sure the resolve step does not change the version
                    # of this package
  end
  resolve()  # add dependencies

  build(repo_name)  # build this package

  if branch_ish != ""
    free(repo_name)
  end

  return nothing
end

"""
  Clones a Julia package (and nothing more).  This function does not build
  the package nor install its dependencies

  **Inputs**

   * repo: the repository URL
   * branch_ish: a branch or tag to checkout after cloning.  This argument can
                 be anything accepted by `git checkout`.  This argument is
                 optional.

  **Implementation Notes**

  `branch_ish` is passed to [`PkgFix.checkout`](@ref), see that function
               for details

"""
function clone(repo::AbstractString, branch_ish="")

  repo_name = getRepoName(repo)
  repo_name_ext = repo_name*".jl"
  start_dir = pwd()

  pth1 = joinpath(Pkg.dir(), repo_name)
  pth2 = joinpath(Pkg.dir(), repo_name_ext)
  println("pth1 = ", pth1)
  println("pth2 = ", pth2)
  println("isdir(pth1) = ", isdir(pth1))
  println("isdir(pth2) = ", isdir(pth2))

  if isdir(joinpath( Pkg.dir(), repo_name)) || isdir(joinpath( Pkg.dir(), repo_name_ext))
    println(STDERR, "Package $repo_name already exists, not cloning...")
    return nothing
  end

  try
    cd(Pkg.dir())

    run(`git clone $repo`)

    # rename ot remove the .jl
    run(`mv $repo_name_ext $repo_name`)

    if branch_ish != ""
      checkout(repo_name, branch_ish)
    end

    # add the new repo to the global REQUIRE file
    # so it doesn't get removed the next time Pkg.add is run
    addRequirement(repo_name)
  catch x
    cd(start_dir)
    rethrow(x)
  end

  cd(start_dir)

  return nothing
end

"""
  Checks out a particular branch/tag/commit of a package.

  **Inputs**

   * pkg: the package name
   * branch_ish: the name of the thing to checkout (note that this argumetn
                 can be anything accepted by `git checkout`

  **Implementation Notes**

  If `branch_ish` is a commit, this function will checkout the commit to a new branch
  named `detached_from_branch_ish`, otherwise the `branch_ish` will be checkout
  out directly.

"""
function checkout(pkg::AbstractString, branch_ish)

  start_dir = pwd()

  try
    cd(joinpath(Pkg.dir(), pkg))

    if is_commit(pwd(), branch_ish)
      run(`git checkout -b detatched_from_$branch_ish`)
    else
      run(`git checkout $branch_ish`)
    end
  catch x
    cd(start_dir)
    rethrow(x)
  end

  cd(start_dir)

  return nothing
end

"""
  Given a git url, extracts the package name.  This assumes the url is of the
  format blah/blah/blah/Pkgname.jl.git
"""
function getRepoName(url::AbstractString)

  println("url = ", url)
  if !endswith(url, ".git")
    error("url $url does not end in .git")
  end

  url_stub = url[1:end-4]
  println("url_stub = ", url_stub)

  if !endswith(url_stub, ".jl")
    error("repository name does not end in .jl")
  end

  url_stub = url_stub[1:end-3]
  println("url_stub = ", url_stub)

  # find the final slash (which preceeds the start of the package name)
  idx = 0
  for i=length(url_stub):-1:1
    if url_stub[i] == '/'
      idx = i + 1  # record start of repo name
      break
    end
  end

  repo_name = url_stub[idx:end]

  return repo_name
end

"""
  Add a package to the global REQUIRE file (if it is not already present).

  This function does not check if the directory for the package exists, it
  blindly adds to the REQUIRE file.

  This function does not support adding version requirements

  **Inputs**

   * repo_name: name of the package (not including the .jl extension)
"""
function addRequirement(repo_name::AbstractString)

  f = Pkg.Reqs.read(joinpath(Pkg.dir(), "REQUIRE"))
  already_required = false
  for i=1:length(f)
    if f[i].package == repo_name
      already_required = true
      break
    end
  end

  if !already_required
    line = Reqs.Requirement(repo_name, VersionSet())
    line_arr = Reqs.Line[line]
    f2 = open("REQUIRE", "a")
    Reqs.write(f2, line_arr)
    close(f2)
  end

  return nothing
end

"""
  This function determines if a git identifier is a commit or something else
  (ie. a branch or tag)

  **Inputs**

   * pkg_dir: path to directory containing the repo
   * hash: git object identifier

  **Outputs**

   * bool, true if hash represents a commit, false otherwise
"""
function is_commit(pkg_dir::AbstractString, hash::AbstractString)

  start_dir = pwd()
  cd(pkg_dir)

  iscommit = false
  try
    vals = readall(`git show-ref | grep $hash`)
    # make sure this is an exact match
    name = split(vals, '/')[end]
    if name[1:end-1] == hash
      iscommit = false
    else
      iscommit = true
    end
  catch x
    iscommit = true
  end

  cd(start_dir)
  return iscommit
end

"""
  This function pins a package so that the Julia package manager will not change
  its version.  Note that packages pinned by this function must be freed by
  PkgFix.free and not by Pkg.free.  Similarly, PkgFix.free cannot free packakages
  pinned by Pkg.pin.  You can tell if a package was pinned by PkgFix by inspecting
  the branch name shown by Pkg.status for `PKGFIX`.

  This function does not change the version of any other package.

  **Inputs**

   * pkg: the package name
"""
function pin(pkg::AbstractString)

  start_dir=pwd()
  try
    cd(joinpath(Pkg.dir(), pkg))
    # it appears as long as the branch names starts with pinned., then the
    # package manager considers it pinned

    head = getHeadIdentifier(pkg)
    println("head = ", head)
    run(`git checkout -b pinned.PKGFIX$head.tmp`)
  catch x
    cd(start_dir)
    rethrow(x)
  end

  cd(start_dir)

  return nothing
end

"""
  This function frees a package pinned by PkgFree.pin.  It cannot free a
  a package pinned by the Julia package manager (it prints a warning and exits
  without making changes in this case).

  This function will restore the package to its state previous to pinning
  (if it was on a branch, it will return to that branch).

  This function does not change the version of any other packages.

  **Inputs**

   * pkg: the package name
"""
function free(pkg::AbstractString)
  start_dir = pwd()
  try
    cd(joinpath(Pkg.dir(), pkg))
    head = getHeadIdentifier(pkg)

    # make sure the format of head is the format written by pin
    if !startswith(head, "pinned")
      println(STDERR, "Warning: package $pkg is not pinned, not freeing...")
      cd(start_dir)
      return nothing
    end

    # remove the prefix pinned.
    head = head[8:end]

    # remove the suffix .tmp
    head = head[1:end-4]

    if !startswith(head, "PKGFIX")
      println(STDERR, "Warning: package was not pinned by PkgFix, not freeing...")
      cd(start_dir)
      return nothing
    end

    # if we get here, this package was pinned by PkgFee, so restore it
    # to the previous state
    old_head = head[7:end]
    run(`git checkout $old_head`)

  catch x
    cd(start_dir)
    rethrow(x)
  end

  cd(start_dir)

  return nothing
end

function getHeadIdentifier(pkg::AbstractString)

  start_dir = pwd()

  str = ""
  try
    cd(joinpath(Pkg.dir(), pkg))
    str = readchomp(`git rev-parse --abbrev-ref HEAD`)
    println("str = ", str)

    if str == "HEAD"  # the current head is not a branch
      str = readchomp(`git rev-parse HEAD`)
      println("str = ", str)
    end
  catch x
    cd(start_dir)
    rethrow(x)
  end

  cd(start_dir)
  return str
end

"""
  Checks if a package is installed, returning Void if not and
  the version number if true

  **Inputs**

   * pkg: the package name
"""
function installed(pkg::AbstractString)
  if isdir(joinpath(Pkg.dir(), pkg))
    return Pkg.installed(pkg)
  else
    return nothing
  end
end

# functions that delegate to things in Pkg
function status()
  Pkg.status()
end

function rm(pkg::AbstractString)
  Pkg.rm(pkg)
end

function installed()
  Pkg.installed()
end

function resolve()
  Pkg.resolve()
end

function build(pkg::AbstractString)
  Pkg.build(pkg)
end

function dir()
  Pkg.dir()
end

function dir(pkg::AbstractString)
  Pkg.dir(pkg)
end


text:
"""
 Holds the tableau of a stochastic implicit partitioned Runge-Kutta method.
 qdrift, pdrift hold the RK coefficients for the drift part,
 and qdiff, pdiff hold the RK coefficients for the diffusion part of the SDE.
"""
struct TableauSIPRK{T} <: AbstractTableauIRK{T}
    name::Symbol
    s::Int
    qdrift::CoefficientsRK{T}
    qdiff::CoefficientsRK{T}
    pdrift::CoefficientsRK{T}
    pdiff::CoefficientsRK{T}

    # Order of the tableau is not included, because unlike in the deterministic
    # setting, it depends on the properties of the noise (e.g., the dimension of
    # the Wiener process and the commutativity properties of the diffusion matrix)
    #
    # Orders stored in qdrift and qdiff are understood as the classical orders of these methods.

    function TableauSIPRK{T}(name, s, qdrift, qdiff, pdrift, pdiff) where {T}
        # THE COMMENTED OUT PART WAS FOR TableauFIRK. MAY IMPLEMENT SOMETHING
        # SIMILAR FOR TableauSFIRK LATER.

        # if (q.s > 1 && istrilstrict(q.a)) || (q.s==1 && q.a[1,1] == 0)
        #     warn("Initializing TableauFIRK with explicit tableau ", q.name, ".\n",
        #          "You might want to use TableauERK instead.")
        # elseif q.s > 1 && istril(q.a)
        #     warn("Initializing TableauFIRK with diagonally implicit tableau ", q.name, ".\n",
        #          "You might want to use TableauDIRK instead.")
        # end

        @assert s == qdrift.s == qdiff.s == pdrift.s == pdiff.s

        new(name, s, qdrift, qdiff, pdrift, pdiff)
    end
end

function TableauSIPRK(name::Symbol, qdrift::CoefficientsRK{T}, qdiff::CoefficientsRK{T}, pdrift::CoefficientsRK{T}, pdiff::CoefficientsRK{T}) where {T}
    TableauSIPRK{T}(name, qdrift.s, qdrift, qdiff, pdrift, pdiff)
end

function TableauSIPRK(name::Symbol, qorder_drift::Int, qa_drift::Matrix{T}, qb_drift::Vector{T}, qc_drift::Vector{T}, qorder_diff::Int, qa_diff::Matrix{T}, qb_diff::Vector{T}, qc_diff::Vector{T},
                                    porder_drift::Int, pa_drift::Matrix{T}, pb_drift::Vector{T}, pc_drift::Vector{T}, porder_diff::Int, pa_diff::Matrix{T}, pb_diff::Vector{T}, pc_diff::Vector{T}) where {T}
    TableauSIPRK{T}(name, length(qc_drift), CoefficientsRK(name, qorder_drift, qa_drift, qb_drift, qc_drift), CoefficientsRK(name, qorder_diff, qa_diff, qb_diff, qc_diff), CoefficientsRK(name, porder_drift, pa_drift, pb_drift, pc_drift), CoefficientsRK(name, porder_diff, pa_diff, pb_diff, pc_diff))
end

# TODO function readTableauSFIRKFromFile(dir::AbstractString, name::AbstractString)

"""
Parameters for right-hand side function of implicit Runge-Kutta methods.
  A - if positive, the upper bound of the Wiener process increments; if A=0.0, no truncation
"""
mutable struct ParametersSIPRK{DT, TT, ET <: PSDE{DT,TT}, D, M, S} <: Parameters{DT,TT}
    equ::ET
    tab::TableauSIPRK{TT}
    Δt::TT
    ΔW::Vector{DT}
    ΔZ::Vector{DT}
    A::DT

    t::TT
    q::Vector{DT}
    p::Vector{DT}
end

function ParametersSIPRK(equ::ET, tab::TableauSIPRK{TT}, Δt::TT, ΔW::Vector{DT}, ΔZ::Vector{DT}, A::DT) where {DT, TT, ET <: PSDE{DT,TT}}
    @assert equ.m == length(ΔW) == length(ΔZ)
    ParametersSIPRK{DT, TT, ET, equ.d, equ.m, tab.s}(equ, tab, Δt, ΔW, ΔZ, A, 0, zeros(DT, equ.d), zeros(DT, equ.d))
end

"""
Structure for holding the internal stages Q, the values of the drift vector
and the diffusion matrix evaluated at the internal stages VQ=v(Q), BQ=B(Q),
and the increments Y = Δt*a_drift*v(Q) + a_diff*B(Q)*ΔW
"""
struct NonlinearFunctionCacheSIPRK{DT}
    Q::Vector{Vector{DT}}
    P::Vector{Vector{DT}}
    V::Vector{Vector{DT}}
    F::Vector{Vector{DT}}
    B::Vector{Matrix{DT}}
    G::Vector{Matrix{DT}}
    Y::Vector{Vector{DT}}
    Z::Vector{Vector{DT}}

    v::Vector{DT}
    f::Vector{DT}
    b::Matrix{DT}
    g::Matrix{DT}
    y::Vector{DT}
    z::Vector{DT}

    function NonlinearFunctionCacheSIPRK{DT}(d, m, s) where {DT}

        # create internal stage vectors
        Q = create_internal_stage_vector(DT, d, s)
        P = create_internal_stage_vector(DT, d, s)
        V = create_internal_stage_vector(DT, d, s)
        F = create_internal_stage_vector(DT, d, s)
        B = create_internal_stage_vector(DT, d, m, s)
        G = create_internal_stage_vector(DT, d, m, s)
        Y = create_internal_stage_vector(DT, d, s)
        Z = create_internal_stage_vector(DT, d, s)

        # create velocity and update vector
        v = zeros(DT,d)
        f = zeros(DT,d)
        b = zeros(DT,d,m)
        g = zeros(DT,d,m)
        y = zeros(DT,d)
        z = zeros(DT,d)

        new(Q, P, V, F, B, G, Y, Z, v, f, b, g, y, z)
    end
end

"""
Unpacks the data stored in x = (Y[1][1], Y[1][2], ... Y[1][D], Y[2][1], ..., Z[1][1], Z[1][2], ... Z[1][D], Z[2][1], ...)
into Y, Z::Vector{Vector}, calculates the internal stages Q, P, the values of the RHS
of the SDE ( v(Q,P), f(Q,P), B(Q,P) and G(Q,P) ), and assigns them to V, F, B and G.
Unlike for FIRK, here
Y = Δt a_drift v(Q,P) + a_diff B(Q,P) ΔW,
Z = Δt ̃a_drift v(Q,P) + ̃a_diff B(Q,P) ΔW.
"""
function compute_stages!(x::Vector{ST}, Q::Vector{Vector{ST}}, P::Vector{Vector{ST}},
                                        V::Vector{Vector{ST}}, F::Vector{Vector{ST}},
                                        B::Vector{Matrix{ST}}, G::Vector{Matrix{ST}},
                                        Y::Vector{Vector{ST}}, Z::Vector{Vector{ST}},
                                        params::ParametersSIPRK{DT,TT,ET,D,M,S}) where {ST,DT,TT,ET,D,M,S}

    local tqᵢ::TT       #times for the q internal stages
    local tpᵢ::TT       #times for the p internal stages

    @assert S == length(Q) == length(V) == length(B) == length(P) == length(F) == length(G)

    # copy x to Y, Z and calculate Q, P
    for i in 1:S
        @assert D == length(Q[i]) == length(V[i]) == size(B[i],1) == length(P[i]) == length(F[i]) == size(G[i],1)
        @assert M == size(B[i],2) == size(G[i],2)

        for k in 1:D
            Y[i][k] = x[D*(  i-1)+k]
            Z[i][k] = x[D*(S+i-1)+k]
        end

        Q[i] .= params.q .+ Y[i]
        P[i] .= params.p .+ Z[i]
    end

    # compute V=v(t,Q,P), F=f(t,Q,P), B=B(t,Q,P) and G=G(t,Q,P)
    for i in 1:S
        tqᵢ = params.t + params.Δt * params.tab.qdrift.c[i]
        tpᵢ = params.t + params.Δt * params.tab.pdrift.c[i]
        # calculates v(t,Q,P) and f(t,Q,P) and assigns to the i-th column of V and F
        params.equ.v(tqᵢ, Q[i], P[i], V[i])
        params.equ.f(tpᵢ, Q[i], P[i], F[i])
        # calculates B(t,Q,P) and assigns to the matrices B[i] and G[i]
        params.equ.B(tqᵢ, Q[i], P[i], B[i])
        params.equ.G(tpᵢ, Q[i], P[i], G[i])
    end
end

"Compute stages of implicit Runge-Kutta methods."
@generated function function_stages!(x::Vector{ST}, b::Vector{ST}, params::ParametersSIPRK{DT,TT,ET,D,M,S}) where {ST,DT,TT,ET,D,M,S}

    cache = NonlinearFunctionCacheSIPRK{ST}(D, M, S)

    quote
        compute_stages!(x, $cache.Q, $cache.P, $cache.V, $cache.F, $cache.B, $cache.G, $cache.Y, $cache.Z, params)

        local y1::ST
        local y2::ST
        local z1::ST
        local z2::ST

        # compute b = - (Y-AV)
        for i in 1:S
            for k in 1:D
                y1 = 0
                y2 = 0
                z1 = 0
                z2 = 0
                for j in 1:S
                    y1 += params.tab.qdrift.a[i,j] * $cache.V[j][k] * params.Δt + params.tab.qdiff.a[i,j] * dot($cache.B[j][k,:], params.ΔW)
                    y2 += params.tab.qdrift.â[i,j] * $cache.V[j][k] * params.Δt + params.tab.qdiff.â[i,j] * dot($cache.B[j][k,:], params.ΔW)
                    z1 += params.tab.pdrift.a[i,j] * $cache.F[j][k] * params.Δt + params.tab.pdiff.a[i,j] * dot($cache.G[j][k,:], params.ΔW)
                    z2 += params.tab.pdrift.â[i,j] * $cache.F[j][k] * params.Δt + params.tab.pdiff.â[i,j] * dot($cache.G[j][k,:], params.ΔW)
                end
                b[D*(  i-1)+k] = - $cache.Y[i][k] + (y1 + y2)
                b[D*(S+i-1)+k] = - $cache.Z[i][k] + (z1 + z2)
            end
        end
    end
end

"Stochastic implicit partitioned Runge-Kutta integrator."
struct IntegratorSIPRK{DT, TT, PT <: ParametersSIPRK{DT,TT},
                               ST <: NonlinearSolver{DT}, N} <: StochasticIntegrator{DT,TT}
    params::PT
    solver::ST
    # InitialGuessPSDE not implemented for SIPRK
    #iguess::IT
    fcache::NonlinearFunctionCacheSIPRK{DT}

    q::Matrix{Vector{TwicePrecision{DT}}}
    p::Matrix{Vector{TwicePrecision{DT}}}
end

# K - the integer in the bound A = √(2 K Δt |log Δt|) due to Milstein & Tretyakov; K=0 no truncation
function IntegratorSIPRK(equation::PSDE{DT,TT,VT,BT,N}, tableau::TableauSIPRK{TT}, Δt::TT; K::Int=0) where {DT,TT,VT,BT,N}
    D = equation.d
    M = equation.m
    NS= equation.ns
    NI= equation.n
    S = tableau.s

    # create params
    K==0 ? A = 0.0 : A = sqrt( 2*K*Δt*abs(log(Δt)) )
    params = ParametersSIPRK(equation, tableau, Δt, zeros(DT,M), zeros(DT,M), A)

    # create solver
    solver = create_nonlinear_solver(DT, 2*D*S, params)

    # Not implementing InitialGuessSDE
    # create initial guess
    #iguess = InitialGuessODE(get_config(:ig_interpolation), equation, Δt)

    # create cache for internal stage vectors and update vectors
    fcache = NonlinearFunctionCacheSIPRK{DT}(D, M, S)

    # create solution vectors
    q = create_solution_vector(DT, D, NS, NI)
    p = create_solution_vector(DT, D, NS, NI)

    # create integrator
    IntegratorSIPRK{DT, TT, typeof(params), typeof(solver), N}(params, solver, fcache, q, p)
end

equation(integrator::IntegratorSIPRK) = integrator.params.equ
timestep(integrator::IntegratorSIPRK) = integrator.params.Δt
tableau(integrator::IntegratorSIPRK) = integrator.params.tab
Base.eltype(integrator::IntegratorSIPRK{DT, TT, PT, ST, N}) where {DT, TT, PT, ST, N} = DT

function initialize!(int::IntegratorSIPRK{DT,TT}, sol::SolutionPSDE, k::Int, m::Int) where {DT,TT}
    check_solution_dimension_asserts(sol, k, m)

    # copy the m-th initial condition for the k-th sample path
    get_initial_conditions!(sol, int.q[k,m], int.p[k,m], k, m)

    # Not implementing InitialGuessSDE
    # # initialise initial guess
    # initialize!(int.iguess, m, sol.t[0], int.q[m])
end

"""
This function computes initial guesses for Y, Z and assigns them to int.solver.x
The prediction is calculated using an explicit integrator.
"""
function initial_guess!(int::IntegratorSIPRK{DT,TT}) where {DT,TT}

    # NOT IMPLEMENTING InitialGuessSDE

    # SIMPLE SOLUTION
    # The simplest initial guess for Y, Z is 0
    # int.solver.x .= zeros(eltype(int), 2*int.params.tab.s*ndims(int))

    # USING AN EXPLICIT INTEGRATOR TO COMPUTE AN INITIAL GUESS
    # Below we use the R2 method of Burrage & Burrage to calculate
    # the internal stages at the times c[1]...c[s].
    # This approach seems to give very good approximations if the time step
    # and magnitude of noise are not too large. If the noise intensity is too big,
    # one may have to perform a few iterations of the explicit method with a smaller
    # time step, use a higher-order explicit method (e.g. CL or G5), or use
    # the simple solution above.

    local tV1 = zeros(DT,int.params.equ.d)
    local tV2 = zeros(DT,int.params.equ.d)
    local tF1 = zeros(DT,int.params.equ.d)
    local tF2 = zeros(DT,int.params.equ.d)
    local tB1 = zeros(DT,int.params.equ.d, int.params.equ.m)
    local tB2 = zeros(DT,int.params.equ.d, int.params.equ.m)
    local tG1 = zeros(DT,int.params.equ.d, int.params.equ.m)
    local tG2 = zeros(DT,int.params.equ.d, int.params.equ.m)
    local Q   = zeros(DT,int.params.equ.d)
    local P   = zeros(DT,int.params.equ.d)
    local t2::TT
    local Δt_local::TT
    local ΔW_local = zeros(DT,int.params.equ.m)

    # When calling this function, int.params should contain the data:
    # int.params.q - the q solution at the previous time step
    # int.params.p - the p solution at the previous time step
    # int.params.t - the time of the previous step
    # int.params.ΔW- the increment of the Brownian motion for the current step

    #Evaluating the functions v and B at t,q - same for all stages
    int.params.equ.v(int.params.t, int.params.q, int.params.p, tV1)
    int.params.equ.B(int.params.t, int.params.q, int.params.p, tB1)
    int.params.equ.f(int.params.t, int.params.q, int.params.p, tF1)
    int.params.equ.G(int.params.t, int.params.q, int.params.p, tG1)

    # Calculating the positions q at the points qdrift.c[i]
    # if qdrift.c==pdrift.c, then also calculating the momenta
    for i in 1:int.params.tab.s

        # Taking the c[i] from the qdrift tableau.
        Δt_local  = int.params.tab.qdrift.c[i]*int.params.Δt
        ΔW_local .= int.params.tab.qdrift.c[i]*int.params.ΔW

        Q = int.params.q + 2. / 3. * Δt_local * tV1 + 2. / 3. * tB1 * ΔW_local
        P = int.params.p + 2. / 3. * Δt_local * tF1 + 2. / 3. * tG1 * ΔW_local

        t2 = int.params.t + 2. / 3. * Δt_local

        int.params.equ.v(t2, Q, P, tV2)
        int.params.equ.B(t2, Q, P, tB2)
        int.params.equ.f(t2, Q, P, tF2)
        int.params.equ.G(t2, Q, P, tG2)

        #Calculating the Y's and assigning them to the array int.solver.x as initial guesses
        for j in 1:int.params.equ.d
            int.solver.x[(i-1)*int.params.equ.d+j] =  Δt_local*(1. / 4. * tV1[j] + 3. / 4. * tV2[j]) + dot( (1. / 4. * tB1[j,:] + 3. / 4. * tB2[j,:]), ΔW_local )
        end

        # if the collocation points are the same for both q and p parts
        if int.params.tab.qdrift.c==int.params.tab.pdrift.c
            for j in 1:int.params.equ.d
                int.solver.x[(int.params.tab.s+i-1)*int.params.equ.d+j] =  Δt_local*(1. / 4. * tF1[j] + 3. / 4. * tF2[j]) + dot( (1. / 4. * tG1[j,:] + 3. / 4. * tG2[j,:]), ΔW_local )
            end
        end
    end

    # If qdrift.c != pdrift.c, then calculating the momenta p at the points pdrift.c[i]
    if int.params.tab.qdrift.c != int.params.tab.pdrift.c

        for i in 1:int.params.tab.s

            # Taking the c[i] from the pdrift tableau.
            Δt_local  = int.params.tab.pdrift.c[i]*int.params.Δt
            ΔW_local .= int.params.tab.pdrift.c[i]*int.params.ΔW

            Q = int.params.q + 2. / 3. * Δt_local * tV1 + 2. / 3. * tB1 * ΔW_local
            P = int.params.p + 2. / 3. * Δt_local * tF1 + 2. / 3. * tG1 * ΔW_local

            t2 = int.params.t + 2. / 3. * Δt_local

            int.params.equ.v(t2, Q, P, tV2)
            int.params.equ.B(t2, Q, P, tB2)
            int.params.equ.f(t2, Q, P, tF2)
            int.params.equ.G(t2, Q, P, tG2)

            # Calculating the Z's and assigning them to the array int.solver.x as initial guesses
            # The guesses for the Y's have already been written to x above
            for j in 1:int.params.equ.d
                int.solver.x[(int.params.tab.s+i-1)*int.params.equ.d+j] =  Δt_local*(1. / 4. * tF1[j] + 3. / 4. * tF2[j]) + dot( (1. / 4. * tG1[j,:] + 3. / 4. * tG2[j,:]), ΔW_local )
            end
        end
    end

end

"""
Integrate PSDE with a stochastic implicit partitioned Runge-Kutta integrator.
Integrating the k-th sample path for the m-th initial condition
"""
function integrate_step!(int::IntegratorSIPRK{DT,TT}, sol::SolutionPSDE{DT,TT,NQ,NW}, k::Int, m::Int, n::Int) where {DT,TT,NQ,NW}
    check_solution_dimension_asserts(sol, k, m, n)

    # set time for nonlinear solver
    int.params.t  = sol.t[0] + (n-1)*int.params.Δt
    int.params.q .= int.q[k, m]
    int.params.p .= int.p[k, m]

    # copy the increments of the Brownian Process
    if NW==1
        #1D Brownian motion, 1 sample path
        int.params.ΔW[1] = sol.W.ΔW[n-1]
        int.params.ΔZ[1] = sol.W.ΔZ[n-1]
    elseif NW==2
        #Multidimensional Brownian motion, 1 sample path
        int.params.ΔW .= sol.W.ΔW[n-1]
        int.params.ΔZ .= sol.W.ΔZ[n-1]
    elseif NW==3
        #1D or Multidimensional Brownian motion, k-th sample path
        int.params.ΔW .= sol.W.ΔW[n-1,k]
        int.params.ΔZ .= sol.W.ΔZ[n-1,k]
    end

    # truncate the increments ΔW with A
    if int.params.A>0
        for i in eachindex(int.params.ΔW)
            if int.params.ΔW[i]<-int.params.A
                int.params.ΔW[i] = -int.params.A
            elseif int.params.ΔW[i]>int.params.A
                int.params.ΔW[i] = int.params.A
            end
        end
    end

    # compute initial guess and assign to int.solver.x
    initial_guess!(int)

    # call nonlinear solver
    solve!(int.solver)

    # print solver status
    print_solver_status(int.solver.status, int.solver.params)

    # check if solution contains NaNs or error bounds are violated
    check_solver_status(int.solver.status, int.solver.params)

    # compute the drift vector field and the diffusion matrix at internal stages
    compute_stages!(int.solver.x, int.fcache.Q, int.fcache.P, int.fcache.V, int.fcache.F, int.fcache.B, int.fcache.G, int.fcache.Y, int.fcache.Z, int.params)

    # compute final update
    update_solution!(int.q[k,m], int.p[k,m], int.fcache.V, int.fcache.F, int.fcache.B, int.fcache.G,
                    int.params.tab.qdrift.b, int.params.tab.qdrift.b̂, int.params.tab.qdiff.b, int.params.tab.qdiff.b̂,
                    int.params.tab.pdrift.b, int.params.tab.pdrift.b̂, int.params.tab.pdiff.b, int.params.tab.pdiff.b̂,
                    int.params.Δt, int.params.ΔW)

    # # NOT IMPLEMENTING InitialGuessSDE
    # # # copy solution to initial guess
    # # update!(int.iguess, m, sol.t[0] + n*int.params.Δt, int.q[m])

    # take care of periodic solutions
    cut_periodic_solution!(int.q[k,m], int.params.equ.periodicity)

    # # copy to solution
    set_solution!(sol, int.q[k,m], int.p[k,m], n, k, m)
end


text:
using Trilinos
using Test
using MPI

if Test.get_testset_depth() == 0
  MPI.Init()
end

comm = Teuchos.MpiComm(MPI.CComm(MPI.COMM_WORLD))
const n = 20
rowmap = Tpetra.Map(n, 0, comm)
num_my_elements = Int(Tpetra.getNodeNumElements(rowmap))

my_rank = Teuchos.getRank(comm)
if my_rank != 0
  redirect_stdout(open("/dev/null", "w"))
  redirect_stderr(open("/dev/null", "w"))
end

# Number of columns for the multivector
mv_cols = 3

v = Tpetra.Vector(rowmap)
mv = Tpetra.MultiVector(rowmap, mv_cols)

Tpetra.putScalar(v, 1)
Tpetra.putScalar(mv, 2)

@show device_vt = Tpetra.device_view_type(v)
@show host_vt = Tpetra.host_view_type(v)

@show device_mvt = Tpetra.device_view_type(mv)
@show Kokkos.value_type(device_mvt)
@show Kokkos.array_layout(device_mvt)
@show host_mvt = Tpetra.host_view_type(mv)

host_v = Tpetra.getLocalView(host_vt, v)
host_mv = Tpetra.getLocalView(host_mvt, mv)
@test Kokkos.dimension(host_v, 0) == num_my_elements
@test Int(Kokkos.dimension(host_v, 1)) == 1
@test unsafe_load(host_v(0,0)) == 1.0
@test unsafe_load(host_mv(0,2)) == 2.0

# high-level interface, views are 0-based for consistency with the Trilinos API
hv1 = Tpetra.host_view(v)
@test length(LinearIndices(hv1)) == num_my_elements
if my_rank == 0
  hv1[0] = 3
end
if my_rank == (Teuchos.getSize(comm)-1)
  hv1[last(axes(hv1,1))] = 4
end
@show hv1
@test Tpetra.dot(v,v) == (n-2)+9.0+16.0

dv2 = Tpetra.device_view(mv)
@test map(length,axes(dv2)) == (num_my_elements, mv_cols)
if my_rank == 0
  dv2[0,0] = 3
end
if my_rank == (Teuchos.getSize(comm)-1)
  dv2[last(axes(dv2,1)),2] = 4
end

sim_mat = similar(dv2)
@test typeof(sim_mat) == Array{Float64, 2}
@test size(sim_mat) == (n,3)
copy!(sim_mat, dv2)
println("sim_mat:")
println(sim_mat)

dots_check = [(n-1)*4.0+9.0, n*4.0, (n-1)*4.0+16.0]
mv_dots = zeros(mv_cols)
Tpetra.dot(mv,mv,Teuchos.ArrayView(mv_dots))
@test mv_dots == dots_check

if my_rank == 0
  println("hv1:\n",hv1)
  println("dv2:\n",dv2)
end

function benchmark_fill_lowlevel(rowmap, a)
  num_my_elements = Int(Tpetra.getNodeNumElements(rowmap))
  av = Tpetra.getLocalView(Tpetra.device_view_type(a), a)
  for i in 0:(num_my_elements-1)
    unsafe_store!(av(i,0), i)
  end
end

function benchmark_fill_abstractarray(rowmap, a)
  num_my_elements = Int(Tpetra.getNodeNumElements(rowmap))
  av = Tpetra.device_view(a)
  for i in LinearIndices(av)
    av[i] = i
  end
end

function benchmark_fill_abstractarray2(rowmap, a)
  num_my_elements = Int(Tpetra.getNodeNumElements(rowmap))
  av = Tpetra.device_view(a)
  for j in axes(av,2)
    for i in axes(av,1)
      av[i,j] = i*(j+1)
    end
  end
end

const bench_size = 1000000
benchmap = Tpetra.Map(bench_size, 0, comm)
benchvec = Tpetra.Vector(benchmap)
benchview = Tpetra.device_view(benchvec)
n_my_elms = Int(Tpetra.getNodeNumElements(benchmap))
v_end = n_my_elms-1

Tpetra.putScalar(benchvec, 5.0)
println("lowlevel timings:")
@time benchmark_fill_lowlevel(benchmap, benchvec)
@time benchmark_fill_lowlevel(benchmap, benchvec)
@time benchmark_fill_lowlevel(benchmap, benchvec)

@test benchview[0] == 0
@test benchview[v_end] == n_my_elms-1

Tpetra.putScalar(benchvec, 5.0)
println("abstractarray timings:")
@time benchmark_fill_abstractarray(benchmap, benchvec)
@time benchmark_fill_abstractarray(benchmap, benchvec)
@time benchmark_fill_abstractarray(benchmap, benchvec)

@test benchview[0] == 0
@test benchview[v_end] == n_my_elms-1

Tpetra.putScalar(benchvec, 5.0)
println("C++ timings:")
@time Trilinos.Benchmark.vector_fill(benchmap, benchvec)
@time Trilinos.Benchmark.vector_fill(benchmap, benchvec)
@time Trilinos.Benchmark.vector_fill(benchmap, benchvec)

@test benchview[0] == 0
@test benchview[v_end] == n_my_elms-1

benchmv = Tpetra.MultiVector(benchmap, mv_cols)
benchmvview = Tpetra.device_view(benchmv)

Tpetra.putScalar(benchmv, 5.0)
println("abstractarray timings, MultiVector with linear indexing:")
@time benchmark_fill_abstractarray(benchmap, benchmv)
@time benchmark_fill_abstractarray(benchmap, benchmv)
@time benchmark_fill_abstractarray(benchmap, benchmv)

@test benchmvview[0,0] == 1
@test benchmvview[v_end,mv_cols-1] == n_my_elms*mv_cols

Tpetra.putScalar(benchmv, 5.0)
println("abstractarray timings, MultiVector with [i,j] indexing:")
@time benchmark_fill_abstractarray2(benchmap, benchmv)
@time benchmark_fill_abstractarray2(benchmap, benchmv)
@time benchmark_fill_abstractarray2(benchmap, benchmv)

@test benchmvview[0,0] == 0
@test benchmvview[v_end,mv_cols-1] == (n_my_elms-1)*mv_cols

Tpetra.putScalar(benchmv, 5.0)
println("C++ timings, MultiVector:")
@time Trilinos.Benchmark.multivector_fill(benchmap, benchmv)
@time Trilinos.Benchmark.multivector_fill(benchmap, benchmv)
@time Trilinos.Benchmark.multivector_fill(benchmap, benchmv)

@test benchmvview[0,0] == 0
@test benchmvview[v_end,mv_cols-1] == (n_my_elms-1)*mv_cols

if Test.get_testset_depth() == 0
  MPI.Finalize()
end


text:
# ---
# title: 1615. Maximal Network Rank
# id: problem1615
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Graph
# link: <https://leetcode.com/problems/maximal-network-rank/description/>
# hidden: true
# ---
#
# There is an infrastructure of `n` cities with some number of `roads`
# connecting these cities. Each `roads[i] = [ai, bi]` indicates that there is a
# bidirectional road between cities `ai` and `bi`.
#
# The **network rank** __ of **two different cities** is defined as the total
# number of  **directly** connected roads to **either** city. If a road is
# directly connected to both cities, it is only counted **once**.
#
# The **maximal network rank** of the infrastructure is the **maximum network
# rank** of all pairs of different cities.
#
# Given the integer `n` and the array `roads`, return _the **maximal network
# rank** of the entire infrastructure_.
#
#
#
# **Example 1:**
#
# **![](https://assets.leetcode.com/uploads/2020/09/21/ex1.png)**
#
#
#
#     Input: n = 4, roads = [[0,1],[0,3],[1,2],[1,3]]
#     Output: 4
#     Explanation: The network rank of cities 0 and 1 is 4 as there are 4 roads that are connected to either 0 or 1. The road between 0 and 1 is only counted once.
#
#
# **Example 2:**
#
# **![](https://assets.leetcode.com/uploads/2020/09/21/ex2.png)**
#
#
#
#     Input: n = 5, roads = [[0,1],[0,3],[1,2],[1,3],[2,3],[2,4]]
#     Output: 5
#     Explanation: There are 5 roads that are connected to cities 1 or 2.
#
#
# **Example 3:**
#
#
#
#     Input: n = 8, roads = [[0,1],[1,2],[2,3],[2,4],[5,6],[5,7]]
#     Output: 5
#     Explanation: The network rank of 2 and 5 is 5. Notice that all the cities do not have to be connected.
#
#
#
#
# **Constraints:**
#
#   * `2 <= n <= 100`
#   * `0 <= roads.length <= n * (n - 1) / 2`
#   * `roads[i].length == 2`
#   * `0 <= ai, bi <= n-1`
#   * `ai != bi`
#   * Each pair of cities has **at most one** road connecting them.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
function extend!(
    crr::CorefResolver,
    params::ExtendParams,
    )
    source_pair, target_pair = crr |> innerspaces
    rel_pair = crr.crs.relps[params.rel_pair_ind]

    for (source_space, target_space, joiner) in zip(
            [source_pair.erspace1, source_pair.erspace2],
            [target_pair.erspace1, target_pair.erspace2],
            [rel_pair.sj1, rel_pair.sj2]
            )
        extend!(source_space, target_space, joiner, params.root_side)
        target_space.entity_types = source_space.entity_types
        m = target_space.width
        target_space.entity_types[m] = rel_pair.colmap[other(params.root_side)]
    end
end

function extend!(
    source_space::EntityRefSpace,
    target_space::EntityRefSpace,
    joiner::SortedJoiner,
    side::Int,
    )
    opp = other(side)
    source_arr = source_space.space
    target_arr = target_space.space

    n_source = source_space.height
    m_source = source_space.width

    i_target = 0
    m_target = m_source + 1
    # @inbounds
    for i_source in 1:n_source
        for j in searchsorted(joiner.sorted_views[side], source_arr[i_source, m_source])
            i_target += 1
            # inbounds problem - planning needed
            for k in 1:m_source target_arr[i_target, k] = source_arr[i_source, k] end
            target_arr[i_target, m_source + 1] = joiner.arr[joiner.sorted_inds[side][j], opp]
        end
    end
    target_space.width = m_target
    target_space.height = i_target
end


text:
using Test

using PointCloudRasterizers
using LazIO
using GeoArrays
using Statistics

lazfn = joinpath(dirname(pathof(LazIO)), "..", "test/libLAS_1.2.laz")
pointcloud = LazIO.open(lazfn)

cellsizes = (1.,1.)
idx = index(pointcloud, cellsizes)
raster = reduce(idx, field=:Z, reducer=median)

# indexing and reducing tests
@test typeof(idx) == PointCloudRasterizers.PointCloudIndex
@test count(ismissing,raster) == 24503457
@test count(!ismissing,raster) == 496543
@test isapprox(mean(skipmissing(raster.A)), 861.5422515270361)
@test maximum(idx.counts) == 2

# filtering tests
last_return(p) = LazIO.return_number(p) == LazIO.number_of_returns(p)
filter!(idx, last_return)
@test sum(idx.counts) == 497534

# file IO tests
GeoArrays.write!("last_return_median.tif", raster)
@test isfile("last_return_median.tif")
rm("last_return_median.tif")


text:
using Zygote
using Zygote: @adjoint, @nograd

@adjoint function free_energy(K::Matrix{T}, model::AbstractSampler, samples) where T <: Real
    free_energy(K, model, samples), function (adjy)
        adjmodel = grad_model(K, model, samples) .* adjy
        return (nothing, adjmodel, nothing)
    end
end

@adjoint function PSAModel(nbit::Int)
    PSAmodel(nbit) , _ -> nothing
end

@adjoint function AutoRegressiveModel(nbit::Int, nhiddens)
    AutoRegressivemodel(nbit) , _ -> nothing
end

@nograd gen_samples, createmasks


text:
export Linear

const WeightType{T} = AbstractMatrix{T}
const BiasType{T} = Union{Compat.Nothing, AbstractVector{T}}

"""
    Linear{T, WeightType, BiasType, InputType} <: AbstractBlock

Linear layer.

``z^{l+1} = W^l * z^{l} + b^l``

expression of gradients:

``\\delta^{l} = (W^l)^T\\delta^{l+1}``
``\\delta^{l}_w = \\delta^{l+1}(z^l)^T``
``\\delta^{l}_b = \\delta^{l+1}``
"""
mutable struct Linear{
                        T, # Numeric Type
                        WType<:WeightType{T},
                        BType<:BiasType{T},
                        IType<:AbstractVecOrMat{T}
                    } <: AbstractBlock

    weight::WType
    bias::BType

    grad_weight::WType
    grad_bias::BType

    input::IType

    function Linear{T, W, B, I}(
        weight::W,
        bias::B,
        grad_weight::W,
        grad_bias::B
        ) where {T, W, B, I}
        new{T, W, B, I}(weight, bias, grad_weight, grad_bias)
    end
end

include("Utils.jl")
include("Constructor.jl")
include("Forward.jl")
include("Backward.jl")

import Base: show

function show(io::IO, op::Linear)
    dname = string(device(op.weight))
    print(io, "Linear{$(dname)} ($(size(op, 1)) -> $(size(op, 2)))")
end


text:
module ITensorMPOCompression

using ITensors
using ITensors.NDTensors

import ITensors.BlockSparseTensor,ITensors.DenseTensor,ITensors.tensor
# For some reason we import these in order for the mere existence of
# the ITensorMPOCompression versions to be acknowledged!?!?
import ITensors.orthogonalize!,ITensors.truncate!

export block_qx,ql,lq,rq,set_scale! #qx related
export slice,assign!,redim #operator handling
export getV,setV,growRL,V_offsets #blocking related
# lots of characterization functions
export reg_form,orth_type,matrix_state,upper,lower,left,right,mirror
export parse_links,is_regular_form,getM,grow,detect_regular_form
export is_lower_regular_form,is_upper_regular_form,is_canonical,is_orthogonal
export detect_upper_lower,is_upper_lower
# Hamiltonian related
export make_transIsing_MPO,make_Heisenberg_AutoMPO,make_transIsing_AutoMPO,to_openbc,fast_GS,fix_autoMPO!
# MPO and bond spectrum
export get_Dw,min,max
export bond_spectrum,bond_spectrums

export add_or_replace # Handle inpenetrable kwargs

export orthogonalize!,truncate,truncate! #the punchline

default_eps=1e-14 #for characterization routines, floats abs()<default_eps are considered to be zero.

"""
    @enum reg_form  upper lower

Indicates that an MPO or operator-valued matrix has either an `upper` or `lower` regular form.
This becomes non-trival for rectangular matrices.
See also [`detect_regular_form`](@ref) and related functions
"""
@enum reg_form  upper lower

"""
    @enum orth_type left right

Indicates that an MPO matrix satisfies the conditions for `left` or `right` canonical form
"""
@enum orth_type left right

# """
#     mirror(lr::orth_type)::orth_type
#     returns this mirror of lr.  `left`->`right` and `right`->`left`
# """
function mirror(lr::orth_type)::orth_type
    if lr==left
        ret=right
    else #must be right
        ret=left
    end
    return ret
end

 """
     Indicates both the `orth_type` and `reg_form` of an MPO
     See also [`orth_type`](@ref) [`reg_form`](@ref)
 """
struct matrix_state
    ul::reg_form
    lr::orth_type
end

mirror(ms::matrix_state)::matrix_state=matrix_state(ms.ul,mirror(ms.lr))

"""
    V_offsets

A simple struct for encapsulating offsets for V-blocks

"""
struct V_offsets
    o1::Int64 #currently o1=o2 for std. Parker compression.
    o2::Int64 #Leave them distinct for now until we know more
    #
    # The purpose of this struct is to ensure the asserts below
    #
    V_offsets(o1_::Int64, o2_::Int64) = begin
        @assert o1_==0 || o1_==1
        @assert o2_==0 || o2_==1
        new(o1_,o2_)
    end
end

"""
    V_offsets(ms::matrix_state)

Derives the correct V-block offsets for the given `matrix_state`
"""
V_offsets(ms::matrix_state) = begin
    if ms.lr==left
        if ms.ul ==lower
            o1_=1
            o2_=1
        else #upper
            o1_=0
            o2_=0
        end
    else #right
        if ms.ul ==lower
            o1_=0
            o2_=0
        else #upper
            o1_=1
            o2_=1
        end
    end
    V_offsets(o1_,o2_)
end

struct bond_spectrum
    spectrum::Vector{Float64}
    link_number::Int64
    bond_spectrum(s::ITensor,link::Int64) = begin
        @assert link>0
        @assert order(s)==2
        new(diag(array(s)),link)
    end
end

bond_spectrums = Vector{bond_spectrum}

function max(s::bond_spectrum)::Float64 s.spectrum[1] end
function min(s::bond_spectrum)::Float64 s.spectrum[end] end
function max(ss::bond_spectrums)::Float64
    ret=max(ss[1])
    for n in 2:length(ss)
        ms=max(ss[n])
        if ms>ret ret=ms end
    end
    return ret
end

function min(ss::bond_spectrums)::Float64
    ret=min(ss[1])
    for n in 2:length(ss)
        ms=min(ss[n])
        if ms<ret ret=ms end
    end
    return ret
end

function slice(A::ITensor,iv::IndexVal...)::ITensor
    iv_dagger=[dag(x.first)=>x.second for x in iv]
    return A*onehot(iv_dagger...)
end

"""
    assign!(W::ITensor,i1::IndexVal,i2::IndexVal,op::ITensor)

Assign an operator to an element of the operator valued matrix W
    W[i1,i2]=op
"""
function assign!(W::ITensor,op::ITensor,ivs::IndexVal...)
    assign!(W,tensor(op),ivs...)
end

function assign!(W::ITensor,op::DenseTensor{ElT,N},ivs::IndexVal...) where {ElT,N}
    iss=inds(op)
    for s in eachindval(iss)
        s2=[x.second for x in s]
        W[ivs...,s...]=op[s2...]
    end
end

function assign!(W::ITensor,op::BlockSparseTensor{ElT,N},ivs::IndexVal...) where {ElT,N}
    iss=inds(op)
    for b in eachnzblock(op)
        isv=[iss[i]=>b[i] for i in 1:length(b)]
        W[ivs...,isv...]=op[b][1] #not sure why we need [1] here
    end
end

"""
    function redim(i::Index,Dw::Int64)::Index

    Create an index with the same tags ans plev, but different dimension(s) and and id
"""
function redim(i::Index,Dw::Int64...)::Index
    @assert length(Dw)==nblocks(i)
    if hasqns(i)
        j=0
        new_qns=[(j+=1;q.first=>Dw[j]) for q in space(i)]
        return Index(new_qns...;dir=dir(i),tags=tags(i),plev=plev(i))
    else
        return Index(Dw[1];tags=tags(i),plev=plev(i))
    end
end

"""
    set_scale!(RL::ITensor,Q::ITensor,off::V_offsets)

Fix the gauge freedom between QR/QL after a block respecting QR/QL decomposition. The
gauge fix is to ensure that either the top left or bottom right corner of `R` is 1.0.
"""
function set_scale!(RL::ITensor,Q::ITensor,off::V_offsets)
    @assert order(RL)==2
    is=inds(RL)
    Dw1,Dw2=map(dim,is)
    i1= off.o1==0 ? 1 : Dw1
    i2= off.o2==0 ? 1 : Dw2
    scale=RL[is[1]=>i1,is[2]=>i2]
    @assert abs(scale)>1e-12
    RL./=scale
    Q.*=scale
end

include("util.jl")
include("hamiltonians.jl")
include("qx.jl")
include("characterization.jl")
include("blocking.jl")
include("orthogonalize.jl")
include("truncate.jl")

end


text:
# Approximations
include("./approximations/gradient_descent.jl")
include("./approximations/SR1.jl")
include("./approximations/DFP.jl")
include("./approximations/BFGS.jl")
include("./approximations/DBFGS.jl")
include("./approximations/newton.jl")
# include("LBFGS.jl")

include("update_qn.jl")
include("find_direction.jl")