text:
"""
    AbstractEstimationMethod

Abstract type for defining statistical estimation methods.
"""
abstract type AbstractEstimationMethod end

"""
    GradientDescentEstimation <: AbstractEstimationMethod

Method for estimation using gradient descent.
"""
struct GradientDescentEstimation <: AbstractEstimationMethod end

"""
    CyclicProximalPointEstimation <: AbstractEstimationMethod

Method for estimation using the cyclic proximal point technique.
"""
struct CyclicProximalPointEstimation <: AbstractEstimationMethod end

"""
    ExtrinsicEstimation <: AbstractEstimationMethod

Method for estimation in the ambient space and projecting to the manifold.

For [`mean`](@ref) estimation, [`GeodesicInterpolation`](@ref) is used for mean estimation
in the ambient space.
"""
struct ExtrinsicEstimation <: AbstractEstimationMethod end

"""
    WeiszfeldEstimation <: AbstractEstimationMethod

Method for estimation using the Weiszfeld algorithm for the [`median`](@ref)
"""
struct WeiszfeldEstimation <: AbstractEstimationMethod end

_unit_weights(n::Int) = StatsBase.UnitWeights{Float64}(n)

@doc raw"""
    GeodesicInterpolation <: AbstractEstimationMethod

Repeated weighted geodesic interpolation method for estimating the Riemannian
center of mass.

The algorithm proceeds with the following simple online update:

```math
\begin{aligned}
μ_1 &= x_1\\
t_k &= \frac{w_k}{\sum_{i=1}^k w_i}\\
μ_{k} &= γ_{μ_{k-1}}(x_k; t_k),
\end{aligned}
```

where $x_k$ are points, $w_k$ are weights, $μ_k$ is the $k$th estimate of the
mean, and $γ_x(y; t)$ is the point at time $t$ along the
[`shortest_geodesic`](@ref shortest_geodesic(::AbstractManifold, ::Any, ::Any, ::Real))
between points $x,y ∈ \mathcal M$. The algorithm
terminates when all $x_k$ have been considered. In the [`Euclidean`](@ref) case,
this exactly computes the weighted mean.

The algorithm has been shown to converge asymptotically with the sample size for
the following manifolds equipped with their default metrics when all sampled
points are in an open geodesic ball about the mean with corresponding radius
(see [`GeodesicInterpolationWithinRadius`](@ref)):

* All simply connected complete Riemannian manifolds with non-positive sectional
  curvature at radius $∞$ [^Cheng2016], in particular:
    + [`Euclidean`](@ref)
    + [`SymmetricPositiveDefinite`](@ref) [^Ho2013]
* Other manifolds:
    + [`Sphere`](@ref): $\frac{π}{2}$ [^Salehian2015]
    + [`Grassmann`](@ref): $\frac{π}{4}$ [^Chakraborty2015]
    + [`Stiefel`](@ref)/[`Rotations`](@ref): $\frac{π}{2 \sqrt 2}$ [^Chakraborty2019]

For online variance computation, the algorithm additionally uses an analogous
recursion to the weighted Welford algorithm [^West1979].

[^Ho2013]:
    > Ho J.; Cheng G.; Salehian H.; Vemuri B. C.; Recursive Karcher expectation
    > estimators and geometric law of large numbers.
    > Proceedings of the 16th International Conference on Artificial Intelligence
    > and Statistics (2013), pp. 325–332.
    > [pdf](http://proceedings.mlr.press/v31/ho13a.pdf).
[^Salehian2015]:
    > Salehian H.; Chakraborty R.; Ofori E.; Vaillancourt D.; An efficient
    > recursive estimator of the Fréchet mean on a hypersphere with applications
    > to Medical Image Analysis.
    > Mathematical Foundations of Computational Anatomy (2015).
    > [pdf](https://www-sop.inria.fr/asclepios/events/MFCA15/Papers/MFCA15_4_2.pdf).
[^Chakraborty2015]:
    > Chakraborty R.; Vemuri B. C.; Recursive Fréchet Mean Computation on the
    > Grassmannian and Its Applications to Computer Vision.
    > Proceedings of the IEEE International Conference on Computer Vision (ICCV) (2015),
    > pp. 4229-4237.
    > doi: [10.1109/ICCV.2015.481](https://doi.org/10.1109/ICCV.2015.481),
    > [link](http://openaccess.thecvf.com/content_iccv_2015/html/Chakraborty_Recursive_Frechet_Mean_ICCV_2015_paper.html).
[^Cheng2016]:
    > Cheng G.; Ho J.; Salehian H.; Vemuri B. C.; Recursive Computation of the
    > Fréchet Mean on Non-positively Curved Riemannian Manifolds with Applications.
    > Riemannian Computing in Computer Vision. Springer, Cham (2016), pp. 21-43.
    > doi: [10.1007/978-3-319-22957-7_2](https://doi.org/10.1007/978-3-319-22957-7_2),
    > [pdf](https://www.cise.ufl.edu/~vemuri/paperphp/article.php?y=2016&i=5).
[^Chakraborty2019]:
    > Chakraborty R.; Vemuri B. C.; Statistics on the (compact) Stiefel manifold:
    > Theory and Applications.
    > The Annals of Statistics (2019), 47(1), pp. 415-438.
    > doi: [10.1214/18-AOS1692](https://doi.org/10.1214/18-AOS1692),
    > arxiv: [1708.00045](https://arxiv.org/abs/1708.00045).
[^West1979]:
    > West D. H. D.; Updating Mean and Variance Estimates: An Improved Method.
    > Communications of the ACM (1979), 22(9), pp. 532–535.
    > doi: [10.1145/359146.359153](https://doi.org/10.1145/359146.359153).
"""
struct GeodesicInterpolation <: AbstractEstimationMethod end

"""
    GeodesicInterpolationWithinRadius{T} <: AbstractEstimationMethod

Estimation of Riemannian center of mass using [`GeodesicInterpolation`](@ref)
with fallback to [`GradientDescentEstimation`](@ref) if any points are outside of a
geodesic ball of specified `radius` around the mean.

# Constructor

    GeodesicInterpolationWithinRadius(radius)
"""
struct GeodesicInterpolationWithinRadius{T} <: AbstractEstimationMethod
    radius::T

    function GeodesicInterpolationWithinRadius(radius::T) where {T}
        radius > 0 && return new{T}(radius)
        return throw(
            DomainError("The radius must be strictly postive, received $(radius)."),
        )
    end
end

function Base.show(io::IO, method::GeodesicInterpolationWithinRadius)
    return print(io, "GeodesicInterpolationWithinRadius($(method.radius))")
end

"""
    Statistics.cov(
        M::AbstractManifold,
        x::AbstractVector;
        basis::AbstractBasis=DefaultOrthonormalBasis(),
        tangent_space_covariance_estimator::CovarianceEstimator=SimpleCovariance(;
            corrected=true,
        ),
        mean_estimation_method::AbstractEstimationMethod=GradientDescentEstimation(),
        inverse_retraction_method::AbstractInverseRetractionMethod=default_inverse_retraction_method(
            M,
        ),
    )

Estimate the covariance matrix of a set of points `x` on manifold `M`. Since the covariance matrix
on a manifold is a rank 2 tensor, the function returns its coefficients in basis induced by
the given tangent space basis. See Section 5 of [^Pennec2006] for details.

The mean is calculated using the specified `mean_estimation_method` using
[mean](@ref Statistics.mean(::AbstractManifold, ::AbstractVector, ::AbstractEstimationMethod),
and tangent vectors at this mean are calculated using the provided `inverse_retraction_method`.
Finally, the covariance matrix in the tangent plane is estimated using the Euclidean space
 estimator `tangent_space_covariance_estimator`. The type `CovarianceEstimator` is defined
 in [`StatsBase.jl`](https://juliastats.org/StatsBase.jl/stable/cov/#StatsBase.CovarianceEstimator)
 and examples of covariance estimation methods can be found in
 [`CovarianceEstimation.jl`](https://github.com/mateuszbaran/CovarianceEstimation.jl/).

[^Pennec2006]:
    > X. Pennec, “Intrinsic Statistics on Riemannian Manifolds: Basic Tools for Geometric
    > Measurements,” J Math Imaging Vis, vol. 25, no. 1, p. 127, Jul. 2006,
    > doi: 10.1007/s10851-006-6228-4.
"""
function Statistics.cov(
    M::AbstractManifold,
    x::AbstractVector;
    basis::AbstractBasis=DefaultOrthonormalBasis(),
    tangent_space_covariance_estimator::CovarianceEstimator=SimpleCovariance(;
        corrected=true,
    ),
    mean_estimation_method::AbstractEstimationMethod=GradientDescentEstimation(),
    inverse_retraction_method::AbstractInverseRetractionMethod=default_inverse_retraction_method(
        M,
    ),
)
    μ = mean(M, x, mean_estimation_method)
    Xs_at_μ = map(p -> inverse_retract(M, μ, p, inverse_retraction_method), x)
    Xcs = map(X -> get_coordinates(M, μ, X, basis), Xs_at_μ)
    return cov(
        tangent_space_covariance_estimator,
        reduce(hcat, Xcs);
        mean=zero(Xcs[1]),
        dims=2,
    )
end

@doc raw"""
    mean(M::AbstractManifold, x::AbstractVector[, w::AbstractWeights]; kwargs...)

Compute the (optionally weighted) Riemannian center of mass also known as
Karcher mean of the vector `x` of points on the [`AbstractManifold`](@ref) `M`, defined
as the point that satisfies the minimizer
````math
\argmin_{y ∈ \mathcal M} \frac{1}{2 \sum_{i=1}^n w_i} \sum_{i=1}^n w_i\mathrm{d}_{\mathcal M}^2(y,x_i),
````
where $\mathrm{d}_{\mathcal M}$ denotes the Riemannian [`distance`](@ref).

In the general case, the [`GradientDescentEstimation`](@ref) is used to compute the mean.
    mean(
        M::AbstractManifold,
        x::AbstractVector,
        [w::AbstractWeights,]
        method::AbstractEstimationMethod;
        kwargs...,
    )

Compute the mean using the specified `method`.

    mean(
        M::AbstractManifold,
        x::AbstractVector,
        [w::AbstractWeights,]
        method::GradientDescentEstimation;
        p0=x[1],
        stop_iter=100,
        retraction::AbstractRetractionMethod = ExponentialRetraction(),
        inverse_retraction::AbstractInverseRetractionMethod = LogarithmicInverseRetraction(),
        kwargs...,
    )

Compute the mean using the gradient descent scheme [`GradientDescentEstimation`](@ref).

Optionally, provide `p0`, the starting point (by default set to the first data
point). `stop_iter` denotes the maximal number of iterations to perform and the
`kwargs...` are passed to [`isapprox`](@ref) to stop, when the minimal change
between two iterates is small. For more stopping criteria check the
[`Manopt.jl`](https://manoptjl.org) package and use a solver therefrom.

Optionally, pass `retraction` and `inverse_retraction` method types to specify
the (inverse) retraction.

The Theory stems from[^Karcher1977] and is also described in[^PennecArsigny2013]
as the exponential barycenter.
The algorithm is further described in[^Afsari2013].

[^Afsari2013]:
    > Afsari, B; Tron, R.; Vidal, R.: On the Convergence of Gradient
    > Descent for Finding the Riemannian Center of Mass,
    > SIAM Journal on Control and Optimization (2013), 51(3), pp. 2230–2260,
    > doi: [10.1137/12086282X](https://doi.org/10.1137/12086282X),
    > arxiv: [1201.0925](https://arxiv.org/abs/1201.0925)
[^PennecArsigny2013]:
    > Pennec X., Arsigny V.: Exponential Barycenters of the Canonical Cartan Connection and
    > Invariant Means on Lie Groups.
    > In: Nielsen F., Bhatia R. (eds) Matrix Information Geometry, (2013), pp. 123-166.
    > doi: [10.1007/978-3-642-30232-9_7](https://doi.org/10.1007/978-3-642-30232-9_7),
    > hal: [https://hal.inria.fr/hal-00699361/document](https://hal.inria.fr/hal-00699361/document)
[^Karcher1977]:
    > Karcher, H.: Riemannian center of mass and mollifier smoothing.
    > Communications on Pure Applied Mathematics (1977), 30, pp. 509–541.
    > doi [10.1002/cpa.3160300502](https://doi.org/10.1002/cpa.3160300502)
"""
mean(::AbstractManifold, ::Any...)
function Statistics.mean(
    M::AbstractManifold,
    x::AbstractVector,
    method::AbstractEstimationMethod...;
    kwargs...,
)
    y = allocate_result(M, mean, x[1])
    return mean!(M, y, x, method...; kwargs...)
end
function Statistics.mean(
    M::AbstractManifold,
    x::AbstractVector,
    w::AbstractVector,
    method::AbstractEstimationMethod...;
    kwargs...,
)
    y = allocate_result(M, mean, x[1])
    return mean!(M, y, x, w, method...; kwargs...)
end

@doc raw"""
    mean!(M::AbstractManifold, y, x::AbstractVector[, w::AbstractWeights]; kwargs...)
    mean!(
        M::AbstractManifold,
        y,
        x::AbstractVector,
        [w::AbstractWeights,]
        method::AbstractEstimationMethod;
        kwargs...,
    )

Compute the [`mean`](@ref mean(::AbstractManifold, args...)) in-place in `y`.
"""
mean!(::AbstractManifold, ::Any...)
function Statistics.mean!(
    M::AbstractManifold,
    y,
    x::AbstractVector,
    method::AbstractEstimationMethod...;
    kwargs...,
)
    w = _unit_weights(length(x))
    return mean!(M, y, x, w, method...; kwargs...)
end
function Statistics.mean!(
    M::AbstractManifold,
    y,
    x::AbstractVector,
    w::AbstractVector;
    kwargs...,
)
    return mean!(M, y, x, w, GradientDescentEstimation(); kwargs...)
end
function Statistics.mean!(
    M::AbstractManifold,
    y,
    x::AbstractVector,
    w::AbstractVector,
    ::GradientDescentEstimation;
    p0=x[1],
    stop_iter=100,
    retraction::AbstractRetractionMethod=ExponentialRetraction(),
    inverse_retraction::AbstractInverseRetractionMethod=LogarithmicInverseRetraction(),
    kwargs...,
)
    n = length(x)
    if length(w) != n
        throw(
            DimensionMismatch(
                "The number of weights ($(length(w))) does not match the number of points for the mean ($(n)).",
            ),
        )
    end
    copyto!(y, p0)
    yold = allocate_result(M, mean, y)
    v = zero_vector(M, y)
    vtmp = copy(v)
    α = w ./ cumsum(w)
    for i in 1:stop_iter
        copyto!(yold, y)
        # Online weighted mean
        @inbounds inverse_retract!(M, v, yold, x[1], inverse_retraction)
        @inbounds for j in 2:n
            inverse_retract!(M, vtmp, yold, x[j], inverse_retraction)
            v .+= α[j] .* (vtmp .- v)
        end
        retract!(M, y, yold, v, 0.5, retraction)
        isapprox(M, y, yold; kwargs...) && break
    end
    return y
end

"""
    mean(
        M::AbstractManifold,
        x::AbstractVector,
        [w::AbstractWeights,]
        method::GeodesicInterpolation;
        shuffle_rng=nothing,
        retraction::AbstractRetractionMethod = ExponentialRetraction(),
        inverse_retraction::AbstractInverseRetractionMethod = LogarithmicInverseRetraction(),
        kwargs...,
    )

Estimate the Riemannian center of mass of `x` in an online fashion using
repeated weighted geodesic interpolation. See
[`GeodesicInterpolation`](@ref) for details.

If `shuffle_rng` is provided, it is used to shuffle the order in which the
points are considered for computing the mean.

Optionally, pass `retraction` and `inverse_retraction` method types to specify
the (inverse) retraction.
"""
mean(::AbstractManifold, ::AbstractVector, ::AbstractVector, ::GeodesicInterpolation)

function Statistics.mean!(
    M::AbstractManifold,
    q,
    x::AbstractVector,
    w::AbstractVector,
    ::GeodesicInterpolation;
    shuffle_rng::Union{AbstractRNG,Nothing}=nothing,
    retraction::AbstractRetractionMethod=ExponentialRetraction(),
    inverse_retraction::AbstractInverseRetractionMethod=LogarithmicInverseRetraction(),
    kwargs...,
)
    n = length(x)
    if length(w) != n
        throw(
            DimensionMismatch(
                "The number of weights ($(length(w))) does not match the number of points for the mean ($(n)).",
            ),
        )
    end
    order = shuffle_rng === nothing ? (1:n) : shuffle(shuffle_rng, 1:n)
    @inbounds begin
        j = order[1]
        s = w[j]
        copyto!(q, x[j])
    end
    v = zero_vector(M, q)
    ytmp = allocate_result(M, mean, q)
    @inbounds for i in 2:n
        j = order[i]
        s += w[j]
        t = w[j] / s
        inverse_retract!(M, v, q, x[j], inverse_retraction)
        retract!(M, ytmp, q, v, t, retraction)
        copyto!(q, ytmp)
    end
    return q
end

"""
    mean(
        M::AbstractManifold,
        x::AbstractVector,
        [w::AbstractWeights,]
        method::GeodesicInterpolationWithinRadius;
        kwargs...,
    )

Estimate the Riemannian center of mass of `x` using
[`GeodesicInterpolationWithinRadius`](@ref).

See [`mean`](@ref mean(::AbstractManifold, ::AbstractVector, ::AbstractVector, ::GeodesicInterpolation))
for a description of `kwargs`.
"""
Statistics.mean(
    ::AbstractManifold,
    ::AbstractVector,
    ::AbstractVector,
    ::GeodesicInterpolationWithinRadius,
)

function Statistics.mean!(
    M::AbstractManifold,
    q,
    x::AbstractVector,
    w::AbstractVector,
    method::GeodesicInterpolationWithinRadius;
    shuffle_rng=nothing,
    kwargs...,
)
    mean!(M, q, x, w, GeodesicInterpolation(); shuffle_rng=shuffle_rng, kwargs...)
    radius = method.radius
    injectivity_radius(M, q) ≤ radius && return q
    for i in eachindex(x)
        @inbounds if distance(M, q, x[i]) ≥ radius
            return mean!(M, q, x, w, GradientDescentEstimation(); p0=q, kwargs...)
        end
    end
    return q
end
function Statistics.mean!(
    M::AbstractManifold,
    q,
    x::AbstractVector,
    w::AbstractVector,
    ::CyclicProximalPointEstimation;
    p0=x[1],
    stop_iter=1000000,
    retraction::AbstractRetractionMethod=ExponentialRetraction(),
    inverse_retraction::AbstractInverseRetractionMethod=LogarithmicInverseRetraction(),
    kwargs...,
)
    n = length(x)
    if length(w) != n
        throw(
            DimensionMismatch(
                "The number of weights ($(length(w))) does not match the number of points for the mean ($(n)).",
            ),
        )
    end
    copyto!(q, p0)
    yold = allocate_result(M, mean, q)
    ytmp = copy(yold)
    X = zero_vector(M, q)
    wv = convert(AbstractVector, w) ./ sum(w)
    for i in 1:stop_iter
        λ = 0.5 / i
        copyto!(yold, q)
        for j in 1:n
            @inbounds t = (2 * λ * wv[j]) / (1 + 2 * λ * wv[j])
            @inbounds inverse_retract!(M, X, q, x[j], inverse_retraction)
            retract!(M, ytmp, q, X, t, retraction)
            copyto!(q, ytmp)
        end
        isapprox(M, q, yold; kwargs...) && break
    end
    return q
end

@decorator_transparent_signature Statistics.mean!(
    M::AbstractDecoratorManifold,
    y,
    x::AbstractVector,
    w::AbstractVector;
    kwargs...,
)

function decorator_transparent_dispatch(
    ::typeof(mean!),
    ::AbstractEmbeddedManifold,
    args...;
    kwargs...,
)
    return Val(:parent)
end
function decorator_transparent_dispatch(
    ::typeof(mean!),
    ::AbstractEmbeddedManifold{𝔽,<:TransparentIsometricEmbedding},
    args...;
    kwargs...,
) where {𝔽}
    return Val(:transparent)
end

"""
    mean(
        M::AbstractManifold,
        x::AbstractVector,
        [w::AbstractWeights,]
        method::ExtrinsicEstimation;
        kwargs...,
    )

Estimate the Riemannian center of mass of `x` using
[`ExtrinsicEstimation`](@ref), i.e. by computing the mean in the embedding and projecting
the result back.
You can specify an `extrinsic_method` to specify which mean estimation method to use in the embedding,
which defaults to [`GeodesicInterpolation`](@ref).

See [`mean`](@ref mean(::AbstractManifold, ::AbstractVector, ::AbstractVector, ::GeodesicInterpolation))
for a description of the remaining `kwargs`.
"""
Statistics.mean(
    ::AbstractManifold,
    ::AbstractVector,
    ::AbstractVector,
    ::ExtrinsicEstimation,
)

function Statistics.mean!(
    M::AbstractManifold,
    y,
    x::AbstractVector,
    w::AbstractVector,
    ::ExtrinsicEstimation;
    extrinsic_method::AbstractEstimationMethod=GeodesicInterpolation(),
    kwargs...,
)
    embedded_x = map(p -> embed(M, p), x)
    embedded_y = mean(get_embedding(M), embedded_x, w, extrinsic_method; kwargs...)
    project!(M, y, embedded_y)
    return y
end

@doc raw"""
    median(M::AbstractManifold, x::AbstractVector[, w::AbstractWeights]; kwargs...)
    median(
        M::AbstractManifold,
        x::AbstractVector,
        [w::AbstractWeights,]
        method::AbstractEstimationMethod;
        kwargs...,
    )

Compute the (optionally weighted) Riemannian median of the vector `x` of points on the
[`AbstractManifold`](@ref) `M`, defined as the point that satisfies the minimizer
````math
\argmin_{y ∈ \mathcal M} \frac{1}{\sum_{i=1}^n w_i} \sum_{i=1}^n w_i\mathrm{d}_{\mathcal M}(y,x_i),
````
where $\mathrm{d}_{\mathcal M}$ denotes the Riemannian [`distance`](@ref).
This function is nonsmooth (i.e nondifferentiable).

In the general case, the [`CyclicProximalPointEstimation`](@ref) is used to compute the
median. However, this default may be overloaded for specific manifolds.

Compute the median using the specified `method`.
"""
Statistics.median(::AbstractManifold, ::Any...)

"""
    median(
        M::AbstractManifold,
        x::AbstractVector,
        [w::AbstractWeights,]
        method::CyclicProximalPointEstimation;
        p0=x[1],
        stop_iter=1000000,
        retraction::AbstractRetractionMethod = ExponentialRetraction(),
        inverse_retraction::AbstractInverseRetractionMethod = LogarithmicInverseRetraction(),
        kwargs...,
    )

Compute the median using [`CyclicProximalPointEstimation`](@ref).

Optionally, provide `p0`, the starting point (by default set to the first
data point). `stop_iter` denotes the maximal number of iterations to perform
and the `kwargs...` are passed to [`isapprox`](@ref) to stop, when the minimal
change between two iterates is small. For more stopping criteria check the
[`Manopt.jl`](https://manoptjl.org) package and use a solver therefrom.

Optionally, pass `retraction` and `inverse_retraction` method types to specify
the (inverse) retraction.

The algorithm is further described in [^Bačák2014].

[^Bačák2014]:
    > Bačák, M: Computing Medians and Means in Hadamard Spaces.
    > SIAM Journal on Optimization (2014), 24(3), pp. 1542–1566,
    > doi: [10.1137/140953393](https://doi.org/10.1137/140953393),
    > arxiv: [1210.2145](https://arxiv.org/abs/1210.2145)
"""
Statistics.median(
    ::AbstractManifold,
    ::AbstractVector,
    ::AbstractVector,
    ::CyclicProximalPointEstimation,
)

"""
    median(
        M::AbstractManifold,
        x::AbstractVector,
        [w::AbstractWeights,]
        method::ExtrinsicEstimation;
        extrinsic_method = CyclicProximalPointEstimation(),
        kwargs...,
    )

Estimate the median of `x` using [`ExtrinsicEstimation`](@ref), i.e. by computing the median
in the embedding and projecting the result back.
You can specify an `extrinsic_method` to specify which median estimation method to use in
the embedding, which defaults to [`CyclicProximalPointEstimation`](@ref).

See [`median`](@ref median(::AbstractManifold, ::AbstractVector, ::AbstractVector, ::CyclicProximalPointEstimation))
for a description of `kwargs`.
"""
Statistics.median(
    ::AbstractManifold,
    ::AbstractVector,
    ::AbstractVector,
    ::ExtrinsicEstimation,
)

@doc raw"""
    median(
        M::AbstractManifold,
        x::AbstractVector,
        [w::AbstractWeights,]
        method::WeiszfeldEstimation;
        α = 1.0,
        p0=x[1],
        stop_iter=2000,
        retraction::AbstractRetractionMethod = ExponentialRetraction(),
        inverse_retraction::AbstractInverseRetractionMethod = LogarithmicInverseRetraction(),
        kwargs...,
    )

Compute the median using [`WeiszfeldEstimation`](@ref).

Optionally, provide `p0`, the starting point (by default set to the first
data point). `stop_iter` denotes the maximal number of iterations to perform
and the `kwargs...` are passed to [`isapprox`](@ref) to stop, when the minimal
change between two iterates is small. For more stopping criteria check the
[`Manopt.jl`](https://manoptjl.org) package and use a solver therefrom.

The parameter ``α\in (0,2]`` is a step size.

The algorithm is further described in [^FletcherVenkatasubramanianJoshi2008],
especially the update rule in Eq. (6), i.e. Let ``q_{k}`` denote the current
iterate, $n$ the number of points ``x_1,\ldots,x_n``, and

```math
I_k = \bigl\{ i \in \{1,\ldots,n\} \big| x_i \neq q_k \bigr\}
```

all indices of points that are not equal to the current iterate. Then the update reads
``q_{k+1} = \exp_{q_k}(αX)``, where

```math
X = \frac{1}{s}\sum_{i\in I_k} \frac{w_i}{d_{\mathcal M}(q_k,x_i)}\log_{q_k}x_i
\quad
\text{ with }
\quad
s = \sum_{i\in I_k} \frac{w_i}{d_{\mathcal M}(q_k,x_i)},
```

and where $\mathrm{d}_{\mathcal M}$ denotes the Riemannian [`distance`](@ref).

Optionally, pass `retraction` and `inverse_retraction` method types to specify
the (inverse) retraction, which by default use the exponential and logarithmic map,
respectively.

[^FletcherVenkatasubramanianJoshi2008]:
    > Fletcher, T., Venkatasubramanian, S., Joshi, S:
    > Robust statistics on Riemannian manifolds via the geometric median
    > 2008 IEEE Conference on Computer Vision and Pattern Recognition,
    > doi: [10.1109/CVPR.2008.4587747](https://doi.org/10.1109/CVPR.2008.4587747),
"""
Statistics.median(
    ::AbstractManifold,
    ::AbstractVector,
    ::AbstractVector,
    ::WeiszfeldEstimation,
)

function Statistics.median(
    M::AbstractManifold,
    x::AbstractVector,
    method::AbstractEstimationMethod...;
    kwargs...,
)
    y = allocate_result(M, median, x[1])
    return median!(M, y, x, method...; kwargs...)
end
function Statistics.median(
    M::AbstractManifold,
    x::AbstractVector,
    w::AbstractVector,
    method::AbstractEstimationMethod...;
    kwargs...,
)
    y = allocate_result(M, median, x[1])
    return median!(M, y, x, w, method...; kwargs...)
end

@doc raw"""
    median!(M::AbstractManifold, y, x::AbstractVector[, w::AbstractWeights]; kwargs...)
    median!(
        M::AbstractManifold,
        y,
        x::AbstractVector,
        [w::AbstractWeights,]
        method::AbstractEstimationMethod;
        kwargs...,
    )

computes the [`median`](@ref) in-place in `y`.
"""
median!(::AbstractManifold, ::Any...)
function Statistics.median!(
    M::AbstractManifold,
    q,
    x::AbstractVector,
    method::AbstractEstimationMethod...;
    kwargs...,
)
    w = _unit_weights(length(x))
    return median!(M, q, x, w, method...; kwargs...)
end
function Statistics.median!(
    M::AbstractManifold,
    y,
    x::AbstractVector,
    w::AbstractVector;
    kwargs...,
)
    return median!(M, y, x, w, CyclicProximalPointEstimation(); kwargs...)
end
function Statistics.median!(
    M::AbstractManifold,
    q,
    x::AbstractVector,
    w::AbstractVector,
    ::CyclicProximalPointEstimation;
    p0=x[1],
    stop_iter=1000000,
    retraction::AbstractRetractionMethod=ExponentialRetraction(),
    inverse_retraction::AbstractInverseRetractionMethod=LogarithmicInverseRetraction(),
    kwargs...,
)
    n = length(x)
    if length(w) != n
        throw(
            DimensionMismatch(
                "The number of weights ($(length(w))) does not match the number of points for the median ($(n)).",
            ),
        )
    end
    copyto!(q, p0)
    yold = allocate_result(M, median, q)
    ytmp = copy(yold)
    v = zero_vector(M, q)
    wv = convert(AbstractVector, w) ./ sum(w)
    for i in 1:stop_iter
        λ = 0.5 / i
        copyto!(yold, q)
        for j in 1:n
            @inbounds t = min(λ * wv[j] / distance(M, q, x[j]), 1.0)
            @inbounds inverse_retract!(M, v, q, x[j], inverse_retraction)
            retract!(M, ytmp, q, v, t, retraction)
            copyto!(q, ytmp)
        end
        isapprox(M, q, yold; kwargs...) && break
    end
    return q
end

function Statistics.median!(
    M::AbstractManifold,
    y,
    x::AbstractVector,
    w::AbstractVector,
    ::ExtrinsicEstimation;
    extrinsic_method::AbstractEstimationMethod=CyclicProximalPointEstimation(),
    kwargs...,
)
    embedded_x = map(p -> embed(M, p), x)
    embedded_y = median(get_embedding(M), embedded_x, w, extrinsic_method; kwargs...)
    project!(M, y, embedded_y)
    return y
end

function Statistics.median!(
    M::AbstractManifold,
    q,
    x::AbstractVector,
    w::AbstractVector,
    ::WeiszfeldEstimation;
    p0=x[1],
    stop_iter=2000,
    α=1.0,
    retraction::AbstractRetractionMethod=ExponentialRetraction(),
    inverse_retraction::AbstractInverseRetractionMethod=LogarithmicInverseRetraction(),
    kwargs...,
)
    n = length(x)
    if length(w) != n
        throw(
            DimensionMismatch(
                "The number of weights ($(length(w))) does not match the number of points for the median ($(n)).",
            ),
        )
    end
    copyto!(q, p0)
    yold = allocate_result(M, median, q)
    ytmp = copy(yold)
    d = zeros(n)
    v = zero_vector(M, q)
    wv = convert(AbstractVector, w) ./ sum(w)
    for i in 1:stop_iter
        d .= [distance(M, q, xi) for xi in x] # compute distances
        # compute new weights / exclude points xi=q
        d .= [di > 0 ? wi / di : zero(typeof(wi / di)) for (di, wi) in zip(d, w)]
        copyto!(yold, q)
        zero_vector!(M, v, q)
        for j in 1:n
            @inbounds v .+= d[j] * inverse_retract(M, q, x[j], inverse_retraction)
        end
        retract!(M, ytmp, q, v, α / sum(d), retraction)
        copyto!(q, ytmp)
        isapprox(M, q, yold; kwargs...) && break
    end
    return q
end

@decorator_transparent_signature Statistics.median!(
    M::AbstractDecoratorManifold,
    y,
    x::AbstractVector,
    w::AbstractVector;
    kwargs...,
)

function decorator_transparent_dispatch(
    ::typeof(median!),
    ::AbstractEmbeddedManifold,
    args...;
    kwargs...,
)
    return Val(:parent)
end
function decorator_transparent_dispatch(
    ::typeof(median!),
    ::AbstractEmbeddedManifold{𝔽,<:TransparentIsometricEmbedding},
    args...;
    kwargs...,
) where {𝔽}
    return Val(:transparent)
end

@doc raw"""
    var(M, x, m=mean(M, x); corrected=true)
    var(M, x, w::AbstractWeights, m=mean(M, x, w); corrected=false)

compute the (optionally weighted) variance of a `Vector` `x` of `n` data points
on the [`AbstractManifold`](@ref) `M`, i.e.

````math
\frac{1}{c} \sum_{i=1}^n w_i d_{\mathcal M}^2 (x_i,m),
````
where `c` is a correction term, see
[Statistics.var](https://juliastats.org/StatsBase.jl/stable/scalarstats/#Statistics.var).
The mean of `x` can be specified as `m`, and the corrected variance
can be activated by setting `corrected=true`. All further `kwargs...` are passed
to the computation of the mean (if that is not provided).
"""
var(::AbstractManifold, ::Any)
function Statistics.var(
    M::AbstractManifold,
    x::AbstractVector,
    w::AbstractWeights,
    m;
    corrected::Bool=false,
)
    wv = convert(AbstractVector, w)
    s = sum(eachindex(x, w)) do i
        return @inbounds w[i] * distance(M, m, x[i])^2
    end
    c = StatsBase.varcorrection(w, corrected)
    return c * s
end
function Statistics.var(M::AbstractManifold, x::AbstractVector, m; corrected::Bool=true)
    n = length(x)
    w = _unit_weights(n)
    return var(M, x, w, m; corrected=corrected)
end
function Statistics.var(
    M::AbstractManifold,
    x::AbstractVector,
    w::AbstractWeights;
    kwargs...,
)
    return mean_and_var(M, x, w; kwargs...)[2]
end
function Statistics.var(M::AbstractManifold, x::AbstractVector; kwargs...)
    return mean_and_var(M, x; kwargs...)[2]
end

@doc raw"""
    std(M, x, m=mean(M, x); corrected=true, kwargs...)
    std(M, x, w::AbstractWeights, m=mean(M, x, w); corrected=false, kwargs...)

compute the optionally weighted standard deviation of a `Vector` `x` of `n` data
points on the [`AbstractManifold`](@ref) `M`, i.e.

````math
\sqrt{\frac{1}{c} \sum_{i=1}^n w_i d_{\mathcal M}^2 (x_i,m)},
````
where `c` is a correction term, see
[Statistics.std](https://juliastats.org/StatsBase.jl/stable/scalarstats/#Statistics.std).
The mean of `x` can be specified as `m`, and the corrected variance
can be activated by setting `corrected=true`.
"""
Statistics.std(M::AbstractManifold, args...; kwargs...) = sqrt(var(M, args...; kwargs...))

@doc raw"""
    mean_and_var(M::AbstractManifold, x::AbstractVector[, w::AbstractWeights]; kwargs...) -> (mean, var)

Compute the [`mean`](@ref mean(::AbstractManifold, args...)) and the [`var`](@ref)iance
simultaneously. See those functions for a description of the arguments.

    mean_and_var(
        M::AbstractManifold,
        x::AbstractVector
        [w::AbstractWeights,]
        method::AbstractEstimationMethod;
        kwargs...,
    ) -> (mean, var)

Use the `method` for simultaneously computing the mean and variance. To use
a mean-specific method, call [`mean`](@ref mean(::AbstractManifold, args...)) and then
[`var`](@ref).
"""
mean_and_var(M::AbstractManifold, ::Any...)
function StatsBase.mean_and_var(
    M::AbstractManifold,
    x::AbstractVector,
    w::AbstractWeights,
    method::AbstractEstimationMethod...;
    corrected=false,
    kwargs...,
)
    m = mean(M, x, w, method...; kwargs...)
    v = var(M, x, w, m; corrected=corrected)
    return m, v
end
function StatsBase.mean_and_var(
    M::AbstractManifold,
    x::AbstractVector,
    method::AbstractEstimationMethod...;
    corrected=true,
    kwargs...,
)
    n = length(x)
    w = _unit_weights(n)
    return mean_and_var(M, x, w, method...; corrected=corrected, kwargs...)
end

@doc raw"""
    mean_and_var(
        M::AbstractManifold,
        x::AbstractVector
        [w::AbstractWeights,]
        method::GeodesicInterpolation;
        shuffle_rng::Union{AbstractRNG,Nothing} = nothing,
        retraction::AbstractRetractionMethod = ExponentialRetraction(),
        inverse_retraction::AbstractInverseRetractionMethod = LogarithmicInverseRetraction(),
        kwargs...,
    ) -> (mean, var)

Use the repeated weighted geodesic interpolation to estimate the mean.
Simultaneously, use a Welford-like recursion to estimate the variance.

If `shuffle_rng` is provided, it is used to shuffle the order in which the
points are considered. Optionally, pass `retraction` and `inverse_retraction`
method types to specify the (inverse) retraction.

See [`GeodesicInterpolation`](@ref) for details on the geodesic
interpolation method.

!!! note
    The Welford algorithm for the variance is experimental and is not guaranteed
    to give accurate results except on [`Euclidean`](@ref).
"""
function StatsBase.mean_and_var(
    M::AbstractManifold,
    x::AbstractVector,
    w::AbstractWeights,
    ::GeodesicInterpolation;
    shuffle_rng::Union{AbstractRNG,Nothing}=nothing,
    corrected=false,
    retraction::AbstractRetractionMethod=ExponentialRetraction(),
    inverse_retraction::AbstractInverseRetractionMethod=LogarithmicInverseRetraction(),
    kwargs...,
)
    n = length(x)
    if length(w) != n
        throw(
            DimensionMismatch(
                "The number of weights ($(length(w))) does not match the number of points for the mean ($(n)).",
            ),
        )
    end
    order = shuffle_rng === nothing ? (1:n) : shuffle(shuffle_rng, 1:n)
    @inbounds begin
        j = order[1]
        s = w[j]
        y = copy(x[j])
    end
    v = zero_vector(M, y)
    M₂ = zero(number_eltype(v))
    ytmp = allocate_result(M, mean, y)
    @inbounds for i in 2:n
        j = order[i]
        snew = s + w[j]
        t = w[j] / snew
        inverse_retract!(M, v, y, x[j], inverse_retraction)
        retract!(M, ytmp, y, v, t, retraction)
        d = norm(M, y, v)
        copyto!(y, ytmp)
        M₂ += t * s * d^2
        s = snew
    end
    c = StatsBase.varcorrection(w, corrected)
    σ² = c * M₂
    return y, σ²
end

"""
    mean_and_var(
        M::AbstractManifold,
        x::AbstractVector
        [w::AbstractWeights,]
        method::GeodesicInterpolationWithinRadius;
        kwargs...,
    ) -> (mean, var)

Use repeated weighted geodesic interpolation to estimate the mean.
Simultaneously, use a Welford-like recursion to estimate the variance.

See [`GeodesicInterpolationWithinRadius`](@ref) and
[`mean_and_var`](@ref mean_and_var(::AbstractManifold, ::AbstractVector, ::AbstractWeights, ::GeodesicInterpolation))
for more information.
"""
function StatsBase.mean_and_var(
    M::AbstractManifold,
    x::AbstractVector,
    w::AbstractWeights,
    method::GeodesicInterpolationWithinRadius;
    shuffle_rng=nothing,
    corrected=false,
    kwargs...,
)
    y, v = mean_and_var(
        M,
        x,
        w,
        GeodesicInterpolation();
        shuffle_rng=shuffle_rng,
        corrected=corrected,
        kwargs...,
    )
    radius = method.radius
    injectivity_radius(M, y) ≤ radius && return y, v
    for i in eachindex(x)
        @inbounds if distance(M, y, x[i]) ≥ radius
            mean!(M, y, x, w, GradientDescentEstimation(); p0=y, kwargs...)
            v = var(M, x, w, y; corrected=corrected)
            return y, v
        end
    end
    return y, v
end

@doc raw"""
    mean_and_std(M::AbstractManifold, x::AbstractVector[, w::AbstractWeights]; kwargs...) -> (mean, std)

Compute the [`mean`](@ref mean(::AbstractManifold, args...)) and the standard deviation
[`std`](@ref) simultaneously.

    mean_and_std(
        M::AbstractManifold,
        x::AbstractVector
        [w::AbstractWeights,]
        method::AbstractEstimationMethod;
        kwargs...,
    ) -> (mean, var)

Use the `method` for simultaneously computing the mean and standard deviation.
To use a mean-specific method, call [`mean`](@ref mean(::AbstractManifold, args...)) and
then [`std`](@ref).
"""
function StatsBase.mean_and_std(M::AbstractManifold, args...; kwargs...)
    m, v = mean_and_var(M, args...; kwargs...)
    return m, sqrt(v)
end

"""
    moment(M::AbstractManifold, x::AbstractVector, k::Int[, w::AbstractWeights], m=mean(M, x[, w]))

Compute the `k`th central moment of points in `x` on manifold `M`. Optionally
provide weights `w` and/or a precomputed
[`mean`](@ref mean(::AbstractManifold, args...)).
"""
function StatsBase.moment(
    M::AbstractManifold,
    x::AbstractVector,
    k::Int,
    w::AbstractWeights,
    m=mean(M, x, w),
)
    s = sum(eachindex(x, w)) do i
        return @inbounds w[i] * distance(M, m, x[i])^k
    end
    return s / sum(w)
end
function StatsBase.moment(M::AbstractManifold, x::AbstractVector, k::Int, m=mean(M, x))
    w = _unit_weights(length(x))
    return moment(M, x, k, w, m)
end

"""
    skewness(M::AbstractManifold, x::AbstractVector, k::Int[, w::AbstractWeights], m=mean(M, x[, w]))

Compute the standardized skewness of points in `x` on manifold `M`. Optionally
provide weights `w` and/or a precomputed
[`mean`](@ref mean(::AbstractManifold, args...)) `m`.
"""
function StatsBase.skewness(M::AbstractManifold, x::AbstractVector, w::AbstractWeights)
    m, s = mean_and_std(M, x, w; corrected=false)
    return moment(M, x, 3, w, m) / s^3
end
function StatsBase.skewness(M::AbstractManifold, x::AbstractVector, w::AbstractWeights, m)
    return moment(M, x, 3, w, m) / std(M, x, w, m; corrected=false)^3
end
function StatsBase.skewness(M::AbstractManifold, x::AbstractVector, args...)
    w = _unit_weights(length(x))
    return skewness(M, x, w, args...)
end

"""
    kurtosis(M::AbstractManifold, x::AbstractVector, k::Int[, w::AbstractWeights], m=mean(M, x[, w]))

Compute the excess kurtosis of points in `x` on manifold `M`. Optionally
provide weights `w` and/or a precomputed
[`mean`](@ref mean(::AbstractManifold, args...)) `m`.
"""
function StatsBase.kurtosis(M::AbstractManifold, x::AbstractVector, w::AbstractWeights)
    m, v = mean_and_var(M, x, w; corrected=false)
    return moment(M, x, 4, w, m) / v^2 - 3
end
function StatsBase.kurtosis(M::AbstractManifold, x::AbstractVector, w::AbstractWeights, m)
    return moment(M, x, 4, w, m) / var(M, x, w, m; corrected=false)^2 - 3
end
function StatsBase.kurtosis(M::AbstractManifold, x::AbstractVector, args...)
    w = _unit_weights(length(x))
    return kurtosis(M, x, w, args...)
end


text:
"""
stores functions needed to create bind and update OpenGl textues
"""
module TextureManag
using Base: Float16
using  ModernGL ,   ..OpenGLDisplayUtils,  ..ForDisplayStructs
using   ..Uniforms, Logging,Setfield, Logging,  ..CustomFragShad,  ..DataStructs,  ..DisplayWords
export activateTextures,addTextToTexture,initializeTextures,createTexture, getProperGL_TEXTURE,updateImagesDisplayed, updateTexture, assignUniformsAndTypesToMasks

"""
uploading data to given texture; of given types associated - specified in TextureSpec
if we want to update only part of the texture we need to specify  offset and size of texture we use
Just for reference openGL function definition
    void glTextureSubImage2D(
        GLuint texture,
         GLint level,
         GLint xoffset,
         GLint yoffset,
         GLsizei width,
         GLsizei height,
         GLenum format,
         GLenum type,
         const void *pixels);
"""
function updateTexture(::Type{Tt}
                    ,data::AbstractArray
                    ,textSpec::TextureSpec
                    ,xoffset::Int
                    ,yoffset::Int
                    ,widthh::Int32
                    ,heightt::Int32) where{Tt}

    glActiveTexture(textSpec.actTextrureNumb); # active proper texture unit before binding
    glBindTexture(GL_TEXTURE_2D, textSpec.ID[]);

    if((parameter_type(textSpec)== Float16) || (parameter_type(textSpec)== Float32))
        glTexSubImage2D(GL_TEXTURE_2D,0,xoffset,yoffset, widthh, heightt, GL_RED, textSpec.OpGlType, collect(data))
    else
	    glTexSubImage2D(GL_TEXTURE_2D,0,xoffset,yoffset, widthh, heightt, GL_RED_INTEGER, textSpec.OpGlType, collect(data))
	   # glTexSubImage2D(GL_TEXTURE_2D,0,xoffset,yoffset, widthh, heightt, GL_RED_INTEGER, textSpec.OpGlType, reduce(vcat,data))

    end

end

"""
creating texture that is storing values like integer, uint, float values that are representing
main image or  mask data and which will be used by a shader to draw appropriate colors
juliaDataType- data type defined as a Julia datatype of data we are dealing with
width,height - dimensions of the texture that we need
GL_RType,OpGlType - Open Gl types needed to properly specify the  texture they need to be compatible with juliaDataType
"""
function createTexture(juliaDataType::Type{juliaDataTyp}
                        , width::Int32
                        , height::Int32
                        ,GL_RType::UInt32 =GL_R8UI
                        ,OpGlType= GL_UNSIGNED_BYTE) where {juliaDataTyp}

#The texture we're going to render to
    texture= Ref(GLuint(0));
    glGenTextures(1, texture);
    glBindTexture(GL_TEXTURE_2D, texture[]);

    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_BASE_LEVEL, 0);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAX_LEVEL, 0);
#we just assign storage
    glTexStorage2D(GL_TEXTURE_2D, 1, GL_RType, width, height);

return texture
end

"""
initializing textures  - so we basically execute specification for configuration and bind all to OpenGL context
listOfTextSpecs - list of TextureSpec structs that  holds data needed to
calcDimStruct - struct holding necessery data about taxture dimensions, quad propertiess etc.
"""
function initializeTextures(listOfTextSpecs::Vector{TextureSpec}
                            ,calcDimStruct ::CalcDimsStruct)::Vector{TextureSpec}

    res = Vector{TextureSpec}()

    for (ind, textSpec ) in enumerate(listOfTextSpecs)
        index=ind-1
        textUreId= createTexture(parameter_type(textSpec),calcDimStruct.imageTextureWidth,calcDimStruct.imageTextureHeight,textSpec.GL_Rtype,textSpec.OpGlType )#binding texture and populating with data

        actTextrureNumb = getProperGL_TEXTURE(index)
        glActiveTexture(actTextrureNumb)
        glUniform1i(textSpec.uniforms.samplerRef,index);# we first look for uniform sampler in shader
        # we set ..Uniforms of visibility and colors according to specified in configuration
        if(!textSpec.isMainImage) setMaskColor(textSpec.color ,textSpec.uniforms)   end

        setTextureVisibility(textSpec.isVisible ,textSpec.uniforms)

        push!(res,setproperties(textSpec, (ID=textUreId, actTextrureNumb=actTextrureNumb,associatedActiveNumer=index )))

    end # for

    return res
end #initializeAndDrawTextures

"""
activating textures that were already initialized in order to be able to use them with diffrent shader program
shader_program- regference to OpenGL program so we will be able to activate  textures
listOfTextSpecs - list of TextureSpec structs that  holds data needed to bind textures to shader program (Hovewer this new shader program have to keep the same ..Uniforms)
return unmodified textures
"""
function activateTextures(listOfTextSpecs::Vector{TextureSpec} )::Vector{TextureSpec}

    for (ind, textSpec ) in enumerate(listOfTextSpecs)
        glBindTexture(GL_TEXTURE_2D, textSpec.ID[]);
        glUniform1i(textSpec.uniforms.samplerRef,textSpec.associatedActiveNumer);# we first look for uniform sampler in shader
        # we set ..Uniforms of visibility and colors according to specified in configuration
        if(!textSpec.isMainImage) setMaskColor(textSpec.color ,textSpec.uniforms)   end
        setTextureVisibility(textSpec.isVisible ,textSpec.uniforms)
    end # for

    return listOfTextSpecs

end#activateTextures

"""
associates GL_TEXTURE UInt32 to given index
"""
function getProperGL_TEXTURE(index::Int)::UInt32
    return eval(Meta.parse("GL_TEXTURE$(index)"))
end#getProperGL_TEXTURE

"""
coordinating updating all of the images, masks...
singleSliceDat - holds data we want to use for update
forDisplayObjects - stores all needed constants that holds reference to GLFW and OpenGL
"""
function updateImagesDisplayed(singleSliceDat::SingleSliceDat
                            ,forDisplayConstants::forDisplayObjects
                            ,wordsDispObj::ForWordsDispStruct
                            ,calcDimStruct::CalcDimsStruct
                            ,valueForMaskToSett::valueForMasToSetStruct )

        forDisplayConstants.stopListening[]=true# freeing GLFW from listening for event
             modulelistOfTextSpecs=forDisplayConstants.listOfTextSpecifications

            for updateDat in singleSliceDat.listOfDataAndImageNames
                findList= findall( (texSpec)-> texSpec.name == updateDat.name, modulelistOfTextSpecs)
                texSpec = !isempty(findList) ? modulelistOfTextSpecs[findList[1]] : throw(DomainError(findList, "no such name specified in start configuration - $( updateDat[1])"))
                updateTexture(updateDat.type,updateDat.dat,texSpec,0,0,calcDimStruct.imageTextureWidth,calcDimStruct.imageTextureHeight )
            end #for
            #render text associated with this slice
            activateForTextDisp(
                wordsDispObj.shader_program_words
                ,wordsDispObj.vbo_words
                ,calcDimStruct)

            matr= addTextToTexture(wordsDispObj
                            ,[singleSliceDat.textToDisp...,valueForMaskToSett.text]
                            ,calcDimStruct )

            glDrawElements(GL_TRIANGLES, 6, GL_UNSIGNED_INT, C_NULL)

            reactivateMainObj(forDisplayConstants.shader_program
                            ,forDisplayConstants.vbo
                            ,calcDimStruct )

            #render onto the screen
             OpenGLDisplayUtils.basicRender(forDisplayConstants.window)
            glFinish()
            forDisplayConstants.stopListening[]=false
end

"""
on the basis of the type supplied in texture characteristic
it supplies given set of ..Uniforms to it
It would also assign proper openGl types to given julia data type, and pass data from texture specification to opengl context
textSpecs - list of texture specificaton that we want to enrich by adding information about ..Uniforms
return list of texture specifications enriched by information about ..Uniforms
"""
function assignUniformsAndTypesToMasks(textSpecs::Vector{TextureSpec}
                ,shader_program::UInt32
                ,windowControlStruct::WindowControlStruct)
    mainTexture,notMainTextures=   divideTexteuresToMainAndRest(textSpecs)
#main texture ..Uniforms
n= mainTexture.name

mainUnifs =MainImageUniforms(
    samplerName= n
    ,samplerRef = glGetUniformLocation(shader_program, n)
    ,isVisibleRef = glGetUniformLocation(shader_program, "$(n)isVisible")
    ,min_shown_white= glGetUniformLocation(shader_program, "min_shown_white")
   , max_shown_black= glGetUniformLocation(shader_program, "max_shown_black")
    ,displayRange=glGetUniformLocation(shader_program, "displayRange")
    ,isMaskDiffrenceVis=glGetUniformLocation(shader_program, "isMaskDiffrenceVis")
    ,mainImageContribution=glGetUniformLocation(shader_program, "mainImageContribution")
)
setCTWindow(windowControlStruct.min_shown_white,windowControlStruct.max_shown_black,mainUnifs)

maintext = setproperties(mainTexture, (uniforms= mainUnifs ))
# joining main and not main textures data
mapped= [ map(x-> setuniforms(x,shader_program) , notMainTextures)..., maintext]

return  (mainUnifs, map(x->setProperOpenGlTypes(x),mapped))

end#assign..UniformsToMasks

 """
helper for assign..UniformsToMasks
On the basis of the name of the Texture it will assign the informs referencs to it
- ..Uniforms for main image will be set separately

"""
function setuniforms(textSpec:: TextureSpec,shader_program::UInt32):: TextureSpec

    n= textSpec.name
    unifs = MaskTextureUniforms(samplerName= n
                        ,samplerRef= glGetUniformLocation(shader_program, n)
                        ,colorsMaskRef=glGetUniformLocation(shader_program, "$(n)ColorMask")
                        ,maskMinValue=glGetUniformLocation(shader_program, "$(n)minValue")
                        ,maskMAxValue=glGetUniformLocation(shader_program, "$(n)maxValue")
                        ,maskRangeValue=glGetUniformLocation(shader_program, "$(n)ValueRange")
                        ,isVisibleRef=glGetUniformLocation(shader_program, "$(n)isVisible")
                        ,maskContribution=glGetUniformLocation(shader_program, "$(n)maskContribution"))

    return setproperties(textSpec, (uniforms= unifs ))

end#assign..UniformsToMasks

 """
On the basis of the type associated to texture we set proper open Gl types associated
based on https://www.khronos.org/registry/OpenGL-Refpages/gl4/html/glTexImage2D.xhtml
and https://www.khronos.org/opengl/wiki/OpenGL_Type
"""
function setProperOpenGlTypes(textSpec:: TextureSpec):: TextureSpec
    if(parameter_type(textSpec)== Float16 ) return  setproperties(textSpec, (GL_Rtype= GL_R16F ,OpGlType= GL_HALF_FLOAT ))     end
    if(parameter_type(textSpec)== Float32 ) return  setproperties(textSpec, (GL_Rtype= GL_R32F ,OpGlType= GL_FLOAT))     end
    if(parameter_type(textSpec)== Int8 ) return  setproperties(textSpec, (GL_Rtype= GL_R8I ,OpGlType= GL_BYTE ))     end
    if(parameter_type(textSpec)== UInt8 ) return  setproperties(textSpec, (GL_Rtype= GL_R8UI,OpGlType= GL_UNSIGNED_BYTE ))     end
    if(parameter_type(textSpec)== Int16 ) return  setproperties(textSpec, (GL_Rtype= GL_R16I,OpGlType= GL_SHORT ))     end
    if(parameter_type(textSpec)== UInt16) return  setproperties(textSpec, (GL_Rtype=GL_R16UI ,OpGlType=GL_UNSIGNED_SHORT))     end
    if(parameter_type(textSpec)== Int32 ) return  setproperties(textSpec, (GL_Rtype= GL_R32I,OpGlType= GL_INT ))     end
    if(parameter_type(textSpec)== UInt32) return  setproperties(textSpec, (GL_Rtype=GL_R32UI ,OpGlType= GL_UNSIGNED_INT))     end

    throw(DomainError(textSpec, "type  of texture is not supported - supported types - Int8,16,32 UInt 8,16,32 float16,32"))
end#

"""
Given  vector of SimpleLineTextStructs it will return matrix of data that will be used
to display text
lines - data about text to be displayed
calcDimStruct - struct holding important data about size of textures etc.
wordsDispObj - object wit needed constants to display text
"""
function addTextToTexture(wordsDispObj::ForWordsDispStruct
                          ,lines::Vector{SimpleLineTextStruct}
                          ,calcDimStruct::CalcDimsStruct)
    textureWidth = calcDimStruct.textTexturewidthh
    fontFace= wordsDispObj.fontFace

    matrPrim=  map(x-> renderSingleLineOfText(x,textureWidth,fontFace) ,reverse(lines)) |>
    (xl)-> reduce( hcat  ,xl)

    sz= size(matrPrim)
    # below just to clear any data from texture uploaded before
    matr= hcat(calcDimStruct.textTextureZeros[:,sz[2]:size(calcDimStruct.textTextureZeros)[2]-1] ,matrPrim  )

    updateTexture(UInt8
                ,matr
                ,wordsDispObj.textureSpec
                ,0
                ,0
                ,calcDimStruct.textTexturewidthh
                ,calcDimStruct.textTextureheightt) #  ,Int32(10000),Int32(1000)
    return matr
end #addTextToTexture

########## puts bytes of image into PBO as fas as I get it  copy an image data to texture buffer

# preparePixelBufferStr="""
# width -width of the image in  number of pixels
# height - height of the image in  number of pixels
# pboNumber - just states which PBO it is
# return reference to the pixel buffer object that we use to upload this texture and data size calculated for this texture
# """
# @doc preparePixelBufferStr
# function preparePixelBuffer(juliaDataTyp::Type{juliaDataType},width,height,pboNumber)where{juliaDataType}
#     DATA_SIZE = 8 * sizeof(juliaDataTyp) *width * height  # number of bytes our image will have so in 2D it will be width times height times number of bytes needed for used datatype we need to multiply by 8 becouse sizeof() return bytes instead of bits
#     pbo = Ref(GLuint(pboNumber))
#     glGenBuffers(1, pbo)
#     return (pbo,DATA_SIZE)
# end

# usePixelBuferAndUploadDataStr = """
# adapted from http://www.songho.ca/opengl/gl_pbo.html
# creates single pixel buffer of given type
# pboID - id of the pixel buffer object that was prepared for some particular texture
# textureId - reference to id of a texture that we want to bind to this PBO
# juliaDataType -julia type that is representing datatype in 2 dimensional array representing ima
# width -width of the image in  number of pixels
# height - height of the image in  number of pixels
# subImageDataType - variable used in glTexSubImage2D to tell open Glo what type of data is in texture
# data one dimensional array o julia type and width*height length
# DATA_SIZE - size of texture in bytes
# """
# @doc usePixelBuferAndUploadDataStr
# function usePixelBuferAndUploadData(
#     juliaDataTyp::Type{juliaDataType}
#                     ,pboID
#                     ,width
#                     ,height
#                     ,data
#                     ,textureId
#                     ,DATA_SIZE
#                     ,subImageDataType = GL_SHORT

#                     )where{juliaDataType}

#     glBindTexture(GL_TEXTURE_2D,textureId[]);
#     # copy pixels from PBO to texture object
#     # Use offset instead of pointer.
#    # glTexSubImage2D(GL_TEXTURE_2D_ARRAY, 0, 0, 0, GLsizei(width), GLsizei(height),  GL_RED_INTEGER, GL_SHORT, Ptr{juliaDataTyp}());

#     glTexSubImage2D(GL_TEXTURE_2D,0,0,0, width, height, GL_RED_INTEGER, subImageDataType, Ptr{juliaDataType}());

#     # bind the PBO
#     glBindBuffer(GL_PIXEL_UNPACK_BUFFER, pboID[]);

#     # Note that glMapBuffer() causes sync issue.
#     # If GPU is working with this buffer, glMapBuffer() will wait(stall)
#     # until GPU to finish its job. To avoid waiting (idle), you can call
#     # first glBufferData() with NULL pointer before glMapBuffer().
#     # If you do that, the previous data in PBO will be discarded and
#     # glMapBuffer() returns a new allocated pointer immediately
#     # even if GPU is still working with the previous data.
#     glBufferData(GL_PIXEL_UNPACK_BUFFER, DATA_SIZE, Ptr{juliaDataType}(), GL_STREAM_DRAW);

#     # map the buffer object into client's memory
#     glMapBuffer(GL_ARRAY_BUFFER, GL_WRITE_ONLY)

#     ptr = Ptr{juliaDataType}(glMapBuffer(GL_PIXEL_UNPACK_BUFFER, GL_WRITE_ONLY))
#     # update data directly on the mapped buffer - this is internal function implemented below

#     updatePixels(ptr,data,length(data));

#     glUnmapBuffer(GL_PIXEL_UNPACK_BUFFER); # release the mapped buffer

#     # it is good idea to release PBOs with ID 0 after use.
#     # Once bound with 0, all pixel operations are back to normal ways.
#     glBindBuffer(GL_PIXEL_UNPACK_BUFFER, 0);

# end

# updatePixelsStr = """
# adapted from https://github.com/JuliaPlots/GLMakie.jl/blob/2717d812fdc66b283f63d5d97237e8d69e2c1f25/src/GLAbstraction/GLBuffer.jl from unsafe copy
# """
# @doc updatePixelsStr
# function updatePixels(ptr, data,length)
#     for i=1:length
#         unsafe_store!(ptr,data[i], i)
#     end
# end

end #..TextureManag

ext: jl


text:
module VPRK

    using Documenter: @doc
    using LinearAlgebra: mul!
    using RungeKutta
    using SimpleSolvers

    import GeometricEquations: IODE, LODE, LDAE, get_functions, hassecondary

    using ..GeometricBase
    using ..Config
    using ..Utils

    import ..Solutions: AtomicSolutionPODE, AtomicSolutionPDAE, SolutionPDAE, SolutionVector
    import ..Solutions: update!

    import ..Integrators

    import ..Integrators: IODEIntegrator, IODEIntegratorCache, InitialGuessIODE,
                          AbstractIntegratorIRK, AbstractIntegratorPRK
    import ..Integrators: IntegratorCache, CacheDict, CacheType, Parameters
    import ..Integrators: AbstractTableau, AbstractCoefficients,
                          CoefficientsPGLRK,
                          @CoefficientsRK, @HeaderTableau, @HeaderCoefficientsRK,
                          get_symplectic_conjugate_coefficients
    import ..Integrators: create_internal_stage_vector, create_nonlinear_solver,
                          update_vector_fields!, update_solution!, update_multiplier!,
                          initialize!
    import ..Integrators: equation, equations, tableau, timestep,
                          eachdim, eachstage, nstages

    export IntegratorVPRK, IntegratorVPRKpNone, TableauVPRK

    export IntegratorVPRKdegenerate

    export IntegratorVPRKpStandard, IntegratorVPRKpSymplectic,
         IntegratorVPRKpInternal, IntegratorVPRKpMidpoint,
         IntegratorVPRKpSymmetric, IntegratorVPRKpTableau,
         IntegratorVPRKpSecondary, IntegratorVPRKpVariational,
         IntegratorVPRKpVariationalQ, IntegratorVPRKpVariationalP

    export IntegratorVPRKpLegendre#, TableauVPRKpLegendre

    include("vprk/integrators_vprk_abstract.jl")
    include("vprk/integrators_vprk_cache.jl")
    include("vprk/integrators_vprk_tableau.jl")
    include("vprk/integrators_vprk_parameters.jl")
    include("vprk/integrators_vprk_common.jl")
    include("vprk/integrators_vprk.jl")
    include("vprk/integrators_vprk_degenerate.jl")
    include("vprk/integrators_vprk_pinternal.jl")
    include("vprk/integrators_vprk_pmidpoint.jl")
    include("vprk/integrators_vprk_pstandard.jl")
    include("vprk/integrators_vprk_psecondary.jl")
    include("vprk/integrators_vprk_psymmetric.jl")
    include("vprk/integrators_vprk_ptableau.jl")
    include("vprk/integrators_vprk_pvariational.jl")
    include("vprk/integrators_vprk_plegendre.jl")

end


text:
using Base.Filesystem

println(isfile("pretrained.bson"))

println(isfile("hoge"))

include("util.jl")
include("./dataset.jl")
include("./iterator.jl")

using Printf
using Flux
using Statistics
using Flux: onehotbatch, onecold, crossentropy
using Base.Iterators: partition
using Base.Filesystem
using Metalhead:VGG19
using BSON: @load, @save
using CuArrays

function define_model()
    vgg = VGG19()
    model = Chain(vgg.layers[1:end-2],
                  Dense(4096, 101),
                  softmax)
    Flux.testmode!(model, false)
    return model
end

function get_dataset(datasetdir)
    train_paris, val_pairs = make_trainval_pairs(datasetdir)
    train_dataset = Dataset(train_paris)
    val_dataset = Dataset(val_pairs, train=false)
    return train_dataset, val_dataset
end

function main()
    datasetdir = expanduser("~/dataSSD120GB/food-101")
    batchsize = 32
    epochs = 100
    cache_multiplier = 10
    train_dataset, val_dataset = get_dataset(datasetdir)

    model = define_model() |> gpu
    if isfile("checkpoint_weights.bson")
        println("loading checkpoin file")
        @load "checkpoint_weights.bson" checkpoint_weights
        Flux.loadparams!(model, weights)

    loss(x,y)= crossentropy(model(x), y)
    accuracy(x, y) = mean(onecold(model(x)) .== onecold(y))
    optimizer = ADAM(params(model))
    for e in 1:epochs
        println("train loop ", e," / ",epochs)
        train_iter = SerialIterator(train_dataset, cache_multiplier * batchsize)
        val_iter = SerialIterator(val_dataset, cache_multiplier * batchsize, shuffle=false)
        for (i, batch) in enumerate(train_iter)
            println("progress ", i," / ", floor(Int, length(train_dataset) / batchsize / cache_multiplier))
            xs = [img for (img, _) in batch]
            ys = [label for (_, label) in batch]
            X = cat(xs..., dims=4)
            Y = onehotbatch(ys, 1:101)
            data = [(X[:,:,:,b] |> gpu, Y[:,b] |> gpu) for b in partition(1: cache_multiplier * batchsize, batchsize)]
            Flux.train!(loss, data, optimizer)
        end

        println("check accuracy")
        total_loss=0
        total_acc=0
        Flux.testmode!(model)
        cnt = 0
        for (i, batch) in enumerate(val_iter)
            println("progress ", i," / ", floor(Int, length(val_dataset) / batchsize / cache_multiplier))
            xs = [img for (img, _) in batch]
            ys = [label for (_, label) in batch]
            X = cat(xs..., dims=4)
            Y = onehotbatch(ys, 1:101)
            data = [(X[:,:,:,b] |> gpu, Y[:,b] |> gpu) for b in partition(1: cache_multiplier * batchsize, batchsize)]
            for (X,Y) in data
                total_loss += loss(X,Y)
                total_acc += accuracy(X,Y)
                cnt+=1
            end
        end
        Flux.testmode!(model, false)
        total_loss /= cnt
        total_acc /= cnt
        @printf("loss = %.3f\n", total_loss)
        @printf("acc = %.3f\n", total_acc)
        checkpoint = model |> cpu
        checkpoint_weights = Tracker.data.(params(checkpoint))
        @save "checkpoint.bson" checkpoint
        @save "checkpoint_weights.bson" checkpoint_weights

    end
    pretrained = model |> cpu
    weights = Tracker.data.(params(pretrained))
    @save "pretrained.bson" pretrained
    @save "weights.bson" weights
    println("Finished to train")
end

main()

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
# Autogenerated wrapper script for MbedTLS_jll for armv6l-linux-musleabihf
export libmbedcrypto, libmbedtls, libmbedx509

JLLWrappers.@generate_wrapper_header("MbedTLS")
JLLWrappers.@declare_library_product(libmbedcrypto, "libmbedcrypto.so.7")
JLLWrappers.@declare_library_product(libmbedtls, "libmbedtls.so.13")
JLLWrappers.@declare_library_product(libmbedx509, "libmbedx509.so.1")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libmbedcrypto,
        "lib/libmbedcrypto.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmbedtls,
        "lib/libmbedtls.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmbedx509,
        "lib/libmbedx509.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Unrolled
using Test
using StaticArrays
using MacroTools
using QuickTypes: type_parameters

# Check that we can extract argument names
@capture(:(function foo(a, b::Int, c=2; d::Int=4) end),
         (function foo(args__; kwargs__) end))
@test map(Unrolled.function_argument_name, vcat(args, kwargs)) == [:a, :b, :c, :d]

@unroll function my_sum(ss)
    total = zero(eltype(ss))
    @unroll for x in ss
        total += x
    end
    return total
end
@unroll function my_sum(; ss::Tuple=1)  # test kwargs
    total = zero(eltype(ss))
    @unroll for x in ss
        total += x
    end
    return total
end

@test my_sum((1,2,3)) == 6
@test my_sum([1,2,3]) == 6
@test my_sum(SVector(1,2,3)) == 6
@test my_sum(; ss=(1,2,3)) == 6

# Had to disable this test on 0.7. @test_throws looks broken?
# @test_throws AssertionError @eval @unroll function my_sum(ss)
#     total = zero(eltype(ss))
#     @unroll for x in ss[1:end-1]
#         total += x
#     end
#     return total
# end

@unroll function _do_sum(sub_seq) # helper for my_sum_but_last
    total = zero(eltype(sub_seq))
    @unroll for x in sub_seq
        total += x
    end
    return total
end

# Sum every number in seq except the last one
function my_sum_but_last(seq)
    return _do_sum(seq[1:end-1])
end

@test my_sum_but_last((1,20,3)) == 21

# Check that unrolled_union & co are correct and type-stable
struct MyVal{T} end
# Otherwise it's not type-stable. Perhaps it should defined in Base for all singleton
# types.
@generated Base.in(val::MyVal, tup::Tuple) = val in type_parameters(tup)

@test (@inferred(unrolled_union((MyVal{1}(), MyVal{2}()), (MyVal{2}(), MyVal{0}()))) ==
       (MyVal{1}(), MyVal{2}(), MyVal{0}()))
@test (@inferred(unrolled_intersect((MyVal{1}(), MyVal{2}()), (MyVal{2}(), MyVal{0}())))==
       (MyVal{2}(),))
@test (@inferred(unrolled_setdiff((MyVal{1}(), MyVal{2}()), (MyVal{2}(), MyVal{0}()))) ==
       (MyVal{1}(),))
@test (@inferred(unrolled_union((MyVal{1}(), MyVal{2}()),
                                (MyVal{2}(), MyVal{0}()),
                                (MyVal{10}(),))) ==
       (MyVal{10}(), MyVal{2}(), MyVal{0}(), MyVal{1}()))
@test @inferred(unrolled_reduce((+), 0, unrolled_map(abs, (1,2,-3,7)))) == 13

const tupl = (1,2,3,4,5.0,6,7...)
@test @inferred(getindex(tupl, FixedRange{4, 5}())) == (4, 5.0)
@test @inferred(getindex(tupl, FixedRange{4, FixedEnd{1}()}())) == (4, 5.0, 6)

f(tupl) = @fixed_range(tupl[4:end-1])
@test @inferred(f(tupl)) == (4, 5.0, 6)

@unroll function summer(tup, t::Type{T}) where T
    s = zero(T)
    @unroll for x in tup
        s += x
    end
    s
end

@test summer((1,2,3), Float64) === 6.0

x = [0.0]
unrolled_foreach((1,2,3, 1.0)) do y
    x[1] += y
end
@test x[1] == 7.0


text:
function BOBYQAOptimizer(arg0::jint)
    return BOBYQAOptimizer((jint,), arg0)
end

function BOBYQAOptimizer(arg0::jint, arg1::jdouble, arg2::jdouble)
    return BOBYQAOptimizer((jint, jdouble, jdouble), arg0, arg1, arg2)
end

function compute_objective_value(obj::MultivariateOptimizer, arg0::Vector{jdouble})
    return jcall(obj, "computeObjectiveValue", jdouble, (Vector{jdouble},), arg0)
end

function get_convergence_checker(obj::BaseOptimizer)
    return jcall(obj, "getConvergenceChecker", ConvergenceChecker, ())
end

function get_evaluations(obj::BaseOptimizer)
    return jcall(obj, "getEvaluations", jint, ())
end

function get_goal_type(obj::MultivariateOptimizer)
    return jcall(obj, "getGoalType", GoalType, ())
end

function get_iterations(obj::BaseOptimizer)
    return jcall(obj, "getIterations", jint, ())
end

function get_lower_bound(obj::BaseMultivariateOptimizer)
    return jcall(obj, "getLowerBound", Vector{jdouble}, ())
end

function get_max_evaluations(obj::BaseOptimizer)
    return jcall(obj, "getMaxEvaluations", jint, ())
end

function get_max_iterations(obj::BaseOptimizer)
    return jcall(obj, "getMaxIterations", jint, ())
end

function get_start_point(obj::BaseMultivariateOptimizer)
    return jcall(obj, "getStartPoint", Vector{jdouble}, ())
end

function get_upper_bound(obj::BaseMultivariateOptimizer)
    return jcall(obj, "getUpperBound", Vector{jdouble}, ())
end

function optimize(obj::BaseOptimizer)
    return jcall(obj, "optimize", Object, ())
end

function optimize(obj::MultivariateOptimizer, arg0::Vector{OptimizationData})
    return jcall(obj, "optimize", PointValuePair, (Vector{OptimizationData},), arg0)
end


text:
"""
$(SIGNATURES)
Extract important information from the SNOPT output file(s)
"""
function parse_snopt_summary(file=joinpath(pwd(),"snopt.out"))
    props = Dict()
    obj = Vector{Float64}()
    c_max = Vector{Float64}()
    next_line = false

    function stash_prop(ln::String,prop::String,prop_name::String=prop,vartype::Type=Float64)
        if occursin(prop,ln)
            loc = findfirst(prop,ln)
            if vartype <: Real
                val = split(ln[loc[end]+1:end])[1]
                val = convert(vartype,parse(Float64,val))
                props[prop_name] = val
                return true
            end
        end
        return false
    end

    function store_itervals(ln::String)
        vals = split(ln)

        if vals[6][1] == "("[1]

            idx = -1
            for i = 2:length(vals[6])
                if vals[6][i] == ")"[1] && idx == -1
                    idx = i-1
                    break
                end
            end
            _c = vals[6][2:idx]

        else
            _c = vals[6]
        end
        push!(obj, parse(Float64,vals[8])) # Merit function
        push!(c_max, parse(Float64,_c))

        next_line = false
    end

    open(file) do f
        for ln in eachline(f)
            stash_prop(ln,"No. of iterations","iterations",Int64)
            stash_prop(ln,"No. of major iterations","major iterations",Int64)
            stash_prop(ln,"No. of calls to funobj","objective calls",Int64)
            stash_prop(ln,"Nonlinear constraint violn","c_max")
            if next_line
                store_itervals(ln)
            end

            if length(ln) > 0 && split(ln)[1] == "Itns" && split(ln)[2] == "Major"
                next_line = true
            end
        end
    end
    props[:cost] = obj
    props[:c_max] = c_max
    return props
end

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


text:
using GPUifyLoops

function kernel3!(A)
    s1 = @shmem eltype(A) (1024,)
    s2 = @shmem eltype(A) (1024,)

    @loop for i in (1:size(A,1); threadIdx().x)
        s1[i] = 2*A[i]
        s2[i] = 3*A[i]
    end
    @synchronize
    @loop for i in (1:size(A,1); threadIdx().x)
        A[i] = s1[i]
    end
    nothing
end

data = rand(Float32, 1024)
cpudata = copy(data)

@launch CPU() kernel3!(cpudata)
@assert cpudata ≈ 2 .* data

@static if Base.find_package("CuArrays") !== nothing
  using CuArrays
  using CUDAnative

  cudata = CuArray(data)
  @launch CUDA() threads=length(cudata) kernel3!(cudata)
  @assert Array(cudata) ≈ 2 .* data
end


text:
using LinearAlgebra
using Statistics
using StaticArrays

@testset "Geometry" begin

for type_under_test in [Vector3, Point3, Normal3]
    typename = string(type_under_test)
    @testset "Common tests: $typename" begin
        X = type_under_test(1.0, 0.0, 0.0)
        Y = type_under_test(0.0, 1.0, 0.0)
        Z = type_under_test(0.0, 0.0, 1.0)
        @testset "Creating" begin
            v = type_under_test(1.0, 2.0, 3.0)
            @test v.x ≈ 1.0
            @test v.y ≈ 2.0
            @test v.z ≈ 3.0
            v = type_under_test(1,2,3)
            @test v.x ≈ 1.0
            @test v.y ≈ 2.0
            @test v.z ≈ 3.0
            v = type_under_test(1//2, 3//4, -2//1)
            @test v.x ≈ 0.5
            @test v.y ≈ 0.75
            @test v.z ≈ -2.0
            v = type_under_test(1.0)
            @test v.x ≈ 1.0
            @test v.y ≈ 1.0
            @test v.x ≈ 1.0
            v = type_under_test(Vector3(1, 2, 3))
            @test v.x ≈ 1.0
            @test v.y ≈ 2.0
            @test v.z ≈ 3.0
        end

        @testset "Arithmetic" begin
            v = type_under_test(1, 2, 3)
            @test 2.0 * v ≈ type_under_test(2, 4, 6)
            @test v * 2.0 ≈ type_under_test(2, 4, 6)
            @test v/2 ≈ type_under_test(0.5, 1.0, 1.5)
        end

        @testset "Normalizing" begin
            v = type_under_test(2, 0, 0)
            @test normalize(v) ≈ type_under_test(1, 0, 0)
            v = type_under_test(0, -2, 0)
            @test normalize(v) ≈ type_under_test(0, -1, 0)
            v = type_under_test(3, -4, 0)
            @test normalize(v) ≈ type_under_test(3//5, -4//5, 0)
        end

        @testset "Dot product" begin
            v = type_under_test(0.5, 1.0, -2.0)
            @test dot(v, X) ≈ 0.5
            @test dot(v, Y) ≈ 1.0
            @test dot(v, Z) ≈ -2.0

            u = type_under_test(2, 1, 1)
            @test dot(v, u) ≈ 0.0
        end

        @testset "Min/max" begin
            v = type_under_test(-1, 0, 2)
            u = type_under_test(0, -1, 0)
            @test min(v,u) ≈ type_under_test(-1, -1, 0)
            @test max(u,v) ≈ type_under_test(0, 0, 2)
        end

        @testset "Iteration" begin
            V = type_under_test(1, 2, 3)
            A = [x for x in V]
            @test A ≈ [1.0, 2.0, 3.0]
        end

        @testset "Cross product" begin
            X = type_under_test(1,0,0)
            Y = type_under_test(0,1,0)
            Z = type_under_test(0,0,1)
            @test cross(X, Y) ≈ Vector3(0, 0, 1)
            @test cross(X, Z) ≈ Vector3(0, -1, 0)
            @test cross(Y, X) ≈ Vector3(0, 0, -1)
            @test cross(Y, Z) ≈ Vector3(1, 0, 0)
            @test cross(Z, X) ≈ Vector3(0, 1, 0)
            @test cross(Z, Y) ≈ Vector3(-1, 0, 0)
        end

        @testset "Various functions" begin
            X = type_under_test(1, 2, 3)
            @test size(X) == 3
            @test length(X) == 3
            @test !isnan(X)
            @test isnan(type_under_test(NaN, 2, 3))
            @test isnan(type_under_test(1, NaN, 3))
            @test isnan(type_under_test(1, 2, NaN))
        end

        @testset "Conversions" begin
            T = type_under_test(1, 2, 3)
            @test convert(Array{Float64, 1}, T) ≈ [1.0, 2.0, 3.0]

            V = Vector3(1, 2, 3)
            P = Point3(1, 4, 9)
            N = Normal3(1, 8, 27)
            @test convert(type_under_test, V) ≈ type_under_test(1, 2, 3)
            @test convert(type_under_test, P) ≈ type_under_test(1, 4, 9)
            @test convert(type_under_test, N) ≈ type_under_test(1, 8, 27)
        end
    end
end

@testset "Type-specific functions" begin
    P1 = Point3(2, 0, -4)
    P2 = Point3(0, 0, -2)
    @test mean(P1, P2) ≈ Point3(1, 0, -3)

    V1 = Vector3(1, 2, 3)
    V2 = Vector3(1, 0, -1)
    @test V1 + V2 ≈ Vector3(2, 2, 2)
    @test V1 - V2 ≈ Vector3(0, 2, 4)
end

@testset "Transformations" begin
    @testset "Identity transformation" begin
        T = Transformation()
        @test T.M ≈ Matrix{Float64}(I, 4, 4)
        @test T.MInv ≈ Matrix{Float64}(I, 4, 4)
    end
    @testset "Translation" begin
        point = Point3(0)
        T = translation(1, 0, 0)
        @test T(point) ≈ Point3(1.0, 0.0, 0.0)
        U = inv(T)
        @test U(T(point)) ≈ point

        point = Point3(-1.0, 0.0, 2.0)
        T = translation(1, 2, 5)
        @test T(point) ≈ Point3(0.0, 2.0, 7.0)
        U = inv(T)
        @test U(T(point)) ≈ point
    end
    @testset "Rotation" begin
        X = Vector3(1, 0, 0)
        Y = Vector3(0, 1, 0)
        Z = Vector3(0, 0, 1)

        point = Point3(1,0,0)
        T = rotation(X, pi/2)
        @test T(point) ≈ Point3(1.0, 0.0, 0.0)
        U = inv(T)
        @test U(T(point)) ≈ point

        point = Point3(1.0, 0.0, 0.0)
        T = rotation(Z, pi/2)
        @test isapprox(T(point), Point3(0.0, 1.0, 0.0) ; atol=1e-16)
        U = inv(T)
        @test U(T(point)) ≈ point
        T = rotation(Y, pi/2)
        @test isapprox(T(point), Point3(0.0, 0.0, -1.0) ; atol=1e-16)
        U = inv(T)
        @test U(T(point)) ≈ point
    end
    @testset "Scaling" begin
        P = Point3(-1, 1, 2)
        T = scaling(2, 1, 1)
        @test T(P) ≈ Point3(-2, 1, 2)
        U = inv(T)
        @test U(T(P)) ≈ P
        T = scaling(1, 2, 1)
        @test T(P) ≈ Point3(-1, 2, 2)
        U = inv(T)
        @test U(T(P)) ≈ P
        T = scaling(1, 1, 2)
        @test T(P) ≈ Point3(-1, 1, 4)
        U = inv(T)
        @test U(T(P)) ≈ P

    end
    @testset "Miscellaneous" begin
        v = Vector3(1/sqrt(3), 1/sqrt(3), 1/sqrt(3))
        for direction in [X_AXIS, Y_AXIS, Z_AXIS, v]
            T = GRAYTR.rotate_z_to(direction)
            @test isapprox(T(Z_AXIS), direction ; atol=1e-16)
        end
    end
end

end # testset Geometry

max_stars_repo_name: dronir/GRAYTR

max_issues_repo_name: dronir/GRAYTR

max_forks_repo_name: dronir/GRAYTR