Move tutorials to Quarto

docs/make.jl

@@ -87,7 +87,11 @@ makedocs(
     sitename="Manifolds.jl",
     pages=[
         "Home" => "index.md",
-        "How to..." => ["🚀 Get Started with `Manifolds.jl`" => "tutorials/getstarted.md"],
+        "How to..." => [
+            "🚀 Get Started with `Manifolds.jl`" => "tutorials/getstarted.md",
+            "work in charts" => "tutorials/working-in-charts.md",
+            "perform Hand gesture analysis" => "tutorials/hand-gestures.md",
+        ],
         "Manifolds" => [
             "Basic manifolds" => [
                 "Centered matrices" => "manifolds/centeredmatrices.md",

tutorials/Project.toml

@@ -1,10 +1,21 @@
 [deps]
+BoundaryValueDiffEq = "764a87c0-6b3e-53db-9096-fe964310641d"
+CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
 IJulia = "7073ff75-c697-5162-941a-fcdaad2a7d2a"
 Manifolds = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e"
 ManifoldsBase = "3362f125-f0bb-47a3-aa74-596ffd7ef2fb"
 Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
+MultivariateStats = "6f286f6a-111f-5878-ab1e-185364afe411"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 
 [compat]
+CSV = "0.10"
+DataFrames = "1"
 IJulia = "1"
 Manifolds = "0.8.46"
 ManifoldsBase = "0.14.5"
+MultivariateStats = "0.10"

tutorials/hand-gestures.qmd

@@ -0,0 +1,110 @@
+---
+title: Hand gesture analysis
+---
+
+In this tutorial we will learn how to use Kendall's shape space to analyze hand gesture data.
+
+```{julia}
+#| echo: false
+#| code-fold: true
+#| output: false
+using Pkg;
+cd(@__DIR__)
+Pkg.activate("."); # for reproducibility use the local tutorial environment.
+using Markdown
+```
+
+Let's start by loading libraries required for our work.
+
+```{julia}
+using Manifolds, CSV, DataFrames, Plots, MultivariateStats
+```
+
+Our first function loads dataset of hand gestures, described [here](https://geomstats.github.io/notebooks/14_real_world_applications__hand_poses_analysis_in_kendall_shape_space.html).
+
+```{julia}
+
+function load_hands()
+    hands_url = "https://raw.githubusercontent.com/geomstats/geomstats/master/geomstats/datasets/data/hands/hands.txt"
+    hand_labels_url = "https://raw.githubusercontent.com/geomstats/geomstats/master/geomstats/datasets/data/hands/labels.txt"
+
+    hands = Matrix(CSV.read(download(hands_url), DataFrame, header=false))
+    hands = reshape(hands, size(hands, 1), 3, 22)
+    hand_labels = CSV.read(download(hand_labels_url), DataFrame, header=false).Column1
+    return hands, hand_labels
+end
+```
+
+The following code plots a sample gesture as a 3D scatter plot of points.
+```{julia}
+hands, hand_labels = load_hands()
+scatter3d(hands[1, 1, :], hands[1, 2, :], hands[1, 3, :])
+```
+
+Each gesture is represented by 22 landmarks in $ℝ³$, so we use the appropriate Kendall's shape space
+```{julia}
+Mshape = KendallsShapeSpace(3, 22)
+```
+
+Hands read from the dataset are projected to the shape space to remove translation
+and scaling variability. Rotational variability is then handled using the quotient
+structure of ``[`KendallsShapeSpace`](@ref)``{=commonmark}
+```{julia}
+#| output: false
+hands_projected = [project(Mshape, hands[i, :, :]) for i in axes(hands, 1)]
+```
+
+In the next part let's do tangent space PCA. This starts with computing a mean point and computing
+logithmic maps at mean to each point in the dataset.
+```{julia}
+#| output: false
+mean_hand = mean(Mshape, hands_projected)
+hand_logs = [log(Mshape, mean_hand, p) for p in hands_projected]
+```
+
+For a tangent PCA, we need coordinates in a basis.
+Some libraries skip this step because the representation of tangent vectors
+forms a linear subspace of an Euclidean space so PCA automatically detects
+which directions have no variance but this is a more generic way to solve
+this issue.
+```{julia}
+#| output: false
+B = get_basis(Mshape, mean_hand, ProjectedOrthonormalBasis(:svd))
+hand_log_coordinates = [get_coordinates(Mshape, mean_hand, X, B) for X in hand_logs]
+```
+
+This code prepares data for MultivariateStats -- `mean=0` is set because we've centered
+the data geometrically to `mean_hand` in the code above.
+```{julia}
+red_coords = reduce(hcat, hand_log_coordinates)
+fp = fit(PCA, red_coords; mean=0)
+```
+
+Now let's show explained variance of each principal component.
+```{julia}
+plot(principalvars(fp), title="explained variance", label="Tangent PCA")
+```
+
+The next plot shows how projections on the first two pricipal components look like.
+```{julia}
+fig = plot(; title="coordinates per gesture of the first two principal components")
+for label_num in [0, 1]
+    mask = hand_labels .== label_num
+    cur_hand_logs = red_coords[:, mask]
+    cur_t = MultivariateStats.transform(fp, cur_hand_logs)
+    scatter!(fig, cur_t[1, :], cur_t[2, :], label="gesture " * string(label_num))
+end
+xlabel!(fig, "principal component 1")
+ylabel!(fig, "principal component 2")
+fig
+```
+
+The following heatmap displays pairwise distances between gestures.
+We can use them for clustering, classification, etc.
+```{julia}
+hand_distances = [
+    distance(Mshape, hands_projected[i], hands_projected[j]) for
+    i in eachindex(hands_projected), j in eachindex(hands_projected)
+]
+heatmap(hand_distances, aspect_ratio=:equal)
+```

tutorials/working-in-charts.qmd

@@ -0,0 +1,180 @@
+---
+title: Working in charts
+---
+
+In this tutorial we will learn how to use charts for basic geometric operations like exponential map, logarithmic map and parallel transport.
+
+```{julia}
+#| echo: false
+#| code-fold: true
+#| output: false
+using Pkg;
+cd(@__DIR__)
+Pkg.activate("."); # for reproducibility use the local tutorial environment.
+using Markdown
+```
+
+
+There are two conceptually different approaches to working on a manifold: working in charts and chart-free representations.
+
+The first one, widespread in differential geometry textbooks, is based on defining an atlas on the manifold and performing computations in selected charts. This approach, while generic, is not ideally suitable in all circumstances. For example, working in charts that do not cover the entire manifold causes issues with having to switch charts when operating on a manifold.
+
+The second one is beneficital, if there exist a representation of points and tangent vectors for a manifold, which allow for efficient closed-form formulas for standard functions like the exponential map or Riemannian distance in this representation. These computations are then chart-free. `Manifolds.jl` supports both approaches, although the chart-free approach is the main focus of the library.
+
+In this tutorial we focus on chart-based computation.
+
+```{julia}
+#| output: false
+using Manifolds, RecursiveArrayTools, OrdinaryDiffEq, DiffEqCallbacks, BoundaryValueDiffEq
+```
+
+The manifold we consider is the `M` is the torus in form of the [`EmbeddedTorus`](https://juliamanifolds.github.io/Manifolds.jl/latest/manifolds/torus.html#Manifolds.EmbeddedTorus), that is the representation defined as a surface of revolution of a circle of radius 2 around a circle of radius 3.
+The atlas we will perform computations in is its [`DefaultTorusAtlas`](https://juliamanifolds.github.io/Manifolds.jl/latest/manifolds/torus.html#Manifolds.DefaultTorusAtlas) `A`, consistting of a family of charts indexed by two angles, that specify the base point of the chart.
+
+We will draw geodesics time between `0` and `t_end`, and then sample the solution at multiples of `dt` and draw a line connecting sampled points.
+
+```{julia}
+M = Manifolds.EmbeddedTorus(3, 2)
+A = Manifolds.DefaultTorusAtlas()
+```
+
+## Setup
+
+We will first set up our plot with an empty torus.
+`param_points` are points on the surface of the torus that will be used for basic surface shape in `Makie.jl`.
+The torus will be colored according to its Gaussian curvature stored in `gcs`. We later want to have a color scale that has negative curvature blue, zero curvature white and positive curvature red so `gcs_mm` is the largest absolute value of the curvature that will be needed to properly set range of curvature values.
+
+In the documentation this tutorial represents a static situation (without interactivity). `Makie.jl` rendering is turned off.
+
+```{julia}
+# using GLMakie, Makie
+# GLMakie.activate!()
+
+"""
+	torus_figure()
+
+This function generates a simple plot of a torus and returns the new figure containing the plot.
+"""
+function torus_figure()
+    fig = Figure(resolution=(1400, 1000), fontsize=16)
+    ax = LScene(fig[1, 1], show_axis=true)
+    ϴs, φs = LinRange(-π, π, 50), LinRange(-π, π, 50)
+    param_points = [Manifolds._torus_param(M, θ, φ) for θ in ϴs, φ in φs]
+    X1, Y1, Z1 = [[p[i] for p in param_points] for i in 1:3]
+    gcs = [gaussian_curvature(M, p) for p in param_points]
+    gcs_mm = max(abs(minimum(gcs)), abs(maximum(gcs)))
+    pltobj = surface!(
+        ax,
+        X1,
+        Y1,
+        Z1;
+        shading=true,
+        ambient=Vec3f(0.65, 0.65, 0.65),
+        backlight=1.0f0,
+        color=gcs,
+        colormap=Reverse(:RdBu),
+        colorrange=(-gcs_mm, gcs_mm),
+        transparency=true,
+    )
+    wireframe!(ax, X1, Y1, Z1; transparency=true, color=:gray, linewidth=0.5)
+    zoom!(ax.scene, cameracontrols(ax.scene), 0.98)
+    Colorbar(fig[1, 2], pltobj, height=Relative(0.5), label="Gaussian curvature")
+    return ax, fig
+end
+```
+
+
+## Values for the geodesic
+
+`solve_for` is a helper function that solves a parallel transport along geodesic problem on the torus `M`.
+`p0x` is the $(\theta, \varphi)$ parametrization of the point from which we will transport the vector.
+We first calculate the coordinates in the embedding of `p0x` and store it as `p`, and then get the initial chart from atlas `A` appropriate for starting working at point `p`.
+The vector we transport has coordinates `Y_transp` in the induced tangent space basis of chart `i_p0x`.
+The function returns the full solution to the parallel transport problem, containing the sequence of charts that was used and solutions of differential equations computed using `OrdinaryDiffEq`.
+
+`bvp_i` is needed later for a different purpose, it is the chart index we will use for solving the logarithmic map boundary value problem in.
+
+Next we solve the vector transport problem `solve_for([θₚ, φₚ], [θₓ, φₓ], [θy, φy])`, sample the result at the selected time steps and store the result in `geo`. The solution includes the geodesic which we extract and convert to a sequence of points digestible by `Makie.jl`, `geo_ps`.
+`[θₚ, φₚ]` is the parametrization in chart (0, 0) of the starting point of the geodesic.
+The direction of the geodesic is determined by `[θₓ, φₓ]`, coordinates of the tangent vector at the starting point expressed in the induced basis of chart `i_p0x` (which depends on the initial point).
+Finally, `[θy, φy]` are the coordinates of the tangent vector that will be transported along the geodesic, which are also expressed in same basis as `[θₓ, φₓ]`.
+
+We won't draw the transported vector at every point as there would be too many arrows, which is why we select every 100th point only for that purpose with `pt_indices`. Then, `geo_ps_pt` contains points at which the transported vector is tangent to and `geo_Ys` the transported vector at that point, represented in the embedding.
+
+The logarithmic map will be solved between points with parametrization `bvp_a1` and `bvp_a2` in chart `bvp_i`.
+The result is assigned to variable `bvp_sol` and then sampled with time step 0.05. The result of this sampling is converted from parameters in chart `bvp_i` to point in the embedding and stored in `geo_r`.
+
+
+```{julia}
+function solve_for(p0x, X_p0x, Y_transp, T)
+    p = [Manifolds._torus_param(M, p0x...)...]
+    i_p0x = Manifolds.get_chart_index(M, A, p)
+    p_exp = Manifolds.solve_chart_parallel_transport_ode(
+        M,
+        [0.0, 0.0],
+        X_p0x,
+        A,
+        i_p0x,
+        Y_transp;
+        final_time=T,
+    )
+    return p_exp
+end;
+```
+
+
+### Solving parallel Transport ODE
+
+We set the end time `t_end` and time step `dt`.
+
+```{julia}
+t_end = 2.0
+dt = 1e-1
+```
+
+We also parametrise the start point and direction.
+
+```{julia}
+θₚ = π/10
+φₚ = -π/4
+θₓ = π/2
+φₓ = 0.7
+θy = 0.2
+φy = -0.1
+
+geo = solve_for([θₚ, φₚ], [θₓ, φₓ], [θy, φy], t_end)(0.0:dt:t_end);
+# geo_ps = [Point3f(s[1]) for s in geo]
+# pt_indices = 1:div(length(geo), 10):length(geo)
+# geo_ps_pt = [Point3f(s[1]) for s in geo[pt_indices]]
+# geo_Ys = [Point3f(s[3]) for s in geo[pt_indices]]
+
+# ax1, fig1 = torus_figure()
+# arrows!(ax1, geo_ps_pt, geo_Ys, linewidth=0.05, color=:blue)
+# lines!(geo_ps; linewidth=4.0, color=:green)
+# fig1
+```
+
+![fig-pt](working-in-charts/working-in-charts-transport.png)
+
+### Solving the logairthmic map ODE
+
+```{julia}
+θ₁=π/2
+φ₁=-1.0
+θ₂=-π/8
+φ₂=π/2
+
+bvp_i = (0, 0)
+bvp_a1 = [θ₁, φ₁]
+bvp_a2 = [θ₂, φ₂]
+bvp_sol = Manifolds.solve_chart_log_bvp(M, bvp_a1, bvp_a2, A, bvp_i);
+# geo_r = [Point3f(get_point(M, A, bvp_i, p[1:2])) for p in bvp_sol(0.0:0.05:1.0)]
+
+# ax2, fig2 = torus_figure()
+# lines!(geo_r; linewidth=4.0, color=:green)
+# fig2
+```
+
+![fig-geodesic](working-in-charts/working-in-charts-geodesic.png)
+
+An interactive Pluto version of this tutorial is available in file [`tutorials/working-in-charts.jl`](https://github.com/JuliaManifolds/Manifolds.jl/blob/616855447996fb1ee7dfb2a779341b962a1323f8/tutorials/working-in-charts.jl).