text:
### A Pluto.jl notebook ###
# v0.19.0

using Markdown
using InteractiveUtils

# ╔═╡ 70c1055f-6fd7-4dce-afca-32b7ab8237c0
using Manopt, Manifolds, Colors, Random, PlutoUI

# ╔═╡ 1c3a0022-43f1-4f2f-a4c5-9dfbe0d9437e
using LinearAlgebra: svd

# ╔═╡ 09eb0cd6-8c08-11ec-3c61-a1da0a6e4937
md"""
# Geodesic Regression

Geodesic regression generalizes [linear regression](https://en.wikipedia.org/wiki/Linear_regression)
to Riemannian manifolds. Let's first phrase it informally as follows:

> For given data points ``d_1,\ldots,d_n`` on a Riemannian manifold ``\mathcal M`` find
> the geodesic that “best explains” the data.

The meaning of “best explain” has still to be clarified. We distinguish two cases: time labelled data and unlabelled data
"""

# ╔═╡ c867bd9c-da61-49a0-b6e6-c5fd12536b55
md"We define some colors from [Paul Tol](https://personal.sron.nl/~pault/)."

# ╔═╡ 39ac3269-1cc9-4eaa-9340-2e7e37648dc8
begin
    black = RGBA{Float64}(colorant"#000000")
    TolVibrantOrange = RGBA{Float64}(colorant"#EE7733")
    TolVibrantBlue = RGBA{Float64}(colorant"#0077BB")
    TolVibrantTeal = RGBA{Float64}(colorant"#009988")
    TolVibrantMagenta = RGBA{Float64}(colorant"#EE3377")
    TolVibrantCyan = RGBA{Float64}(colorant"#33BBEE")
    Random.seed!(42)
end;

# ╔═╡ c901b5ab-fdc7-43e9-b728-5dbd3b21104a
md"""
We use the following data, where you can choose one data point, that you wakt to highlight.
"""

# ╔═╡ 501378b8-ea0b-467f-a67c-d804de9b8d01
begin
    n = 7
    highlighted = 4
    (highlighted > n - 1) && error(
        "Please choose a highlighted point from {1,...,$(n-1)} – you set it to 		highlighted.",
    )
    σ = π / 8
    S = Sphere(2)
    base = 1 / sqrt(2) * [1.0, 0.0, 1.0]
    dir = [-0.75, 0.5, 0.75]
    data_orig = [exp(S, base, dir, t) for t in range(-0.5, 0.5; length=n)]
    # add noise to the points on the geodesic
    data = map(x -> exp(S, x, random_tangent(S, x, :Gaussian, σ)), data_orig)
    localpath = join(splitpath(@__FILE__)[1:(end - 1)], "/") # files folder
    image_prefix = localpath * "/regression"
    @info image_prefix
    render_asy = false # on CI or when you do not have asymptote, this should be false
end;

# ╔═╡ 75b3aff8-842e-407d-84e4-6373129911f3
md"""
This notebook loads data if `render_asymptote` is set to false, and renders the images (using [Asymptote](https://asymptote.sourceforge.io)) otherwise.
"""

# ╔═╡ 115d380c-2729-46f1-82bf-32a318c4d8aa
render_asy && asymptote_export_S2_signals(
    image_prefix * "/regression_data.asy";
    points=[data],
    colors=Dict(:points => [TolVibrantBlue]),
    dot_size=3.5,
    camera_position=(1.0, 0.5, 0.5),
);

# ╔═╡ a8184bf9-133b-4952-b11d-abaaf0ee39d8
render_asy && render_asymptote(image_prefix * "/regression_data.asy"; render=2);

# ╔═╡ 24337e9c-b5c3-4aed-9043-2e7eaca8f984
PlutoUI.LocalResource(image_prefix * "/regression_data.png")

# ╔═╡ dd4a7653-edad-4400-9e29-59b9cca72679
md"""
## Time labeled data
if for each data item $d_i$ we are also given a time point $t_i\in\mathbb R$, which are pairwise different.

Then we can use the least squares error to state the objetive function as [^Fletcher2013]

```math
F(p,X) = \frac{1}{2}\sum_{i=1}^n d_{\mathcal M}^2(γ_{p,X}(t_i), d_i),
```

where ``d_{\mathcal M}`` is the Riemannian distance and ``γ_{p,X}`` is the geodesic
with ``γ(0) = p`` and ``\dot\gamma(0) = X``.

For the real-valued case ``\mathcal M = \mathbb R^m`` the solution ``(p^*, X^*)`` is given in closed form
as follows: with ``d^* = \frac{1}{n}\displaystyle\sum_{i=1}^{n}d_i`` and ``t^* = \frac{1}{n}\displaystyle\sum_{i=1}^n t_i``
we get

```math
 X^* = \frac{\sum_{i=1}^n (d_i-d^*)(t-t^*)}{\sum_{i=1}^n (t_i-t^*)^2}
\quad\text{ and }\quad
p^* = d^* - t^*X^*
```

and hence the linear regression result is the line ``γ_{p^*,x^*}(t) = p^* + tX^*``.

On a Riemannian manifold we can phrase this as an optimization problem on the [tangent bundle](https://en.wikipedia.org/wiki/Tangent_bundle),
i.e. the disjoiint union of all tangent spaces, as

```math
\operatorname*{arg\,min}_{(p,X) \in \mathrm{T}\mathcal M} F(p,X)
```

Due to linearity, the gradient of ``F(p,X)`` is the sum of the single gradients of

```math
 \frac{1}{2}d_{\mathcal M}^2\bigl(γ_{p,X}(t_i),d_i\bigr)
 = \frac{1}{2}d_{\mathcal M}^2\bigl(\exp_p(t_iX),d_i\bigr)
 ,\quad i∈\{1,\ldots,n\}
```

which can be computed using a chain rule of the squared distance and the exponential map,
see for example [^BergmannGousenbourger2018] for details or Equations (7) and (8) of [^Fletcher2013]:
"""

# ╔═╡ 3d4bd7d0-876b-4915-88a4-49b188df16b5
M = TangentBundle(S)

# ╔═╡ 92b18619-f806-412f-a125-52a7731a1389
begin
    struct RegressionCost{T,S}
        data::T
        times::S
    end
    RegressionCost(data::T, times::S) where {T,S} = RegressionCost{T,S}(data, times)
    function (a::RegressionCost)(M, x)
        pts = [geodesic(M.manifold, x[M, :point], x[M, :vector], ti) for ti in a.times]
        return 1 / 2 * sum(distance.(Ref(M.manifold), pts, a.data) .^ 2)
    end
end;

# ╔═╡ 5bd7dbf5-90cb-4e7b-8a12-eb1876cda2f2
begin
    struct RegressionGradient!{T,S}
        data::T
        times::S
    end
    function RegressionGradient!(data::T, times::S) where {T,S}
        return RegressionGradient!{T,S}(data, times)
    end
    function (a::RegressionGradient!)(M, Y, x)
        pts = [geodesic(M.manifold, x[M, :point], x[M, :vector], ti) for ti in a.times]
        gradients = grad_distance.(Ref(M.manifold), a.data, pts)
        Y[M, :point] .= sum(
            adjoint_differential_exp_basepoint.(
                Ref(M.manifold),
                Ref(x[M, :point]),
                [ti * x[M, :vector] for ti in a.times],
                gradients,
            ),
        )
        Y[M, :vector] .= sum(
            adjoint_differential_exp_argument.(
                Ref(M.manifold),
                Ref(x[M, :point]),
                [ti * x[M, :vector] for ti in a.times],
                gradients,
            ),
        )
        return Y
    end
end;

# ╔═╡ 876e011d-fc77-4fe4-a61b-1071b29d7d81
md"""
Now we need just a start point.

For the Euclidean case, the result is given by the first principal component of a principal component analysis,
see [PCR](https://en.wikipedia.org/wiki/Principal_component_regression), i.e. with ``p^* = \frac{1}{n}\displaystyle\sum_{i=1}^n d_i``
the direction ``X^*`` is obtained by defining the zero mean data matrix

```math
D = \bigl(d_1-p^*, \ldots, d_n-p^*\bigr) \in \mathbb R^{m,n}
```

and taking ``X^*`` as an eigenvector to the larges eigenvalue of ``D^{\mathrm{T}}D``.

We can do something similar, when considering the tangent space at the (Riemannian) mean
of the data and then do a PCA on the coordinate coefficients with respect to a basis.
"""

# ╔═╡ 39015e8f-1fcd-4dc1-87d4-15be5ea9f8fa
begin
    m = mean(S, data)
    A = hcat(
        map(x -> get_coordinates(S, m, log(S, m, x), DefaultOrthonormalBasis()), data)...
    )
    pca1 = get_vector(S, m, svd(A).U[:, 1], DefaultOrthonormalBasis())
    x0 = ProductRepr(m, pca1)
end

# ╔═╡ 9302c3ff-f7ac-4104-b80e-9678a0d68ab8
md"""
The optimal “time labels” are then just the projections ``t_i = ⟨d_i,X^*⟩``, ``i=1,\ldots,n``.
"""

# ╔═╡ ff15718a-5edb-4ccc-8782-597b261be970
t = map(d -> inner(S, m, pca1, log(S, m, d)), data)

# ╔═╡ f6deecdd-31b9-493d-a5ed-d9e62cfea9f7
md"""
And we can call the gradient descent. Note that since `gradF!` works in place of `Y`, we have to set the
`evalutation` type accordingly.
"""

# ╔═╡ bada1b3c-6a62-4896-aa4a-7f1f52694682
with_terminal() do
    global y = gradient_descent(
        M,
        RegressionCost(data, t),
        RegressionGradient!(data, t),
        x0;
        evaluation=MutatingEvaluation(),
        stepsize=ArmijoLinesearch(1.0, ExponentialRetraction(), 0.99, 0.1),
        stopping_criterion=StopAfterIteration(100) | StopWhenGradientNormLess(1e-8),
        debug=[:Iteration, " | ", :Cost, "\n", :Stop, 50],
    )
end

# ╔═╡ 72f5ffef-6133-440b-842a-566c0eb09377
md"The result looks like"

# ╔═╡ f58e540f-e83e-4a6c-9e21-74bae7b3d69f
begin
    dense_t = range(-0.5, 0.5; length=100)
    geo = geodesic(S, y[M, :point], y[M, :vector], dense_t)
    init_geo = geodesic(S, x0[M, :point], x0[M, :vector], dense_t)
    geo_pts = geodesic(S, y[M, :point], y[M, :vector], t)
    geo_conn_highlighted = shortest_geodesic(
        S, data[highlighted], geo_pts[highlighted], 0.5 .+ dense_t
    )
end

# ╔═╡ 06c520d1-720c-4eb4-97fb-96870aba73a0
render_asy && asymptote_export_S2_signals(
    image_prefix * "/regression_result1.asy";
    points=[data, [y[M, :point]], geo_pts],
    curves=[init_geo, geo, geo_conn_highlighted],
    tangent_vectors=[[Tuple([y[M, :point], y[M, :vector]])]],
    colors=Dict(
        :curves => [black, TolVibrantTeal, TolVibrantBlue],
        :points => [TolVibrantBlue, TolVibrantOrange, TolVibrantTeal],
        :tvectors => [TolVibrantOrange],
    ),
    dot_sizes=[3.5, 3.5, 2],
    line_widths=[0.33, 0.66, 0.33, 1.0],
    camera_position=(1.0, 0.5, 0.5),
);

# ╔═╡ ab26bdda-9c18-4df3-9a39-c81423a84c16
render_asy && render_asymptote(image_prefix * "/regression_result1.asy"; render=2);

# ╔═╡ 7c491e8b-44d6-420b-89a1-cbf947a98ffc
PlutoUI.LocalResource(image_prefix * "/regression_result1.png")

# ╔═╡ c1f18ab3-c804-43bd-8aba-b83446c355f0
md"""
In this image, together with the blue data points, you see the geodesic of the initialization in black
(evaluated on ``[-\frac{1}{2},\frac{1}{2}]``),
the final point on the tangent bundle in orange, as well as the resulting regression geodesic in teal,
(on the same interval as the start) as well as small teal points indicating the time points on the geodesic corresponding to the data.
Additionally, a thin blue line indicates the geodesic between a data point and its corresponding data point on the geodesic.
While this would be the closest point in Euclidean space and hence the two directions (along the geodesic vs. to the data point) orthogonal, here we have
"""

# ╔═╡ 1d408cc3-7fff-4b85-a264-0841372b3b2d
inner(
    S,
    geo_pts[highlighted],
    log(S, geo_pts[highlighted], geo_pts[highlighted + 1]),
    log(S, geo_pts[highlighted], data[highlighted]),
)

# ╔═╡ 8a371cb2-101f-448b-b219-61e298761026
md"""
But we also started with one of the best scenarios, i.e. equally spaced points on a geodesic obstructed by noise

this gets worse if you start with less even distributed data
"""

# ╔═╡ 0bf7e944-1d9f-490f-8a86-d9dffc564256
begin
    data2 = [exp(S, base, dir, t) for t in [-0.5, -0.49, -0.48, 0.1, 0.48, 0.49, 0.5]]
    data2 = map(x -> exp(S, x, random_tangent(S, x, :Gaussian, σ / 2)), data2)
    m2 = mean(S, data2)
    A2 = hcat(
        map(x -> get_coordinates(S, m, log(S, m, x), DefaultOrthonormalBasis()), data2)...
    )
    pca2 = get_vector(S, m, svd(A2).U[:, 1], DefaultOrthonormalBasis())
    x1 = ProductRepr(m, pca2)
    t2 = map(d -> inner(S, m2, pca2, log(S, m2, d)), data2)
end

# ╔═╡ 3f584132-46c7-42dd-9e2c-0bf2bbf38eb4
with_terminal() do
    global y2 = gradient_descent(
        M,
        RegressionCost(data2, t2),
        RegressionGradient!(data2, t2),
        x1;
        evaluation=MutatingEvaluation(),
        stepsize=ArmijoLinesearch(1.0, ExponentialRetraction(), 0.990, 0.066),
        stopping_criterion=StopAfterIteration(200) | StopWhenGradientNormLess(1e-8),
        debug=[:Iteration, " | ", :Cost, "\n", :Stop, 50],
    )
end

# ╔═╡ c8badbb1-b999-4398-9967-6966a81bcbc1
begin
    geo2 = geodesic(S, y2[M, :point], y2[M, :vector], dense_t)
    init_geo2 = geodesic(S, x1[M, :point], x1[M, :vector], dense_t)
    geo_pts2 = geodesic(S, y2[M, :point], y2[M, :vector], t2)
    geo_conn_highlighted2 = shortest_geodesic(
        S, data2[highlighted], geo_pts2[highlighted], 0.5 .+ dense_t
    )
end

# ╔═╡ e4d16b82-3d18-4c1a-bb85-db400c5c65f7
render_asy && asymptote_export_S2_signals(
    image_prefix * "/regression_result2.asy";
    points=[data2, [y2[M, :point]], geo_pts2],
    curves=[init_geo2, geo2, geo_conn_highlighted2],
    tangent_vectors=[[Tuple([y2[M, :point], y2[M, :vector]])]],
    colors=Dict(
        :curves => [black, TolVibrantTeal, TolVibrantCyan],
        :points => [TolVibrantBlue, TolVibrantOrange, TolVibrantTeal],
        :tvectors => [TolVibrantOrange],
    ),
    dot_sizes=[3.5, 3.5, 2.5],
    line_widths=[0.33, 0.66, 0.33, 1.0],
    camera_position=(1.0, 0.5, 0.5),
);

# ╔═╡ 837f2bc0-a1c1-4d39-adb3-8a7e262dda7d
render_asy && render_asymptote(image_prefix * "/regression_result2.asy"; render=2);

# ╔═╡ ea577515-04d1-4d80-af32-aa941d5a2efb
PlutoUI.LocalResource(image_prefix * "/regression_result2.png")

# ╔═╡ 4bb8db61-2f35-4c41-af1d-4b053536b291
md"""
## Unlabeled data

If we are not given time points $t_i$, then the optimization problem extends – informally speaking –
to also finding the “best fitting” (in the sense of smallest error).
To formalize, the objective function here reads

```math
F(p, X, t) = \frac{1}{2}\sum_{i=1}^n d_{\mathcal M}^2(γ_{p,X}(t_i), d_i),
```

where ``t = (t_1,\ldots,t_n) \in \mathbb R^n`` is now an additional parameter of the objective function.
We write ``F_1(p, X)`` to refer to the function on the tangent bundle for fixed values of ``t`` (as the one in the last part)
and ``F_2(t)`` for the function ``F(p, X, t)`` as a function in ``t`` with fixed values ``(p, X)``.

For the Euclidean case, there is no neccessity to optimize with respect to ``t``, as we saw
above for the initialisation of the fixed time points.

On a Riemannian manifold this can be stated as a problem on the product manifold ``\mathcal N = \mathrm{T}\mathcal M \times \mathbb R^n``, i.e.

"""

# ╔═╡ d762b0e0-d722-4fab-a46c-647a07853354
N = M × Euclidean(length(t2))

# ╔═╡ 03c99340-7927-4404-ba23-e0226e3f32a4
md"""
```math
  \operatorname*{arg\,min}_{\bigl((p,X),t\bigr)\in\mathcal N} F(p, X, t).
```

In this tutorial we present an approach to solve this using an alternating gradient descent scheme.
To be precise, we define the cost funcion now on the product manifold
"""

# ╔═╡ 10b85d9f-e1ec-46e2-99ce-7e25d1dc7ac0
begin
    struct RegressionCost2{T}
        data::T
    end
    RegressionCost2(data::T) where {T} = RegressionCost2{T}(data)
    function (a::RegressionCost2)(N, x)
        TM = N[1]
        pts = [
            geodesic(TM.manifold, x[N, 1][TM, :point], x[N, 1][TM, :vector], ti) for
            ti in x[N, 2]
        ]
        return 1 / 2 * sum(distance.(Ref(TM.manifold), pts, a.data) .^ 2)
    end
end

# ╔═╡ 498dde67-324c-431e-970d-b200f9835f83
md"""
The gradient in two parts, namely (a) the same gradient as before w.r.t. ``(p,X) ∈ T\mathcal M``
just now with a fixed `t` in mind for the second component of the product manifold ``\mathcal N``
"""

# ╔═╡ 52ee5cdd-4a45-47c6-ba25-482d81cb1b6a
begin
    struct RegressionGradient2a!{T}
        data::T
    end
    RegressionGradient2a!(data::T) where {T} = RegressionGradient2a!{T}(data)
    function (a::RegressionGradient2a!)(N, Y, x)
        TM = N[1]
        p = x[N, 1]
        pts = [geodesic(TM.manifold, p[TM, :point], p[TM, :vector], ti) for ti in x[N, 2]]
        gradients = grad_distance.(Ref(TM.manifold), a.data, pts)
        Y[TM, :point] .= sum(
            adjoint_differential_exp_basepoint.(
                Ref(TM.manifold),
                Ref(p[TM, :point]),
                [ti * p[TM, :vector] for ti in x[N, 2]],
                gradients,
            ),
        )
        Y[TM, :vector] .= sum(
            adjoint_differential_exp_argument.(
                Ref(TM.manifold),
                Ref(p[TM, :point]),
                [ti * p[TM, :vector] for ti in x[N, 2]],
                gradients,
            ),
        )
        return Y
    end
end

# ╔═╡ 9da6b6b0-303e-4865-afdd-aad397495517
md"""
Finally we addionally look for a fixed point ``x=(p,X) ∈ \mathrm{T}\mathcal M`` at
the gradient with respect to ``t∈\mathbb R^n``, i.e. the second component, which is given by

```math
  (\operatorname{grad}F_2(t))_i
  = - ⟨\dot γ_{p,X}(t_i), \log_{γ_{p,X}(t_i)}d_i⟩_{γ_{p,X}(t_i)}, i = 1, \ldots, n.
```
"""

# ╔═╡ 96907f2c-e7e4-4bbe-8a2e-91295fac2a8c
begin
    struct RegressionGradient2b!{T}
        data::T
    end
    RegressionGradient2b!(data::T) where {T} = RegressionGradient2b!{T}(data)
    function (a::RegressionGradient2b!)(N, Y, x)
        TM = N[1]
        p = x[N, 1]
        pts = [geodesic(TM.manifold, p[TM, :point], p[TM, :vector], ti) for ti in x[N, 2]]
        logs = log.(Ref(TM.manifold), pts, a.data)
        pt = map(
            d -> vector_transport_to(TM.manifold, p[TM, :point], p[TM, :vector], d), pts
        )
        Y .= -inner.(Ref(TM.manifold), pts, logs, pt)
        return Y
    end
end

# ╔═╡ 52138cac-404a-41e0-ad6e-a7cb51dbec3e
md"""
We can reuse the computed initial values from before, just that now we are on a product manifold
"""

# ╔═╡ a64b78bb-c226-4c68-af5a-02659a4a2dd9
begin
    x2 = ProductRepr(x1, t2)
    F3 = RegressionCost2(data2)
    gradF3_vector = [RegressionGradient2a!(data2), RegressionGradient2b!(data2)]
end;

# ╔═╡ 27548918-eb14-404c-97a6-d39a1fc87929
with_terminal() do
    global y3 = alternating_gradient_descent(
        N,
        F3,
        gradF3_vector,
        x2;
        evaluation=MutatingEvaluation(),
        debug=[:Iteration, " | ", :Cost, "\n", :Stop, 50],
        stepsize=ArmijoLinesearch(1.0, ExponentialRetraction(), 0.990, 0.066),
        inner_iterations=1,
    )
end
# ╔═╡ 1099f1b0-6304-491b-84c9-6f39e08968a9
begin
    geo3 = geodesic(S, y3[N, 1][M, :point], y3[N, 1][M, :vector], dense_t)
    init_geo3 = geodesic(S, x1[M, :point], x1[M, :vector], dense_t)
    geo_pts3 = geodesic(S, y3[N, 1][M, :point], y3[N, 1][M, :vector], y3[N, 2])
    t3 = y3[N, 2]
    geo_conns = shortest_geodesic.(Ref(S), data2, geo_pts3, Ref(0.5 .+ dense_t))
end;

# ╔═╡ 2d8d4f2f-85cf-49fc-a204-02e28566b582
render_asy && asymptote_export_S2_signals(
    image_prefix * "/regression_result3.asy";
    points=[data2, [y3[N, 1][M, :point]], geo_pts3],
    curves=[init_geo3, geo3, geo_conns...],
    tangent_vectors=[[Tuple([y3[N, 1][M, :point], y3[N, 1][M, :vector]])]],
    colors=Dict(
        :curves =>
            [black, TolVibrantTeal, repeat([TolVibrantCyan], length(geo_conns))...],
        :points => [TolVibrantBlue, TolVibrantOrange, TolVibrantTeal],
        :tvectors => [TolVibrantOrange],
    ),
    dot_sizes=[3.5, 3.5, 2.5],
    line_widths=[0.33, 0.66, repeat([0.33], length(geo_conns))..., 1.0],
    camera_position=(1.0, 0.5, 0.5),
);

# ╔═╡ 695478c0-e8e5-4621-bce4-e76eb1f00835
render_asy && render_asymptote(image_prefix * "/regression_result3.asy"; render=2);

# ╔═╡ 6f72b449-e51d-4b6f-bc2a-cd3631343286
PlutoUI.LocalResource(image_prefix * "/regression_result3.png")

# ╔═╡ 426d965f-04e9-4d8e-b44b-bf8d1b85d84b
md"""
Note that the geodesics from the data to the regression geodesic meet at an nearly orthogonal angle.

**Acknowledgement.** Parts of this tutorial are based on the bachelor thesis of
[Jeremias Arf](https://orcid.org/0000-0003-3765-0130).
"""

# ╔═╡ c54d2401-023f-424d-a9a0-db5683de50de
md"""
## References

[^BergmannGousenbourger2018]:
    > Bergmann, R. and Gousenbourger, P.-Y.: _A variational model for data fitting on manifolds
    > by minimizing the acceleration of a Bézier curve_.
    > Frontiers in Applied Mathematics and Statistics, 2018.
    > doi: [10.3389/fams.2018.00059](https://dx.doi.org/10.3389/fams.2018.00059),
    > arXiv: [1807.10090](https://arxiv.org/abs/1807.10090)
[^Fletcher2013]:
    > Fletcher, P. T., _Geodesic regression and the theory of least squares on Riemannian manifolds_,
    > International Journal of Computer Vision(105), 2, pp. 171–185, 2013.
    > doi: [10.1007/s11263-012-0591-y](https://doi.org/10.1007/s11263-012-0591-y)
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
Manifolds = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e"
Manopt = "0fc0a36d-df90-57f3-8f93-d78a9fc72bb5"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[compat]
Colors = "~0.12.8"
Manifolds = "~0.7.7"
Manopt = "~0.3.20"
PlutoUI = "~0.7.37"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.2"
manifest_format = "2.0"

[[deps.AbstractFFTs]]
deps = ["ChainRulesCore", "LinearAlgebra"]
git-tree-sha1 = "6f1d9bc1c08f9f4a8fa92e3ea3cb50153a1b40d4"
uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c"
version = "1.1.0"

[[deps.AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "8eaf9f1b4921132a4cff3f36a1d9ba923b14a481"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.4"

[[deps.Adapt]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "af92965fb30777147966f58acb05da51c5616b5f"
uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
version = "3.3.3"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[deps.ArnoldiMethod]]
deps = ["LinearAlgebra", "Random", "StaticArrays"]
git-tree-sha1 = "62e51b39331de8911e4a7ff6f5aaf38a5f4cc0ae"
uuid = "ec485272-7323-5ecc-a04f-4719b315124d"
version = "0.2.0"

[[deps.ArrayInterface]]
deps = ["Compat", "IfElse", "LinearAlgebra", "Requires", "SparseArrays", "Static"]
git-tree-sha1 = "6e8fada11bb015ecf9263f64b156f98b546918c7"
uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
version = "5.0.5"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[deps.Calculus]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "f641eb0a4f00c343bbc32346e1217b86f3ce9dad"
uuid = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9"
version = "0.5.1"

[[deps.ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "9950387274246d08af38f6eef8cb5480862a435f"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.14.0"

[[deps.ChangesOfVariables]]
deps = ["ChainRulesCore", "LinearAlgebra", "Test"]
git-tree-sha1 = "bf98fa45a0a4cee295de98d4c1462be26345b9a1"
uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
version = "0.1.2"

[[deps.ColorSchemes]]
deps = ["ColorTypes", "Colors", "FixedPointNumbers", "Random"]
git-tree-sha1 = "12fc73e5e0af68ad3137b886e3f7c1eacfca2640"
uuid = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
version = "3.17.1"

[[deps.ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[deps.Colors]]
deps = ["ColorTypes", "FixedPointNumbers", "Reexport"]
git-tree-sha1 = "417b0ed7b8b838aa6ca0a87aadf1bb9eb111ce40"
uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
version = "0.12.8"

[[deps.Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "96b0bc6c52df76506efc8a441c6cf1adcb1babc4"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.42.0"

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[deps.CovarianceEstimation]]
deps = ["LinearAlgebra", "Statistics", "StatsBase"]
git-tree-sha1 = "a3e070133acab996660d31dcf479ea42849e368f"
uuid = "587fd27a-f159-11e8-2dae-1979310e6154"
version = "0.2.7"

[[deps.DataAPI]]
git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.9.0"

[[deps.DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "3daef5523dd2e769dad2365274f760ff5f282c7d"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.11"

[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[deps.DensityInterface]]
deps = ["InverseFunctions", "Test"]
git-tree-sha1 = "80c3e8639e3353e5d2912fb3a1916b8455e2494b"
uuid = "b429d917-457f-4dbc-8f4c-0cc954292b1d"
version = "0.4.0"

[[deps.Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[deps.Distributions]]
deps = ["ChainRulesCore", "DensityInterface", "FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SparseArrays", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns", "Test"]
git-tree-sha1 = "c43e992f186abaf9965cc45e372f4693b7754b22"
uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
version = "0.25.52"

[[deps.DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "b19534d1895d702889b219c382a6e18010797f0b"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.6"

[[deps.Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[deps.DualNumbers]]
deps = ["Calculus", "NaNMath", "SpecialFunctions"]
git-tree-sha1 = "5837a837389fccf076445fce071c8ddaea35a566"
uuid = "fa6b7ba4-c1ee-5f82-b5fc-ecf0adba8f74"
version = "0.6.8"

[[deps.Einsum]]
deps = ["Compat"]
git-tree-sha1 = "4a6b3eee0161c89700b6c1949feae8b851da5494"
uuid = "b7d42ee7-0b51-5a75-98ca-779d3107e4c0"
version = "0.4.1"

[[deps.FillArrays]]
deps = ["LinearAlgebra", "Random", "SparseArrays", "Statistics"]
git-tree-sha1 = "246621d23d1f43e3b9c368bf3b72b2331a27c286"
uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
version = "0.13.2"

[[deps.FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[deps.Graphs]]
deps = ["ArnoldiMethod", "Compat", "DataStructures", "Distributed", "Inflate", "LinearAlgebra", "Random", "SharedArrays", "SimpleTraits", "SparseArrays", "Statistics"]
git-tree-sha1 = "57c021de207e234108a6f1454003120a1bf350c4"
uuid = "86223c79-3864-5bf0-83f7-82e725a168b6"
version = "1.6.0"

[[deps.HybridArrays]]
deps = ["LinearAlgebra", "Requires", "StaticArrays"]
git-tree-sha1 = "eb6b23460f5544c5d09efae0818b86736cefcd3d"
uuid = "1baab800-613f-4b0a-84e4-9cd3431bfbb9"
version = "0.4.10"

[[deps.HypergeometricFunctions]]
deps = ["DualNumbers", "LinearAlgebra", "SpecialFunctions", "Test"]
git-tree-sha1 = "65e4589030ef3c44d3b90bdc5aac462b4bb05567"
uuid = "34004b35-14d8-5ef3-9330-4cdb6864b03a"
version = "0.3.8"

[[deps.Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[deps.HypertextLiteral]]
git-tree-sha1 = "2b078b5a615c6c0396c77810d92ee8c6f470d238"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.3"

[[deps.IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[deps.IfElse]]
git-tree-sha1 = "debdd00ffef04665ccbb3e150747a77560e8fad1"
uuid = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
version = "0.1.1"

[[deps.Inflate]]
git-tree-sha1 = "f5fc07d4e706b84f72d54eedcc1c13d92fb0871c"
uuid = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9"
version = "0.1.2"

[[deps.InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[deps.InverseFunctions]]
deps = ["Test"]
git-tree-sha1 = "91b5dcf362c5add98049e6c29ee756910b03051d"
uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
version = "0.1.3"

[[deps.IrrationalConstants]]
git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151"
uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
version = "0.1.1"

[[deps.JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.4.1"

[[deps.JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "3c837543ddb02250ef42f4738347454f95079d4e"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.3"

[[deps.Kronecker]]
deps = ["LinearAlgebra", "NamedDims", "SparseArrays", "StatsBase"]
git-tree-sha1 = "a51f46415c844dee694cb8b20a3fcbe6dba342c2"
uuid = "2c470bb0-bcc8-11e8-3dad-c9649493f05e"
version = "0.5.1"

[[deps.LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[deps.LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[deps.LogExpFunctions]]
deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"]
git-tree-sha1 = "58f25e56b706f95125dcb796f39e1fb01d913a71"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.3.10"

[[deps.Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[deps.MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "3d3e902b31198a27340d0bf00d6ac452866021cf"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.9"

[[deps.Manifolds]]
deps = ["Colors", "Distributions", "Einsum", "Graphs", "HybridArrays", "Kronecker", "LinearAlgebra", "ManifoldsBase", "Markdown", "Random", "RecipesBase", "RecursiveArrayTools", "Requires", "SimpleWeightedGraphs", "SpecialFunctions", "StaticArrays", "Statistics", "StatsBase"]
git-tree-sha1 = "e8fb7c69d6e67d377152fd882a20334535db050f"
uuid = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e"
version = "0.7.7"

[[deps.ManifoldsBase]]
deps = ["LinearAlgebra", "Markdown"]
git-tree-sha1 = "058ba95cf4a41d4c7b88879f5b961352880ec919"
uuid = "3362f125-f0bb-47a3-aa74-596ffd7ef2fb"
version = "0.12.12"

[[deps.Manopt]]
deps = ["ColorSchemes", "ColorTypes", "Colors", "DataStructures", "Dates", "LinearAlgebra", "ManifoldsBase", "Markdown", "Printf", "Random", "Requires", "SparseArrays", "StaticArrays", "Statistics", "Test"]
git-tree-sha1 = "bcd147170ad2699518155ba9b12ad33f835c5a4d"
uuid = "0fc0a36d-df90-57f3-8f93-d78a9fc72bb5"
version = "0.3.20"

[[deps.Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[deps.Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "bf210ce90b6c9eed32d25dbcae1ebc565df2687f"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.2"

[[deps.Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[deps.NaNMath]]
git-tree-sha1 = "737a5957f387b17e74d4ad2f440eb330b39a62c5"
uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
version = "1.0.0"

[[deps.NamedDims]]
deps = ["AbstractFFTs", "ChainRulesCore", "CovarianceEstimation", "LinearAlgebra", "Pkg", "Requires", "Statistics"]
git-tree-sha1 = "0856b62716585eb90cc1dada226ac9eab5f69aa5"
uuid = "356022a1-0364-5f58-8944-0da4b18d706f"
version = "0.2.47"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

[[deps.OpenLibm_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "05823500-19ac-5b8b-9628-191a04bc5112"

[[deps.OpenSpecFun_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
version = "0.5.5+0"

[[deps.OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[deps.PDMats]]
deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse"]
git-tree-sha1 = "e8185b83b9fc56eb6456200e873ce598ebc7f262"
uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
version = "0.11.7"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "85b5da0fa43588c75bb1ff986493443f821c70b7"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.2.3"

[[deps.Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[deps.PlutoUI]]
deps = ["AbstractPlutoDingetjes", "Base64", "ColorTypes", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "bf0a1121af131d9974241ba53f601211e9303a9e"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.37"

[[deps.Preferences]]
deps = ["TOML"]
git-tree-sha1 = "d3538e7f8a790dc8903519090857ef8e1283eecd"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.5"

[[deps.Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[deps.QuadGK]]
deps = ["DataStructures", "LinearAlgebra"]
git-tree-sha1 = "78aadffb3efd2155af139781b8a8df1ef279ea39"
uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
version = "2.4.2"

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[deps.Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[deps.RecipesBase]]
git-tree-sha1 = "6bf3f380ff52ce0832ddd3a2a7b9538ed1bcca7d"
uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
version = "1.2.1"

[[deps.RecursiveArrayTools]]
deps = ["Adapt", "ArrayInterface", "ChainRulesCore", "DocStringExtensions", "FillArrays", "LinearAlgebra", "RecipesBase", "Requires", "StaticArrays", "Statistics", "ZygoteRules"]
git-tree-sha1 = "f5dd036acee4462949cc10c55544cc2bee2545d6"
uuid = "731186ca-8d62-57ce-b412-fbd966d074cd"
version = "2.25.1"

[[deps.Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[deps.Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.3.0"

[[deps.Rmath]]
deps = ["Random", "Rmath_jll"]
git-tree-sha1 = "bf3188feca147ce108c76ad82c2792c57abe7b1f"
uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa"
version = "0.7.0"

[[deps.Rmath_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "68db32dff12bb6127bac73c209881191bf0efbb7"
uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f"
version = "0.3.0+0"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[deps.SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[deps.SimpleTraits]]
deps = ["InteractiveUtils", "MacroTools"]
git-tree-sha1 = "5d7e3f4e11935503d3ecaf7186eac40602e7d231"
uuid = "699a6c99-e7fa-54fc-8d76-47d257e15c1d"
version = "0.9.4"

[[deps.SimpleWeightedGraphs]]
deps = ["Graphs", "LinearAlgebra", "Markdown", "SparseArrays", "Test"]
git-tree-sha1 = "a6f404cc44d3d3b28c793ec0eb59af709d827e4e"
uuid = "47aef6b3-ad0c-573a-a1e2-d07658019622"
version = "1.2.1"

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

[[deps.SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[deps.SpecialFunctions]]
deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
git-tree-sha1 = "5ba658aeecaaf96923dce0da9e703bd1fe7666f9"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "2.1.4"

[[deps.Static]]
deps = ["IfElse"]
git-tree-sha1 = "87e9954dfa33fd145694e42337bdd3d5b07021a6"
uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
version = "0.6.0"

[[deps.StaticArrays]]
deps = ["LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "4f6ec5d99a28e1a749559ef7dd518663c5eca3d5"
uuid = "90137ffa-7385-5640-81b9-e52037218182"
version = "1.4.3"

[[deps.Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[deps.StatsAPI]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "c3d8ba7f3fa0625b062b82853a7d5229cb728b6b"
uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
version = "1.2.1"

[[deps.StatsBase]]
deps = ["DataAPI", "DataStructures", "LinearAlgebra", "LogExpFunctions", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"]
git-tree-sha1 = "8977b17906b0a1cc74ab2e3a05faa16cf08a8291"
uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
version = "0.33.16"

[[deps.StatsFuns]]
deps = ["ChainRulesCore", "HypergeometricFunctions", "InverseFunctions", "IrrationalConstants", "LogExpFunctions", "Reexport", "Rmath", "SpecialFunctions"]
git-tree-sha1 = "25405d7016a47cf2bd6cd91e66f4de437fd54a07"
uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
version = "0.9.16"

[[deps.SuiteSparse]]
deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"

[[deps.TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[deps.Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[deps.Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[deps.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[deps.Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[deps.Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[deps.ZygoteRules]]
deps = ["MacroTools"]
git-tree-sha1 = "8c1a8e4dfacb1fd631745552c8db35d0deb09ea0"
uuid = "700de1a5-db45-46bc-99cf-38207098b444"
version = "0.2.2"

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─09eb0cd6-8c08-11ec-3c61-a1da0a6e4937
# ╠═70c1055f-6fd7-4dce-afca-32b7ab8237c0
# ╠═1c3a0022-43f1-4f2f-a4c5-9dfbe0d9437e
# ╟─c867bd9c-da61-49a0-b6e6-c5fd12536b55
# ╠═39ac3269-1cc9-4eaa-9340-2e7e37648dc8
# ╟─c901b5ab-fdc7-43e9-b728-5dbd3b21104a
# ╠═501378b8-ea0b-467f-a67c-d804de9b8d01
# ╟─75b3aff8-842e-407d-84e4-6373129911f3
# ╠═115d380c-2729-46f1-82bf-32a318c4d8aa
# ╠═a8184bf9-133b-4952-b11d-abaaf0ee39d8
# ╠═24337e9c-b5c3-4aed-9043-2e7eaca8f984
# ╟─dd4a7653-edad-4400-9e29-59b9cca72679
# ╠═3d4bd7d0-876b-4915-88a4-49b188df16b5
# ╠═92b18619-f806-412f-a125-52a7731a1389
# ╠═5bd7dbf5-90cb-4e7b-8a12-eb1876cda2f2
# ╟─876e011d-fc77-4fe4-a61b-1071b29d7d81
# ╠═39015e8f-1fcd-4dc1-87d4-15be5ea9f8fa
# ╟─9302c3ff-f7ac-4104-b80e-9678a0d68ab8
# ╠═ff15718a-5edb-4ccc-8782-597b261be970
# ╟─f6deecdd-31b9-493d-a5ed-d9e62cfea9f7
# ╠═bada1b3c-6a62-4896-aa4a-7f1f52694682
# ╟─72f5ffef-6133-440b-842a-566c0eb09377
# ╠═f58e540f-e83e-4a6c-9e21-74bae7b3d69f
# ╠═06c520d1-720c-4eb4-97fb-96870aba73a0
# ╠═ab26bdda-9c18-4df3-9a39-c81423a84c16
# ╠═7c491e8b-44d6-420b-89a1-cbf947a98ffc
# ╟─c1f18ab3-c804-43bd-8aba-b83446c355f0
# ╠═1d408cc3-7fff-4b85-a264-0841372b3b2d
# ╟─8a371cb2-101f-448b-b219-61e298761026
# ╠═0bf7e944-1d9f-490f-8a86-d9dffc564256
# ╠═3f584132-46c7-42dd-9e2c-0bf2bbf38eb4
# ╠═c8badbb1-b999-4398-9967-6966a81bcbc1
# ╠═e4d16b82-3d18-4c1a-bb85-db400c5c65f7
# ╠═837f2bc0-a1c1-4d39-adb3-8a7e262dda7d
# ╠═ea577515-04d1-4d80-af32-aa941d5a2efb
# ╟─4bb8db61-2f35-4c41-af1d-4b053536b291
# ╠═d762b0e0-d722-4fab-a46c-647a07853354
# ╟─03c99340-7927-4404-ba23-e0226e3f32a4
# ╠═10b85d9f-e1ec-46e2-99ce-7e25d1dc7ac0
# ╟─498dde67-324c-431e-970d-b200f9835f83
# ╠═52ee5cdd-4a45-47c6-ba25-482d81cb1b6a
# ╟─9da6b6b0-303e-4865-afdd-aad397495517
# ╠═96907f2c-e7e4-4bbe-8a2e-91295fac2a8c
# ╟─52138cac-404a-41e0-ad6e-a7cb51dbec3e
# ╠═a64b78bb-c226-4c68-af5a-02659a4a2dd9
# ╠═27548918-eb14-404c-97a6-d39a1fc87929
# ╠═1099f1b0-6304-491b-84c9-6f39e08968a9
# ╠═2d8d4f2f-85cf-49fc-a204-02e28566b582
# ╠═695478c0-e8e5-4621-bce4-e76eb1f00835
# ╠═6f72b449-e51d-4b6f-bc2a-cd3631343286
# ╟─426d965f-04e9-4d8e-b44b-bf8d1b85d84b
# ╟─c54d2401-023f-424d-a9a0-db5683de50de
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_path: pluto/GeodesicRegression backup 1.jl

max_issues_repo_path: pluto/GeodesicRegression backup 1.jl

max_forks_repo_path: pluto/GeodesicRegression backup 1.jl


text:
@testset "shrink small array" begin
    x = [1, 2, 3, 4]
    @test x[1] == 1
    @test x[2] == 2
    @test x[3] == 3
    @test x[4] == 4
    @test ccall(:jl_array_size, Int, (Any, UInt), x, 0) == 4
    @test ccall(:jl_array_size, Int, (Any, UInt), x, 1) == 4
    sizehint!(x, 10000)
    @test x[1] == 1
    @test x[2] == 2
    @test x[3] == 3
    @test x[4] == 4
    @test ccall(:jl_array_size, Int, (Any, UInt), x, 0) == 4
    @test ccall(:jl_array_size, Int, (Any, UInt), x, 1) == 10000
    sizehint!(x, 4)
    @test x[1] == 1
    @test x[2] == 2
    @test x[3] == 3
    @test x[4] == 4
    @test ccall(:jl_array_size, Int, (Any, UInt), x, 0) == 4
    @test ccall(:jl_array_size, Int, (Any, UInt), x, 1) == 4

    x = [1, 2, 3, 4]
    @test x[1] == 1
    @test x[2] == 2
    @test x[3] == 3
    @test x[4] == 4
    @test ccall(:jl_array_size, Int, (Any, UInt), x, 0) == 4
    @test ccall(:jl_array_size, Int, (Any, UInt), x, 1) == 4
    sizehint!(x, 1000000)
    @test x[1] == 1
    @test x[2] == 2
    @test x[3] == 3
    @test x[4] == 4
    @test ccall(:jl_array_size, Int, (Any, UInt), x, 0) == 4
    @test ccall(:jl_array_size, Int, (Any, UInt), x, 1) == 1000000
    sizehint!(x, 4)
    @test x[1] == 1
    @test x[2] == 2
    @test x[3] == 3
    @test x[4] == 4
    @test ccall(:jl_array_size, Int, (Any, UInt), x, 0) == 4
    @test ccall(:jl_array_size, Int, (Any, UInt), x, 1) == 4
end

ext: jl

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
# Noncentral hypergeometric distribution
# TODO: this distribution needs clean-up and testing

abstract NoncentralHypergeometric <: DiscreteUnivariateDistribution

### handling support

function insupport(d::NoncentralHypergeometric, x::Real)
    isinteger(x) && minimum(d) <= x <= maximum(d)
end

@distr_support NoncentralHypergeometric max(d.n - d.nf, 0) min(d.ns, d.n)

# Functions

function quantile(d::NoncentralHypergeometric, q::Float64)
    if !(zero(q) <= q <= one(q))
        NaN
    else
        range = support(d)
        if q > 0.5
            q = 1-q
            range = reverse(range)
        end

        qsum, i = 0.0, 0
        while qsum < q
            i += 1
            qsum += pdf(d, range[i])
        end
        range[i]
    end
end

params(d::NoncentralHypergeometric) = (d.ns, d.nf, d.n, d.ω)

## Fisher's noncentral hypergeometric distribution

immutable FisherNoncentralHypergeometric <: NoncentralHypergeometric
    ns::Int    # number of successes in population
    nf::Int    # number of failures in population
    n::Int     # sample size
    ω::Float64 # odds ratio

    function FisherNoncentralHypergeometric(ns::Real, nf::Real, n::Real, ω::Float64)
        @check_args(FisherNoncentralHypergeometric, ns >= zero(ns) && nf >= zero(nf))
        @check_args(FisherNoncentralHypergeometric, zero(n) < n < ns + nf)
        @check_args(FisherNoncentralHypergeometric, ω > zero(ω))
        new(ns, nf, n, ω)
    end
end

# Properties
binomial{T<:Integer}(n::T, k::T) = Base.binomial(BigInt(n), BigInt(k))
binomial{T<:Integer}(n::T, kR::OrdinalRange{T}) = [binomial(n,k) for k in kR]
function _P(d::FisherNoncentralHypergeometric, k::Int)
    y = support(d)
    p = binomial(d.ns, y) .* binomial(d.nf, d.n-y) .* d.ω.^y .* y.^k
    sum(p)
end

function _mode(d::FisherNoncentralHypergeometric)
    A = d.ω - 1
    B = d.n - d.nf - (d.ns + d.n + 2)*d.ω
    C = (d.ns + 1)*(d.n + 1)*d.ω
    -2*C / (B - sqrt(B^2-4*A*C))
end

mean(d::FisherNoncentralHypergeometric) =_P(d,1) / _P(d,0)
var(d::FisherNoncentralHypergeometric) = _P(d,2)/_P(d,0) - (_P(d,1) / _P(d,0))^2
mode(d::FisherNoncentralHypergeometric) = floor(Int, _mode(d))

logpdf(d::FisherNoncentralHypergeometric, k::Int) =
    Float64(log(binomial(d.ns, k)) + log(binomial(d.nf, d.n-k)) + k*log(d.ω) - log(_P(d,0)))

pdf(d::FisherNoncentralHypergeometric, k::Int) = exp(logpdf(d, k))

## Wallenius' noncentral hypergeometric distribution

immutable WalleniusNoncentralHypergeometric <: NoncentralHypergeometric
    ns::Int    # number of successes in population
    nf::Int    # number of failures in population
    n::Int     # sample size
    ω::Float64 # odds ratio
    function WalleniusNoncentralHypergeometric(ns::Real, nf::Real, n::Real, ω::Float64)
        @check_args(WalleniusNoncentralHypergeometric, ns >= zero(ns) && nf >= zero(nf))
        @check_args(WalleniusNoncentralHypergeometric, zero(n) < n < ns + nf)
        @check_args(WalleniusNoncentralHypergeometric, ω > zero(ω))
        new(ns, nf, n, ω)
    end
end

# Properties
mean(d::WalleniusNoncentralHypergeometric)=sum(support(d).*pdf(d,support(d)))
var(d::WalleniusNoncentralHypergeometric) = sum((support(d)-mean(d)).^2.*pdf(d,support(d)))
mode(d::WalleniusNoncentralHypergeometric) = support(d)[indmax(pdf(d,support(d)))]

function pdf(d::WalleniusNoncentralHypergeometric, k::Int)
    D = d.ω*(d.ns-k)+(d.nf-d.n+k)
    f(t) = (1-t^(d.ω/D))^k * (1-t^(1/D))^(d.n-k)
    I,_ = quadgk(f,0,1)
    Float64(binomial(d.ns,k)*binomial(d.nf,d.n-k)*I)
end

logpdf(d::WalleniusNoncentralHypergeometric, k::Int) = log(pdf(d, k))


text:
get_degree(offsets, i::Int) = offsets[i+1] - offsets[i]
get_degree(offsets, is::Vector{Int}) = [get_degree(offsets, i) for i in is]

function get_max_degree(offsets)
    allnode = collect(1:(length(offsets)-1))
    degrees = get_degree(offsets, allnode)
    return max(degrees...)
end

function get_neighbor(adjacency_list, offsets, nodeid)
    degree = get_degree(offsets, nodeid)
    adjacency_list[offsets[nodeid]:offsets[nodeid]+degree-1]
end

function onewalk(adjacency_list, offsets, query_nodeid)
    degree = get_degree(offsets, query_nodeid)
    if degree == 0
        # isolated node
        return -1
    end
    rand_position = floor(Int, rand(Uniform(0, degree)))
    return adjacency_list[offsets[query_nodeid]+rand_position]
end

function randomwalk(
    adjacency_list,
    offsets,
    query_nodeid,
    count_same_nodetype,
    terminate_prob,
    total_walk_length,
    min_high_visited_candidates,
    high_visited_count_threshold,
)::Dict{Int,Int}
    visited_count = Dict{Int,Int}()
    current_length = 0
    n_high_visited_candidates = 0
    while (current_length < total_walk_length) &&
        (n_high_visited_candidates < min_high_visited_candidates)
        walk_length = rand(Geometric(terminate_prob))
        current_nodeid = query_nodeid

        for _ = 1:walk_length
            current_nodeid = onewalk(adjacency_list, offsets, current_nodeid)
            if count_same_nodetype
                current_nodeid = onewalk(adjacency_list, offsets, current_nodeid)
                # TODO: what if -1 returned?
            end

            if current_nodeid in keys(visited_count)
                visited_count[current_nodeid] += 1
            else
                visited_count[current_nodeid] = 1
            end

            if visited_count[current_nodeid] == high_visited_count_threshold
                n_high_visited_candidates += 1
            end

            if !count_same_nodetype
                current_nodeid = onewalk(adjacency_list, offsets, current_nodeid)
                # TODO: what if -1 returned?
            end
        end

        current_length += walk_length
    end

    return visited_count
end

function pixie_multi_hit_boost(visited_counts::Vector{Dict{Int,Int}})::Dict{Int,Real}
    if length(visited_counts) == 1
        allnodes = keys(visited_counts[1])
    else
        allnodes = union(keys.(visited_counts)...)
    end
    boosted_visited_count = Dict{Int,Real}()
    for nodeid in allnodes
        counts = get.(visited_counts, nodeid, 0)
        boosted_visited_count[nodeid] = (sum(sqrt.(counts)))^2
    end
    return boosted_visited_count
end

function aggregate_multi_randomwalk(
    visited_counts::Vector{Dict{Int,Int}},
    aggregate_function = sum,
)::Dict{Int,Int}
    if length(visited_counts) == 1
        allnodes = keys(visited_counts[1])
    else
        allnodes = union(keys.(visited_counts)...)
    end
    total_visited_count = Dict{Int,Int}()
    for nodeid in allnodes
        total_visited_count[nodeid] = aggregate_function(get.(visited_counts, nodeid, 0))
    end
    return total_visited_count
end

function randomwalk_multiple(
    adjacency_list,
    offsets,
    query_nodeids,
    count_same_nodetype,
    terminate_prob,
    total_walk_length,
    min_high_visited_candidates,
    high_visited_count_threshold,
    pixie_walk_length_scaling = false,
    pixie_multi_hit_boosting = false,
    max_degree = nothing,
    aggregate_function = sum,
)::Dict{Int,Real}

    if max_degree === nothing
        max_degree = get_max_degree(offsets)
    end

    degrees = get_degree(offsets, query_nodeids)

    if pixie_walk_length_scaling
        scaling_factors = degrees .* (max_degree .- log.(degrees))
        scaling_factors = scaling_factors / sum(scaling_factors)
    else
        scaling_factors = ones(length(query_nodeids)) / length(query_nodeids)
    end

    @assert length(query_nodeids) == length(scaling_factors)

    visited_counts = Vector{Dict{Int,Int}}(undef, length(query_nodeids))

    for (i, query_nodeid) in enumerate(query_nodeids)
        this_total_walk_length = round(Int, total_walk_length * scaling_factors[i])
        visited_counts[i] = randomwalk(
            adjacency_list,
            offsets,
            query_nodeid,
            count_same_nodetype,
            terminate_prob,
            this_total_walk_length,
            min_high_visited_candidates,
            high_visited_count_threshold,
        )
    end

    if pixie_multi_hit_boosting
        return pixie_multi_hit_boost(visited_counts)
    else
        return aggregate_multi_randomwalk(visited_counts, aggregate_function)
    end
end

function build_graph(table; col_user = :userid, col_item = :itemid)
    table, user2uidx, item2iidx, _ =
        make_idmap(table, col_user = col_user, col_item = col_item)

    n_user = length(user2uidx)
    n_item = length(item2iidx)

    df = DataFrame(table)

    # nodeid = 1 - n_item: item node
    # nodeid = (n_item + 1) - (n_item + n_user): user node
    df[!, col_user] = df[!, col_user] .+ n_item
    for user in keys(user2uidx)
        user2uidx[user] = user2uidx[user] + n_item
    end

    adjacency_list = sort(df, [col_item])[!, col_user]
    append!(adjacency_list, sort(df, [col_user])[!, col_item])

    item_degrees = sort(combine(groupby(df, col_item), nrow), [col_item])[!, :nrow]
    user_degrees = sort(combine(groupby(df, col_user), nrow), [col_user])[!, :nrow]

    @assert length(adjacency_list) == sum(item_degrees) + sum(user_degrees)

    offsets = Vector{Int}(undef, n_item + n_user + 1)
    item_cumsum = cumsum(item_degrees) .+ 1
    user_cumsum = cumsum(user_degrees) .+ item_cumsum[end]
    offsets[1] = 1
    offsets[2:(n_item+1)] = item_cumsum
    offsets[(n_item+2):end] = user_cumsum

    return adjacency_list, offsets, user2uidx, item2iidx
end


text:
using Expr2LaTeX
using Base.Test

include("julia2math.jl")


text:
# deprecates related to the version up from WORLD v0.1.4 to v0.2.0

# Since `libworld.0.2.0` doesn't provide an interface to Platinum,
# I decided `platinum` simply throws an error for deprecation.
function platinum(x, fs, timeaxis, f0, spectrogram)
    error("Platinum was removed since WORLD v0.2.0. Please use d4c for aperiodicity (not excitation) estimation.")
end

@deprecate star cheaptrick
@deprecate aperiodicityratio(x, fs, f0, timeaxis) d4c(x, fs, timeaxis, f0)
@deprecate synthesis_from_aperiodicity synthesis
@deprecate get_fftsize_for_star(fs) get_fftsize_for_cheaptrick(fs)


text:
function logmsg(msg...; doflush=false)
    println("[",now(), "]", msg...)
    doflush && flush(stdout)
end

function path_prefix(workspace)
    try
        abspath(joinpath(readchomp(`ws_find $workspace`)))
    catch e
        path_prefix()
    end
end

path_prefix() = abspath(joinpath(pwd(), "data"))

const GEOMETRIES = [:box, :box_pbc, :noisy_chain, :noisy_chain_pbc]

# The OBC Volume coeffs should probably be a bit bigger to accommodate for the extra space outside
# Right now results won't really be comparable, but we will not need this (right now that is)
const VOLUME_COEFFS = (; sphere=[2.0, π, 4/3*π], box=[1.0,1.0,1.0], box_pbc=[1.0,1.0,1.0])

parse_geometry(s::AbstractString) = Symbol(lowercase(s))
parse_geometry(s::Symbol) = s

## General idea:
## set r_bl=1 and scale volume via the density rho
## rho is defined as the ratio (blockaded volume) / (total volume)
## total volume = VOLUME_FACTOR[geom][d]*L**d
## blockaded volume = N*VOLUME_FACTOR[sphere][d]*1**d
length_from_density(geometry, ρ, N, dim) = (1/ρ * N * VOLUME_COEFFS.sphere[dim] / VOLUME_COEFFS[geometry][dim])^(1/dim)

geometry_from_density(geom::Symbol, ρ, N, dim) = geometry_from_density(Val(geom), ρ, N, dim)
geometry_from_density(::Val{:box}, ρ, N, dim) = Box(ones(dim)*length_from_density(:box, ρ, N, dim))
geometry_from_density(::Val{:box_pbc}, ρ, N, dim) = BoxPBC(ones(dim)*length_from_density(:box_pbc, ρ, N, dim))
function geometry_from_density(::Val{:noisy_chain}, ρ, N, dim)
    dim == 1 || error("Chain is only 1D!")
    # know that rho = 2/spacing, since r_bl = 1, and V = N*spacing
    spacing = 2/ρ
    σ = 1.5*(spacing-1) # this ensures ~90% success rate when generating positions -> see demo/rho_scaling.ipynb
    NoisyChain(N, spacing, σ)
end

function geometry_from_density(::Val{:noisy_chain_pbc}, ρ, N, dim)
    dim == 1 || error("Chain is only 1D!")
    # know that rho = 2/spacing, since r_bl = 1, and V = N*spacing
    spacing = 2/ρ
    σ = 1.5*(spacing-1) # this ensures ~90% success rate when generating positions -> see demo/rho_scaling.ipynb
    NoisyChainPBC(N, spacing, σ)
end

# Determined by testing - within these ranges the position sampling converges almost certainly
const SAFE_RHO_RANGES = (;box     = [(0.1, 1.25), (0.1, 1.7), (0.1,2.05)],
                          box_pbc = [(0.1, 1.25), (0.1, 1.7), (0.1,2.05)],
                          noisy_chain_pbc = [(1.0, 1.99)],
                          noisy_chain     = [(1.0, 1.99)])

meandrop(arr; dims) = dropdims(mean(arr; dims); dims)
stddrop(arr; dims) = dropdims(std(arr; dims); dims)


text:
function circleShape(x, y, r)
    θ = range(0, 2π, length=100)
    x .+ r*sin.(θ), y .+ r*cos.(θ)
end

@recipe function f(con::CircleConstraint; circle_color=:red)
    p = length(con)
    for i = 1:p
        circ = circleShape(con.x[i], con.y[i], con.radius[i])
        @series begin
            aspect_ratio --> :equal
            seriestype := :shape
            label --> :none
            linecolor --> circle_color
            linewidth --> 0.5
            fillalpha --> 1
            fillcolor --> circle_color
            primary := false
            circ
        end
    end
    primary := false
    ()
end

@recipe function f(con::GoalConstraint; xind=1, yind=2)
    @series begin
        seriestype := :scatter
        label --> "goal"
        markersize --> 10
        seriescolor --> :green
        [con.xf[xind]],[con.xf[yind]]
    end
end

@recipe function f(con::T) where T <: TrajectoryOptimization.AbstractConstraint
    primary := false
    ()
end

@recipe f(::Type{T}, conval::T) where T <: TrajectoryOptimization.ConVal = conval.con

@recipe function f(conSet::TrajectoryOptimization.AbstractConstraintSet)
    for con in conSet
        @series begin
            con
        end
    end
    aspect_ratio --> :equal
    primary := false
    ()
end

## Plotting trajectories for specific models
@recipe function f(prob::Problem)
    @series begin
        prob.constraints
    end
    (prob.model, prob.Z)
end

@recipe function f(car::RobotZoo.DubinsCar, Z::Traj)
    linecolor --> :black
    linewidth --> 2
    aspect_ratio --> :equal
    xind := 1
    yind := 2
    RobotDynamics.Traj2((Z,))
end

@recipe function f(model::RobotZoo.DoubleIntegrator, Z::Traj)
    label --> ["position" "velocity"]
    linewidth --> 2
    (RobotDynamics.get_times(Z), states(Z))
end


text:
mutable struct qp_controller_input_t <: LCMType
    be_silent::Bool
    timestamp::Int64
    zmp_data::zmp_data_t

    num_support_data::Int32
    support_data::Vector{support_data_t}

    num_tracked_bodies::Int32
    body_motion_data::Vector{body_motion_data_t}

    num_external_wrenches::Int32
    body_wrench_data::Vector{body_wrench_data_t}

    whole_body_data::whole_body_data_t

    num_joint_pd_overrides::Int32
    joint_pd_override::Vector{joint_pd_override_t}

    param_set_name::String
    torque_alpha_filter::Float64
end

@lcmtypesetup(qp_controller_input_t,
    support_data => (num_support_data,),
    body_motion_data => (num_tracked_bodies,),
    body_wrench_data => (num_external_wrenches,),
    joint_pd_override => (num_joint_pd_overrides,),
)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct AdminAppsApi <: SwaggerApi
    client::Swagger.Client
end

function _swaggerinternal_adminAppsApprove(_api::AdminAppsApi; app_id=nothing, request_id=nothing, team_id=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/admin.apps.approve", ["slackAuth"])
    Swagger.set_param(_ctx.form, "app_id", app_id)  # type String
    Swagger.set_param(_ctx.form, "request_id", request_id)  # type String
    Swagger.set_param(_ctx.form, "team_id", team_id)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Approve an app for installation on a workspace.
Param: app_id::String
Param: request_id::String
Param: team_id::String
Return: Dict{String, Any}
"""
function adminAppsApprove(_api::AdminAppsApi; app_id=nothing, request_id=nothing, team_id=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_adminAppsApprove(_api; app_id=app_id, request_id=request_id, team_id=team_id, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function adminAppsApprove(_api::AdminAppsApi, response_stream::Channel; app_id=nothing, request_id=nothing, team_id=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_adminAppsApprove(_api; app_id=app_id, request_id=request_id, team_id=team_id, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_adminAppsRestrict(_api::AdminAppsApi; app_id=nothing, request_id=nothing, team_id=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/admin.apps.restrict", ["slackAuth"])
    Swagger.set_param(_ctx.form, "app_id", app_id)  # type String
    Swagger.set_param(_ctx.form, "request_id", request_id)  # type String
    Swagger.set_param(_ctx.form, "team_id", team_id)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Restrict an app for installation on a workspace.
Param: app_id::String
Param: request_id::String
Param: team_id::String
Return: Dict{String, Any}
"""
function adminAppsRestrict(_api::AdminAppsApi; app_id=nothing, request_id=nothing, team_id=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_adminAppsRestrict(_api; app_id=app_id, request_id=request_id, team_id=team_id, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function adminAppsRestrict(_api::AdminAppsApi, response_stream::Channel; app_id=nothing, request_id=nothing, team_id=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_adminAppsRestrict(_api; app_id=app_id, request_id=request_id, team_id=team_id, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export adminAppsApprove, adminAppsRestrict


text:
# The following errors may be thrown by the parser
const E_EXPECTED_EOF    = "Expected end of input"
const E_UNEXPECTED_EOF  = "Unexpected end of input"
const E_UNEXPECTED_CHAR = "Unexpected character"
const E_BAD_KEY         = "Invalid object key"
const E_BAD_ESCAPE      = "Invalid escape sequence"
const E_BAD_CONTROL     = "ASCII control character in string"
const E_LEADING_ZERO    = "Invalid leading zero in number"
const E_BAD_NUMBER      = "Invalid number"

export E_EXPECTED_EOF, E_UNEXPECTED_EOF, E_UNEXPECTED_CHAR, E_BAD_KEY,
       E_BAD_ESCAPE, E_BAD_CONTROL, E_LEADING_ZERO, E_BAD_NUMBER


text:
# Test transforms.jl
using Test, LinearAlgebra

println("testing transforms")

tol = 1e-10

# rotx, roty, rotz, invrpy
roll = 0.3
pitch = -1.0
yaw = 0.8
T = rotz(yaw)*roty(pitch)*rotx(roll)

(rpy1, rpy2) = invrpy(T)

# Check 1st solution matches roll pitch and yaw
@test abs(rpy1[1]-yaw) < tol
@test abs(rpy1[2]-pitch) < tol
@test abs(rpy1[3]-roll) < tol

# Check 2nd solution produces the same transform
T2 = rotz(rpy2[1])*roty(rpy2[2])*rotx(rpy2[3])
@test maximum(abs.(T2-T)) < tol

# inveuler
phi = 0.5
theta = -0.2
psi = 0.1
T = rotz(phi) * roty(theta) * rotz(psi)
(euler1, euler2) = inveuler(T)

# Check both solutions produce the same transform
T1 = rotz(euler1[1])*roty(euler1[2])*rotz(euler1[3])
@test maximum(abs.(T1-T)) < tol

T2 = rotz(euler2[1])*roty(euler2[2])*rotz(euler2[3])
@test maximum(abs.(T2-T)) < tol

# invht
v = [1,-2,3]
T = rotz(yaw)*roty(pitch)*rotx(roll) * trans(v)

@test maximum(abs.(invht(T)*T - I)) < tol
@test maximum(abs.(trans(v) - trans(v[1],v[2],v[3]))) < tol

# angleaxis, angleaxis2matrix, matrix2angleaxis
# matrix2quaternion, quaternion2matrix
T1 = rotz(yaw)*roty(pitch)*rotx(roll)
ax = matrix2angleaxis(T1)
T2 = angleaxis2matrix(ax)
@test maximum(abs.(T2-T1)) < tol

theta = 0.35
axis = T1[1:3,1]
ax = angleaxis(theta, axis)
Q = quaternion(theta,axis)

T1 = quaternion2matrix(Q)
T2 = angleaxis2matrix(ax)
@test maximum(abs.(T2-T1)) < tol

Q2 = matrix2quaternion(T1)
@test maximum(abs.(Q-Q2)) < tol

# angleaxisrotate, quaternionrotate
ax = angleaxis(theta, [1,0,0])
Q = quaternion(theta, [1,0,0])
v = [3,-2,7,1]
vnew1 = quaternionrotate(Q,v)
vnew2 = angleaxisrotate(ax,v)
vnew3 = rotx(theta)*v

@test maximum(abs.(vnew3-vnew1)) < tol
@test maximum(abs.(vnew3-vnew2)) < tol

# quaternionconjugate, quaternionproduct,
# vector2quaternion


text:
import Distributions
import SpecialFunctions

"""    gaussianMechConstant(ϵ::Real,δ::Real)
Compute the proportionality constant κ(ϵ,δ) for the Gaussian mechanism.
"""
function gaussianMechConstant(ϵ::Real, δ::Real)
    K = sqrt(2) * SpecialFunctions.erfinv(1-2δ)
    return (K+sqrt(K^2+2ϵ))/(2ϵ)
end

"""    (k₁, k₂) = gaussianMechConstant2(ϵ::Real,δ::Real)
Compute the proportionality constant κ(ϵ,δ) for the Gaussian mechanism,
both with the formula using the inverse Q-function (k₁), and according to the
approximation of Theorem A.1 in Dwork and Roth's book (k₂, higher value).
"""
function gaussianMechConstant2(ϵ::Real, δ::Real)
    K = sqrt(2) * SpecialFunctions.erfinv(1-2δ)
    return ((K+sqrt(K^2+2ϵ))/(2ϵ), sqrt(2log(1.25/δ))/ϵ)
end

"""    mean_dp(x::Vector,epsilon::Float64,delta::Float64=0.0)
Compute the average of a set of numbers contained in the vector x,
in an (ϵ,δ)-differentially private way. It is assumed that each entry
x_i of the vector is in [0,1], and the adjacency relation looks at
arbitrary variations in one entry x_i, within this interval.
"""
function mean_dp(x::Vector, ϵ::Real, δ::Real=0.0)
    if δ == 0.0
        d = Distributions.Laplace(0, 1/(ϵ*length(x)))
    else
        d = Distributions.Normal(0, gaussianMechConstant(ϵ,δ)/length(x))
    end
    return Distributions.mean(x) + Distributions.rand(d)
end

"""    laplaceMech(x,f,l1sens::Real,ϵ::Real)
Computes a randomized version to f(x) according to the
Laplace mechanism, for ϵ-differential privacy. l1sens is
the l1-sensitivity of f, which must be computed and provided by the user.
f must take values in R^k, for some k, i.e., return an array of k real values.
"""
function laplaceMech(x, f, l1sens::Real, ϵ::Real)
    t = f(x)
    d = Distributions.Laplace(0, l1sens/ϵ)
    return (t .+ Distributions.rand(d, length(t)))
end

"""    gaussianMech(x,f,l2sens::Real,ϵ::Real,δ::Real)
Computes a randomized version to f(x) according to the
Gaussian mechanism, for (ϵ,δ)-differential privacy. l2sens is
the l2-sensitivity of f, which must be computed by the user.
f must take values in R^k, from some k, i.e., return an array
of k real values.
"""
function gaussianMech(x, f, l2sens::Real, ϵ::Real, δ::Real)
    t = f(x)
    d = Distributions.Normal(0, gaussianMechConstant(ϵ, δ) * l2sens)
    return (t .+ Distributions.rand(d, length(t)))
end

"""    truncatedLaplaceMech(x,f,l1sens::Real,ϵ::Real,δ::Real)
Computes a randomized version to f(x) according to the
truncated Laplace mechanism, for (ϵ,δ)-differential privacy.
This scheme adds bounded noise.
l1sens is the l1-sensitivity of f, which must be computed and provided by the
user. f must take values in R^k, for some k, i.e., return an array of k real
values.

Returns: (r, a) where r is a noisy version of f(x), and a defines the
support of the noise distribution (in the interval [f(x)-a, f(x)+a])
"""
function truncatedLaplaceMech(x, f, l1sens::Real, ϵ::Real, δ::Real)
    t = f(x)
    k = length(t)
    λ = l1sens/ϵ
    a = λ * log( 1 + exp(ϵ) * (k*(1-exp(-ϵ/k)))/(2*δ) )
    # alternative, more conservative but indpt of k
    # a = λ * log(1 + ϵ * exp(ϵ) / (2*δ))
    d1 = Distributions.Laplace(0, λ)
    d = Distributions.truncated(d1, -a, a)
    return (t .+ Distributions.rand(d, length(t)), a)
end


text:
"""
    Eddy-Diffusivity-Mass-Flux (EDMF)

Module for solving the sub-grid-scale equations
using the Eddy-Diffusivity-Mass-Flux (EDMF) model.
"""
module EDMF

using CLIMAParameters
using CLIMAParameters.Planet
using Thermodynamics
using UnPack
using LinearAlgebra
const TD = Thermodynamics

using ..FiniteDifferenceGrids
using ..StateVecs
using ..DomainDecompositions

include("cases.jl")
include("turb_conv_models.jl")
include("aux_funcs.jl")
include("init_params.jl")
include("updraft_vars.jl")
include("update_forcing.jl")
include("ref_states.jl")
include("initial_conditions.jl")
include("update_aux.jl")
include("to_be_depricated.jl")
include("sum_tendencies.jl")
include("apply_bcs.jl")
include("edmf_funcs.jl")
include("update_surface.jl")
include("solve.jl")
include("main.jl")

end # module


text:
using ReachabilityAnalysis, LinearAlgebra, SparseArrays

function tline(;η=3, R=1.00, Rd=10.0, L=1e-10, C=1e-13 * 4.00)
    A₁₁ = zeros(η, η)
    A₁₂ = Bidiagonal(fill(-1/C, η), fill(1/C, η-1), :U)
    A₂₁ = Bidiagonal(fill(1/L, η), fill(-1/L, η-1), :L)
    A₂₂ = Diagonal(vcat(-Rd/L, fill(-R/L, η-1)))
    A  = [A₁₁ A₁₂; A₂₁ A₂₂]
    B = sparse([η+1], [1], 1/L, 2η, 1)
    return A, B
end

using Plots

A, _ = tline(η=20)
spy(A, legend=nothing, markersize=2.0, title="Sparsity pattern of A", xlab="columns", ylab="rows")

function scale!(s, α=1.0)
    s.A .*= α
    s.B .*= α
    return s
end

η = 20 # order
n = 2η # state-space dimension
A, B = tline(η=η);

Uin_ss = Interval(-0.2, 0.2)
□(ϵ) = BallInf(zeros(n), ϵ)
X0 = -inv(Matrix(A)) * B * Uin_ss ⊕ □(0.001);

Uin = Interval(0.99, 1.01)
s = @system(x'= A*x + B*u, x ∈ Universe(n), u ∈ Uin)
α = 1e-9 # scaling factor
scale!(s, α);

T = 0.7 * 1e-11 # time horizon
P = InitialValueProblem(s, X0);

sol = solve(P, T=0.7, alg=BOX(δ=1e-3));

Uout_vs_t = @. (-1.0) * project(sol, η);

plot(Uout_vs_t, vars=(0, η), color=:blue, xlab="t", ylab="Uout", alpha=.5, lw=0.5)


text:
module Nanomsg

using Sockets

export Socket, CSymbols, jl_nn_errno_check, poll
export _nn_errno, _nn_strerror, _nn_symbol_info, _NNSymbolProperties

const LIB = @static Sys.iswindows() ? "nanomsg.dll" : "libnanomsg"

mutable struct NanomsgError <: Exception
    context::AbstractString
    errno::Cint

    NanomsgError(m) = new(m,0)
    NanomsgError(m,e) = new(m,e)
end

function Base.show(io::IO, err::NanomsgError)
	print(io, (err.errno == 0 ? "Nanomsg.jl error, " : "Nanomsg.lib error, errno=$(err.errno), message=$(nn_strerror(err.errno)), "), err.context)
end

mutable struct Socket
    s::Cint

    function Socket(domain::Integer, protocol::Integer)
        p = _nn_socket(domain, protocol)
        if p == -1
            throw(NanomsgError("Socket creation failed", _nn_errno()))
        end

        # 0.8-beta max recv size
        maxSize = convert(Ptr{Nothing}, -1)
        size = convert(Csize_t, sizeof(maxSize))
        _nn_setsockopt(p, CSymbols.NN_SOL_SOCKET, CSymbols.NN_RCVMAXSIZE, maxSize, size)

        socket = new(p)
        finalizer(close, socket)
        return socket
    end
end

function Base.close(socket::Socket)
    if socket.s != -1
        rc = _nn_close(socket.s)
        if rc == -1
            throw(NanomsgError("Socket close failed", _nn_errno()))
        end
        socket.s = -1
    end
	nothing
end

function Base.bind(socket::Socket, endpoint::AbstractString)
    rc = _nn_bind(socket.s, pointer(endpoint))
    if rc == -1
        throw(NanomsgError("Socket bind failed", _nn_errno()))
    end
    nothing
end

function Sockets.connect(socket::Socket, endpoint::AbstractString)
    rc = _nn_connect(socket.s, pointer(endpoint))
    if rc == -1
        throw(NanomsgError("Socket connect failed", _nn_errno()))
    end
    nothing
end

# TODO: Use a in memory IO in addition to the helper string methods
function Sockets.send(socket::Socket, msg::AbstractString, flags::Integer = CSymbols.NN_DONTWAIT)
	size = convert(Csize_t, length(msg.data))
	_send(socket, pointer(msg), size, flags)
end

function Sockets.send(socket::Socket, msg::Array{UInt8}, flags::Integer = CSymbols.NN_DONTWAIT)
	size = convert(Csize_t, length(msg))
	_send(socket, pointer(msg), size, flags)
end

function _send(socket::Socket, msg::Ptr{UInt8}, size::Csize_t, flags::Integer = CSymbols.NN_DONTWAIT)
    rc = _nn_send(socket.s, convert(Ptr{Nothing}, msg), size, flags)
    if rc == -1
    	err = _nn_errno()
    	if err == CSymbols.EAGAIN
    		return nothing
    	end
        throw(NanomsgError("Socket send failed", err))
    end

    if size != rc
    	throw(NanomsgError("Socket sent bytes $rc != $size failed"))
    end

    # println("sent:", rc)
    return rc
end

function Sockets.recv(socket::Socket, ::Type{AbstractString}, flags::Integer = CSymbols.NN_DONTWAIT)
    buf = Array{Ptr{Cchar}}(undef, 1)
    rc = _nn_recv(socket.s, convert(Ptr{Nothing}, pointer(buf)), CSymbols.NN_MSG, flags)
    if rc == -1
    	err = _nn_errno()
    	if err == CSymbols.EAGAIN
    		return nothing
    	end
        throw(NanomsgError("Socket recv failed", err))
    end
    str = bytestring(buf[1], rc)
    _nn_freemsg(convert(Ptr{Nothing}, buf[1]))

    # println("recv(s):", rc, str)
    return str
end

function Sockets.recv(socket::Socket, flags::Integer = CSymbols.NN_DONTWAIT)
    buf = Array{Ptr{Cuchar}}(undef, 1)
    rc = _nn_recv(socket.s, convert(Ptr{Nothing}, pointer(buf)), CSymbols.NN_MSG, flags)
    if rc == -1
    	err = _nn_errno()
    	if err == CSymbols.EAGAIN
    		return nothing
    	end
        throw(NanomsgError("Socket recv failed", err))
    end

    result::Ptr{Cuchar} = buf[1]
    arr = pointer_to_array(result, rc)
    finalizer((val) -> @async _nn_freemsg(convert(Ptr{Nothing}, result)), arr)

    #println("recv(a):", length(arr))
    return arr
end

function poll(sockets::Array{Socket}, pollIn::Bool = true, pollOut::Bool = true, timeout::Integer = -1)
	ct = length(sockets)
	go = ct > 0
	f = 0
	if pollIn
		f = CSymbols.NN_POLLIN
	end
	if pollOut
		f |= CSymbols.NN_POLLOUT
	end

	watch = map(s -> _NNPollFD(s.s, f, 0), sockets)
	ptr = convert(Ptr{Nothing}, pointer(watch))

	@task while go
		rc = _nn_poll(ptr, convert(Cint, ct), convert(Cint, timeout))
		if rc == -1
			throw(NanomsgError("Socket recv failed", _nn_errno()))
		end

		if rc == 0
			# timeout
			go = false
			continue
		end

		i = 1
		for t in watch
			rIn = ((t.revents & CSymbols.NN_POLLIN) != 0)
			rOut = ((t.revents & CSymbols.NN_POLLOUT) != 0)

			if rIn || rOut
				produce((sockets[i], i, rIn, rOut))
			end
			i = i + 1
		end
	end
end

struct _NNPollFD
    fd::Cint
    events::Cshort
    revents::Cshort
end

# Bit type for direct mapping to C struct
struct _NNSymbolProperties
    value::Cint
    name::Ptr{UInt8}
    ns::Cint
    typ::Cint
    unit::Cint

    _NNSymbolProperties() = new(0, 0, 0, 0, 0)
end

function Base.show(io::IO, prop::_NNSymbolProperties)
	name = bytestring(prop.name)

	print(io, "_NNSymbolProperties: value=", prop.value, " name=", name, ", namespace=", prop.ns, ", type=", prop.typ, ", unit=", prop.unit)
end

# Create an SP socket
_nn_socket(domain::Cint, protocol::Cint) = ccall((:nn_socket, LIB), Cint, (Cint, Cint), domain, protocol)

# Close an SP socket
_nn_close(s::Cint) = ccall((:nn_close, LIB), Cint, (Cint,), s)

# Set a socket option
_nn_setsockopt(s::Cint, level::Cint, option::Cint, optval::Ptr{Nothing}, optvallen::Csize_t) = ccall((:nn_setsockopt, LIB), Cint, (Cint, Cint, Cint, Ptr{Nothing}, Csize_t), s, level, option, optval, optvallen)

# Retrieve a socket option
_nn_getsockopt(s::Cint, level::Cint, option::Cint, optval::Ptr{Nothing}, optvallen::Csize_t) = ccall((:nn_getsockopt, LIB), Cint, (Cint, Cint, Cint, Ptr{Nothing}, Csize_t), s, level, option, optval, optvallen)

# Add a local endpoint to the socket
_nn_bind(s::Cint, addr::Ptr{UInt8}) = ccall((:nn_bind, LIB), Cint, (Cint,Ptr{UInt8}), s, addr)

# Add a remote endpoint to the socket
_nn_connect(s::Cint, addr::Ptr{UInt8}) = ccall((:nn_connect, LIB), Cint, (Cint,Ptr{UInt8}), s, addr)

# Remove an endpoint from the socket
_nn_shutdown(s::Cint, how::Cint) = ccall((:nn_shutdown, LIB), Cint, (Cint,Cint), s, how)

# Send a message
_nn_send(s::Cint, buf::Ptr{Nothing}, len::Csize_t, flags::Cint) = ccall((:nn_send, LIB), Cint, (Cint,Ptr{Nothing},Csize_t,Cint), s, buf, len, flags)

# Receive a message
_nn_recv(s::Cint, buf::Ptr{Nothing}, len::Csize_t, flags::Cint) = ccall((:nn_recv, LIB), Cint, (Cint,Ptr{Nothing},Csize_t,Cint), s, buf, len, flags)

# Fine-grained alternative to nn_send
####### nn_sendmsg(3)

# Fine-grained alternative to nn_recv
####### nn_recvmsg(3)

# Allocation of messages
####### nn_allocmsg(3)
####### nn_reallocmsg(3)

_nn_freemsg(buf::Ptr{Nothing}) = ccall((:nn_freemsg, LIB), Cint, (Ptr{Nothing},), buf)

# Manipulation of message control data
####### nn_cmsg(3)

# Multiplexing
_nn_poll(fds::Ptr{Nothing}, nfds::Cint, timeout::Cint) = ccall((:nn_poll, LIB), Cint, (Ptr{Nothing},Cint,Cint), fds, nfds, timeout)

# Retrieve the current errno
_nn_errno() = ccall((:nn_errno, LIB), Cint, ())

# Convert an error number into human-readable string
_nn_strerror(errno::Cint) = ccall((:nn_strerror, LIB), Ptr{UInt8}, (Cint,), errno)
function nn_strerror(errno::Cint)
	if errno == 0
		return nothing
	else
		cstr = _nn_strerror(errno)
		if cstr != C_NULL
			return bytestring(cstr)
		else
			return "Unknown errno"
		end

	end
end

# Check errno and throw an exception if set
function jl_nn_errno_check(context::AbstractString)
	errno = _nn_errno()
	if errno != 0
		throw(NanomsgError(context, errno))
	end
end

# Query the names and values of nanomsg symbols
_nn_symbol(i::Cint, value::Ptr{UInt8}) = ccall((:nn_symbol, LIB), Ptr{UInt8}, (Cint,Ptr{Cint}), i, value)

# Query properties of nanomsg symbols
function _nn_symbol_info(i::Cint)
	buflen = sizeof(_NNSymbolProperties)
	buf = Array{UInt8}(undef, buflen)

	r = ccall((:nn_symbol_info, LIB), Cint, (Cint,Ptr{UInt8},Csize_t), i, buf, buflen)

	if r == 0
		return nothing
	end

	if r != buflen
		throw(NanomsgError("Failed to query symbol info: " * r))
	end

	ptr = convert(Ptr{_NNSymbolProperties}, pointer(buf))
	unsafe_load(ptr)
end

# Start a device
_nn_device(s1::Cint, s2::Cint) = ccall((:nn_device, LIB), Cint, (Cint,Cint), s1, s2)

# Notify all sockets about process termination
_nn_term() = ccall((:nn_term, LIB), Nothing, ())

function j_nn_load_symbols()
	symbols = Dict{Cint, Dict{Cint, AbstractString}}()
	index::Cint = 0
	while true
		value = _nn_symbol_info(index)
		if value == nothing
			break
		end

		entry = get(symbols, value.ns, Dict{Cint, AbstractString}())
		entry[value.value] = unsafe_string(value.name)
		symbols[value.ns] = entry
		index = index + 1
	end

	symbols
end

macro load_symbols()
	blk = quote
	end

	symbols = CSymbols.J_NN_MAP
	ns = symbols[0]
	rstr = string("")
	for (v,k) in ns
		rstr = string(rstr, string("$k=$v\n"))
		ks = Symbol(k)
		push!(blk.args, Expr(:module, false, esc(ks::Symbol), Expr(:begin)))

		if v == 0
			# Namespace itself, skip as the constants will appear later
			continue
		end

		for (ev, ek) in symbols[v]
			rstr = string(rstr, string("\t$ek=$ev\n"))
			eks = Symbol(ek)
			push!(blk.args, :(const $(esc(eks)) = $ev))
		end
	end

	push!(blk.args, :(const $(Expr(:escape, :J_NN_ALL)) = $rstr))
	#println(blk)
	push!(blk.args, :nothing)
	blk.head = :toplevel
	return blk
end

# This module is dynamically loaded with constants from the
# nanomg library. Refer to the official docs for a list
baremodule CSymbols
	using Base
	using Nanomsg: j_nn_load_symbols
	using Nanomsg: @load_symbols

	const J_NN_MAP = j_nn_load_symbols()

	const NN_MSG = ~convert(Csize_t, 0)
	const NN_NO_FLAG = convert(Cint, 0)
	const NN_DONTWAIT = convert(Cint, 1)
	@load_symbols()
end

end


text:
# Credit to https://github.com/bischtob/SEC/blob/main/src/Sample/mcmc.jl

using Distributions
import Distributions: Uniform

"""
acceptance criteria for Metropolis-Hastings
# Definition
accept(Δ) = log(rand(Uniform(0, 1))) < Δ
# Arguments
- `Δ`: (scalar): change in negative log-likelihood
# Output
- true or false (bool)
# Notes
Always go downhill but sometimes go uphill
"""
accept(Δ) = log(rand(Uniform(0, 1))) < Δ

"""
markov_link(nll, proposal, current_X, current_nll)
# Description
- Takes a single step in the random walk markov chain monte carlo algorithm and outputs proposal parameters,
  new parameters, and the evaluate of the loss function
# Arguments
- `nll`: The negative log-likelihood function. In the absence of priors this becomes a loss function
- `proposal`: (function), determines the proposal step
- `current_X`: (array), current parameter
- `current_nll`: (scalar), proposal_nll = nll(X). The value of negative log-likelihood of the current parameter
# Return
- `new_X`: The value of the accepted X
- `new_nll`: value of nll(new_X)
- `proposal_X`: The X from the "proposal step". Was either rejected or accepted.
- `proposal_nll`: value of nll(proposal_X)
"""
function markov_link(proposal_X, proposal_nll::AbstractFloat, current_X, current_nll::AbstractFloat)
    Δ = (current_nll - proposal_nll)

    if accept(Δ)
        new_X, new_nll = proposal_X, proposal_nll
    else
        new_X, new_nll = current_X, current_nll
    end

    return new_X, new_nll
end

function markov_link(proposal_X, proposal_nll::Vector, current_X, current_nll::Vector)

    # Vector of NamedTuples
    result = markov_link.(proposal_X, proposal_nll, current_X, current_nll)
    # result = NamedTuple(key => getproperty.(result, key) for key in keys(result))

    return getindex.(result, 1), getindex.(result, 2)
end

function proposal_X_nll(nll, proposal, current_X, current_nll::AbstractFloat)
    proposal_X = proposal(current_X)
    proposal_nll = nll(proposal_X)

    return proposal_X, proposal_nll
end

function proposal_X_nll(nll, proposal, current_X, current_nll::Vector)
    proposal_X = proposal.(current_X)
    proposal_nll = nll(proposal_X)

    return proposal_X, proposal_nll
end

"""
markov_chain(nll, proposal, seed_X, chain_length; random_seed = 1234)
# Description
- A random walk that computes the posterior distribution
# Arguments
- `nll`: The negative log-likelihood function. In the absence of priors this becomes a loss function
- `proposal`: (function), proposal function for MCMC
- `seed_X`: (Array), initial parameter values
- `chain_length`: (Int) number of markov chain monte carlo steps
- `perturb`: a function that performs a perturbation of X
# Keyword Arguments
- `burn_in`: (Int) number of initial steps to throw away
- `n_chains`: (Int) number of parallel chains of MCMC to run, in the case that `nll` is parallelizable.
If `n_chains == 1`, it is assumed that `nll` maps arrays of size (n_param,) to scalar values.
If `n_chains > 1`, it is assumed that `nll` maps vectors of arrays of size `(n_param,)` to arrays of size (1, n_chains)`,
and the `chain_length` and `burn_in` period will be divided amongst the `n_chain` chains.
# Return
- `chain_X`: The matrix of accepted parameters in the random walk
- `chain_nll`: The array of errors associated with each step in param chain
"""
function markov_chain(nll, proposal, seed_X, chain_length::Int; burn_in=0, n_chains=1)

    current_X = seed_X # [n_param, n_chains]
    current_nll = nll(seed_X) # [1, n_chains]
    chain_X = typeof(current_X)[]
    chain_nll = typeof(current_nll)[]
    push!(chain_X, current_X)
    push!(chain_nll, current_nll)

    for i = 1:chain_length-1
        for chain in n_chains
            proposal_X, proposal_nll = proposal_X_nll(nll, proposal, current_X, current_nll)
            new_X, new_nll = markov_link(proposal_X, proposal_nll, current_X, current_nll)
            current_X, current_nll = new_X, new_nll # mcmc update
            push!(chain_X, new_X)
            push!(chain_nll, new_nll)
        end
    end

    if n_chains > 1
        chain_X = vcat(chain_X...)
        chain_nll = vcat(chain_nll...)
    end

    total_burn_in = burn_in * n_chains

    n_forward_runs = n_chains * (chain_length - burn_in)
    @info "$n_forward_runs samples were generated by MCMC."

    return (chain_X[total_burn_in:end], chain_nll[total_burn_in:end])
end


text:
@testset "cmp_search_results" begin
    idxs = shuffle!(randsubseq(1:100, 0.5))
    ds = rand(Float64, length(idxs))

    p = randperm(length(idxs))
    @test cmp_search_results((idxs, ds), (idxs[p], ds[p]))

    @test !cmp_search_results((idxs, ds), (idxs[2:end], ds[2:end]))
    @test !cmp_search_results((idxs, ds), (idxs[p], ds))
end


text:
struct TrackedNeuralODE{R,Z,M,P,RE,T,A,K} <: DiffEqFlux.NeuralDELayer
    model::M
    p::P
    re::RE
    tspan::T
    args::A
    kwargs::K
    time_dep::Bool

    function TrackedNeuralODE(model, tspan, time_dep, regularize, args...; kwargs...)
        return_multiple = get(kwargs, :save_everystep, false) || haskey(kwargs, :saveat)
        p, re = Flux.destructure(model)
        kwargs = Dict(kwargs)
        new{
            regularize,
            return_multiple,
            typeof(model),
            typeof(p),
            typeof(re),
            typeof(tspan),
            typeof(args),
            typeof(kwargs),
        }(
            model,
            p,
            re,
            tspan,
            args,
            kwargs,
            time_dep,
        )
    end
end

function update_saveat!(kwargs, saveat)
    if isnothing(saveat)
        return kwargs, () -> nothing
    end
    original_saveat = get(kwargs, :saveat, nothing)
    function restore!()
        kwargs[:saveat] = original_saveat
        return kwargs
    end
    kwargs[:saveat] = saveat
    return kwargs, restore!
end

@fastmath function (n::TrackedNeuralODE{false,false})(
    x,
    p = n.p;
    func = (u, t, int) -> 0,
    tspan = nothing,
    saveat = nothing,
)
    dudt_(u, p, t) = n.time_dep ? n.re(p)(u, t) : n.re(p)(u)

    tspan = _convert_tspan(isnothing(tspan) ? n.tspan : tspan, p)

    kwargs, restore_kwargs! = update_saveat!(n.kwargs, saveat)

    ff = ODEFunction{false}(dudt_, tgrad = DiffEqFlux.basic_tgrad)
    prob = ODEProblem{false}(ff, x, tspan, p)

    sol = solve(
        prob,
        n.args...;
        sensealg = SensitivityADPassThrough(),
        callback = nothing,
        kwargs...,
    )
    res = diffeqsol_to_trackedarray(sol)::typeof(x)
    nfe = sol.destats.nf::Int

    restore_kwargs!()

    return res, nfe, nothing
end

@fastmath function (n::TrackedNeuralODE{false,true})(
    x,
    p = n.p;
    func = (u, t, int) -> 0,
    tspan = nothing,
    saveat = nothing,
)
    dudt_(u, p, t) = n.time_dep ? n.re(p)(u, t) : n.re(p)(u)

    tspan = _convert_tspan(isnothing(tspan) ? n.tspan : tspan, p)

    kwargs, restore_kwargs! = update_saveat!(n.kwargs, saveat)

    ff = ODEFunction{false}(dudt_, tgrad = DiffEqFlux.basic_tgrad)
    prob = ODEProblem{false}(ff, x, tspan, p)

    sol = solve(
        prob,
        n.args...;
        sensealg = SensitivityADPassThrough(),
        callback = nothing,
        kwargs...,
    )
    res = diffeqsol_to_3dtrackedarray(sol)::TrackedArray{Float32,3,CuArray{Float32,3}}
    nfe = sol.destats.nf::Int

    restore_kwargs!()

    return res, nfe, nothing
end

@fastmath function (n::TrackedNeuralODE{true,false})(
    x,
    p = n.p;
    # Default is to regularize using Error Estimates. Alternative tested
    # strategy is using Stiffness Estimate:
    # (u, t, integrator) -> integrator.eigen_est * integrator.dt
    func = (u, t, integrator) -> integrator.EEst * integrator.dt,
    tspan = nothing,
    saveat = nothing,
)
    dudt_(u, p, t) = n.time_dep ? n.re(p)(u, t) : n.re(p)(u)

    tspan = _convert_tspan(isnothing(tspan) ? n.tspan : tspan, p)

    kwargs, restore_kwargs! = update_saveat!(n.kwargs, saveat)

    sv = SavedValues(eltype(tspan), eltype(p))
    svcb = SavingCallback(func, sv)
    ff = ODEFunction{false}(dudt_, tgrad = DiffEqFlux.basic_tgrad)
    prob = ODEProblem{false}(ff, x, tspan, p)

    sol = solve(
        prob,
        n.args...;
        sensealg = SensitivityADPassThrough(),
        callback = svcb,
        kwargs...,
    )
    res = diffeqsol_to_trackedarray(sol)::typeof(x)

    restore_kwargs!()

    nfe = sol.destats.nf::Int
    return res, nfe, sv
end

@fastmath function (n::TrackedNeuralODE{true,true})(
    x,
    p = n.p;
    # Default is to regularize using Error Estimates. Alternative tested
    # strategy is using Stiffness Estimate:
    # (u, t, integrator) -> integrator.eigen_est * integrator.dt
    func = (u, t, integrator) -> integrator.EEst * integrator.dt,
    tspan = nothing,
    saveat = nothing,
)
    dudt_(u, p, t) = n.time_dep ? n.re(p)(u, t) : n.re(p)(u)

    tspan = _convert_tspan(isnothing(tspan) ? n.tspan : tspan, p)

    kwargs, restore_kwargs! = update_saveat!(n.kwargs, saveat)

    sv = SavedValues(eltype(tspan), eltype(p))
    svcb = SavingCallback(func, sv)
    ff = ODEFunction{false}(dudt_, tgrad = n.time_dep ? nothing : DiffEqFlux.basic_tgrad)
    prob = ODEProblem{false}(ff, x, tspan, p)

    sol = solve(
        prob,
        n.args...;
        sensealg = SensitivityADPassThrough(),
        callback = svcb,
        kwargs...,
    )
    res = diffeqsol_to_3dtrackedarray(sol)::TrackedArray{Float32,3,CuArray{Float32,3}}

    restore_kwargs!()

    nfe = sol.destats.nf::Int
    return res, nfe, sv
end

function solution(
    n::TrackedNeuralODE,
    x,
    p = n.p;
    solver = nothing,
    tspan = nothing,
    saveat = nothing,
)
    solver = isnothing(solver_override) ? n.args[1] : solver
    dudt_(u, p, t) = n.time_dep ? n.re(p)(u, t) : n.re(p)(u)

    tspan = _convert_tspan(isnothing(tspan) ? n.tspan : tspan, p)

    kwargs, restore_kwargs! = update_saveat!(n.kwargs, saveat)

    ff = ODEFunction{false}(dudt_, tgrad = n.time_dep ? nothing : DiffEqFlux.basic_tgrad)
    prob = ODEProblem{false}(ff, x, tspan, p)

    sol = solve(
        prob,
        solver;
        sensealg = SensitivityADPassThrough(),
        callback = nothing,
        kwargs...,
    )

    restore_kwargs!()

    return sol
end


text:
struct LiftUpdater{T, TV<:AbstractVector{T}, TPD}
    pd::TPD
    ρ::TV
    g::TV
    ng_approx::TV
    n_j::Int
end

function (lu::LiftUpdater{T})() where T
    @unpack ρ, pd = lu
    @unpack x, x1, σ  = pd
    n_i = length(ρ)
    w = computew(x, x1, σ)
    w /= 2
    lift = false
    for i in 1:n_i
        lift = lift || lu(w, i)
    end
    return lift
end
function computew(x::Vector{T}, x1, σ) where {T}
    tmapreduce(+, 1:length(x), init = zero(T)) do j
        (x[j] - x1[j])^2 / (σ[j]^2 - (x[j] - x1[j])^2)
    end
end

function computew(x::CuVector{T}, x1, σ, ::Val{blocksize} = Val(80), ::Val{threads} = Val(256)) where {T, blocksize, threads}
    result = similar(x, T, (blocksize,))
    args = (x, x1, σ, result, Val(threads))
    @cuda blocks = blocksize threads = threads computew_kernel(args...)
    CUDAnative.synchronize()
    w = sum(Array(result))
end
function computew_kernel(x::AbstractVector{T}, x1, σ, result, ::Val{LMEM}) where {T, LMEM}
    @mapreduce_block(j, length(x), +, T, LMEM, result, begin
        (x[j] - x1[j])^2 / (σ[j]^2 - (x[j] - x1[j])^2)
    end)

    return
end

function (lu::LiftUpdater{T})(w, i::Int) where T
    @unpack ρ, g, ng_approx = lu
    gi, g_approxi = g[i], -ng_approx[i]
    lift = (gi > g_approxi)
    δi = (gi - g_approxi)/w
    ρi = ρ[i]
    ρ[i] = ifelse(lift, min(10ρi, T(1.1)*(ρi+δi)), ρi)
    return lift
end

struct LiftResetter{T, TV}
    ρ::TV
    ρmin::T
end
function (lr::LiftResetter)(k::Iteration)
    @unpack ρ = lr
    if k.i > 1
        map!(lr, ρ, 1:length(ρ))
    end
    return
end
function (lr::LiftResetter{T})(i::Int) where T
    @unpack ρmin, ρ = lr
    return max(T(0.1)*ρ[i], ρmin)
end


text:
#####
##### Broadcasting utilities
#####

using Base.Broadcast: DefaultArrayStyle
using Base.Broadcast: Broadcasted
using CUDA

using Oceananigans.Architectures: device_event

struct FieldBroadcastStyle <: Broadcast.AbstractArrayStyle{3} end

Base.Broadcast.BroadcastStyle(::Type{<:AbstractField}) = FieldBroadcastStyle()

# Precedence rule
Base.Broadcast.BroadcastStyle(::FieldBroadcastStyle, ::DefaultArrayStyle{N}) where N = FieldBroadcastStyle()
Base.Broadcast.BroadcastStyle(::FieldBroadcastStyle, ::CUDA.CuArrayStyle{N}) where N = FieldBroadcastStyle()

# For use in Base.copy when broadcasting with numbers and arrays (useful for comparisons like f::AbstractField .== 0)
Base.similar(bc::Broadcasted{FieldBroadcastStyle}, ::Type{ElType}) where ElType = similar(Array{ElType}, axes(bc))

# Bypass style combining for in-place broadcasting with arrays / scalars to use built-in broadcasting machinery
const BroadcastedArrayOrCuArray = Union{Broadcasted{<:DefaultArrayStyle},
                                        Broadcasted{<:CUDA.CuArrayStyle}}

@inline function Base.Broadcast.materialize!(dest::Field, bc::BroadcastedArrayOrCuArray)
    if any(a isa OffsetArray for a in bc.args)
        return Base.Broadcast.materialize!(dest.data, bc)
    else
        return Base.Broadcast.materialize!(interior(dest), bc)
    end
end

# TODO: make this support Field that are windowed in _only_ 1 or 2 dimensions.
# Right now, this may only produce expected behavior (re: dimensionality) for
# WindowedField that are windowed in three-dimensions. Of course, broadcasting with
# scalar `bc` is no issue.
@inline Base.Broadcast.materialize!(dest::WindowedField, bc::BroadcastedArrayOrCuArray) =
    Base.Broadcast.materialize!(parent(dest), bc)

#####
##### Kernels
#####

@inline offset_compute_index(::Colon, i) = i
@inline offset_compute_index(range::UnitRange, i) = range[1] + i - 1

@kernel function broadcast_kernel!(dest, bc, index_ranges)
    i, j, k = @index(Global, NTuple)

    i′ = offset_compute_index(index_ranges[1], i)
    j′ = offset_compute_index(index_ranges[2], j)
    k′ = offset_compute_index(index_ranges[3], k)

    @inbounds dest[i′, j′, k′] = bc[i′, j′, k′]
end

# Interface for getting AbstractOperation right
broadcasted_to_abstract_operation(loc, grid, a) = a

# Broadcasting with interpolation breaks Base's default rules,
# so we bypass the infrastructure for checking axes compatibility,
# and head straight to copyto! from materialize!.
@inline Base.Broadcast.materialize!(::Base.Broadcast.BroadcastStyle,
                                    dest::Field,
                                    bc::Broadcasted{<:FieldBroadcastStyle}) = copyto!(dest, convert(Broadcasted{Nothing}, bc))

@inline function Base.copyto!(dest::Field, bc::Broadcasted{Nothing})

    grid = dest.grid
    arch = architecture(dest)

    bc′ = broadcasted_to_abstract_operation(location(dest), grid, bc)

    event = launch!(arch, grid, size(dest), broadcast_kernel!, dest, bc′, dest.indices,
                    dependencies = device_event(arch))

    wait(device(arch), event)

    return dest
end


text:
function tr4!(r::Robot)#главная функция
    move_corner(r)#перемещаем робота в угол
    mark!(r)
end
function mark!(r::Robot)
    for i in 1:11
        for i in i:11
            if isborder(r, Ost) == false
                putmarker!(r)
                move!(r, Ost)
            end
        end
        putmarker!(r)
        for i in i:11
            if isborder(r, West) == false
                move!(r, West)
            end
        end
        if isborder(r,Nord) == false
            move!(r,Nord)
        end
        i=i+1
    end
end
function move_corner(r::Robot) #перемещаем робота в угол
    while isborder(r, West) == false
        move!(r, West)
    end
    while isborder(r, Sud) == false
        move!(r, Sud)
    end
end


text:
@testset "xy06" begin
    x, y = xy06(2400000.5, 53736.0)
    @test x ≈ 0.5791308486706010975e-3 atol=1e-16
    @test y ≈ 0.4020579816732958141e-4 atol=1e-17
end

@testset "xys00a" begin
    x, y, s = xys00a(2400000.5, 53736.0)
    @test x ≈ 0.5791308472168152904e-3 atol=1e-16
    @test y ≈ 0.4020595661591500259e-4 atol=3e-17  # originally eps=1e-17; relaxed for 32-bit windows
    @test s ≈ -0.1220040848471549623e-7 atol=1e-20
end

@testset "xys00b" begin
    x, y, s = xys00b(2400000.5, 53736.0)
    @test x ≈ 0.5791301929950208873e-3 atol=1e-16
    @test y ≈ 0.4020553681373720832e-4 atol=1e-16
    @test s ≈ -0.1220027377285083189e-7 atol=1e-19
end

@testset "xys06a" begin
    x, y, s = xys06a(2400000.5, 53736.0)
    @test x ≈ 0.5791308482835292617e-3 atol=1e-16
    @test y ≈ 0.4020580099454020310e-4 atol=1e-15
    @test s ≈ -0.1220032294164579896e-7 atol=1e-19
end


text:
#!/cvmfs/soft.computecanada.ca/easybuild/software/2017/avx2/Compiler/gcc7.3/julia/1.3.0/bin/julia
#SBATCH --array=1-200
#SBATCH -o reproduce_test.%A_%a.out # Standard output
#SBATCH -e reproduce_test.%A_%a.err # Standard error
#SBATCH --mem-per-cpu=512M # Memory request of 2 GB
#SBATCH --time=0:01:00 # Running time of 12 hours
#SBATCH --account=rrg-whitem

using Pkg
Pkg.activate(".")

using Reproduce
const config_file = "./examples/arg_iter_config.toml"

experiment = Experiment(config_file)
create_experiment_dir(experiment; tldr="hello tldr")
ret = job(experiment)


text:
using MyPkgDemo
using Test

@testset "Tests for cat_str" begin
  @test cat_str("a", "b") == "ab"
  @test cat_str("a", "c") == "ac"
  @test cat_str("a", "b", "c") == "abc"
  @test cat_str("a", "c", "e") == "ace"
end


text:
using StatsBase
using Plots

Parameters = [[[0 0 ; 0 0.16], [0 ; 0]],
        [[0.85 0.04 ; -0.04 0.85], [0 ; 1.6]],
        [[0.20 -0.26 ; 0.23 0.22], [0 ; 1.6]],
        [[-0.15 0.28 ; 0.26 0.24], [0 ; 0.44]]]

f1(vector, Pars) = Pars[1][1] * vector + Pars[1][2]
f2(vector, Pars) = Pars[2][1] * vector + Pars[2][2]
f3(vector, Pars) = Pars[3][1] * vector + Pars[3][2]
f4(vector, Pars) = Pars[4][1] * vector + Pars[4][2]

function RandomBarnsleyFern(P, PointNums)
    Functions = [f1, f2, f3, f4]
    X = [0.0 for i = 1:PointNums]
    Y = [0.0 for i = 1:PointNums]
    for i in 1:PointNums
        point = [rand() ; rand()]
        for j in 1:P
            s = sample([1,2,3,4], Weights([0.01,0.85,0.07,0.07]))
            point = Functions[s](point, Parameters)
        end
        X[i]= point[1]
        Y[i]= point[2]
    end
    return X, Y
end

for i in [10, 100]
    for j in [10000,1000000]
        x,y = RandomBarnsleyFern(i, j)

        scatter(x,y,markersize = 10^4/j, legend = false, border=:none, color =:green,dpi=200, label = "Points number = $j Operations number = $i")
        savefig("C:\\Users\\Yaghoub\\Documents\\GitHub\\ComputationalPhysics-Fall2021\\ProblemSet1\\Figs\\Q6\\RandomBarnsleyFern-p$(i)num$(j).png")
    end
end


text:
"""
Reference
---------
- `run.m`
- `gradient_gen.m`
"""

# gradient_func(dt, X)
gradient_func = gradient_gen(sp.ODE!, sp.p, 2, 3)

@testset "Gradient Results" begin
@test gradient_func(1., ones(6,1)) == sp.gradient(1., ones(6,1))
@test gradient_func(1., 2*ones(6,1)) == sp.gradient(1., 2*ones(6,1))
@test gradient_func(1., zeros(6,1)) == sp.gradient(1., zeros(6,1))
@test gradient_func(2., ones(6,1)) == sp.gradient(2., ones(6,1))
end

@testset "Gradient Results: Screening test" for v in ParameterGrid( [1:2, 1:2, 1:2, 1:2, 1:2, 1:2] )

# Convet Num to Float
grad = map(y->y, gradient_func(0.1, v));

@test grad ≈ sp.gradient(0.1, v);

end;

# Gradient function : mutated
a = ones(6,1)
a_ = deepcopy(a)
x = deepcopy(a)
g! = LeastAction.grad!(gradient_func, 0.1, 2, 3)

g!(a, x)

@test a != a_


text:
# Autogenerated wrapper script for ALPS_jll for armv7l-linux-gnueabihf-cxx11
export libalps

using CoinUtils_jll
using Osi_jll
using Clp_jll
using Cgl_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("ALPS")
JLLWrappers.@declare_library_product(libalps, "libAlps.so.0")
function __init__()
    JLLWrappers.@generate_init_header(CoinUtils_jll, Osi_jll, Clp_jll, Cgl_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libalps,
        "lib/libAlps.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Pkg
Pkg.add("JSON3")

import JSON3
using Sockets

struct Coordinate
    x::Float64
    y::Float64
    z::Float64
end

function calc(text)
    jobj = JSON3.read(text)
    coordinates = jobj["coordinates"]
    len = length(coordinates)
    x = y = z = 0

    for coord in coordinates
        x += coord["x"]
        y += coord["y"]
        z += coord["z"]
    end

    Coordinate(x / len, y / len, z / len)
end

function notify(msg)
    try
        socket = connect("localhost", 9001)
        write(socket, msg)
        close(socket)
    catch
        # standalone usage
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    left = calc("""{"coordinates":[{"x":1.1,"y":2.2,"z":3.3}]}""")
    right = Coordinate(1.1, 2.2, 3.3)
    if left != right
        println(stderr, "$(left) != $(right)")
        exit(1)
    end

    text = open("/tmp/1.json") do file
        read(file, String)
    end

    notify("Julia (JSON3)\t$(getpid())")

    println(calc(text))

    notify("stop")
end

max_stars_repo_name: edubart/benchmarks

max_issues_repo_name: edubart/benchmarks

max_forks_repo_name: edubart/benchmarks