diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 24fa896..5042ffd 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -3,6 +3,7 @@ name: CI on: push: branches: [main] + tags: ['v*'] # a release tag builds the versioned docs pull_request: jobs: @@ -23,3 +24,27 @@ jobs: # pyparity/gen_fixtures.py against the Python reference package), so the # test job needs no Python toolchain. - uses: julia-actions/julia-runtest@v1 + + docs: + runs-on: ubuntu-latest + # `deploydocs` pushes the built site to the gh-pages branch; on a fork PR the + # token is read-only and it skips the push and just builds (and doctests). + permissions: + contents: write + pull-requests: read + statuses: write + steps: + - uses: actions/checkout@v4 + - uses: julia-actions/setup-julia@v2 + with: + version: '1' + - uses: julia-actions/cache@v2 + - name: Build docs (doctests included) + run: | + julia --project=docs -e ' + using Pkg + Pkg.develop(PackageSpec(path=pwd())) + Pkg.instantiate()' + julia --project=docs docs/make.jl + env: + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} diff --git a/.gitignore b/.gitignore index 5505293..206cb78 100644 --- a/.gitignore +++ b/.gitignore @@ -6,3 +6,4 @@ Manifest.toml .DS_Store __pycache__/ *.pyc +docs/Manifest.toml diff --git a/PORTING_PLAN.md b/PORTING_PLAN.md index dcae933..5b03147 100644 --- a/PORTING_PLAN.md +++ b/PORTING_PLAN.md @@ -412,7 +412,7 @@ sets `self.xp = numpy`) is excluded. **Upstream bugs.** Four defects found in the Python reference over the course of the port are catalogued — with reproductions, root causes and patches — in -[`docs/upstream-bugs.md`](docs/upstream-bugs.md). +[`docs/src/upstream-bugs.md`](docs/src/upstream-bugs.md). **Harness conventions established** (also in `README.md`): 1-based indexing with the `v .+ 1` parity convention; `NPZ.jl` cannot read zero-element arrays, so the diff --git a/Project.toml b/Project.toml index d471639..8a53c3a 100644 --- a/Project.toml +++ b/Project.toml @@ -30,8 +30,9 @@ julia = "1.10" [extras] CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" NPZ = "15e1cf62-19b3-5cfa-8e77-841668bca605" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Test", "NPZ", "CairoMakie"] +test = ["Test", "NPZ", "CairoMakie", "Documenter"] diff --git a/README.md b/README.md index 3fe0742..3e78193 100644 --- a/README.md +++ b/README.md @@ -32,7 +32,7 @@ the heat diffusion of the data — no mesh, no charts, no prescribed metric. ```julia using Pkg -Pkg.add(url="https://github.com/wardjm/DiffusionGeometryJ.jl") +Pkg.add(url="https://github.com/wardjm/DiffusionGeometryJ") ``` Plotting is provided through a package extension that loads only when a @@ -119,6 +119,25 @@ wedge(α, β) # wedge product, degree k₁ + k₂ — also spelled α ^ f * α # pointwise product with a function (so is wedge(f, α)) ``` +Read one back with the accessors — a tensor stores *coefficients* in the diffusion +basis, not pointwise values, and `to_pointwise_basis` is what evaluates it: + +```julia +coeffs(X) # the coefficient array +space(X) # the tensor space it lives in +geometry(X) # the DiffusionGeometry behind it +to_pointwise_basis(X) # values at the sample points, (batch…, n·C) — the inverse + # of the dg_* factories, up to the basis truncation +``` + +Tensors, spaces and the geometry itself print as a one-line summary rather than +dumping the eigenbasis: + +``` +VectorField(space=VectorFieldSpace(dim=48), shape=(48,), batch_shape=()) +DiffusionGeometry(n=120, ambient_dim=3, n_function_basis=16, n_coefficients=16) +``` + The learned metric gives you geometry: ```julia @@ -128,6 +147,24 @@ l2_norm(X) # global L² norm (or l2_norm(dg, X)) pointwise_norm(X) # pointwise norm (n,) (or pointwise_norm(dg, X)) ``` +### Batches + +Leading axes are *batch* dimensions — many fields carried in one tensor: + +```julia +F = dg_function(dg, ones(3, 60)) # three functions at once +batch_shape(F) # (3,) +``` + +A plain numeric array multiplying a tensor is a bag of **batch-wise scalars**, never a +coefficient vector, and `.*` / `./` are exact synonyms for `*` / `/` (no other +broadcast is defined). That is what makes a spectral filter a one-liner: + +```julia +evals, evecs = spectrum(laplacian(dg, 0)) +smoothed = exp.(-evals) .* evecs # damp each mode by e^(-λ) +``` + ## Differential operators Operator accessors take the geometry and return a `LinearOperator` or @@ -257,7 +294,7 @@ dgplot(T) # ellipsoids for a (0,2)-tensor `dgplot!` adds to an existing axis. Specialised recipes (`dgscatter`, `dgquiver`, `dg2form`, `dg3form`, `dgellipsoids`, `dgeiglines`, `dgtangentplanes`) are available directly. Animate a time-evolving field with `dganimate`. See -[`docs/plotting.md`](docs/plotting.md). +[`docs/src/plotting.md`](docs/src/plotting.md). ## Conventions @@ -265,6 +302,12 @@ available directly. Animate a time-evolving field with `dganimate`. See in keeping with Julia. Point-cloud data is `n × d` (rows are points). - Arrays keep a leading point axis `p`, so a field's coefficients are laid out per point. +- **Batch axes broadcast from the right** (numpy's rule, *not* Julia's): they are + matched from the trailing axis, with length-1 axes stretching. A `(3, 4)` batch and + a `(4,)` batch combine to `(3, 4)`. +- A tensor holds **coefficients in the diffusion basis, not pointwise values** — + `to_pointwise_basis` evaluates it, and `dg_function` and friends project the other + way. This is why a numeric array scaling a tensor is read as batch-wise scalars. - **`*` on two `Form`s is the tensor product**, following Python; the wedge is `wedge(α, β)` or `α ^ β`. Julia's `^` binds *tighter* than `+` and `*` where Python's binds looser than both, so a Python expression that leans on that @@ -302,6 +345,28 @@ One caveat: a fresh Julia process pays a one-time JIT compilation cost (~10 s) o the first `from_point_cloud` call, which the interpreted Python does not; every call after that is warm. +## Documentation + +Every exported function carries a docstring with runnable examples. Build the +[Documenter](https://documenter.juliadocs.org) site with: + +```bash +julia --project=docs -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()' +julia --project=docs docs/make.jl # writes docs/build, and runs every doctest +``` + +The examples are doctests: they assert real values (the Laplacian spectrum of the +circle, the adjointness of `δ` against `d`, the antisymmetry of the wedge) and run as +part of the test suite — both the docstring examples and those in the guide pages — so +they cannot drift from the code. They share a preamble — `dg` (60 points on a circle), +`dg3` (200 on a sphere), and `f = cos θ` — defined in `docs/doctest_setup.jl`. + +The site also carries a [Conventions](docs/src/conventions.md) page (coefficients vs +values, batch broadcasting, weak vs strong forms) and a catalogue of the +[upstream bugs](docs/src/upstream-bugs.md) found in the Python reference during the +port. CI builds the site with `checkdocs=:exports`, so an undocumented export or a +broken example fails the build. + ## Testing ```bash diff --git a/docs/Project.toml b/docs/Project.toml new file mode 100644 index 0000000..5ec6819 --- /dev/null +++ b/docs/Project.toml @@ -0,0 +1,6 @@ +[deps] +DiffusionGeometryJ = "f8140dce-05b1-4e22-80fe-cf341a3061cb" +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" + +[compat] +Documenter = "1" diff --git a/docs/doctest_setup.jl b/docs/doctest_setup.jl new file mode 100644 index 0000000..ae4a53b --- /dev/null +++ b/docs/doctest_setup.jl @@ -0,0 +1,36 @@ +# The doctest preamble, shared by `docs/make.jl` and `test/doctests.jl`. +# +# Every `jldoctest` block in the package runs with this code already evaluated, so +# examples can reach for the running example geometries without rebuilding them. +# Keep it small: whatever appears here is invisible in the rendered docs, so an +# example that leans on more than `dg`, `dg3`, `f` and their point clouds becomes +# hard to read. It is documented for the reader in `docs/src/index.md`. + +const DOCTEST_SETUP = quote + using DiffusionGeometryJ + using LinearAlgebra + + # Nothing unexported is imported here on purpose. The internal helpers documented on + # the "Utilities" page (`np_reshape`, the batch helpers, `_from_pointwise_basis`) are + # *not* exported, so their examples spell them `DiffusionGeometryJ.np_reshape(...)` — + # as a reader would have to. + + # 60 points on the unit circle in ℝ² — the running example. Its intrinsic + # dimension is 1, but the carré du champ lives in the d = 2 ambient basis. + θ = range(0, 2π; length=61)[1:60] + circle = [cos.(θ) sin.(θ)] + dg = from_point_cloud(circle; knn_kernel=16, n_function_basis=8) + + # The coordinate function x₁ = cos θ on that circle. + f = dg_function(dg, cos.(θ)) + + # 200 points on the unit 2-sphere in ℝ³ (Fibonacci lattice), for the examples + # that need a surface: 2-forms, curvature, the Hodge star. + sphere = let n = 200, ϕ = (1 + sqrt(5)) / 2 + z = [1 - 2 * (i - 0.5) / n for i in 1:n] + r = sqrt.(max.(1 .- z .^ 2, 0.0)) + ψ = [2π * (i - 1) / ϕ for i in 1:n] + [r .* cos.(ψ) r .* sin.(ψ) z] + end + dg3 = from_point_cloud(sphere; knn_kernel=20, n_function_basis=12) +end diff --git a/docs/make.jl b/docs/make.jl new file mode 100644 index 0000000..1c0a361 --- /dev/null +++ b/docs/make.jl @@ -0,0 +1,35 @@ +using Documenter +using DiffusionGeometryJ + +include(joinpath(@__DIR__, "doctest_setup.jl")) +DocMeta.setdocmeta!(DiffusionGeometryJ, :DocTestSetup, DOCTEST_SETUP; recursive=true) + +makedocs(; + sitename="DiffusionGeometryJ", + authors="Jeffrey Ward", + modules=[DiffusionGeometryJ], + doctest=true, + checkdocs=:exports, # every exported name must appear in an @docs block + # No `warnonly`: an undocumented export, a broken doctest or a dead cross-reference + # must fail the build, not warn. + format=Documenter.HTML(; + canonical="https://wardjm.github.io/DiffusionGeometryJ", + prettyurls=get(ENV, "CI", "false") == "true", + ), + pages=[ + "Home" => "index.md", + "Conventions" => "conventions.md", + "API" => [ + "Building a geometry" => "api/geometry.md", + "Tensor fields" => "api/tensors.md", + "Operators" => "api/operators.md", + "Numerical methods" => "api/methods.md", + "Diffusion core" => "api/diffusion.md", + "Utilities" => "api/utils.md", + ], + "Plotting" => "plotting.md", + "Upstream bugs" => "upstream-bugs.md", + ], +) + +deploydocs(; repo="github.com/wardjm/DiffusionGeometryJ", devbranch="main") diff --git a/docs/src/api/diffusion.md b/docs/src/api/diffusion.md new file mode 100644 index 0000000..af09161 --- /dev/null +++ b/docs/src/api/diffusion.md @@ -0,0 +1,50 @@ +```@meta +CurrentModule = DiffusionGeometryJ +``` + +# Diffusion core + +The pipeline that turns a point cloud into a diffusion: kNN graph → Markov chain → +symmetric kernel → eigenfunction basis → carré du champ. [`from_point_cloud`](@ref) +runs all of it; these are the stages, for when you want to intervene. + +## The diffusion process + +```@docs +knn_graph +compute_local_bandwidths +tune_kernel +markov_chain +build_symmetric_kernel_matrix +compute_eigenfunction_basis +``` + +## Carré du champ + +```@docs +carre_du_champ_knn +carre_du_champ_graph +gamma_compound +gamma_02 +gamma_02_sym +``` + +## Markov triples + +```@docs +MarkovTriple +ImmersedMarkovTriple +cdc +regularise +immersed_triple_from_point_cloud +immersed_triple_from_knn_kernel +immersed_triple_from_graph_kernel +immersed_triple_from_edges +``` + +## Regularisation + +```@docs +regularise_diffusion +regularise_bandlimit +``` diff --git a/docs/src/api/geometry.md b/docs/src/api/geometry.md new file mode 100644 index 0000000..43401d4 --- /dev/null +++ b/docs/src/api/geometry.md @@ -0,0 +1,91 @@ +```@meta +CurrentModule = DiffusionGeometryJ +DocTestSetup = Main.DOCTEST_SETUP +``` + +# Building a geometry + +A `DiffusionGeometry` is built from a point cloud (or a kernel, or a graph), and from +then on it is the thing everything else hangs off: it owns the eigenfunction basis, the +measure, the carré du champ cache, the tensor spaces, and every differential operator. + +```jldoctest +julia> dg = from_point_cloud(circle; knn_kernel=16, n_function_basis=8) +DiffusionGeometry(n=60, ambient_dim=2, n_function_basis=8, n_coefficients=8) +``` + +## The type + +```@docs +DiffusionGeometry +``` + +## Constructors + +Six ways in, from the most convenient to the most raw. + +```@docs +from_point_cloud +from_knn_graph +from_knn_kernel +from_graph_kernel +from_edges +from_sparse_matrix +``` + +## Accessors + +```@docs +npoints +ambient_dim +n_coefficients +n_function_basis +function_basis +measure +immersion_coords +``` + +## Tensor spaces + +Each is cached on the geometry, so repeated calls return the same object. + +```@docs +function_space +vector_field_space +form_space +tensor02_space +tensor02sym_space +``` + +## Building tensors + +```@docs +dg_function +dg_vector_field +dg_form +dg_tensor02 +dg_tensor02sym +``` + +## The metric + +```@docs +g +pointwise_norm +inner +l2_norm +``` + +## The γ cache + +The carré du champ contractions are the expensive part of the package, so they are +memoised on the geometry as `dg.cache`. + +```@docs +GammaCache +gamma_coords +gamma_functions +gamma_mixed +gamma_ambient +gamma_coords_compound +``` diff --git a/docs/src/api/methods.md b/docs/src/api/methods.md new file mode 100644 index 0000000..11873ed --- /dev/null +++ b/docs/src/api/methods.md @@ -0,0 +1,38 @@ +```@meta +CurrentModule = DiffusionGeometryJ +``` + +# Numerical methods + +## PDEs + +```@docs +solve_differential_operator +``` + +## Geodesic distances + +```@docs +geodesic_distances_function +``` + +## Topology + +Betti numbers come from the spectral gap of a penalised Hodge operator — the discrete +`d` does not square to zero, so there is no exact kernel to count. Build the geometry +with the full coefficient basis, and review the spectrum when the count matters. + +```@docs +BettiSpectrum +betti_spectrum +betti_spectra +betti_number +betti_numbers +betti_gap +``` + +## The Hodge star (for 2-forms) + +```@docs +hodge_star_2_form +``` diff --git a/docs/src/api/operators.md b/docs/src/api/operators.md new file mode 100644 index 0000000..1aa7870 --- /dev/null +++ b/docs/src/api/operators.md @@ -0,0 +1,109 @@ +```@meta +CurrentModule = DiffusionGeometryJ +DocTestSetup = Main.DOCTEST_SETUP +``` + +# Operators + +Every differential operator is a [`LinearOperator`](@ref) (or, for the Lie bracket, a +[`BilinearOperator`](@ref)) between two tensor spaces. They are memoised on the +geometry, compose with `∘`, take adjoints with `'`, and apply by being called. + +```jldoctest +julia> laplacian(dg, 0) +LinearOperator(domain=FunctionSpace(dim=8), codomain=FunctionSpace(dim=8), shape=(8, 8)) + +julia> laplacian(f) +ScalarFunction(space=FunctionSpace(dim=8), shape=(8,), batch_shape=()) +``` + +## Differential operators + +Each takes the geometry (giving the operator) or a tensor (giving the result). + +```@docs +grad +d +codifferential +divergence +up_laplacian +down_laplacian +laplacian +hessian +levi_civita +lie_bracket +``` + +## Curvature + +```@docs +riemann_curvature +sectional_curvature +``` + +## Tensors acting as operators + +```@docs +vf_operator +t02_operator +wedge_operator +``` + +## Hodge theory + +```@docs +hodge_decomposition +``` + +## Ambient representations + +```@docs +to_ambient +vector_field_to_quiver +vector_field_from_reconstruction +``` + +## The operator types + +```@docs +LinearOperator +BilinearOperator +weak +matrix +strong +spectrum +inverse +is_self_adjoint +component_shape +partial_apply +full_apply +transpose_operator +identity_operator +zero_operator +real_if_close +``` + +## Block operators + +```@docs +block +hstack +vstack +``` + +## Weak-form builders + +The assemblers underneath the operators. Each returns the weak matrix the corresponding +operator carries; you only need them to build an operator by hand. + +```@docs +derivative_weak +up_delta_weak +hessian_functions +hessian_coords +hessian_02_weak +hessian_02_sym_weak +levi_civita_02_weak +lie_bracket_weak +metric +``` diff --git a/docs/src/api/tensors.md b/docs/src/api/tensors.md new file mode 100644 index 0000000..6a6a531 --- /dev/null +++ b/docs/src/api/tensors.md @@ -0,0 +1,84 @@ +```@meta +CurrentModule = DiffusionGeometryJ +``` + +# Tensor fields + +A tensor is a coefficient vector plus the space it lives in — see +[Conventions](@ref) for what that means and how batch axes behave. + +## Abstract types + +```@docs +AbstractTensor +AbstractTensorSpace +coeffs +space +geometry +batch_shape +space_degree +degree +``` + +## Spaces + +```@docs +FunctionSpace +VectorFieldSpace +FormSpace +Tensor02Space +Tensor02SymSpace +DirectSumSpace +``` + +### Space properties + +```@docs +component_dim +space_dim +cdc_components +gram +gram_inv +orthonormal_basis +metric_tensor +metric_apply +``` + +## Elements + +```@docs +ScalarFunction +VectorField +Form +Tensor02 +Tensor02Sym +DirectSumElement +``` + +### Building and unpacking + +```@docs +from_pointwise +to_pointwise_basis +wrap +Base.zeros(::AbstractTensorSpace, ::Tuple) +``` + +## Algebra + +```@docs +wedge +flat +sharp +symmetrise +full_tensor +transpose_tensor +``` + +## Direct sums + +```@docs +pack +unpack +split_coeffs +``` diff --git a/docs/src/api/utils.md b/docs/src/api/utils.md new file mode 100644 index 0000000..9e2629e --- /dev/null +++ b/docs/src/api/utils.md @@ -0,0 +1,51 @@ +```@meta +CurrentModule = DiffusionGeometryJ +``` + +# Utilities + +Combinatorics of the wedge and symmetric bases, the numpy-compatible reshape and batch +rules, and the basis conversions. Mostly internal, but documented because the +coefficient layout they define is visible in every array the package returns. + +## Basis combinatorics + +```@docs +get_wedge_basis_indices +get_symmetric_basis_indices +get_wedge_product_indices +kp1_children_and_signs +lex_rank +permutations_with_signs +expand_symmetric_tensor_coeffs +symmetrise_tensor_coeffs +``` + +## Row-major reshape + +```@docs +np_reshape +``` + +## Batch shapes + +```@docs +broadcast_batch_shape +compatible_batches +pad_batch_dims +align_batch_pair +broadcast_batch_to +broadcast_flatten_batch +infer_batch_shape +flatten_batch_dims +restore_batch_dims +``` + +## Basis conversions + +```@docs +_from_pointwise_basis +_to_pointwise_basis +_is_orthonormal +_gram_spectrum +``` diff --git a/docs/src/conventions.md b/docs/src/conventions.md new file mode 100644 index 0000000..15119e8 --- /dev/null +++ b/docs/src/conventions.md @@ -0,0 +1,86 @@ +```@meta +CurrentModule = DiffusionGeometryJ +DocTestSetup = Main.DOCTEST_SETUP +``` + +# Conventions + +Four things to know before the API. Everything else follows from them. + +## Coefficients, not values + +A tensor is *not* an array of values at the sample points. It is a vector of +coefficients in the basis `{φ_i e_a}` — the eigenfunctions `φ_i` of the diffusion, +times the component frame `e_a` of its space. The coefficient axis is always the +**last** axis of the array. + +```jldoctest +julia> f # 8 coefficients, not 60 values +ScalarFunction(space=FunctionSpace(dim=8), shape=(8,), batch_shape=()) + +julia> length(to_pointwise_basis(f)) # the values, when you want them +60 +``` + +This is why elementwise arithmetic on a tensor is *not* offered: `abs.(ω)` would act on +coefficients, which is not the operation it looks like. Round-trip through +[`to_pointwise_basis`](@ref) and [`dg_function`](@ref) if you need to touch values. + +The two truncation levels matter. `n_function_basis` is how many eigenfunctions the +diffusion resolves; `n_coefficients ≤ n_function_basis` is how many a *tensor* is +expanded in. Functions always use the full basis, everything else uses +`n_coefficients`, and the topology functions need them equal. + +## Batch axes broadcast from the right (numpy, not Julia) + +Any leading axes are batch dimensions. They broadcast numpy-style — matched from the +right, with length-1 axes stretching — which is *not* Julia's rule, so the package +carries its own [`broadcast_batch_shape`](@ref). + +```jldoctest +julia> batch_shape(dg_function(dg, ones(3, 60))) # three functions in one tensor +(3,) + +julia> DiffusionGeometryJ.broadcast_batch_shape((3, 4), (4,)) # right-aligned +(3, 4) +``` + +A plain numeric array multiplying a tensor is a bag of *batch-wise scalars*, never a +coefficient vector — that is what makes `exp.(-evals) .* evecs` a spectral filter. +`.*` and `./` are exact synonyms for `*` and `/`; no other broadcast is defined. + +## Everything is 1-based, and stored row-major + +Indices are 1-based throughout: neighbour indices, edge indices, wedge-basis ranks, +the source point of [`geodesic_distances_function`](@ref). The Python reference is +0-based, so a Python array `v` of indices corresponds to `v .+ 1` here. + +Coefficient arrays are flattened **row-major** (C-order, basis index slowest), matching +numpy and the parity fixtures — hence [`np_reshape`](@ref) rather than `Base.reshape` +wherever a multi-index is packed. + +## Weak and strong forms + +Every operator carries two matrices: the *weak* one (the bilinear form `⟨w, Lv⟩`, which +is what the einsums assemble) and the *strong* one (the coefficient map, which is what +applying it multiplies by). They are related by the codomain's [`gram`](@ref) matrix, +and each is derived from the other lazily. + +```jldoctest +julia> L = laplacian(dg, 0); + +julia> weak(L) ≈ gram(function_space(dg)) * matrix(L) +true +``` + +The strong form goes through a *pseudo*-inverse with a spectral cutoff (`dg.rcond`), so +composing operators is not the same computation as assembling the composite directly. +They agree on everything the basis resolves and drift apart on the modes it barely +does — `δd` assembled in one weak form is the trustworthy one, which is why +[`laplacian`](@ref) does not just compose [`codifferential`](@ref) with [`d`](@ref). + +## And one thing that is not true + +`d ∘ d ≠ 0`. The discrete exterior derivative is not a chain complex, so there is no +exact de Rham kernel to read topology from. Betti numbers come from a *spectral gap* in +a penalised Hodge operator instead — see [`betti_spectrum`](@ref). diff --git a/docs/src/index.md b/docs/src/index.md new file mode 100644 index 0000000..7c28df2 --- /dev/null +++ b/docs/src/index.md @@ -0,0 +1,92 @@ +```@meta +CurrentModule = DiffusionGeometryJ +DocTestSetup = Main.DOCTEST_SETUP +``` + +# DiffusionGeometryJ + +Data-driven calculus, geometry, and topology on point clouds — via heat diffusion +and the carré du champ operator. + +This is a Julia port of the +[`DiffusionGeometry`](https://github.com/Iolo-Jones/DiffusionGeometry) Python +package by **Iolo Jones and David Lanners** (*Computing Diffusion Geometry*, 2026). +All of the underlying mathematics, the algorithms, and the reference implementation +are their work; please cite them when using this software. + +Given nothing but a point cloud, the package builds a discrete approximation of the +manifold the points are sampled from and lets you do differential geometry on it: +gradients and Hessians, the exterior derivative and the Hodge Laplacian on +differential forms, lengths and angles under the learned metric, curvature, spectral +PDE solutions, geodesic distances, and Betti numbers read off the harmonic forms. +No mesh, no charts, no prescribed metric. + +## Installation + +```julia +using Pkg +Pkg.add(url="https://github.com/wardjm/DiffusionGeometryJ") +``` + +Plotting comes from a package extension that loads when a +[Makie](https://docs.makie.org) backend is present (`Pkg.add("GLMakie")`); see +[Plotting](@ref). + +## Quick start + +Sixty points on the unit circle, and the calculus that follows from them. The +Laplacian eigenvalues of the unit circle are `0, 1, 1, 4, 4, 9, 9, …`, and the +geometry recovers them from the samples alone: + +```jldoctest +julia> θ = range(0, 2π; length=61)[1:60]; + +julia> circle = [cos.(θ) sin.(θ)]; # 60 × 2 point cloud + +julia> dg = from_point_cloud(circle; knn_kernel=16, n_function_basis=8) +DiffusionGeometry(n=60, ambient_dim=2, n_function_basis=8, n_coefficients=8) + +julia> round.(abs.(spectrum(laplacian(dg, 0); eigvals_only=true)); digits=2) +8-element Vector{Float64}: + 0.0 + 1.0 + 1.0 + 3.78 + 3.78 + 7.76 + 7.76 + 12.21 + +julia> f = dg_function(dg, cos.(θ)); # the coordinate function x₁ + +julia> round(l2_norm(grad(f)); digits=3) # ∫|∇cos θ|² dθ/2π = 1/2 +0.707 +``` + +## The running example + +Every doctest in this documentation runs with the following already defined, so the +examples can get to the point instead of rebuilding a geometry each time: + +```julia +θ = range(0, 2π; length=61)[1:60] +circle = [cos.(θ) sin.(θ)] # 60 points on the unit circle ⊂ ℝ² +dg = from_point_cloud(circle; knn_kernel=16, n_function_basis=8) +f = dg_function(dg, cos.(θ)) # the coordinate function x₁ on it + +sphere = # 200 points on the unit 2-sphere ⊂ ℝ³ (Fibonacci lattice) +dg3 = from_point_cloud(sphere; knn_kernel=20, n_function_basis=12) +``` + +`dg` is the workhorse: a curve, so its forms only run to degree 2 and its geometry is +cheap. `dg3` is there for the examples that need a surface — 2-forms, curvature, the +Hodge star. The exact definitions live in `docs/doctest_setup.jl`. + +## Where to go next + + * [Conventions](@ref) — indexing, batch axes, and how a tensor is stored. Read this + before the API. + * [Building a geometry](@ref) — the constructors, the spaces, the metric. + * [Tensor fields](@ref) — functions, vector fields, forms, (0,2)-tensors. + * [Operators](@ref) — `grad`, `d`, the Laplacians, curvature, and the operator algebra. + * [Numerical methods](@ref) — PDEs, geodesics, Betti numbers. diff --git a/docs/plotting.md b/docs/src/plotting.md similarity index 93% rename from docs/plotting.md rename to docs/src/plotting.md index b207d3a..85b5876 100644 --- a/docs/plotting.md +++ b/docs/src/plotting.md @@ -107,3 +107,34 @@ deliberately absent: - **`overpic_labels`**, which emitted LaTeX `\put` coordinates for subplot centres, is a LaTeX-layout helper with no Makie analogue — Makie lays out figures with `GridLayout` and places text with `Label`. + +## API + +```@meta +CurrentModule = DiffusionGeometryJ +``` + +The plot functions have no methods until a Makie backend is loaded, so their examples +are shown rather than run. + +```@docs +dgplot +dgplot! +dgscatter +dgscatter! +dgquiver +dgquiver! +dg2form +dg2form! +dg3form +dg3form! +dgellipsoids +dgellipsoids! +dgeiglines +dgeiglines! +dgtangentplanes +dgtangentplanes! +dganimate +DIVERGING_COLORMAP +CYCLIC_COLORMAP +``` diff --git a/docs/upstream-bugs.md b/docs/src/upstream-bugs.md similarity index 100% rename from docs/upstream-bugs.md rename to docs/src/upstream-bugs.md diff --git a/src/DiffusionGeometryJ.jl b/src/DiffusionGeometryJ.jl index b78d841..6c6f251 100644 --- a/src/DiffusionGeometryJ.jl +++ b/src/DiffusionGeometryJ.jl @@ -59,6 +59,8 @@ export function_space, vector_field_space, form_space, tensor02_space, tensor02s export dg_function, dg_vector_field, dg_form, dg_tensor02, dg_tensor02sym export g, inner, l2_norm, pointwise_norm export AbstractTensor, AbstractTensorSpace +# The tensor/space accessors: the read side of every tensor a user holds. +export coeffs, space, geometry, batch_shape, space_degree export FunctionSpace, VectorFieldSpace, FormSpace, Tensor02Space, Tensor02SymSpace, DirectSumSpace export ScalarFunction, VectorField, Form, Tensor02, Tensor02Sym, DirectSumElement diff --git a/src/core/diffusion/carre_du_champ.jl b/src/core/diffusion/carre_du_champ.jl index ed3cbb0..efe7743 100644 --- a/src/core/diffusion/carre_du_champ.jl +++ b/src/core/diffusion/carre_du_champ.jl @@ -10,6 +10,36 @@ using LinearAlgebra: det, mul! Γ_p(f, h) = (1/2ρ) Cov_p(f, h) over each point's kNN neighbourhood. `f`, `h` have shape `(n, tail...)`; `diffusion_kernel` and `nbr_indices` are `(n, k)` (1-based indices). Result has shape `(n, f_tail..., h_tail...)`. + +The carré du champ is the object the whole package is built on: `Γ(f, h) = ∇f · ∇h` +in the continuum, recovered from the data as a local covariance under the diffusion. + +# Examples +Three points in a ring, each averaging with its successor. The local variance of +`f = (1, 2, 3)` around point 3 is larger because its neighbourhood straddles the +jump from 3 back to 1: + +```jldoctest +julia> f = [1.0, 2.0, 3.0]; + +julia> carre_du_champ_knn(f, f, fill(0.5, 3, 2), [1 2; 2 3; 3 1]) +3-element Vector{Float64}: + 0.125 + 0.125 + 0.5 +``` + +Applied to the immersion coordinates it gives the metric — `Γ(xᵢ, xⱼ)` is the first +fundamental form, one `d × d` matrix per point: + +```jldoctest +julia> kernel, _ = markov_chain(knn_graph(circle, 16)...); + +julia> Γ = carre_du_champ_knn(circle, circle, kernel, knn_graph(circle, 16)[2]); + +julia> size(Γ) +(60, 2, 2) +``` """ function carre_du_champ_knn(f::AbstractArray, h::AbstractArray, diffusion_kernel::AbstractMatrix, @@ -74,6 +104,26 @@ is a length-`num_edges` vector of (row-stochastic) edge weights; `edge_index` is `(2, num_edges)` matrix of **1-based** node indices with row 1 the source `j` and row 2 the target `i` (mirrors Python's `[source; target]`). `f`, `h` have shape `(n, tail...)`; result has shape `(n, f_tail..., h_tail...)`. + +The graph counterpart of [`carre_du_champ_knn`](@ref): same covariance, but the +neighbourhood of `i` is whatever points have an edge *into* `i`. + +# Examples +A 3-cycle with both directions on every edge, so node 2's neighbours are 1 and 3 — +the values `1` and `3` around it, whose spread is the largest of the three: + +```jldoctest +julia> edges = [1 2 2 3 3 1; # source + 2 1 3 2 1 3]; # target + +julia> f = [1.0, 2.0, 3.0]; + +julia> carre_du_champ_graph(f, f, fill(0.5, 6), edges) +3-element Vector{Float64}: + 0.125 + 0.5 + 0.125 +``` """ function carre_du_champ_graph(f::AbstractArray, h::AbstractArray, diffusion_kernel::AbstractVector, @@ -153,6 +203,27 @@ end k-th compound submatrices and their determinants per point. `gamma_coords` is `(n, d, d)`. `dets[p, a, b] = det(Γ_p[J_a, J_b])` for k-combinations `J_a, J_b` in lex order; `submatrices` is `(n, Ck, Ck, k, k)`. + +The compound determinants are the metric on k-forms — `⟨dx_I, dx_J⟩ = det Γ[I, J]`, +the Gram determinant of the induced inner product on Λᵏ. They are exactly what +[`FormSpace`](@ref) uses as its `cdc_components`. + +# Examples +For `Γ = diag(2, 3)` the only 2-form basis element is `dx₁∧dx₂`, and its squared +length is the determinant `2 · 3 = 6`: + +```jldoctest +julia> Γ = reshape([2.0 0.0; 0.0 3.0], 1, 2, 2); # one point, d = 2 + +julia> gamma_compound(Γ, 2)[2][1, :, :] +1×1 Matrix{Float64}: + 6.0 + +julia> gamma_compound(Γ, 1)[2][1, :, :] # k = 1 is Γ itself +2×2 Matrix{Float64}: + 2.0 0.0 + 0.0 3.0 +``` """ function gamma_compound(gamma_coords::AbstractArray{T,3}, k::Integer) where {T} n, d, _ = size(gamma_coords) @@ -191,6 +262,21 @@ end Expanded Γ for full (0,2)-tensors: `Γ_{02}[p, j⊗k, J⊗K] = Γ_p(x_j,x_J) Γ_p(x_k,x_K)`, flattened row-major. `gamma_coords` is `(n, d, d)`; result is `(n, d², d²)`. + +The metric on `T*M ⊗ T*M` — the Kronecker square of Γ, and the `cdc_components` of +[`Tensor02Space`](@ref). + +# Examples +```jldoctest +julia> Γ = reshape([1.0 2.0; 3.0 4.0], 1, 2, 2); + +julia> gamma_02(Γ)[1, :, :] # Γ ⊗ Γ, in the row-major dx_j ⊗ dx_k basis +4×4 Matrix{Float64}: + 1.0 2.0 2.0 4.0 + 3.0 4.0 6.0 8.0 + 3.0 6.0 4.0 8.0 + 9.0 12.0 12.0 16.0 +``` """ function gamma_02(gamma_coords::AbstractArray{T,3}) where {T} n, d, _ = size(gamma_coords) @@ -213,6 +299,25 @@ end Symmetric-part Γ for symmetric (0,2)-tensors, with off-diagonal basis elements weighted (×2 per off-diagonal index). `gamma_coords` is `(n, d, d)`; result is `(n, d_sym, d_sym)` with `d_sym = d(d+1)/2`. + +The `cdc_components` of [`Tensor02SymSpace`](@ref). The ×2 weights are what make the +`d_sym` stored components carry the norm of all `d²` components of the tensor they +stand for. + +# Examples +With `Γ = diag(1, 2)`, the basis is `(dx₁⊙dx₁, dx₁⊙dx₂, dx₂⊙dx₂)`. The off-diagonal +element gets `½(Γ₁₁Γ₂₂ + Γ₁₂Γ₂₁) · 2 · 2 = 4`, matching the two slots `dx₁⊗dx₂` and +`dx₂⊗dx₁` it represents: + +```jldoctest +julia> Γ = reshape([1.0 0.0; 0.0 2.0], 1, 2, 2); + +julia> gamma_02_sym(Γ)[1, :, :] +3×3 Matrix{Float64}: + 1.0 0.0 0.0 + 0.0 4.0 0.0 + 0.0 0.0 4.0 +``` """ function gamma_02_sym(gamma_coords::AbstractArray{T,3}) where {T} n, d, _ = size(gamma_coords) diff --git a/src/core/diffusion/diffusion_process.jl b/src/core/diffusion/diffusion_process.jl index ff1e113..d7baea4 100644 --- a/src/core/diffusion/diffusion_process.jl +++ b/src/core/diffusion/diffusion_process.jl @@ -15,6 +15,28 @@ using Random: MersenneTwister k-nearest-neighbour graph. `data_matrix` is `(n, d)`. Returns `(n, k)` distances and 1-based neighbour indices, sorted by increasing distance (each point's first neighbour is itself). + +Step one of [`from_point_cloud`](@ref); its output feeds [`markov_chain`](@ref). + +# Examples +Three points on a line. Each is its own nearest neighbour (distance 0), and the +next column is the closest other point: + +```jldoctest +julia> nbr_distances, nbr_indices = knn_graph([0.0 0.0; 1.0 0.0; 3.0 0.0], 2); + +julia> nbr_distances +3×2 Matrix{Float64}: + 0.0 1.0 + 0.0 1.0 + 0.0 2.0 + +julia> nbr_indices # 1-based; column 1 is the point itself +3×2 Matrix{Int64}: + 1 2 + 2 1 + 3 2 +``` """ function knn_graph(data_matrix::AbstractMatrix, knn_kernel::Integer=32) n = size(data_matrix, 1) @@ -36,6 +58,28 @@ end Local bandwidths from neighbour distances. `k_bandwidth` matches the Python (0-based) argument: `"k"` uses the `k_bandwidth`-th neighbour distance; `"l1"` and `"l2"` average over the first `k_bandwidth` neighbours excluding self. + +The bandwidth is the local length scale of the kernel — it is what lets the +diffusion adapt to a point cloud of varying density. + +# Examples +Two points whose sorted neighbour distances are `(0, 1, 2)` and `(0, 1, 3)`. The +`"k"` rule reads off the `k_bandwidth`-th (0-based) column; `"l2"` takes the +root-mean-square over the neighbours before it, excluding self: + +```jldoctest +julia> nbr_distances = [0.0 1.0 2.0; 0.0 1.0 3.0]; + +julia> compute_local_bandwidths(nbr_distances, 2, "k") # the 2nd neighbour +2-element Vector{Float64}: + 2.0 + 3.0 + +julia> compute_local_bandwidths(nbr_distances, 3, "l2") # √((1² + 2²)/2), √((1² + 3²)/2) +2-element Vector{Float64}: + 1.5811388300841898 + 2.23606797749979 +``` """ function compute_local_bandwidths(nbr_distances::AbstractMatrix, k_bandwidth::Integer=8, bandwidth_type::AbstractString="l2") @@ -56,6 +100,25 @@ end Select the kernel scale `epsilon` maximising the log-log slope of the average kernel value, and the corresponding intrinsic dimension estimate `dim`. + +The slope of `log Σ exp(-dᵢⱼ²/ε)` against `log ε` peaks at the scale where the +kernel sees the manifold rather than the noise below it or the whole cloud above it, +and twice that peak slope estimates the *intrinsic* dimension. + +# Examples +The circle is one-dimensional, and the tuner says so — it never sees the ℝ² the +points are embedded in: + +```jldoctest +julia> nbr_distances, _ = knn_graph(circle, 16); + +julia> epsilons = 2 .^ collect(-10:0.25:9.75); + +julia> epsilon, dim = tune_kernel(nbr_distances .^ 2, epsilons); + +julia> round(dim; digits=2) +1.01 +``` """ function tune_kernel(kernel_entries::AbstractMatrix, epsilons::AbstractVector) ne = length(epsilons) @@ -80,6 +143,38 @@ end Row-stochastic diffusion kernel `(n, k)` and local bandwidths `ρ` `(n,)` via the variable-bandwidth ("Diffusion Maps") construction with the α-normalisation. + +`knn_bandwidth` neighbours set the local length scale, so the kNN graph must be +wider than it: `size(nbr_distances, 2) > knn_bandwidth`. + +# Examples +The kernel is a Markov chain — every row is a probability distribution over the +point's neighbours: + +```jldoctest +julia> kernel, bandwidths = markov_chain(knn_graph(circle, 16)...); + +julia> size(kernel) +(60, 16) + +julia> all(≈(1.0), sum(kernel; dims=2)) # row-stochastic +true + +julia> round.(kernel[1, 1:3]; digits=4) # self first, then the two nearest +3-element Vector{Float64}: + 0.2157 + 0.1861 + 0.1861 +``` + +The circle is sampled uniformly, so every point gets the same bandwidth: + +```jldoctest +julia> _, bandwidths = markov_chain(knn_graph(circle, 16)...); + +julia> round(maximum(bandwidths) - minimum(bandwidths); digits=6) +0.0 +``` """ function markov_chain(nbr_distances::AbstractMatrix, nbr_indices::AbstractMatrix{<:Integer}; c::Real=0, bandwidth_variability::Real=-0.5, knn_bandwidth::Integer=8) @@ -127,6 +222,37 @@ end Symmetric sparse kernel `K = (A + Aᵀ)/2` from the row-stochastic kernel, and its column sums (the unnormalised measure μ). `nbr_indices` is 1-based. + +Symmetrising is what makes the eigenproblem in [`compute_eigenfunction_basis`](@ref) +self-adjoint, and the column sums are the measure the whole geometry integrates +against. + +# Examples +Two points, each diffusing half its mass to the other: + +```jldoctest +julia> K, row_sums = build_symmetric_kernel_matrix([0.5 0.5; 0.5 0.5], [1 2; 2 1]); + +julia> Matrix(K) +2×2 Matrix{Float64}: + 0.5 0.5 + 0.5 0.5 + +julia> row_sums +2-element Vector{Float64}: + 1.0 + 1.0 +``` + +On the circle it is a 60×60 sparse matrix, symmetric by construction: + +```jldoctest +julia> K, _ = build_symmetric_kernel_matrix(markov_chain(knn_graph(circle, 16)...)[1], + knn_graph(circle, 16)[2]); + +julia> size(K), Matrix(K) == Matrix(K)' +((60, 60), true) +``` """ function build_symmetric_kernel_matrix(diffusion_kernel::AbstractMatrix, nbr_indices::AbstractMatrix{<:Integer}) @@ -153,6 +279,25 @@ end Diffusion-maps eigenfunction basis `{φ_i}` `(n, n0)`, ordered by increasing complexity with `φ_0 ≡ 1`. Uses Arpack for `n0 < n`, otherwise a dense Hermitian eigendecomposition. + +This is the basis every coefficient in the package is expressed in. The columns are +ordered by decreasing kernel eigenvalue — the low-frequency modes of the data come +first — and the *sign* of each column beyond `φ_0` is whatever the eigensolver +returns, so never depend on it. + +# Examples +```jldoctest +julia> K, row_sums = build_symmetric_kernel_matrix(markov_chain(knn_graph(circle, 16)...)[1], + knn_graph(circle, 16)[2]); + +julia> u = compute_eigenfunction_basis(K, row_sums; n0=4); + +julia> size(u) +(60, 4) + +julia> all(≈(1.0), u[:, 1]) # φ₀ is the constant function, normalised to 1 +true +``` """ function compute_eigenfunction_basis(symmetric_kernel_matrix, row_sums; n0::Integer=40) D = Diagonal(row_sums .^ (-1 / 2)) diff --git a/src/core/diffusion/markov_triples.jl b/src/core/diffusion/markov_triples.jl index 7f4e59d..3c89034 100644 --- a/src/core/diffusion/markov_triples.jl +++ b/src/core/diffusion/markov_triples.jl @@ -10,6 +10,19 @@ using LinearAlgebra: I Coefficient-form Markov triple `(M, μ, Γ)`: coefficient functions `function_basis` `(n, n0)`, measure `measure` `(n,)`, a carré du champ callback `cdc(f, h)`, and an optional `regularise` map (defaults to `identity`). + +This is the whole of the "diffusion" input to the geometry: a basis to expand in, a +measure to integrate against, and a carré du champ to differentiate with. Everything +else in the package is derived from it. + +# Examples +```jldoctest +julia> t = MarkovTriple(ones(3, 1), fill(1/3, 3), (f, h) -> zeros(3)) +MarkovTriple(n=3, n_function_basis=1) + +julia> t.n, t.n_function_basis +(3, 1) +``` """ struct MarkovTriple{C,R} function_basis::Matrix{Float64} @@ -27,7 +40,41 @@ function MarkovTriple(function_basis::AbstractMatrix, measure::AbstractVector, c Matrix{Float64}(function_basis), Vector{Float64}(measure), cdc, regularise, n, n0) end +Base.show(io::IO, t::MarkovTriple) = + print(io, "MarkovTriple(n=", t.n, ", n_function_basis=", t.n_function_basis, ")") + +""" + cdc(triple, f, h) -> Array + +Apply the triple's carré du champ: `Γ(f, h)`, the diffusion-geometric stand-in for +`∇f · ∇h`. `f` and `h` carry a leading point axis; the trailing axes of each appear +in order in the result. + +# Examples +```jldoctest +julia> Γ = cdc(dg.triple, dg.triple.immersion_coords, dg.triple.immersion_coords); + +julia> size(Γ) # the (n, d, d) first fundamental form +(60, 2, 2) +``` +""" cdc(t::MarkovTriple, f, h) = t.cdc(f, h) + +""" + regularise(triple, x) -> Array + +Apply the triple's regularisation map to `x` (leading axis = points). Depending on +how the triple was built this is a diffusion step ([`regularise_diffusion`](@ref)), a +spectral projection ([`regularise_bandlimit`](@ref)), or the identity. + +# Examples +```jldoctest +julia> noisy = gamma_coords(dg.cache); + +julia> size(regularise(dg.triple, noisy)) # smoothed, same shape +(60, 2, 2) +``` +""" regularise(t::MarkovTriple, x) = t.regularise(x) """ @@ -35,6 +82,19 @@ regularise(t::MarkovTriple, x) = t.regularise(x) `MarkovTriple` plus an immersion `x : M → ℝ^d` (`immersion_coords` `(n, d)`), required to generate the tensor algebra as an A-module. + +The immersion is what generates the tensor algebra: the coordinate functions `xⱼ` +give the frame `∇xⱼ` that every vector field and form is expanded in. It is the +extra ingredient a [`DiffusionGeometry`](@ref) needs beyond the diffusion itself. + +# Examples +```jldoctest +julia> dg.triple +ImmersedMarkovTriple(n=60, dim=2, n_function_basis=8) + +julia> size(dg.triple.immersion_coords), size(dg.triple.function_basis) +((60, 2), (60, 8)) +``` """ struct ImmersedMarkovTriple{C,R} base::MarkovTriple{C,R} @@ -59,16 +119,35 @@ Base.getproperty(t::ImmersedMarkovTriple, s::Symbol) = s in (:function_basis, :measure, :n, :n_function_basis) ? getfield(getfield(t, :base), s) : getfield(t, s) +Base.show(io::IO, t::ImmersedMarkovTriple) = + print(io, "ImmersedMarkovTriple(n=", t.n, ", dim=", t.dim, + ", n_function_basis=", t.n_function_basis, ")") + """ immersed_triple_from_knn_kernel(nbr_indices, kernel, immersion_coords; kwargs...) Resolve measure and function basis from the symmetric kernel, wire the cdc and -regularisation closures, and return an [`ImmersedMarkovTriple`]. Mirrors +regularisation closures, and return an [`ImmersedMarkovTriple`](@ref). Mirrors `DiffusionGeometry.from_knn_kernel`. Keyword args: `bandwidths`, `n_function_basis=50`, `regularisation_method="diffusion"` (`"diffusion"`/`"bandlimit"`/`"none"`), `measure`, `function_basis`, `use_mean_centres=true`, `data_matrix`. + +Use this when you already have a kernel — from another diffusion-maps +implementation, say. [`from_knn_kernel`](@ref) wraps it in a +[`DiffusionGeometry`](@ref). + +# Examples +```jldoctest +julia> kernel, bandwidths = markov_chain(knn_graph(circle, 16)...); + +julia> nbrs = knn_graph(circle, 16)[2]; + +julia> immersed_triple_from_knn_kernel(nbrs, kernel, circle; + bandwidths=bandwidths, n_function_basis=8) +ImmersedMarkovTriple(n=60, dim=2, n_function_basis=8) +``` """ function immersed_triple_from_knn_kernel(nbr_indices::AbstractMatrix{<:Integer}, kernel::AbstractMatrix, @@ -121,7 +200,16 @@ end Full pipeline: kNN graph → Markov chain → symmetric kernel → basis → triple. Mirrors `DiffusionGeometry.from_point_cloud`. Keyword args add `immersion_coords` (defaults to `data_matrix`), `knn_kernel=32`, `c=0`, `bandwidth_variability=-0.5`, -`knn_bandwidth=8` on top of [`immersed_triple_from_knn_kernel`]. +`knn_bandwidth=8` on top of [`immersed_triple_from_knn_kernel`](@ref). + +The triple underneath [`from_point_cloud`](@ref) — reach for it only if you want the +triple without the geometry wrapped around it. + +# Examples +```jldoctest +julia> immersed_triple_from_point_cloud(circle; knn_kernel=16, n_function_basis=8) +ImmersedMarkovTriple(n=60, dim=2, n_function_basis=8) +``` """ function immersed_triple_from_point_cloud(data_matrix::AbstractMatrix; immersion_coords=nothing, @@ -150,7 +238,7 @@ end """ immersed_triple_from_graph_kernel(edge_index, kernel, immersion_coords; kwargs...) -Build an [`ImmersedMarkovTriple`] from a precomputed kernel on an arbitrary +Build an [`ImmersedMarkovTriple`](@ref) from a precomputed kernel on an arbitrary directed graph. `edge_index` is `(2, num_edges)` with **1-based** node indices (row 1 source, row 2 target); `kernel` is the length-`num_edges` weight vector. Mirrors `DiffusionGeometry.from_graph_kernel`. @@ -158,6 +246,22 @@ Mirrors `DiffusionGeometry.from_graph_kernel`. Keyword args: `bandwidths=nothing`, `measure=nothing` (defaults to the source-weighted degree, normalised), `function_basis=nothing` (defaults to the `n×n` identity), `use_mean_centres=false`. + +Note the defaults differ from the kNN path: the function basis is the *pointwise* +identity (no spectral truncation) and the covariance is centred on the node's own +value, not the neighbourhood mean. + +# Examples +A directed 3-cycle, with weights on the edges: + +```jldoctest +julia> edges = [1 2 3; 2 3 1]; # 1 → 2 → 3 → 1 + +julia> coords = Matrix{Float64}(I, 3, 3); + +julia> immersed_triple_from_graph_kernel(edges, [1.0, 1.0, 1.0], coords) +ImmersedMarkovTriple(n=3, dim=3, n_function_basis=3) +``` """ function immersed_triple_from_graph_kernel(edge_index::AbstractMatrix{<:Integer}, kernel::AbstractVector, immersion_coords; @@ -181,11 +285,25 @@ end """ immersed_triple_from_edges(edge_index; immersion_coords=nothing, kwargs...) -Build an [`ImmersedMarkovTriple`] from a graph given only its edges, with the +Build an [`ImmersedMarkovTriple`](@ref) from a graph given only its edges, with the row-stochastic kernel `w_{ji} = 1/d(i)` (`d(i)` the in-degree of target `i`) and measure `μ(i) = d(i)`. `edge_index` is `(2, num_edges)`, **1-based** (row 1 source, row 2 target). Without `immersion_coords`, `n = maximum(edge_index)` and the immersion defaults to the `n×n` identity. Mirrors `DiffusionGeometry.from_edges`. + +The lightest way in: hand it a graph and it invents the rest. Combinatorial rather +than metric — this is the path to use for a graph that is not sampled from anything. + +# Examples +A 4-cycle, undirected (each edge listed both ways): + +```jldoctest +julia> edges = [1 2 3 4 2 3 4 1; + 2 3 4 1 1 2 3 4]; + +julia> immersed_triple_from_edges(edges) +ImmersedMarkovTriple(n=4, dim=4, n_function_basis=4) +``` """ function immersed_triple_from_edges(edge_index::AbstractMatrix{<:Integer}; immersion_coords=nothing, kwargs...) diff --git a/src/core/diffusion/regularise.jl b/src/core/diffusion/regularise.jl index 460cff9..a98c198 100644 --- a/src/core/diffusion/regularise.jl +++ b/src/core/diffusion/regularise.jl @@ -12,6 +12,24 @@ using LinearAlgebra: adjoint Regularise `x` by one diffusion step: `x_reg[p] = Σ_j kernel[p,j] * x[nbr_indices[p,j]]`. `kernel` and `nbr_indices` are `(n, k)`; `nbr_indices` is 1-based. `x` has shape `(n, tail...)` and the result has the same shape. + +The default regularisation of a geometry built from a point cloud: it is the local +average that damps the pointwise noise in the carré du champ. + +# Examples +Three points in a ring, each averaging itself with its successor: + +```jldoctest +julia> kernel = fill(0.5, 3, 2); + +julia> nbrs = [1 2; 2 3; 3 1]; # point p, then its successor + +julia> regularise_diffusion([1.0, 2.0, 3.0], kernel, nbrs) +3-element Vector{Float64}: + 1.5 + 2.5 + 2.0 +``` """ function regularise_diffusion(x::AbstractArray, kernel::AbstractMatrix, nbr_indices::AbstractMatrix{<:Integer}) @@ -41,6 +59,23 @@ end Regularise `x` by projecting onto the span of the coefficient functions `u` (shape `(n, n0)`) under the measure `measure` (shape `(n,)`): `x_reg = u * (u' * (measure .* x))`. `x` has shape `(n, tail...)`. + +The alternative to [`regularise_diffusion`](@ref), selected with +`regularisation_method="bandlimit"`: instead of averaging locally, it discards +everything outside the span of the eigenfunction basis. + +# Examples +Band-limiting to the constant function alone leaves the μ-weighted mean: + +```jldoctest +julia> u = ones(3, 1); # φ₀ ≡ 1, the only basis function + +julia> regularise_bandlimit([1.0, 2.0, 3.0], u, fill(1/3, 3)) +3-element Vector{Float64}: + 2.0 + 2.0 + 2.0 +``` """ function regularise_bandlimit(x::AbstractArray, u::AbstractMatrix, measure::AbstractVector) n = size(u, 1) diff --git a/src/core/geometry/cache.jl b/src/core/geometry/cache.jl index 3d43ecc..9607e60 100644 --- a/src/core/geometry/cache.jl +++ b/src/core/geometry/cache.jl @@ -9,9 +9,25 @@ """ GammaCache(triple) -Lazily memoised γ-tensors for an [`ImmersedMarkovTriple`]: `gamma_coords` -`Γ(x_i, x_j)`, `gamma_functions` `Γ(φ_i, φ_j)`, `gamma_mixed` `Γ(x_i, φ_j)`, their -regularised variants, and the per-degree compound matrices. +Lazily memoised γ-tensors for an [`ImmersedMarkovTriple`](@ref): [`gamma_coords`](@ref) +`Γ(x_i, x_j)`, [`gamma_functions`](@ref) `Γ(φ_i, φ_j)`, [`gamma_mixed`](@ref) +`Γ(x_i, φ_j)`, their regularised variants, and the per-degree compound matrices. + +These contractions are the expensive part of the package and every operator wants +them, so they are computed once, on first use, and kept. Each `DiffusionGeometry` +owns one, as `dg.cache`. + +# Examples +```jldoctest +julia> cache = dg.cache +GammaCache(n=60, dim=2, n_function_basis=8) + +julia> size(gamma_coords(cache)) # Γ(xᵢ, xⱼ): the metric, one d×d per point +(60, 2, 2) + +julia> gamma_coords(cache) === gamma_coords(cache) # memoised, not recomputed +true +``` """ mutable struct GammaCache triple::ImmersedMarkovTriple @@ -29,6 +45,37 @@ GammaCache(triple::ImmersedMarkovTriple) = GammaCache(triple, nothing, nothing, nothing, nothing, nothing, nothing, nothing, Dict{Int,Any}()) +Base.show(io::IO, c::GammaCache) = + print(io, "GammaCache(n=", c.triple.n, ", dim=", c.triple.dim, + ", n_function_basis=", c.triple.n_function_basis, ")") + +""" + gamma_coords(cache) -> Array + +`Γ(x_i, x_j)`, shape `(n, dim, dim)` — the carré du champ of the immersion +coordinates against each other. This is the *induced metric* (first fundamental +form) of the data: at each point, the `dim × dim` matrix of the ambient frame's +inner products. + +Its rank is the intrinsic dimension: directions tangent to the manifold have +eigenvalue ≈ 1, normal directions ≈ 0. + +# Examples +The circle is a curve in ℝ², so its metric has one tangent direction and one normal +one — an eigenvalue near 1, and one near 0: + +```jldoctest +julia> Γ = gamma_coords(dg.cache); + +julia> size(Γ) +(60, 2, 2) + +julia> round.(eigvals(Γ[1, :, :]); digits=2) +2-element Vector{Float64}: + 0.02 + 0.94 +``` +""" function gamma_coords(c::GammaCache) if c._gamma_coords === nothing t = c.triple @@ -37,6 +84,24 @@ function gamma_coords(c::GammaCache) return c._gamma_coords end +""" + gamma_functions(cache) -> Array + +`Γ(φ_i, φ_j)`, shape `(n, n0, n0)` — the carré du champ of the eigenfunction basis +against itself, i.e. `∇φ_i · ∇φ_j` at each point. The Dirichlet energies that the +Laplacian's weak form is assembled from. + +# Examples +```jldoctest +julia> G = gamma_functions(dg.cache); + +julia> size(G) # (n, n_function_basis, n_function_basis) +(60, 8, 8) + +julia> round(maximum(abs, G[:, 1, :]); digits=12) # φ₀ ≡ 1, so ∇φ₀ = 0 +0.0 +``` +""" function gamma_functions(c::GammaCache) if c._gamma_functions === nothing t = c.triple @@ -45,6 +110,20 @@ function gamma_functions(c::GammaCache) return c._gamma_functions end +""" + gamma_mixed(cache) -> Array + +`Γ(x_i, φ_j)`, shape `(n, dim, n0)` — the carré du champ of the immersion coordinates +against the eigenfunction basis. This is the bridge between the two: it expresses +`∇φ_j` in the coordinate frame `∇x_i`, which is exactly what [`grad`](@ref) and the +exterior derivative are built from. + +# Examples +```jldoctest +julia> size(gamma_mixed(dg.cache)) # (n, dim, n_function_basis) +(60, 2, 8) +``` +""" function gamma_mixed(c::GammaCache) if c._gamma_mixed === nothing t = c.triple @@ -64,7 +143,18 @@ end gamma_ambient(cache) -> Array `Γ(x^ambient_i, x_j)` `(n, ambient_dim, dim)` when a `data_matrix` is present, -otherwise falls back to `gamma_coords`. +otherwise falls back to [`gamma_coords`](@ref). + +The map that pushes an intrinsic tensor back out into the coordinates the data was +given in — [`to_ambient`](@ref) and the quiver plots are built on it. It differs from +`gamma_coords` only when the immersion is not the raw data (e.g. a regularised or +reduced immersion). + +# Examples +```jldoctest +julia> size(gamma_ambient(dg.cache)) # (n, ambient D, intrinsic frame dim) +(60, 2, 2) +``` """ function gamma_ambient(c::GammaCache) if c._gamma_ambient === nothing @@ -79,6 +169,14 @@ end hessian_functions(cache) -> Array Regularised Hessian tensor `H(φ_k)(∇x_i, ∇x_j)`, shape `(n, dim, dim, n0)`. +The second-order counterpart of [`gamma_mixed`](@ref), and the raw material of the +[`hessian`](@ref) operator. + +# Examples +```jldoctest +julia> size(hessian_functions(dg.cache)) # (n, dim, dim, n_function_basis) +(60, 2, 2, 8) +``` """ function hessian_functions(c::GammaCache) if c._hessian_functions === nothing @@ -97,6 +195,14 @@ end hessian_coords(cache) -> Array Regularised coordinate Hessian `H(x_k)(∇x_i, ∇x_j)`, shape `(n, dim, dim, dim)`. +Geometrically the second fundamental form of the immersion; it supplies the +Christoffel symbols that make [`levi_civita`](@ref) a connection. + +# Examples +```jldoctest +julia> size(hessian_coords(dg.cache)) # (n, dim, dim, dim) +(60, 2, 2, 2) +``` """ function hessian_coords(c::GammaCache) if c._hessian_coords === nothing @@ -114,6 +220,19 @@ end Compound γ submatrices and their determinants for form degree `k` (`1 ≤ k ≤ d`), mirroring the Python `lru_cache`d accessor. `dets` (shape `(n, C, C)`, `C = binomial(d, k)`) are what the k-form space uses as its `cdc_components`. + +Memoised per degree; see [`gamma_compound`](@ref) for the determinants themselves. + +# Examples +```jldoctest +julia> subs, dets = gamma_coords_compound(dg.cache, 2); + +julia> size(dets) # d = 2, so there is one 2-form component per point +(60, 1, 1) + +julia> size(gamma_coords_compound(dg.cache, 1)[2]) # k = 1 is Γ itself +(60, 2, 2) +``` """ function gamma_coords_compound(c::GammaCache, k::Integer) @assert 1 <= k <= c.triple.dim "Form degree k=$k must be between 1 and $(c.triple.dim)" diff --git a/src/core/geometry/diffusion_geometry.jl b/src/core/geometry/diffusion_geometry.jl index 04eed9b..e91ccca 100644 --- a/src/core/geometry/diffusion_geometry.jl +++ b/src/core/geometry/diffusion_geometry.jl @@ -15,9 +15,32 @@ import SparseArrays """ DiffusionGeometry(triple; rcond=1e-5, n_coefficients=nothing) -Global coordinate-free representation built on an [`ImmersedMarkovTriple`]. +Global coordinate-free representation built on an [`ImmersedMarkovTriple`](@ref). `n_coefficients` (defaulting to, and capped at, `n_function_basis`) sets the low-rank truncation used by tensor expansions; `rcond` is the spectral cutoff. + +The central object: it owns the basis, the measure, the γ cache, the tensor spaces, +and (memoised) every differential operator. Build one with [`from_point_cloud`](@ref) +or one of its siblings rather than calling this constructor directly. + +There are two truncation levels, and they are not the same thing. `n_function_basis` +is how many eigenfunctions the *diffusion* resolves; `n_coefficients ≤ +n_function_basis` is how many of them a *tensor* is expanded in. Functions always use +the full basis; vector fields, forms and (0,2)-tensors use `n_coefficients`. Lowering +it is the main performance dial — but the topology functions need the full basis, or +the harmonic forms are truncated away. + +# Examples +```jldoctest +julia> dg +DiffusionGeometry(n=60, ambient_dim=2, n_function_basis=8, n_coefficients=8) + +julia> from_point_cloud(circle; knn_kernel=16, n_function_basis=8, n_coefficients=4) +DiffusionGeometry(n=60, ambient_dim=2, n_function_basis=8, n_coefficients=4) + +julia> space_dim(function_space(dg)), space_dim(vector_field_space(dg)) +(8, 16) +``` """ mutable struct DiffusionGeometry triple::ImmersedMarkovTriple @@ -41,6 +64,11 @@ function DiffusionGeometry(triple::ImmersedMarkovTriple; rcond=1e-5, n_coefficie Dict{Int,AbstractTensorSpace}(), Dict{Any,Any}()) end +Base.show(io::IO, dg::DiffusionGeometry) = + print(io, "DiffusionGeometry(n=", npoints(dg), ", ambient_dim=", ambient_dim(dg), + ", n_function_basis=", n_function_basis(dg), + ", n_coefficients=", n_coefficients(dg), ")") + # ── Constructors from different data sources (mirror the Python classmethods) ── """ from_knn_kernel(nbr_indices, kernel, immersion_coords; rcond=1e-5, @@ -49,7 +77,23 @@ end Build from a precomputed kernel on a neighbour graph. Extra keywords (`n_function_basis`, `regularisation_method`, `bandwidths`, `measure`, `function_basis`, `use_mean_centres`, `data_matrix`) pass through to -[`immersed_triple_from_knn_kernel`]. +[`immersed_triple_from_knn_kernel`](@ref). + +Use this to bring your own diffusion kernel; [`from_point_cloud`](@ref) is the same +thing with the kernel computed for you. + +# Examples +Rebuilding the running example the long way round gives back the same geometry: + +```jldoctest +julia> nbr_distances, nbr_indices = knn_graph(circle, 16); + +julia> kernel, bandwidths = markov_chain(nbr_distances, nbr_indices); + +julia> from_knn_kernel(nbr_indices, kernel, circle; bandwidths=bandwidths, + n_function_basis=8) +DiffusionGeometry(n=60, ambient_dim=2, n_function_basis=8, n_coefficients=8) +``` """ function from_knn_kernel(nbr_indices::AbstractMatrix{<:Integer}, kernel::AbstractMatrix, immersion_coords=nothing; rcond=1e-5, n_coefficients=nothing, kwargs...) @@ -62,7 +106,19 @@ end bandwidth_variability=-0.5, knn_bandwidth=8, rcond=1e-5, n_coefficients=nothing, kwargs...) -> DiffusionGeometry -Build a Markov chain from a neighbour graph, then defer to [`from_knn_kernel`]. +Build a Markov chain from a neighbour graph, then defer to [`from_knn_kernel`](@ref). + +The entry point when you have neighbours and distances but not the point cloud they +came from — an approximate kNN index, or a graph with edge lengths. + +# Examples +```jldoctest +julia> nbr_distances, nbr_indices = knn_graph(circle, 16); + +julia> from_knn_graph(nbr_indices, nbr_distances; immersion_coords=circle, + n_function_basis=8) +DiffusionGeometry(n=60, ambient_dim=2, n_function_basis=8, n_coefficients=8) +``` """ function from_knn_graph(nbr_indices::AbstractMatrix{<:Integer}, nbr_distances::AbstractMatrix; immersion_coords=nothing, c::Real=0, bandwidth_variability::Real=-0.5, @@ -80,6 +136,29 @@ end Full pipeline from raw data: kNN graph → Markov chain → symmetric kernel → eigenbasis → triple → geometry. Mirrors `DiffusionGeometry.from_point_cloud`. + +The usual way in. `data_matrix` is `(n_points, ambient_dim)`. Keywords worth knowing: +`knn_kernel=32` (neighbours per point — must exceed `knn_bandwidth=8`), +`n_function_basis=50` (how many eigenfunctions to resolve), `n_coefficients` +(the tensor truncation, defaults to all of them), and `regularisation_method`. + +# Examples +Sixty points on the unit circle. The geometry knows nothing about circles — it reads +the intrinsic dimension, the metric, and the Laplacian spectrum (`0, 1, 1, 4, 4, …`) +off the samples: + +```jldoctest +julia> θ = range(0, 2π; length=61)[1:60]; + +julia> dg = from_point_cloud([cos.(θ) sin.(θ)]; knn_kernel=16, n_function_basis=8) +DiffusionGeometry(n=60, ambient_dim=2, n_function_basis=8, n_coefficients=8) + +julia> round.(abs.(spectrum(laplacian(dg, 0); eigvals_only=true))[1:3]; digits=3) +3-element Vector{Float64}: + 0.0 + 0.999 + 1.0 +``` """ function from_point_cloud(data_matrix::AbstractMatrix; rcond=1e-5, n_coefficients=nothing, kwargs...) triple = immersed_triple_from_point_cloud(data_matrix; kwargs...) @@ -93,7 +172,16 @@ end Build from a precomputed kernel on an arbitrary directed graph. `edge_index` is `(2, num_edges)`, **1-based** (row 1 source, row 2 target); `kernel` the edge-weight vector. Extra keywords (`bandwidths`, `measure`, `function_basis`, -`use_mean_centres`) pass through to [`immersed_triple_from_graph_kernel`]. +`use_mean_centres`) pass through to [`immersed_triple_from_graph_kernel`](@ref). + +# Examples +```jldoctest +julia> edges = [1 2 3 4 2 3 4 1; # an undirected 4-cycle + 2 3 4 1 1 2 3 4]; + +julia> from_graph_kernel(edges, fill(0.5, 8), Matrix{Float64}(I, 4, 4)) +DiffusionGeometry(n=4, ambient_dim=4, n_function_basis=4, n_coefficients=4) +``` """ function from_graph_kernel(edge_index::AbstractMatrix{<:Integer}, kernel::AbstractVector, immersion_coords; rcond=1e-5, n_coefficients=nothing, kwargs...) @@ -107,7 +195,24 @@ end Build from a graph given only its (1-based) `(2, num_edges)` edges, using the row-stochastic kernel `w_{ji} = 1/d(i)` and measure `μ(i) = d(i)`. Defers to -[`immersed_triple_from_edges`]. +[`immersed_triple_from_edges`](@ref). + +The purely combinatorial entry point: no coordinates, no kernel, just a graph. + +# Examples +A 4-cycle. Its graph Laplacian has the eigenvalues of the discrete circle, and its +Betti numbers are those of a loop: + +```jldoctest +julia> edges = [1 2 3 4 2 3 4 1; + 2 3 4 1 1 2 3 4]; + +julia> dgg = from_edges(edges) +DiffusionGeometry(n=4, ambient_dim=4, n_function_basis=4, n_coefficients=4) + +julia> betti_number(dgg, 0) +1 +``` """ function from_edges(edge_index::AbstractMatrix{<:Integer}; rcond=1e-5, n_coefficients=nothing, kwargs...) @@ -120,7 +225,17 @@ end Bridge for graph libraries emitting sparse kernels `K[i, j]` (point `i`, neighbour `j`). The `(row, col) = (i, j)` entries become edges `(source=j, target=i)`, then -[`from_graph_kernel`] takes over. +[`from_graph_kernel`](@ref) takes over. + +# Examples +```jldoctest +julia> using SparseArrays + +julia> K = sparse([1, 2, 2, 3, 3, 1], [2, 1, 3, 2, 1, 3], fill(0.5, 6), 3, 3); + +julia> from_sparse_matrix(K, Matrix{Float64}(I, 3, 3)) +DiffusionGeometry(n=3, ambient_dim=3, n_function_basis=3, n_coefficients=3) +``` """ function from_sparse_matrix(sparse_matrix::SparseArrays.AbstractSparseMatrix, immersion_coords; kwargs...) @@ -130,38 +245,218 @@ function from_sparse_matrix(sparse_matrix::SparseArrays.AbstractSparseMatrix, end # ── Data accessors (mirror the Python @property forwarders) ──────────────────── +""" + npoints(dg) -> Int + +Number of sample points `n`. + +# Examples +```jldoctest +julia> npoints(dg) +60 +``` +""" npoints(dg::DiffusionGeometry) = dg.triple.n + +""" + ambient_dim(dg) -> Int + +Dimension `d` of the immersion — the number of coordinate functions, hence the number +of components a vector field has per point. This is the dimension of the *ambient* +frame, not the intrinsic dimension of the manifold (the circle below is a 1-manifold +immersed in ℝ², so `ambient_dim` is 2). + +# Examples +```jldoctest +julia> ambient_dim(dg), ambient_dim(dg3) +(2, 3) +``` +""" ambient_dim(dg::DiffusionGeometry) = dg.triple.dim + +""" + n_coefficients(dg) -> Int + +How many eigenfunctions a *tensor* is expanded in — the truncation rank of vector +fields, forms and (0,2)-tensors. At most [`n_function_basis`](@ref); lowering it is +the main way to make the operators cheaper. + +# Examples +```jldoctest +julia> n_coefficients(dg) +8 + +julia> n_coefficients(from_point_cloud(circle; knn_kernel=16, n_function_basis=8, + n_coefficients=4)) +4 +``` +""" n_coefficients(dg::DiffusionGeometry) = dg.n_coefficients + +""" + n_function_basis(dg) -> Int + +How many eigenfunctions the diffusion resolves — the dimension of the function space, +and the ceiling on [`n_coefficients`](@ref). + +# Examples +```jldoctest +julia> n_function_basis(dg) +8 + +julia> space_dim(function_space(dg)) # functions always use the full basis +8 +``` +""" n_function_basis(dg::DiffusionGeometry) = dg.n_function_basis + +""" + function_basis(dg) -> Matrix + +The eigenfunction basis `{φ_i}` evaluated at the sample points, shape `(n, n0)`. +Column 1 is the constant function; later columns are the increasingly oscillatory +modes of the diffusion. Their *signs* come from the eigensolver — don't depend on them. + +# Examples +```jldoctest +julia> size(function_basis(dg)) +(60, 8) + +julia> all(≈(1.0), function_basis(dg)[:, 1]) # φ₀ ≡ 1 +true +``` +""" function_basis(dg::DiffusionGeometry) = dg.triple.function_basis + +""" + measure(dg) -> Vector + +The measure `μ` on the sample points (shape `(n,)`, summing to 1). Every integral in +the package — [`inner`](@ref), [`gram`](@ref), the weak forms — is taken against it. + +# Examples +The circle is sampled uniformly, so the measure is close to flat: + +```jldoctest +julia> sum(measure(dg)) ≈ 1.0 +true + +julia> round(maximum(measure(dg)) - minimum(measure(dg)); digits=6) +0.0 +``` +""" measure(dg::DiffusionGeometry) = dg.triple.measure + +""" + immersion_coords(dg) -> Matrix + +The immersion `x : M → ℝ^d` at the sample points, shape `(n, d)`. Its coordinate +functions generate the tensor algebra (see [`ImmersedMarkovTriple`](@ref)). Usually +the point cloud itself, possibly regularised. + +# Examples +```jldoctest +julia> size(immersion_coords(dg)) +(60, 2) +``` +""" immersion_coords(dg::DiffusionGeometry) = dg.triple.immersion_coords + regularise_fn(dg::DiffusionGeometry) = x -> regularise(dg.triple, x) Base.:(==)(a::DiffusionGeometry, b::DiffusionGeometry) = a === b # ── Tensor spaces (lazily constructed & cached) ──────────────────────────────── +""" + function_space(dg) -> FunctionSpace + +The space of scalar functions on `dg`. Cached, so it is the *same* object every time +— which is what lets operators compare domains and codomains by identity. + +# Examples +```jldoctest +julia> function_space(dg) +FunctionSpace(dim=8) + +julia> function_space(dg) === function_space(dg) +true +``` +""" function function_space(dg::DiffusionGeometry) dg._function_space === nothing && (dg._function_space = FunctionSpace(dg)) return dg._function_space end +""" + vector_field_space(dg) -> VectorFieldSpace + +The space of vector fields on `dg`: `n_coefficients × ambient_dim` coefficients. +Cached like [`function_space`](@ref). + +# Examples +```jldoctest +julia> vector_field_space(dg) # 8 basis functions × 2 components +VectorFieldSpace(dim=16) +``` +""" function vector_field_space(dg::DiffusionGeometry) dg._vector_field_space === nothing && (dg._vector_field_space = VectorFieldSpace(dg)) return dg._vector_field_space end +""" + tensor02_space(dg) -> Tensor02Space + +The space of general (0,2)-tensors on `dg`: `n_coefficients × ambient_dim²` +coefficients. + +# Examples +```jldoctest +julia> tensor02_space(dg) # 8 basis functions × 2² components +Tensor02Space(dim=32) +``` +""" function tensor02_space(dg::DiffusionGeometry) dg._tensor02_space === nothing && (dg._tensor02_space = Tensor02Space(dg)) return dg._tensor02_space end +""" + tensor02sym_space(dg) -> Tensor02SymSpace + +The space of symmetric (0,2)-tensors on `dg`: `n_coefficients × d(d+1)/2` +coefficients. The codomain of [`hessian`](@ref). + +# Examples +```jldoctest +julia> tensor02sym_space(dg) # 8 basis functions × 3 symmetric components +Tensor02SymSpace(dim=24) +``` +""" function tensor02sym_space(dg::DiffusionGeometry) dg._tensor02sym_space === nothing && (dg._tensor02sym_space = Tensor02SymSpace(dg)) return dg._tensor02sym_space end +""" + form_space(dg, degree) -> AbstractTensorSpace + +The space of differential `degree`-forms on `dg`, with `binomial(ambient_dim, degree)` +components per point. Degree 0 is [`function_space`](@ref); degrees run up to +`ambient_dim`. Cached per degree. + +# Examples +```jldoctest +julia> form_space(dg, 1) # 8 × binomial(2, 1) = 16 +FormSpace(degree=1, dim=16) + +julia> form_space(dg, 2) # 8 × binomial(2, 2) = 8 +FormSpace(degree=2, dim=8) + +julia> form_space(dg, 0) === function_space(dg) +true +``` +""" function form_space(dg::DiffusionGeometry, degree::Integer) degree == 0 && return function_space(dg) return get!(dg._form_spaces, Int(degree)) do @@ -170,19 +465,105 @@ function form_space(dg::DiffusionGeometry, degree::Integer) end # ── Factories from pointwise data (mirror dg.function / vector_field / …) ─────── +""" + dg_function(dg, data) -> ScalarFunction + +The function whose values at the sample points are `data` (trailing axis of length +`n`) — or rather, its L²-projection onto the eigenfunction basis, which is the closest +thing the geometry can represent. Leading axes are batch dims. + +# Examples +```jldoctest +julia> f = dg_function(dg, cos.(θ)) +ScalarFunction(space=FunctionSpace(dim=8), shape=(8,), batch_shape=()) + +julia> maximum(abs.(to_pointwise_basis(f) .- cos.(θ))) < 1e-4 # 8 modes suffice +true + +julia> dg_function(dg, ones(3, 60)) # a batch of three functions +ScalarFunction(space=FunctionSpace(dim=8), shape=(3, 8), batch_shape=(3,)) +``` +""" dg_function(dg::DiffusionGeometry, data) = from_pointwise(function_space(dg), data) + """ -`mode=:pullback` projects ambient components onto the diffusion basis; -`mode=:reconstruct` least-squares fits them (the inverse of `to_ambient`). + dg_vector_field(dg, data; mode=:pullback) -> VectorField + +The vector field whose components at the sample points are `data` (trailing shape +`(n, d)`). + +`mode=:pullback` projects the ambient components onto the diffusion basis; +`mode=:reconstruct` least-squares fits them, which is the left inverse of +[`to_ambient`](@ref) (see [`vector_field_from_reconstruction`](@ref)). Use +`:reconstruct` when `data` is a quiver you want to recover exactly; `:pullback` +otherwise. + +# Examples +The tangent field of the circle, `(-sin θ, cos θ)`: + +```jldoctest +julia> T = dg_vector_field(dg, [-sin.(θ) cos.(θ)]) +VectorField(space=VectorFieldSpace(dim=16), shape=(16,), batch_shape=()) + +julia> round(l2_norm(T); digits=3) # ‖T‖_L² = 1, up to the metric's error +0.972 +``` """ function dg_vector_field(dg::DiffusionGeometry, data; mode::Symbol=:pullback) mode === :pullback && return from_pointwise(vector_field_space(dg), data) mode === :reconstruct && return vector_field_from_reconstruction(dg, data) throw(ArgumentError("mode must be :pullback or :reconstruct, got :$mode")) end + +""" + dg_form(dg, data, degree) -> Form + +The `degree`-form whose components at the sample points are `data` (trailing shape +`(n, binomial(d, degree))`, in the lexicographic wedge-basis order of +[`get_wedge_basis_indices`](@ref)). Degree 0 gives a [`ScalarFunction`](@ref). + +# Examples +```jldoctest +julia> dg_form(dg, [-sin.(θ) cos.(θ)], 1) +Form(space=FormSpace(degree=1, dim=16), shape=(16,), batch_shape=()) + +julia> dg_form(dg, ones(60, 1), 2) # one component: dx₁∧dx₂ +Form(space=FormSpace(degree=2, dim=8), shape=(8,), batch_shape=()) +``` +""" dg_form(dg::DiffusionGeometry, data, degree::Integer) = degree == 0 ? dg_function(dg, data) : from_pointwise(form_space(dg, degree), data) + +""" + dg_tensor02(dg, data) -> Tensor02 + +The (0,2)-tensor whose components at the sample points are `data` (trailing shape +`(n, d²)`, row-major in `(j, k)`). + +# Examples +The ambient identity tensor `δ_{jk}`: + +```jldoctest +julia> data = repeat([1.0 0.0 0.0 1.0], 60); + +julia> dg_tensor02(dg, data) +Tensor02(space=Tensor02Space(dim=32), shape=(32,), batch_shape=()) +``` +""" dg_tensor02(dg::DiffusionGeometry, data) = from_pointwise(tensor02_space(dg), data) + +""" + dg_tensor02sym(dg, data) -> Tensor02Sym + +The symmetric (0,2)-tensor whose components at the sample points are `data` (trailing +shape `(n, d(d+1)/2)`, in the order of [`get_symmetric_basis_indices`](@ref)). + +# Examples +```jldoctest +julia> dg_tensor02sym(dg, repeat([1.0 0.0 1.0], 60)) # (T₁₁, T₁₂, T₂₂) +Tensor02Sym(space=Tensor02SymSpace(dim=24), shape=(24,), batch_shape=()) +``` +""" dg_tensor02sym(dg::DiffusionGeometry, data) = from_pointwise(tensor02sym_space(dg), data) # ── Riemannian metric and global inner product ───────────────────────────────── @@ -190,7 +571,23 @@ dg_tensor02sym(dg::DiffusionGeometry, data) = from_pointwise(tensor02sym_space(d g(dg, a, b) -> Array Pointwise Riemannian metric `g_p(a, b)` of two tensors in the same space -(shape `(batch..., n)`). +(shape `(batch..., n)`). The learned inner product — this is where the geometry of +the data enters, and it is a *function on the manifold*, not a number. + +# Examples +On the unit circle, `∇cos θ = -sin θ ∂_θ`, so `g(∇f, ∇f) = sin²θ` point by point: + +```jldoctest +julia> X = grad(f); + +julia> vals = g(dg, X, X); + +julia> size(vals) +(60,) + +julia> maximum(abs.(vals .- sin.(θ) .^ 2)) < 0.1 +true +``` """ function g(dg::DiffusionGeometry, a::AbstractTensor, b::AbstractTensor) @assert a.space == b.space "Spaces $(a.space) vs $(b.space) do not match" @@ -200,8 +597,25 @@ end """ pointwise_norm(dg, a) -> Array + pointwise_norm(a) -> Array + +Pointwise norm `‖a‖_p = √(g_p(a, a))`, shape `(batch..., n)`. The one-argument form +takes the geometry from the tensor. + +# Examples +`‖∇cos θ‖ = |sin θ|`, which peaks at 1 on the sides of the circle. It does not quite +reach 0 at the poles — an 8-mode band-limited gradient cannot, and 0.23 is the size of +that error: -Pointwise norm `‖a‖_p = √(g_p(a, a))`. +```jldoctest +julia> nrm = pointwise_norm(grad(f)); + +julia> round(maximum(nrm); digits=2) +0.97 + +julia> round(minimum(nrm); digits=2) +0.23 +``` """ pointwise_norm(dg::DiffusionGeometry, a::AbstractTensor) = sqrt.(max.(real.(g(dg, a, a)), 0.0)) @@ -210,6 +624,34 @@ pointwise_norm(dg::DiffusionGeometry, a::AbstractTensor) = sqrt.(max.(real.(g(dg Global L² inner product `⟨a, b⟩ = ∫ g(a, b) dμ = coeffsᵀ G coeffs` (`G` the Gram matrix). Returns a scalar for unbatched tensors, otherwise a batched array. + +Integrating [`g`](@ref) against the [`measure`](@ref) — the inner product that makes +each tensor space a Hilbert space, and the one every adjoint in the package is taken +with respect to. + +# Examples +`∫ cos²θ dθ/2π = 1/2`: + +```jldoctest +julia> round(inner(dg, f, f); digits=4) +0.5 + +julia> h = dg_function(dg, sin.(θ)); + +julia> abs(inner(dg, f, h)) < 1e-5 # sin ⟂ cos +true +``` + +Batched tensors give one number per batch element: + +```jldoctest +julia> fb = dg_function(dg, [cos.(θ) sin.(θ)]'); + +julia> round.(inner(dg, fb, fb); digits=4) +2-element Vector{Float64}: + 0.5 + 0.5 +``` """ function inner(dg::DiffusionGeometry, a::AbstractTensor, b::AbstractTensor) @assert a.space == b.space "Spaces $(a.space) vs $(b.space) do not match" @@ -225,8 +667,19 @@ end """ l2_norm(dg, a) -> scalar or Array + l2_norm(a) -> scalar or Array + +Global L² norm `‖a‖ = √⟨a, a⟩`. The one-argument form takes the geometry from the +tensor. + +# Examples +```jldoctest +julia> round(l2_norm(f); digits=4) # √(1/2) +0.7071 -Global L² norm `‖a‖ = √⟨a, a⟩`. +julia> round(l2_norm(grad(f)); digits=4) # ‖∇cos θ‖ = √(∫sin²θ) = √(1/2) +0.7069 +``` """ function l2_norm(dg::DiffusionGeometry, a::AbstractTensor) v = inner(dg, a, a) diff --git a/src/methods/geodesics.jl b/src/methods/geodesics.jl index 87171ac..5b9af4d 100644 --- a/src/methods/geodesics.jl +++ b/src/methods/geodesics.jl @@ -30,6 +30,32 @@ Laplacian. Returns the pointwise distances `dist` (length `n`) and the correctio `index` is 1-based. `num_subsample` (or `nothing` for all points) sets how many random points carry the pointwise Lipschitz constraint. + +Requires a `data_matrix` on the triple (so: a geometry built from a point cloud), and +solves a second-order cone program with SCS — this is much slower than the rest of the +package, and its accuracy is bounded by the resolution of the function basis. + +# Examples +Distance along the circle from point 1. The exact answer is the arc length, `min(θ, +2π - θ)`: zero at the source, and maximal at the antipode (point 31 of 60), where the +true value is `π ≈ 3.14`. The 8-mode basis under-resolves the cone at the source, so +the far side comes out short — the shape is right, the scale is approximate: + +```jldoctest +julia> dist, v = geodesic_distances_function(dg, 1); + +julia> size(dist), v +((60,), ScalarFunction(space=FunctionSpace(dim=8), shape=(8,), batch_shape=())) + +julia> round(dist[1]; digits=6) # d(x, x) = 0, imposed as a constraint +0.0 + +julia> argmax(dist) # farthest point: the antipode +31 + +julia> round(maximum(dist); digits=1) # π = 3.14, under-resolved by the basis +2.7 +``` """ function geodesic_distances_function(dg::DiffusionGeometry, index::Integer; reg::Bool=false, num_subsample=nothing, diff --git a/src/methods/pde.jl b/src/methods/pde.jl index 4e89f5a..6bb72ec 100644 --- a/src/methods/pde.jl +++ b/src/methods/pde.jl @@ -24,6 +24,30 @@ Python original inherits whatever sign LAPACK happens to return, which is not portable across BLAS implementations). We canonicalise each eigenvector's sign so its largest-magnitude entry is positive, making the result deterministic; the parity fixture applies the same canonicalisation. + +# Examples +The heat equation `∂ₜu = -Δu` on the circle, from `u₀ = cos θ`. That is (almost +exactly) the `λ = 1` eigenfunction, so it should decay by `e^{-t}` and keep its shape: + +```jldoctest +julia> u = solve_differential_operator(-laplacian(dg, 0), f, [0.0, 0.5, 1.0]) +ScalarFunction(space=FunctionSpace(dim=8), shape=(3, 8), batch_shape=(3,)) + +julia> u.coeffs[1, :] ≈ f.coeffs # t = 0 returns the initial condition +true + +julia> round.(l2_norm(u); digits=3) # ‖u(t)‖ = e^{-t}‖u₀‖ +3-element Vector{Float64}: + 0.707 + 0.429 + 0.26 + +julia> round(exp(-1) * l2_norm(f); digits=3) # the expected value at t = 1 +0.26 +``` + +Run it backwards (`+Δ`) and the same modes grow instead — which is exactly why the +backward heat equation is unstable, and the solver will happily show you that. """ function solve_differential_operator(operator::LinearOperator, initial_condition::AbstractTensor, diff --git a/src/methods/topology.jl b/src/methods/topology.jl index ae579f7..adbc5b7 100644 --- a/src/methods/topology.jl +++ b/src/methods/topology.jl @@ -31,6 +31,33 @@ weight), `values` (ascending eigenvalues). Get the auto-counted Betti number with [`betti_number`](@ref); print it to review the low cluster and spectral gap by eye. + +# Examples +The circle has one independent loop, and its degree-1 spectrum says so loudly: one +eigenvalue near zero, then a jump of several orders of magnitude. + +```julia +julia> betti_spectrum(dg, 1) +BettiSpectrum: degree k=1, penalty w=1.0e10 + suggested b₁ = 1 (spectral gap 1650.0× — clear) + smallest eigenvalues: + 1: 4.582e-5 ← harmonic + 2: 0.07536 ┃ gap ×1650.0 + 3: 0.6018 + 4: 5.411 + 5: 5.412 + 6: 8.133 + 7: 8.139 + 8: 27.57 + … (16 total) +``` + +!!! warning "The harmonic eigenvalue's *size* is not reproducible" + The penalty `w = 1e10` multiplies the down-Laplacian energy, so it also multiplies + the float noise in it: the near-zero eigenvalues land somewhere around `1e-5`–`1e-3` + and move from one BLAS/machine to the next. Only their *separation* from the rest of + the spectrum is meaningful, which is why [`betti_number`](@ref) counts a gap and + never a threshold. The example above is illustrative, not a fixed value. """ struct BettiSpectrum k::Int @@ -48,6 +75,29 @@ Primary API. Ascending eigenvalues of the penalised Hodge operator for degree `k Build `dg` with the full coefficient basis (`n_coefficients == n_function_basis`) or the harmonic forms are truncated away. Feed the result to [`betti_number`](@ref) to auto-count, or print it to read the gap by eye. + +# Examples +```jldoctest +julia> bs = betti_spectrum(dg, 1); + +julia> (bs.k, bs.w) +(1, 1.0e10) + +julia> length(bs.values) # one per coefficient of Ω¹ +16 + +julia> round(bs.values[2]; sigdigits=3) # the first non-harmonic eigenvalue +0.0754 + +julia> bs.values[1] < bs.values[2] / 100 # the harmonic one is orders of magnitude below +true + +julia> betti_number(bs) # so b₁ = 1: the circle has one loop +1 +``` + +The size of `values[1]` itself is float noise amplified by the penalty — see the warning +on [`BettiSpectrum`](@ref). Count the gap, never a threshold. """ function betti_spectrum(dg::DiffusionGeometry, k::Integer; w::Real=1e10) H = k == 0 ? laplacian(dg, 0) : up_laplacian(dg, k) + w * down_laplacian(dg, k) @@ -59,6 +109,20 @@ end `betti_spectrum` for every degree `0:kmax`, for reviewing the whole complex at once. `result[k+1]` is the degree-`k` spectrum. + +# Examples +```jldoctest +julia> spectra = betti_spectra(dg); + +julia> length(spectra) # degrees 0, 1, 2 +3 + +julia> [betti_number(bs) for bs in spectra] +3-element Vector{Int64}: + 1 + 1 + 0 +``` """ betti_spectra(dg::DiffusionGeometry; kmax::Integer=ambient_dim(dg), w::Real=1e10) = [betti_spectrum(dg, k; w) for k in 0:kmax] @@ -71,6 +135,17 @@ Auto-count `b_k`: the low-eigenvalue cluster cut off by the largest multiplicati gap of at least `gap_ratio`, among eigenvalues below the penalty scale. Use [`betti_gap`](@ref) for the confidence (the gap size) and print the [`BettiSpectrum`](@ref) if the count looks off. + +# Examples +A circle is connected (`b₀ = 1`) with one loop (`b₁ = 1`): + +```jldoctest +julia> betti_number(dg, 0), betti_number(dg, 1) +(1, 1) + +julia> betti_number(betti_spectrum(dg, 1)) # same, from a spectrum you already have +1 +``` """ betti_number(bs::BettiSpectrum; gap_ratio::Real=3.0) = _count_gap(bs.values; gap_ratio, w=bs.w)[1] @@ -83,6 +158,16 @@ betti_number(dg::DiffusionGeometry, k::Integer; w::Real=1e10, gap_ratio::Real=3. The multiplicative spectral gap at the auto-detected cut — a confidence signal for [`betti_number`](@ref). A large gap (≳10×) is a clear Betti number; a gap near `gap_ratio` is marginal and worth an eyeball check. + +# Examples +```jldoctest +julia> betti_gap(betti_spectrum(dg, 1)) > 100 # a gap of ~10³×: nothing marginal here +true +``` + +(The exact number moves between machines — the harmonic eigenvalue below the gap is +noise amplified by the penalty. Compare it against `gap_ratio`, not against a +remembered value.) """ betti_gap(bs::BettiSpectrum; gap_ratio::Real=3.0) = _count_gap(bs.values; gap_ratio, w=bs.w)[2] @@ -93,6 +178,17 @@ betti_gap(bs::BettiSpectrum; gap_ratio::Real=3.0) = Auto-count the Betti numbers `[b_0, …, b_kmax]` (so `result[k+1] == b_k`) over all degrees. Convenience over [`betti_number`](@ref); see it for the method and caveats. +# Examples +The circle: connected, one loop, and no 2-dimensional void. + +```jldoctest +julia> betti_numbers(dg) +3-element Vector{Int64}: + 1 + 1 + 0 +``` + !!! note "It's a spectral-gap heuristic" `b_0` is robust. `b_1` is reliable when the gap is clear (a marginal gap, e.g. a lone circle, can flip). Top-degree counts inflate to the ambient form-space diff --git a/src/operators/differential_operators/derivative.jl b/src/operators/differential_operators/derivative.jl index eb5bacd..00d6014 100644 --- a/src/operators/differential_operators/derivative.jl +++ b/src/operators/differential_operators/derivative.jl @@ -11,6 +11,22 @@ Weak exterior derivative matrix `d^{(k),weak}`, of shape `⟨ d(φ_i dx_J), φ_I dx_{J'} ⟩ = ∫ φ_I det(Γ(dφ_i | dx_J, dx_{J'})) dμ`. +The assembler behind [`grad`](@ref) (`k = 0`) and [`d`](@ref) (`k ≥ 1`) — reach for +those unless you are building an operator by hand. + +# Examples +```jldoctest +julia> W = derivative_weak(function_basis(dg), gamma_mixed(dg.cache), nothing, + measure(dg), 0, n_coefficients(dg)); + +julia> size(W) # (n_coefficients · d, n_function_basis) +(16, 8) + +julia> W ≈ weak(grad(dg)) # exactly what the operator carries +true +``` + +# Arguments - `u` `(n, n0)`: coefficient functions φ. - `gamma_mixed` `(n, d, n_function_basis)`: `Γ(x_j, φ_i)`. - `compound_matrices_k` `(n, Ck, Ck)`: `det(Γ(x_J, x_{J'}))` (needed for `k > 0`). diff --git a/src/operators/differential_operators/geometry_operators.jl b/src/operators/differential_operators/geometry_operators.jl index 06793c9..1e1bc6d 100644 --- a/src/operators/differential_operators/geometry_operators.jl +++ b/src/operators/differential_operators/geometry_operators.jl @@ -11,7 +11,24 @@ _memo(f, dg::DiffusionGeometry, key) = get!(f, dg._op_cache, key) # ── First-order derivatives: d, grad, codifferential, div ────────────────────── -"""Gradient `∇ : A → 𝔛(M)`, mapping functions to vector fields.""" +""" + grad(dg) -> LinearOperator + grad(f::ScalarFunction) -> VectorField + +The gradient `∇ : A → 𝔛(M)`. Memoised on `dg`, so the weak matrix is assembled once. + +# Examples +```jldoctest +julia> grad(dg) +LinearOperator(domain=FunctionSpace(dim=8), codomain=VectorFieldSpace(dim=16), shape=(16, 8)) + +julia> X = grad(f) # ∇cos θ = -sin θ ∂_θ +VectorField(space=VectorFieldSpace(dim=16), shape=(16,), batch_shape=()) + +julia> round(maximum(pointwise_norm(X)); digits=2) # ‖∇cos θ‖ peaks at |sin θ| = 1 +0.97 +``` +""" function grad(dg::DiffusionGeometry) _memo(dg, :grad) do wm = derivative_weak(dg.triple.function_basis, gamma_mixed(dg.cache), nothing, @@ -20,7 +37,34 @@ function grad(dg::DiffusionGeometry) end end -"""Exterior derivative `d : Ωᵏ(M) → Ωᵏ⁺¹(M)`.""" +""" + d(dg, k) -> LinearOperator + d(ω) -> Form + +The exterior derivative `d : Ωᵏ(M) → Ωᵏ⁺¹(M)`. On functions (`k = 0`) it is +[`grad`](@ref) with its index lowered, so `df = (∇f)♭`. + +!!! warning "d² is not zero here" + The discrete exterior derivative does not satisfy `d ∘ d = 0` — the diffusion + complex is not a chain complex. This is expected, and it is why the topology + functions read Betti numbers off a *spectral gap* rather than an exact kernel. + See [`betti_spectrum`](@ref). + +# Examples +```jldoctest +julia> d(dg, 1) +LinearOperator(domain=FormSpace(degree=1, dim=16), codomain=FormSpace(degree=2, dim=8), shape=(8, 16)) + +julia> ω = d(f) # df, a 1-form +Form(space=FormSpace(degree=1, dim=16), shape=(16,), batch_shape=()) + +julia> ω.coeffs ≈ flat(grad(f)).coeffs # df = (∇f)♭ +true + +julia> maximum(abs.(d(ω).coeffs)) < 1e-10 # d²f: small, but not zero +false +``` +""" function d(dg::DiffusionGeometry, k::Integer) @assert 0 <= k < ambient_dim(dg) "Exterior derivative undefined for k=$k in dim $(ambient_dim(dg))." _memo(dg, (:d, Int(k))) do @@ -35,7 +79,24 @@ function d(dg::DiffusionGeometry, k::Integer) end end -"""Codifferential `δ = d* : Ωᵏ(M) → Ωᵏ⁻¹(M)`.""" +""" + codifferential(dg, k) -> LinearOperator + codifferential(ω) -> Form + +The codifferential `δ = d* : Ωᵏ(M) → Ωᵏ⁻¹(M)`, the L²-adjoint of the exterior +derivative: `⟨dα, ω⟩ = ⟨α, δω⟩`. + +# Examples +```jldoctest +julia> codifferential(dg, 1) +LinearOperator(domain=FormSpace(degree=1, dim=16), codomain=FunctionSpace(dim=8), shape=(8, 16)) + +julia> ω = d(dg_function(dg, sin.(θ))); + +julia> isapprox(inner(dg, d(f), ω), inner(dg, f, codifferential(ω)); rtol=1e-8) +true +``` +""" function codifferential(dg::DiffusionGeometry, k::Integer) @assert 1 <= k <= ambient_dim(dg) "Codifferential undefined for degree k=$k." _memo(dg, (:codifferential, Int(k))) do @@ -43,13 +104,52 @@ function codifferential(dg::DiffusionGeometry, k::Integer) end end -"""Divergence `div = -δ ∘ ♭ : 𝔛(M) → A` (here `-∇*`).""" +""" + divergence(dg) -> LinearOperator + divergence(X::VectorField) -> ScalarFunction + +The divergence `div = -∇* : 𝔛(M) → A`, minus the adjoint of the gradient — so +`∫ div(X) f dμ = -∫ g(X, ∇f) dμ`, the integration-by-parts identity, holds by +construction. + +# Examples +`div ∇cos θ = -Δcos θ = -cos θ`: + +```jldoctest +julia> vals = to_pointwise_basis(divergence(grad(f))); + +julia> maximum(abs.(vals .+ cos.(θ))) < 1e-2 +true +``` +""" divergence(dg::DiffusionGeometry) = _memo(dg, :div) do -adjoint(grad(dg)) end # ── Laplacians ───────────────────────────────────────────────────────────────── -"""Up-Laplacian `Δ_up = δ d : Ωᵏ(M) → Ωᵏ(M)`.""" +""" + up_laplacian(dg, k) -> LinearOperator + up_laplacian(ω) -> Form + +The up-Laplacian `Δ_up = δd : Ωᵏ(M) → Ωᵏ(M)`. On functions this is the whole Hodge +Laplacian (there is no `Ω⁻¹` for the down part to use). At top degree it is zero. + +Its kernel is the *closed* forms — half of what a harmonic form must be. That is why +[`betti_spectrum`](@ref) penalises the down-Laplacian on top of it rather than using +either alone. + +# Examples +```jldoctest +julia> up_laplacian(dg, 1) +LinearOperator(domain=FormSpace(degree=1, dim=16), codomain=FormSpace(degree=1, dim=16), shape=(16, 16)) + +julia> up_laplacian(dg, 0)(f).coeffs ≈ laplacian(f).coeffs # on functions, Δ = δd +true + +julia> matrix(up_laplacian(dg, 2)) == zeros(8, 8) # top degree: nothing above +true +``` +""" function up_laplacian(dg::DiffusionGeometry, k::Integer) @assert 0 <= k <= ambient_dim(dg) "Up-Laplacian undefined for degree k=$k." _memo(dg, (:up_laplacian, Int(k))) do @@ -69,7 +169,22 @@ function up_laplacian(dg::DiffusionGeometry, k::Integer) end end -"""Down-Laplacian `Δ_down = d δ : Ωᵏ(M) → Ωᵏ(M)`.""" +""" + down_laplacian(dg, k) -> LinearOperator + down_laplacian(ω) -> Form + +The down-Laplacian `Δ_down = dδ : Ωᵏ(M) → Ωᵏ(M)`. Zero on functions. Its kernel is the +*coclosed* forms. + +# Examples +```jldoctest +julia> down_laplacian(dg, 1) +LinearOperator(domain=FormSpace(degree=1, dim=16), codomain=FormSpace(degree=1, dim=16), shape=(16, 16)) + +julia> matrix(down_laplacian(dg, 0)) == zeros(8, 8) # δ on functions has nowhere to go +true +``` +""" function down_laplacian(dg::DiffusionGeometry, k::Integer) @assert 0 <= k <= ambient_dim(dg) "Down-Laplacian undefined for degree k=$k." _memo(dg, (:down_laplacian, Int(k))) do @@ -78,7 +193,34 @@ function down_laplacian(dg::DiffusionGeometry, k::Integer) end end -"""Hodge Laplacian `Δ = δ d + d δ : Ωᵏ(M) → Ωᵏ(M)`.""" +""" + laplacian(dg, k) -> LinearOperator + laplacian(f::ScalarFunction) -> ScalarFunction + laplacian(ω::Form) -> Form + +The Hodge Laplacian `Δ = δd + dδ : Ωᵏ(M) → Ωᵏ(M)`. Self-adjoint and positive +semi-definite, so its [`spectrum`](@ref) is real and ascending from ~0. + +Sign convention: this is the *geometer's* Laplacian, `Δ = -div grad`, whose eigenvalues +are non-negative. The heat equation is therefore `∂ₜu = -Δu`. + +# Examples +The unit circle has Laplacian eigenvalues `0, 1, 1, 4, 4, …`, and `Δcos θ = cos θ`: + +```jldoctest +julia> round.(abs.(spectrum(laplacian(dg, 0); eigvals_only=true))[1:3]; digits=3) +3-element Vector{Float64}: + 0.0 + 0.999 + 1.0 + +julia> maximum(abs.(to_pointwise_basis(laplacian(f)) .- cos.(θ))) < 1e-2 +true + +julia> is_self_adjoint(laplacian(dg, 1)) +true +``` +""" function laplacian(dg::DiffusionGeometry, k::Integer) @assert 0 <= k <= ambient_dim(dg) "Hodge Laplacian undefined for degree k=$k." _memo(dg, (:laplacian, Int(k))) do @@ -89,7 +231,24 @@ function laplacian(dg::DiffusionGeometry, k::Integer) end # ── Second-order derivatives: Hessian, Lie bracket, Levi-Civita ──────────────── -"""Hessian `Hess : A → Sym²(T*M)`, mapping functions to symmetric (0,2)-tensors.""" +""" + hessian(dg) -> LinearOperator + hessian(f::ScalarFunction) -> Tensor02Sym + +The Hessian `Hess : A → Sym²(T*M)`, the second covariant derivative of a function. + +# Examples +```jldoctest +julia> hessian(dg) +LinearOperator(domain=FunctionSpace(dim=8), codomain=Tensor02SymSpace(dim=24), shape=(24, 8)) + +julia> H = hessian(f) +Tensor02Sym(space=Tensor02SymSpace(dim=24), shape=(24,), batch_shape=()) + +julia> size(H(grad(f), grad(f))) # evaluate it on a pair of vector fields +(60,) +``` +""" function hessian(dg::DiffusionGeometry) _memo(dg, :hessian) do wm = hessian_02_sym_weak(dg.triple.function_basis, hessian_functions(dg.cache), @@ -98,7 +257,26 @@ function hessian(dg::DiffusionGeometry) end end -"""Lie bracket `[·,·] : 𝔛(M) × 𝔛(M) → 𝔛(M)`.""" +""" + lie_bracket(dg) -> BilinearOperator + +The Lie bracket `[·,·] : 𝔛(M) × 𝔛(M) → 𝔛(M)` — the commutator of two vector fields as +derivations. A [`BilinearOperator`](@ref), so it can be called on two fields or +partially applied to one. + +# Examples +```jldoctest +julia> B = lie_bracket(dg); + +julia> X = grad(f); Y = grad(dg_function(dg, sin.(θ))); + +julia> B(X, Y).coeffs ≈ -B(Y, X).coeffs # antisymmetric +true + +julia> maximum(abs.(B(X, X).coeffs)) < 1e-15 # and so [X, X] = 0 +true +``` +""" function lie_bracket(dg::DiffusionGeometry) _memo(dg, :lie_bracket) do wt = lie_bracket_weak(dg.triple.function_basis, dg.triple.immersion_coords, @@ -109,7 +287,26 @@ function lie_bracket(dg::DiffusionGeometry) end end -"""Levi-Civita connection `∇ : 𝔛(M) → Ω⁰²(M)`.""" +""" + levi_civita(dg) -> LinearOperator + levi_civita(X::VectorField) -> Tensor02 + +The Levi-Civita connection `∇ : 𝔛(M) → Ω⁰²(M)`: it sends a vector field `Z` to the +(0,2)-tensor `∇Z`, whose value on `(X, W)` is `g(∇_X Z, W)`. This is the ingredient +[`riemann_curvature`](@ref) is assembled from. + +# Examples +```jldoctest +julia> levi_civita(dg) +LinearOperator(domain=VectorFieldSpace(dim=16), codomain=Tensor02Space(dim=32), shape=(32, 16)) + +julia> Z = grad(f); ∇Z = levi_civita(Z) +Tensor02(space=Tensor02Space(dim=32), shape=(32,), batch_shape=()) + +julia> ∇Z(grad(f)) # ∇_X Z, a vector field again +VectorField(space=VectorFieldSpace(dim=16), shape=(16,), batch_shape=()) +``` +""" function levi_civita(dg::DiffusionGeometry) _memo(dg, :levi_civita) do wm = levi_civita_02_weak(dg.triple.function_basis, gamma_mixed(dg.cache), @@ -125,6 +322,21 @@ end Riemann curvature `R(X, Y, Z, W)` as pointwise function values `(batch..., n)`, following `R(X,Y,Z,W) = g(∇_X∇_Y Z, W) - g(∇_Y∇_X Z, W) - g(∇_[X,Y] Z, W)`. + +# Examples +The tensor's antisymmetry in its first two slots comes out of the construction: + +```jldoctest +julia> X = grad(f); Y = grad(dg_function(dg, sin.(θ))); + +julia> R = riemann_curvature(dg, X, Y, X, Y); + +julia> size(R) # a value per point +(60,) + +julia> R ≈ -riemann_curvature(dg, Y, X, X, Y) # R(X,Y,·,·) = -R(Y,X,·,·) +true +``` """ function riemann_curvature(dg::DiffusionGeometry, X::VectorField, Y::VectorField, Z::VectorField, W::VectorField) @@ -139,7 +351,28 @@ function riemann_curvature(dg::DiffusionGeometry, X::VectorField, Y::VectorField return term1 - term2 - term3 end -"""Sectional curvature `K(X, Y)` for the plane spanned by `X`, `Y`, shape `(..., n)`.""" +""" + sectional_curvature(dg, X, Y) -> Array + +The sectional curvature `K(X, Y) = R(X,Y,X,Y) / (g(X,X)g(Y,Y) - g(X,Y)²)` of the plane +spanned by `X` and `Y`, as pointwise values `(batch..., n)`. + +!!! warning "Noisy where the plane degenerates" + The denominator vanishes as `X` and `Y` become parallel (and it is only floored at + `eps`), so the pointwise values blow up wherever the two fields fail to span. Read + the median over a region, not a single point, and expect the curvature of a coarse + sample to be badly scaled. + +# Examples +```jldoctest +julia> u = dg_function(dg3, sphere[:, 1]); v = dg_function(dg3, sphere[:, 2]); + +julia> K = sectional_curvature(dg3, grad(u), grad(v)); + +julia> size(K) +(200,) +``` +""" function sectional_curvature(dg::DiffusionGeometry, X::VectorField, Y::VectorField) for V in (X, Y) @assert V.space == vector_field_space(dg) "Sectional curvature arguments must be VectorFields" diff --git a/src/operators/differential_operators/hessian.jl b/src/operators/differential_operators/hessian.jl index 160c9cc..a1fa4a7 100644 --- a/src/operators/differential_operators/hessian.jl +++ b/src/operators/differential_operators/hessian.jl @@ -9,6 +9,19 @@ using OMEinsum: @ein_str, @optein_str Pointwise Hessian tensor `H(φ_I)(∇x_i, ∇x_j)`, shape `(n, d, d, n0)`: `H(f)(∇x_i,∇x_j) = ½(Γ(x_i, Γ(x_j, f)) + Γ(x_j, Γ(x_i, f)) − Γ(f, Γ(x_i, x_j)))`. `cdc(f, h)` is the carré du champ callback. + +The raw (unregularised) form; `hessian_functions(cache)` is the memoised, regularised +one the [`hessian`](@ref) operator actually uses. + +# Examples +```jldoctest +julia> H = hessian_functions(function_basis(dg), immersion_coords(dg), + gamma_coords(dg.cache), gamma_mixed(dg.cache), + (a, b) -> cdc(dg.triple, a, b)); + +julia> size(H) # (n, d, d, n_function_basis) +(60, 2, 2, 8) +``` """ function hessian_functions(u::AbstractMatrix, coords::AbstractMatrix, gamma_coords::AbstractArray{<:Any,3}, @@ -21,7 +34,16 @@ end """ hessian_coords(coords, gamma_coords, cdc) -> Array -Pointwise coordinate Hessian `H(x_k)(∇x_i, ∇x_j)`, shape `(n, d, d, d)`. +Pointwise coordinate Hessian `H(x_k)(∇x_i, ∇x_j)`, shape `(n, d, d, d)` — the second +fundamental form of the immersion, and the Christoffel data of +[`levi_civita`](@ref). The memoised, regularised version is `hessian_coords(cache)`. + +# Examples +```jldoctest +julia> size(hessian_coords(immersion_coords(dg), gamma_coords(dg.cache), + (a, b) -> cdc(dg.triple, a, b))) +(60, 2, 2, 2) +``` """ function hessian_coords(coords::AbstractMatrix, gamma_coords::AbstractArray{<:Any,3}, cdc) gc = cdc(coords, gamma_coords) # Γ(x_i, Γ(x_j, x_k)) (n, d, d, d) @@ -33,6 +55,16 @@ end Weak Hessian on functions, `H : A → Ω¹ ⊗ Ω¹`, shape `(n1 * d², n0)`. `hessian_matrix` is `H(φ_i)(∇x_j, ∇x_k)` of shape `(n, d, d, n0)`. + +The full-tensor form. [`hessian`](@ref) uses the symmetric one +([`hessian_02_sym_weak`](@ref)) instead, since the Hessian is symmetric. + +# Examples +```jldoctest +julia> size(hessian_02_weak(function_basis(dg), hessian_functions(dg.cache), + measure(dg), n_coefficients(dg))) +(32, 8) +``` """ function hessian_02_weak(u::AbstractMatrix, hessian_matrix::AbstractArray{<:Any,4}, measure::AbstractVector, n_coefficients::Integer) @@ -47,6 +79,20 @@ end Weak Hessian as a symmetric (0,2)-tensor, shape `(n1 * d_sym, n0)`. Off-diagonal symmetric-basis entries are doubled to match the symmetric basis convention. + +The assembler behind the [`hessian`](@ref) operator. + +# Examples +```jldoctest +julia> W = hessian_02_sym_weak(function_basis(dg), hessian_functions(dg.cache), + measure(dg), n_coefficients(dg)); + +julia> size(W) # (n_coefficients · d_sym, n_function_basis) +(24, 8) + +julia> W ≈ weak(hessian(dg)) +true +``` """ function hessian_02_sym_weak(u::AbstractMatrix, hessian_matrix::AbstractArray{T,4}, measure::AbstractVector, n_coefficients::Integer) where {T} diff --git a/src/operators/differential_operators/laplacian.jl b/src/operators/differential_operators/laplacian.jl index 14d199a..9b16ce2 100644 --- a/src/operators/differential_operators/laplacian.jl +++ b/src/operators/differential_operators/laplacian.jl @@ -12,6 +12,23 @@ Weak up-Laplacian (up-Hodge energy) matrix on `k`-forms, shape `(n1 Ck, n1 Ck)`. `⟨d(φ_i dx_J), d(φ_I dx_{J'})⟩ = ∫ g(dφ_i dx_J, dφ_I dx_{J'}) dμ`, evaluated via the Schur determinant formula `det([[a,b],[c,D]]) = det(D)(a − b D⁻¹ c)`. +The assembler behind [`up_laplacian`](@ref). Building `δd` directly, rather than +composing the two operators, avoids a Gram inverse in the middle — which is why the two +routes do not agree exactly on the modes the basis barely resolves. + +# Examples +```jldoctest +julia> W = up_delta_weak(gamma_functions(dg.cache), gamma_mixed(dg.cache), + gamma_coords(dg.cache), nothing, nothing, measure(dg), 0); + +julia> size(W) +(8, 8) + +julia> W ≈ weak(up_laplacian(dg, 0)) +true +``` + +# Arguments - `gamma_functions` `(n, n0, n0)`: `Γ(φ_i, φ_I)`. - `gamma_mixed` `(n, d, n0)`: `Γ(x_j, φ_i)`. - `gamma_coords` `(n, d, d)`. diff --git a/src/operators/differential_operators/levi_civita.jl b/src/operators/differential_operators/levi_civita.jl index 57b4e62..daa0d5b 100644 --- a/src/operators/differential_operators/levi_civita.jl +++ b/src/operators/differential_operators/levi_civita.jl @@ -10,6 +10,22 @@ Weak Levi-Civita connection matrix `∇^{weak}_{Ikl, ij}`, shape `(n1 d², n1 d) `⟨∇(φ_i ∇x_j), φ_I dx_k ⊗ dx_l⟩ = ∫ φ_I Γ(x_k, φ_i) Γ(x_l, x_j) dμ + ∫ φ_i φ_I H(x_j)(∇x_k, ∇x_l) dμ`. +The assembler behind the [`levi_civita`](@ref) operator. + +# Examples +```jldoctest +julia> W = levi_civita_02_weak(function_basis(dg), gamma_mixed(dg.cache), + gamma_coords(dg.cache), hessian_coords(dg.cache), + measure(dg), n_coefficients(dg)); + +julia> size(W) # (n_coefficients · d², n_coefficients · d) +(32, 16) + +julia> W ≈ weak(levi_civita(dg)) +true +``` + +# Arguments - `gamma_mixed` `(n, d, n0)`, `gamma_coords` `(n, d, d)`. - `hessian_coords` `(n, d, d, d)`: coordinate Hessian `H(x_K)(∇x_i, ∇x_j)`. """ diff --git a/src/operators/differential_operators/lie_bracket.jl b/src/operators/differential_operators/lie_bracket.jl index 3aa8fd1..6881d60 100644 --- a/src/operators/differential_operators/lie_bracket.jl +++ b/src/operators/differential_operators/lie_bracket.jl @@ -8,6 +8,22 @@ using OMEinsum: @ein_str, @optein_str Weak Lie bracket operator tensor `Lie^{weak}_{IJ, i1j1, i2j2}`, shape `(n1 d, n1 d, n1 d)`. `cdc(f, h)` is the carré du champ callback. + +The assembler behind the [`lie_bracket`](@ref) operator. This is the package's one +three-index contraction, and its cost grows fast in `n_coefficients`. + +# Examples +```jldoctest +julia> W = lie_bracket_weak(function_basis(dg), immersion_coords(dg), + gamma_coords(dg.cache), measure(dg), n_coefficients(dg), + (a, b) -> cdc(dg.triple, a, b)); + +julia> size(W) +(16, 16, 16) + +julia> W ≈ weak(lie_bracket(dg)) +true +``` """ function lie_bracket_weak(u::AbstractMatrix, coords::AbstractMatrix, gamma_coords::AbstractArray{<:Any,3}, diff --git a/src/operators/tensor_actions.jl b/src/operators/tensor_actions.jl index 5c95470..becdeda 100644 --- a/src/operators/tensor_actions.jl +++ b/src/operators/tensor_actions.jl @@ -13,7 +13,19 @@ import LinearAlgebra vf_operator(X) -> LinearOperator The vector field as a directional-derivative operator `A → A`, `f ↦ X(f) = g(X, ∇f)` -(unbatched `X`). +(unbatched `X`). A vector field *is* a derivation; this is that derivation as a matrix. +Calling `X(f)` applies it. + +# Examples +```jldoctest +julia> X = grad(f); + +julia> vf_operator(X) +LinearOperator(domain=FunctionSpace(dim=8), codomain=FunctionSpace(dim=8), shape=(8, 8)) + +julia> vf_operator(X)(f).coeffs ≈ X(f).coeffs +true +``` """ function vf_operator(X::VectorField) @assert isempty(batch_shape(X)) "operator only supports unbatched VectorFields." @@ -25,13 +37,50 @@ function vf_operator(X::VectorField) return LinearOperator(function_space(dg), function_space(dg); weak_matrix=weak_matrix) end -"""Apply the vector field as a directional derivative: `X(f) = g(X, ∇f)`.""" +""" + (X::VectorField)(f::ScalarFunction) -> ScalarFunction + +Apply the vector field as a directional derivative: `X(f) = g(X, ∇f)`. + +# Examples +```jldoctest +julia> X = grad(f); + +julia> X(f) # ∇f(f) = ‖∇f‖², up to the basis projection +ScalarFunction(space=FunctionSpace(dim=8), shape=(8,), batch_shape=()) + +julia> Y = grad(dg_function(dg, sin.(θ))); + +julia> isapprox(inner(dg, X(f), f), inner(dg, f, X(f)); rtol=1e-8) +true +``` +""" (X::VectorField)(f::ScalarFunction) = vf_operator(X)(f) -"""Divergence of the vector field, `div X`.""" +""" + divergence(X::VectorField) -> ScalarFunction + +The divergence of the field, `div X` — shorthand for `divergence(geometry(X))(X)`. + +# Examples +```jldoctest +julia> maximum(abs.(to_pointwise_basis(divergence(grad(f))) .+ cos.(θ))) < 1e-2 +true +``` +""" divergence(X::VectorField) = divergence(geometry(X))(X) -"""Levi-Civita covariant derivative `∇X` as a (0,2)-tensor.""" +""" + levi_civita(X::VectorField) -> Tensor02 + +The covariant derivative `∇X` of the field, as a (0,2)-tensor: `∇X(Y, W) = g(∇_Y X, W)`. + +# Examples +```jldoctest +julia> levi_civita(grad(f)) +Tensor02(space=Tensor02Space(dim=32), shape=(32,), batch_shape=()) +``` +""" levi_civita(X::VectorField) = levi_civita(geometry(X))(X) # ── Tensor02 as an operator / bilinear form ──────────────────────────────────── @@ -39,7 +88,18 @@ levi_civita(X::VectorField) = levi_civita(geometry(X))(X) t02_operator(α) -> LinearOperator Operator form `α^{op} : 𝔛(M) → 𝔛(M)` defined by `⟨α^{op}(X), Y⟩ = ∫ α(X, Y) dμ` -(unbatched `α`). +(unbatched `α`). Calling a `Tensor02` on a single vector field applies it. + +# Examples +```jldoctest +julia> α = levi_civita(grad(f)); + +julia> t02_operator(α) +LinearOperator(domain=VectorFieldSpace(dim=16), codomain=VectorFieldSpace(dim=16), shape=(16, 16)) + +julia> t02_operator(α)(grad(f)).coeffs ≈ α(grad(f)).coeffs +true +``` """ function t02_operator(α::Tensor02) @assert isempty(batch_shape(α)) "operator only supports unbatched Tensor02 objects." @@ -60,7 +120,21 @@ end (α::Tensor02)(X, Y) -> Array Evaluate the (0,2)-tensor pointwise on a pair of vector fields, returning -`(batch..., n)` function values (regularised). +`(batch..., n)` function values (regularised). With one argument it acts as an +operator instead (see [`t02_operator`](@ref)). + +# Examples +```jldoctest +julia> α = d(f) * d(dg_function(dg, sin.(θ))); # the tensor product df ⊗ dh + +julia> X = grad(f); Y = grad(dg_function(dg, sin.(θ))); + +julia> size(α(X, Y)) +(60,) + +julia> α(X) # one argument → a vector field +VectorField(space=VectorFieldSpace(dim=16), shape=(16,), batch_shape=()) +``` """ function (α::Tensor02)(X::VectorField, Y::VectorField) dg = geometry(α) @@ -84,7 +158,14 @@ end to_ambient(α::Tensor02) -> Array Ambient-coordinate representation `(n, d, d)` obtained by raising both indices with -the (regularised) coordinate carré du champ. +the (regularised) coordinate carré du champ — the tensor written as a matrix field in +the coordinates the data came in. + +# Examples +```jldoctest +julia> size(to_ambient(full_tensor(hessian(f)))) +(60, 2, 2) +``` """ function to_ambient(α::Tensor02) @assert isempty(batch_shape(α)) "to_ambient only supports unbatched Tensor02 objects." @@ -106,6 +187,17 @@ t02_operator(S::Tensor02Sym) = t02_operator(full_tensor(S)) Interior product `ω(X) = g(ω♯, X)`, as pointwise values of shape `(batch..., n)`. Only defined for 1-forms. + +# Examples +```jldoctest +julia> ω = d(f); X = grad(f); + +julia> ω(X) ≈ g(dg, ω, flat(X)) # by definition +true + +julia> size(ω(X)) +(60,) +``` """ function (ω::Form)(X::VectorField) @assert degree(ω) == 1 "Only 1-forms can act on vector fields." @@ -118,34 +210,142 @@ end # Each is the corresponding `dg` accessor evaluated at the tensor; the operators # are memoised on `dg`, so repeated calls reuse one weak matrix. -"""Gradient `∇f`, a vector field.""" +""" + grad(f::ScalarFunction) -> VectorField + +The gradient `∇f` — shorthand for `grad(geometry(f))(f)`, using the memoised operator. + +# Examples +```jldoctest +julia> grad(f) +VectorField(space=VectorFieldSpace(dim=16), shape=(16,), batch_shape=()) +``` +""" grad(f::ScalarFunction) = grad(geometry(f))(f) -"""Exterior derivative `df`, a 1-form.""" +""" + d(f::ScalarFunction) -> Form + +The exterior derivative `df`, a 1-form — the gradient with its index lowered. + +# Examples +```jldoctest +julia> d(f).coeffs ≈ flat(grad(f)).coeffs +true +``` +""" d(f::ScalarFunction) = d(geometry(f), 0)(f) -"""Up-Laplacian `Δ_up f = δ d f`.""" +""" + up_laplacian(f::ScalarFunction) -> ScalarFunction + +The up-Laplacian `Δ_up f = δdf`. On functions this is the whole Laplacian. + +# Examples +```jldoctest +julia> up_laplacian(f).coeffs ≈ laplacian(f).coeffs +true +``` +""" up_laplacian(f::ScalarFunction) = up_laplacian(geometry(f), 0)(f) -"""Laplacian `Δf`. On functions the down-Laplacian vanishes, so this is `Δ_up`.""" +""" + laplacian(f::ScalarFunction) -> ScalarFunction + +The Laplacian `Δf = δdf` (the down-Laplacian vanishes on functions). + +# Examples +`Δcos θ = cos θ` on the unit circle: + +```jldoctest +julia> maximum(abs.(to_pointwise_basis(laplacian(f)) .- cos.(θ))) < 1e-2 +true +``` +""" laplacian(f::ScalarFunction) = up_laplacian(geometry(f), 0)(f) -"""Hessian `Hess f`, a symmetric (0,2)-tensor.""" +""" + hessian(f::ScalarFunction) -> Tensor02Sym + +The Hessian `Hess f`, a symmetric (0,2)-tensor. + +# Examples +```jldoctest +julia> hessian(f) +Tensor02Sym(space=Tensor02SymSpace(dim=24), shape=(24,), batch_shape=()) +``` +""" hessian(f::ScalarFunction) = hessian(geometry(f))(f) -"""Exterior derivative `dω : Ωᵏ → Ωᵏ⁺¹`.""" +""" + d(ω::Form) -> Form + +The exterior derivative `dω : Ωᵏ → Ωᵏ⁺¹`. + +# Examples +```jldoctest +julia> d(d(f)) # a 2-form on the circle's ambient frame +Form(space=FormSpace(degree=2, dim=8), shape=(8,), batch_shape=()) +``` +""" d(ω::Form) = d(geometry(ω), degree(ω))(ω) -"""Codifferential `δω : Ωᵏ → Ωᵏ⁻¹`.""" +""" + codifferential(ω::Form) -> Form or ScalarFunction + +The codifferential `δω : Ωᵏ → Ωᵏ⁻¹`. On a 1-form it lands in the functions. + +# Examples +`δdf = Δf`: + +```jldoctest +julia> isapprox(codifferential(d(f)).coeffs, laplacian(f).coeffs; atol=1e-5) +true +``` + +(Only to `atol`: `laplacian` assembles `δd` in one weak form, while composing the two +operators passes through a Gram pseudo-inverse. They agree on everything the basis +resolves.) +""" codifferential(ω::Form) = codifferential(geometry(ω), degree(ω))(ω) -"""Up-Laplacian `Δ_up ω = δ d ω`.""" +""" + up_laplacian(ω::Form) -> Form + +The up-Laplacian `Δ_up ω = δdω`. + +# Examples +```jldoctest +julia> up_laplacian(d(f)) +Form(space=FormSpace(degree=1, dim=16), shape=(16,), batch_shape=()) +``` +""" up_laplacian(ω::Form) = up_laplacian(geometry(ω), degree(ω))(ω) -"""Down-Laplacian `Δ_down ω = d δ ω`.""" +""" + down_laplacian(ω::Form) -> Form + +The down-Laplacian `Δ_down ω = dδω`. + +# Examples +```jldoctest +julia> down_laplacian(d(f)) +Form(space=FormSpace(degree=1, dim=16), shape=(16,), batch_shape=()) +``` +""" down_laplacian(ω::Form) = down_laplacian(geometry(ω), degree(ω))(ω) -"""Hodge Laplacian `Δω = (δd + dδ) ω`.""" +""" + laplacian(ω::Form) -> Form + +The Hodge Laplacian `Δω = (δd + dδ)ω`. + +# Examples +```jldoctest +julia> laplacian(d(f)).coeffs ≈ (up_laplacian(d(f)) + down_laplacian(d(f))).coeffs +true +``` +""" laplacian(ω::Form) = laplacian(geometry(ω), degree(ω))(ω) # ── Hodge decomposition ──────────────────────────────────────────────────────── @@ -154,6 +354,14 @@ laplacian(ω::Form) = laplacian(geometry(ω), degree(ω))(ω) Hodge decomposition of a function: `f = δβ + h` with `β = coexact_potential` a 1-form and `h = harmonic_part`. There is no exact part, since `Ω⁻¹` is trivial. + +# Examples +```jldoctest +julia> β, h = hodge_decomposition(f); + +julia> isapprox((codifferential(β) + h).coeffs, f.coeffs; atol=1e-6) +true +``` """ function hodge_decomposition(f::ScalarFunction) dg = geometry(f) @@ -171,6 +379,27 @@ top degree (`k = dim`) there is no coexact part and `coexact_potential` is `noth The potentials come from spectral pseudo-inverses of the Laplacians, so they are determined only up to the respective kernels. + +# Examples +The decomposition reconstructs the form it came from: + +```jldoctest +julia> ω = d(f) + d(dg_function(dg, sin.(θ))); + +julia> α, β, h = hodge_decomposition(ω); + +julia> isapprox((d(α) + codifferential(β) + h).coeffs, ω.coeffs; atol=1e-6) +true +``` + +An exact form is all exact part — its harmonic component vanishes: + +```jldoctest +julia> _, _, h = hodge_decomposition(d(f)); + +julia> l2_norm(h) < 1e-5 +true +``` """ function hodge_decomposition(ω::Form) dg = geometry(ω) @@ -244,7 +473,16 @@ end Ambient polyvector representation `(n, D, …, D)` (k factors of the ambient dimension `D`) of an unbatched k-form, obtained by raising all k indices with the -ambient carré du champ. +ambient carré du champ. This is how a form becomes something you can draw. + +# Examples +```jldoctest +julia> size(to_ambient(d(f))) # a 1-form → one arrow per point +(60, 2) + +julia> size(to_ambient(d(dg_function(dg3, sphere[:, 1])))) +(200, 3) +``` """ function to_ambient(ω::Form) @assert isempty(batch_shape(ω)) "to_ambient only supports unbatched Form objects." @@ -259,13 +497,41 @@ function to_ambient(ω::Form) return expanded end -"""Ambient representation of a scalar function — the pointwise values themselves.""" +""" + to_ambient(f::ScalarFunction) -> Array + +The ambient representation of a function is just its pointwise values. + +# Examples +```jldoctest +julia> to_ambient(f) == to_pointwise_basis(f) +true +``` +""" to_ambient(f::ScalarFunction) = to_pointwise_basis(f) """ to_ambient(X::VectorField) -> Array -Ambient quiver representation `(n, D)`: the 1-form `X♭` pushed to ambient coordinates. +Ambient quiver representation `(n, D)`: the 1-form `X♭` pushed to ambient coordinates — +the arrows a quiver plot draws. Inverted by +[`vector_field_from_reconstruction`](@ref). + +# Examples +The gradient of `cos θ` on the circle points along the tangent, so the arrows are +orthogonal to the radius: + +```jldoctest +julia> arrows = to_ambient(grad(f)); + +julia> size(arrows) +(60, 2) + +julia> radial = sum(arrows .* circle; dims=2); # ⟨arrow, outward normal⟩ + +julia> maximum(abs.(radial)) < 0.02 # tangent to ~2% of the arrow length +true +``` """ to_ambient(X::VectorField) = to_ambient(flat(X)) @@ -278,6 +544,19 @@ The wedge product with a fixed k-form `a`, as an operator `Ωˡ(M) → Ωᵏ⁺ (The Python reference returns the bare coefficient matrix; here it is wrapped as the `LinearOperator` its name promises — that matrix is `matrix(wedge_operator(a, l))`.) + +# Examples +```jldoctest +julia> α = d(f); + +julia> W = wedge_operator(α, 1) +LinearOperator(domain=FormSpace(degree=1, dim=16), codomain=FormSpace(degree=2, dim=8), shape=(8, 16)) + +julia> β = d(dg_function(dg, sin.(θ))); + +julia> W(β).coeffs ≈ wedge(α, β).coeffs # the operator is the product +true +``` """ function wedge_operator(a::Form, l::Integer) @assert isempty(batch_shape(a)) "wedge_operator only supports unbatched Forms." diff --git a/src/operators/types/bilinear.jl b/src/operators/types/bilinear.jl index 47a8c99..10abc77 100644 --- a/src/operators/types/bilinear.jl +++ b/src/operators/types/bilinear.jl @@ -13,6 +13,27 @@ using OMEinsum: @ein_str, @optein_str Bilinear operator `B : domain_a × domain_b → codomain`. Provide at least one of the weak / strong component 3-tensors `(codim, dim_a, dim_b)`. + +[`lie_bracket`](@ref) is the one the package builds. Call it with two tensors to +evaluate the form, or with one to partially apply on the left and get a +[`LinearOperator`](@ref) back. + +# Examples +```jldoctest +julia> B = lie_bracket(dg) +BilinearOperator(domain_a=VectorFieldSpace(dim=16), domain_b=VectorFieldSpace(dim=16), codomain=VectorFieldSpace(dim=16), component_shape=(16, 16, 16)) + +julia> X = grad(f); Y = grad(dg_function(dg, sin.(θ))); + +julia> B(X, Y) # two arguments → a vector field +VectorField(space=VectorFieldSpace(dim=16), shape=(16,), batch_shape=()) + +julia> B(X) # one → the operator [X, ·] +LinearOperator(domain=VectorFieldSpace(dim=16), codomain=VectorFieldSpace(dim=16), shape=(16, 16)) + +julia> B(X, Y).coeffs ≈ -B(Y, X).coeffs # the bracket is antisymmetric +true +``` """ mutable struct BilinearOperator domain_a::AbstractTensorSpace @@ -30,6 +51,17 @@ function BilinearOperator(domain_a::AbstractTensorSpace, domain_b::AbstractTenso return BilinearOperator(domain_a, domain_b, codomain, weak_tensor, strong_tensor) end +""" + component_shape(B::BilinearOperator) -> Tuple + +Shape of the operator's component 3-tensor: `(codomain, domain_a, domain_b)`. + +# Examples +```jldoctest +julia> component_shape(lie_bracket(dg)) +(16, 16, 16) +``` +""" component_shape(B::BilinearOperator) = (space_dim(B.codomain), space_dim(B.domain_a), space_dim(B.domain_b)) @@ -38,14 +70,37 @@ function Base.show(io::IO, B::BilinearOperator) ", codomain=", B.codomain, ", component_shape=", component_shape(B), ")") end -"""Weak-form component 3-tensor.""" +""" + weak(B::BilinearOperator) -> Array + +The bilinear operator's weak-form 3-tensor, `(codim, dim_a, dim_b)` — the form it is +assembled in. See [`strong`](@ref) for the one that evaluates. + +# Examples +```jldoctest +julia> size(weak(lie_bracket(dg))) +(16, 16, 16) +``` +""" function weak(B::BilinearOperator) B._weak !== nothing && return B._weak B._weak = ein"iC,CAB->iAB"(gram(B.codomain), B._strong) return B._weak end -"""Strong-form component 3-tensor.""" +""" + strong(B::BilinearOperator) -> Array + +The bilinear operator's strong-form 3-tensor: the one that contracts against two +coefficient vectors to give the result's coefficients. Derived from [`weak`](@ref) +through the codomain's [`gram_inv`](@ref). + +# Examples +```jldoctest +julia> size(strong(lie_bracket(dg))) +(16, 16, 16) +``` +""" function strong(B::BilinearOperator) B._strong !== nothing && return B._strong B._strong = ein"iC,CAB->iAB"(gram_inv(B.codomain), B._weak) @@ -53,7 +108,26 @@ function strong(B::BilinearOperator) end # ── Application ──────────────────────────────────────────────────────────────── -"""Partial application `B(x, ·) : domain_b → codomain` (unbatched `x`).""" +""" + partial_apply(B, x) -> LinearOperator + B(x) + +Fix the left argument: the linear operator `B(x, ·) : domain_b → codomain`. `x` must be +unbatched. + +# Examples +`[X, ·]` as an operator, which agrees with the two-argument call: + +```jldoctest +julia> B = lie_bracket(dg); X = grad(f); Y = grad(dg_function(dg, sin.(θ))); + +julia> ad_X = partial_apply(B, X) +LinearOperator(domain=VectorFieldSpace(dim=16), codomain=VectorFieldSpace(dim=16), shape=(16, 16)) + +julia> ad_X(Y).coeffs ≈ B(X, Y).coeffs +true +``` +""" function partial_apply(B::BilinearOperator, x::AbstractTensor) @assert x.space == B.domain_a "Domain A mismatch" @assert isempty(batch_shape(x)) "Batched partial application not supported." @@ -61,7 +135,20 @@ function partial_apply(B::BilinearOperator, x::AbstractTensor) return LinearOperator(B.domain_b, B.codomain; weak_matrix=weak_matrix) end -"""Full application `B(x, y) → codomain` (broadcasting over batch axes).""" +""" + full_apply(B, x, y) -> AbstractTensor + B(x, y) + +Evaluate the bilinear form on both arguments, broadcasting over their batch axes. + +# Examples +```jldoctest +julia> B = lie_bracket(dg); X = grad(f); + +julia> maximum(abs.(full_apply(B, X, X).coeffs)) < 1e-15 # [X, X] = 0 +true +``` +""" function full_apply(B::BilinearOperator, x::AbstractTensor, y::AbstractTensor) @assert x.space == B.domain_a "Domain A mismatch" @assert y.space == B.domain_b "Domain B mismatch" @@ -78,7 +165,25 @@ end (B::BilinearOperator)(x::AbstractTensor, y::AbstractTensor) = full_apply(B, x, y) # ── Transpose ────────────────────────────────────────────────────────────────── -"""Transpose `Bᵀ(y, x) = B(x, y)` (swaps the two input slots).""" +""" + transpose_operator(B::BilinearOperator) -> BilinearOperator + B' + +Swap the two input slots: `Bᵀ(y, x) = B(x, y)`. + +# Examples +The Lie bracket is antisymmetric, so its transpose is its negation: + +```jldoctest +julia> B = lie_bracket(dg); X = grad(f); Y = grad(dg_function(dg, sin.(θ))); + +julia> transpose_operator(B)(X, Y).coeffs ≈ B(Y, X).coeffs +true + +julia> transpose_operator(B)(X, Y).coeffs ≈ -B(X, Y).coeffs +true +``` +""" function transpose_operator(B::BilinearOperator) if B._weak !== nothing return BilinearOperator(B.domain_b, B.domain_a, B.codomain; diff --git a/src/operators/types/direct_sum.jl b/src/operators/types/direct_sum.jl index 0fbbeca..1d7d02e 100644 --- a/src/operators/types/direct_sum.jl +++ b/src/operators/types/direct_sum.jl @@ -15,6 +15,26 @@ _sum_spaces(spaces) = reduce(+, spaces) Assemble a rectangular grid of operators, e.g. `block([[A, B], [C, D]])`. Every operator in column `j` must share a domain, and every operator in row `i` a codomain; the result maps `⊕ⱼ domainⱼ → ⊕ᵢ codomainᵢ`. + +The way to build a coupled system — a mixed-degree operator on functions *and* vector +fields is one `LinearOperator` on their direct sum, with a spectrum and an inverse like +any other. + +# Examples +The identity on `A ⊕ 𝔛(M)`, assembled from four blocks: + +```jldoctest +julia> A, V = function_space(dg), vector_field_space(dg); + +julia> M = block([[identity_operator(A) zero_operator(V, A)], + [zero_operator(A, V) identity_operator(V)]]) +LinearOperator(domain=DirectSumSpace(spaces=[FunctionSpace, VectorFieldSpace], dim=24), codomain=DirectSumSpace(spaces=[FunctionSpace, VectorFieldSpace], dim=24), shape=(24, 24)) + +julia> e = pack(A + V, f, grad(f)); + +julia> M(e).coeffs ≈ e.coeffs +true +``` """ function block(block_rows) @assert !isempty(block_rows) "Provide at least one row of operator blocks" @@ -50,13 +70,35 @@ function block(block_rows) return LinearOperator(full_domain, full_codomain; weak_matrix=full_weak) end -"""Horizontally stack operators sharing a codomain: `[A B]`.""" +""" + hstack(blocks) -> LinearOperator + +Stack operators side by side: `[A B] : domain_A ⊕ domain_B → codomain`. They must share +a codomain. + +# Examples +```jldoctest +julia> hstack([grad(dg), identity_operator(vector_field_space(dg))]) +LinearOperator(domain=DirectSumSpace(spaces=[FunctionSpace, VectorFieldSpace], dim=24), codomain=VectorFieldSpace(dim=16), shape=(16, 24)) +``` +""" function hstack(blocks) @assert !isempty(blocks) "hstack requires at least one operator" return block([collect(blocks)]) end -"""Vertically stack operators sharing a domain: `[A; B]`.""" +""" + vstack(blocks) -> LinearOperator + +Stack operators on top of each other: `[A; B] : domain → codomain_A ⊕ codomain_B`. They +must share a domain. + +# Examples +```jldoctest +julia> vstack([grad(dg), grad(dg)]) +LinearOperator(domain=FunctionSpace(dim=8), codomain=DirectSumSpace(spaces=[VectorFieldSpace, VectorFieldSpace], dim=32), shape=(32, 8)) +``` +""" function vstack(blocks) @assert !isempty(blocks) "vstack requires at least one operator" return block([[op] for op in blocks]) diff --git a/src/operators/types/linear.jl b/src/operators/types/linear.jl index add137e..4be2a62 100644 --- a/src/operators/types/linear.jl +++ b/src/operators/types/linear.jl @@ -14,7 +14,19 @@ using LinearAlgebra: Hermitian, eigen, I real_if_close(A; tol=100) -> Array Return the real part of `A` if every imaginary component is within `tol·eps` -(numpy `real_if_close` semantics), otherwise return `A` unchanged. +(numpy `real_if_close` semantics), otherwise return `A` unchanged. Used to clean up the +round-off imaginary parts an eigensolver leaves on a spectrum that is really real. + +# Examples +```jldoctest +julia> real_if_close([1.0 + 1e-16im, 2.0 - 1e-17im]) +2-element Vector{Float64}: + 1.0 + 2.0 + +julia> real_if_close(1.0 + 0.5im) # genuinely complex: left alone +1.0 + 0.5im +``` """ function real_if_close(A::AbstractArray{<:Complex}; tol::Real=100) thresh = tol * eps(real(eltype(A))) @@ -32,6 +44,38 @@ real_if_close(z::Real; tol::Real=100) = z Linear operator `L : domain → codomain`. Provide at least one of the weak (bilinear form ⟨w, Lv⟩) or strong (coefficient map) matrices; the other is derived lazily from the codomain Gram matrix. + +Every differential operator in the package is one of these: [`grad`](@ref), [`d`](@ref), +the Laplacians, [`hessian`](@ref), [`levi_civita`](@ref). They compose with `∘` (or +`*`), add and scale, take adjoints with `'`, and apply to a tensor by calling them. + +# Examples +```jldoctest +julia> L = laplacian(dg, 0) +LinearOperator(domain=FunctionSpace(dim=8), codomain=FunctionSpace(dim=8), shape=(8, 8)) + +julia> L(f) # apply it +ScalarFunction(space=FunctionSpace(dim=8), shape=(8,), batch_shape=()) + +julia> maximum(abs.(to_pointwise_basis(L(f)) .- cos.(θ))) < 1e-2 # Δcos θ = cos θ +true + +julia> grad(dg)' # the adjoint maps the other way +LinearOperator(domain=VectorFieldSpace(dim=16), codomain=FunctionSpace(dim=8), shape=(8, 16)) +``` + +Composition follows the maths, right to left. `div ∘ grad` and `-Δ` agree wherever the +basis actually resolves the function (they part company only on the top modes, where +the truncation bites): + +```jldoctest +julia> Δf = laplacian(f); + +julia> DGf = (divergence(dg) ∘ grad(dg))(f); + +julia> maximum(abs.(to_pointwise_basis(DGf) .+ to_pointwise_basis(Δf))) < 1e-5 +true +``` """ mutable struct LinearOperator domain::AbstractTensorSpace @@ -52,14 +96,55 @@ function Base.show(io::IO, L::LinearOperator) ", shape=(", space_dim(L.codomain), ", ", space_dim(L.domain), "))") end -"""Operator matrix in weak (bilinear-form) representation.""" +""" + weak(L::LinearOperator) -> Matrix + weak(B::BilinearOperator) -> Array + +The operator in weak form: the bilinear form `⟨w, L v⟩` in the bases of the domain and +codomain. This is the form the operators are *built* in (it is what the weak-form +einsums assemble), and the form in which adjoints are transposes and composition needs +no Gram inverse. Derived from [`matrix`](@ref) via the codomain's [`gram`](@ref) if only +the strong form is known. + +# Examples +```jldoctest +julia> L = laplacian(dg, 0); + +julia> size(weak(L)) +(8, 8) + +julia> weak(L) ≈ weak(L)' # Δ is self-adjoint, so its weak form is symmetric +true + +julia> weak(L) ≈ gram(function_space(dg)) * matrix(L) +true +``` +""" function weak(L::LinearOperator) L._weak !== nothing && return L._weak L._weak = gram(L.codomain) * L._strong return L._weak end -"""Operator matrix in strong (coefficient-map) representation.""" +""" + matrix(L::LinearOperator) -> Matrix + +The operator in strong form: the matrix that maps domain coefficients to codomain +coefficients, which is what applying the operator actually multiplies by. Derived from +[`weak`](@ref) through the codomain's [`gram_inv`](@ref), so it inherits the `rcond` +truncation. + +# Examples +```jldoctest +julia> L = laplacian(dg, 0); + +julia> matrix(L) * f.coeffs ≈ L(f).coeffs +true + +julia> size(matrix(grad(dg))) # 16 vector-field coefficients ← 8 function ones +(16, 8) +``` +""" function matrix(L::LinearOperator) L._strong !== nothing && return L._strong L._strong = gram_inv(L.codomain) * L._weak @@ -68,12 +153,52 @@ end Base.size(L::LinearOperator) = size(matrix(L)) -"""Adjoint operator `L* : W → V` (conjugate transpose of the weak matrix).""" +""" + adjoint(L::LinearOperator) -> LinearOperator + L' + +The adjoint `L* : W → V`, defined by `⟨L v, w⟩_W = ⟨v, L* w⟩_V` — the conjugate +transpose of the *weak* matrix. This is how [`codifferential`](@ref) and +[`divergence`](@ref) are defined (`δ = d*`, `div = -∇*`). + +# Examples +The defining identity, checked on the circle: + +```jldoctest +julia> ω = d(dg_function(dg, sin.(θ))); + +julia> lhs = inner(dg, d(f), ω); # ⟨df, ω⟩ + +julia> rhs = inner(dg, f, codifferential(ω)); # ⟨f, δω⟩ + +julia> isapprox(lhs, rhs; rtol=1e-8) +true + +julia> grad(dg)' +LinearOperator(domain=VectorFieldSpace(dim=16), codomain=FunctionSpace(dim=8), shape=(8, 16)) +``` +""" function Base.adjoint(L::LinearOperator) return LinearOperator(L.codomain, L.domain; weak_matrix=Matrix(weak(L)')) end -"""`true` if the (endomorphism) operator equals its adjoint; `nothing` otherwise.""" +""" + is_self_adjoint(L) -> Bool or nothing + +`true` if the operator is an endomorphism equal to its own adjoint, `false` if it is an +endomorphism that is not, and `nothing` if it maps between different spaces (where the +question is meaningless). [`spectrum`](@ref) and [`inverse`](@ref) take the Hermitian +path when this is `true`. + +# Examples +```jldoctest +julia> is_self_adjoint(laplacian(dg, 0)) +true + +julia> is_self_adjoint(grad(dg)) === nothing # A → 𝔛(M): not an endomorphism +true +``` +""" function is_self_adjoint(L::LinearOperator) L.domain == L.codomain || return nothing W = weak(L) @@ -95,7 +220,25 @@ Base.:*(L::LinearOperator, s::Number) = LinearOperator(L.domain, L.codomain; wea Base.:*(s::Number, L::LinearOperator) = L * s Base.:-(L::LinearOperator) = (-1) * L -"""Compose two operators: `(A ∘ B)(x) = A(B(x))`, requiring `B.codomain == A.domain`.""" +""" + ∘(A::LinearOperator, B::LinearOperator) -> LinearOperator + A * B + +Compose two operators: `(A ∘ B)(x) = A(B(x))`. Requires `B.codomain == A.domain`. +Composed in weak form, so no Gram matrix is inverted along the way. + +# Examples +```jldoctest +julia> Δ = codifferential(dg, 1) ∘ d(dg, 0) # δd, the up-Laplacian by hand +LinearOperator(domain=FunctionSpace(dim=8), codomain=FunctionSpace(dim=8), shape=(8, 8)) + +julia> (identity_operator(function_space(dg)) ∘ laplacian(dg, 0))(f).coeffs ≈ laplacian(f).coeffs +true + +julia> grad(dg) ∘ grad(dg) # 𝔛(M) ≠ A +ERROR: AssertionError: Incompatible spaces: VectorFieldSpace → FunctionSpace. +``` +""" function Base.:∘(A::LinearOperator, B::LinearOperator) @assert B.codomain == A.domain "Incompatible spaces: $(typeof(B.codomain)) → $(typeof(A.domain))." # Compose in weak form to avoid inverting Gram matrices unnecessarily. @@ -138,7 +281,32 @@ end Eigenvalues (and, unless `eigvals_only`, the eigenvectors wrapped as a batched tensor of the domain) of an endomorphism, computed in the truncated orthonormal -basis of the domain. +basis of the domain. Eigenvalues come back ascending; the eigenvectors are stacked +along a *batch* axis, so they can be filtered as one batched tensor. + +# Examples +The Laplacian eigenvalues of the unit circle are `0, 1, 1, 4, 4, 9, 9, …`. The first +few are recovered well; the top of the spectrum is where the 8-mode truncation shows: + +```jldoctest +julia> evals = spectrum(laplacian(dg, 0); eigvals_only=true); + +julia> round.(abs.(evals[1:5]); digits=2) # |·| only to pin the sign of the zero +5-element Vector{Float64}: + 0.0 + 1.0 + 1.0 + 3.78 + 3.78 + +julia> evals, evecs = spectrum(laplacian(dg, 0)); + +julia> evecs # one eigenfunction per batch slot +ScalarFunction(space=FunctionSpace(dim=8), shape=(8, 8), batch_shape=(8,)) + +julia> round(l2_norm(evecs)[1]; digits=6) # L²-normalised +1.0 +``` """ function spectrum(L::LinearOperator; eigvals_only::Bool=false) evals, evecs = _spectral_decomposition(L) @@ -159,7 +327,27 @@ end inverse(L; rcond=nothing) -> LinearOperator Spectral (Moore–Penrose) pseudo-inverse of an endomorphism, discarding eigenvalues -with `|λ| ≤ rcond` (default `L.domain.dg.rcond`). +with `|λ| ≤ rcond` (default `L.domain.dg.rcond`). The inverse of a Laplacian is a +Green's function — [`hodge_decomposition`](@ref) is built on exactly this. + +Because the kernel is discarded, `inverse(L) ∘ L` is the projection *off* the kernel, +not the identity. On the circle the kernel of Δ is the constants, so it round-trips any +function with zero mean: + +# Examples +```jldoctest +julia> L = laplacian(dg, 0); + +julia> recovered = inverse(L)(L(f)); # f = cos θ has zero mean + +julia> isapprox(recovered.coeffs, f.coeffs; atol=1e-5) +true + +julia> const_fn = dg_function(dg, ones(60)); # in the kernel of Δ + +julia> maximum(abs.(inverse(L)(L(const_fn)).coeffs)) < 1e-10 +true +``` """ function inverse(L::LinearOperator; rcond::Union{Nothing,Real}=nothing) evals, evecs = _spectral_decomposition(L) @@ -189,7 +377,24 @@ function inverse(L::LinearOperator; rcond::Union{Nothing,Real}=nothing) end # ── Application ──────────────────────────────────────────────────────────────── -"""Apply the operator to a tensor in its domain, returning a tensor in the codomain.""" +""" + (L::LinearOperator)(t) -> AbstractTensor + +Apply the operator to a tensor of its domain, returning one in its codomain. Batch axes +pass straight through. + +# Examples +```jldoctest +julia> laplacian(dg, 0)(f) +ScalarFunction(space=FunctionSpace(dim=8), shape=(8,), batch_shape=()) + +julia> grad(dg)(f) +VectorField(space=VectorFieldSpace(dim=16), shape=(16,), batch_shape=()) + +julia> laplacian(dg, 0)(grad(f)) # wrong space +ERROR: AssertionError: Input tensor belongs to VectorFieldSpace, expected FunctionSpace. +``` +""" function (L::LinearOperator)(t::AbstractTensor) @assert t.space == L.domain "Input tensor belongs to $(typeof(t.space)), expected $(typeof(L.domain))." coeffs_flat, bshape = flatten_batch_dims(t.coeffs) @@ -199,13 +404,37 @@ function (L::LinearOperator)(t::AbstractTensor) end # ── Special operators ────────────────────────────────────────────────────────── -"""Zero operator `domain → codomain` (defaults to an endomorphism).""" +""" + zero_operator(domain, codomain=domain) -> LinearOperator + +The zero operator. Useful as a block in [`block`](@ref), and as the top-degree +up-Laplacian. + +# Examples +```jldoctest +julia> Z = zero_operator(function_space(dg), vector_field_space(dg)) +LinearOperator(domain=FunctionSpace(dim=8), codomain=VectorFieldSpace(dim=16), shape=(16, 8)) + +julia> all(iszero, Z(f).coeffs) +true +``` +""" function zero_operator(domain::AbstractTensorSpace, codomain::AbstractTensorSpace=domain) weak_matrix = zeros(Float64, space_dim(codomain), space_dim(domain)) return LinearOperator(domain, codomain; weak_matrix=weak_matrix) end -"""Identity operator on `space`.""" +""" + identity_operator(space) -> LinearOperator + +The identity on `space`. + +# Examples +```jldoctest +julia> identity_operator(function_space(dg))(f).coeffs ≈ f.coeffs +true +``` +""" function identity_operator(space::AbstractTensorSpace) strong = Matrix{Float64}(I, space_dim(space), space_dim(space)) return LinearOperator(space, space; strong_matrix=strong) diff --git a/src/tensors/base_tensor/abstract.jl b/src/tensors/base_tensor/abstract.jl index e6f32fc..1e839de 100644 --- a/src/tensors/base_tensor/abstract.jl +++ b/src/tensors/base_tensor/abstract.jl @@ -6,21 +6,210 @@ # declared up front so `DiffusionGeometry` can hold spaces as fields and so the # generic arithmetic / metric methods can dispatch on the abstractions. +""" + AbstractTensorSpace + +Supertype of the coefficient spaces: [`FunctionSpace`](@ref), [`VectorFieldSpace`](@ref), +[`FormSpace`](@ref), [`Tensor02Space`](@ref), [`Tensor02SymSpace`](@ref) and +[`DirectSumSpace`](@ref). + +A space is a *representation choice*, not a container: it says how many components a +point carries (`component_dim`), what their pointwise inner product is +(`cdc_components`), and hence what the metric and the [`gram`](@ref) matrix are. Each +is cached on its `DiffusionGeometry`, so spaces of the same kind are `===` and can be +compared by identity. + +# Examples +```jldoctest +julia> function_space(dg), form_space(dg, 1) +(FunctionSpace(dim=8), FormSpace(degree=1, dim=16)) + +julia> function_space(dg) isa AbstractTensorSpace +true + +julia> function_space(dg) + vector_field_space(dg) # `+` builds a direct sum +DirectSumSpace(spaces=[FunctionSpace, VectorFieldSpace], dim=24) +``` +""" abstract type AbstractTensorSpace end + +""" + AbstractTensor + +Supertype of the tensor fields: [`ScalarFunction`](@ref), [`VectorField`](@ref), +[`Form`](@ref), [`Tensor02`](@ref), [`Tensor02Sym`](@ref) and +[`DirectSumElement`](@ref). + +A tensor is a coefficient array plus the space it belongs to. The *last* axis of the +array is the coefficient axis; any leading axes are batch dimensions that broadcast +(numpy-style, right-aligned) through every operation. Arithmetic (`+`, `-`, scaling, +pointwise products with a function) is shared by all of them; see +[`to_pointwise_basis`](@ref) to get back the values at the sample points. + +# Examples +```jldoctest +julia> f +ScalarFunction(space=FunctionSpace(dim=8), shape=(8,), batch_shape=()) + +julia> f isa AbstractTensor +true + +julia> 2f - f == f # tensors are compared by identity, not value +false + +julia> (2f - f).coeffs ≈ f.coeffs +true +``` +""" abstract type AbstractTensor end -"Coefficient array of a tensor (trailing axis = coefficient axis, leading = batch)." +""" + degree(t) -> Int + +The differential-form degree of a tensor: `0` for a [`ScalarFunction`](@ref), `1` for a +[`VectorField`](@ref), and `k` for a degree-`k` [`Form`](@ref). + +# Examples +```jldoctest +julia> degree(f), degree(grad(f)), degree(d(f)), degree(d(d(f))) +(0, 1, 1, 2) +``` +""" +function degree end + +""" + cdc_components(space) -> Array + +The pointwise inner-product matrices of the space's component frame, shape `(n, C, C)` +— the carré du champ evaluated on whatever `e_a` the space is built from. It is the one +thing a concrete space has to supply: the metric, the [`gram`](@ref) matrix, and +[`component_dim`](@ref) are all derived from it. + +| space | components | +|:------|:-----------| +| [`FunctionSpace`](@ref) | `1` (the constant) | +| [`VectorFieldSpace`](@ref) | [`gamma_coords`](@ref), `Γ(xᵢ, xⱼ)` | +| [`FormSpace`](@ref) | the compound determinants, [`gamma_compound`](@ref) | +| [`Tensor02Space`](@ref) | [`gamma_02`](@ref) | +| [`Tensor02SymSpace`](@ref) | [`gamma_02_sym`](@ref) | + +# Examples +```jldoctest +julia> size(cdc_components(vector_field_space(dg))) # (n, d, d) +(60, 2, 2) + +julia> size(cdc_components(tensor02_space(dg))) # (n, d², d²) +(60, 4, 4) + +julia> cdc_components(vector_field_space(dg)) === gamma_coords(dg.cache) +true +``` +""" +function cdc_components end + +""" + wrap(space, coeffs) -> AbstractTensor + +Wrap a raw coefficient array as an element of `space`, checking its trailing dimension +(see [`infer_batch_shape`](@ref)). The low-level constructor — it takes coefficients, +not values, so use [`from_pointwise`](@ref) or the `dg_*` factories for data. + +# Examples +```jldoctest +julia> wrap(function_space(dg), zeros(8)) +ScalarFunction(space=FunctionSpace(dim=8), shape=(8,), batch_shape=()) + +julia> wrap(vector_field_space(dg), zeros(4, 16)) # a batch of four +VectorField(space=VectorFieldSpace(dim=16), shape=(4, 16), batch_shape=(4,)) + +julia> wrap(vector_field_space(dg), zeros(8)) # wrong length +ERROR: AssertionError: VectorField coefficients must have trailing shape (16,), got (8,) +``` +""" +function wrap end + +""" + coeffs(t) -> Array + +Coefficient array of a tensor: trailing axis is the coefficient axis, leading axes are +batch. Same as `t.coeffs`. + +# Examples +```jldoctest +julia> size(coeffs(f)) # 8 basis functions × 1 component +(8,) + +julia> size(coeffs(grad(f))) # 8 basis functions × 2 components +(16,) +``` +""" coeffs(t::AbstractTensor) = t.coeffs -"Tensor space an element lives in." +""" + space(t) -> AbstractTensorSpace + +The tensor space an element lives in. Same as `t.space`. + +# Examples +```jldoctest +julia> space(f) +FunctionSpace(dim=8) + +julia> space(grad(f)) === vector_field_space(dg) +true +``` +""" space(t::AbstractTensor) = t.space -"Parent `DiffusionGeometry` of a tensor / space." +""" + geometry(t) -> DiffusionGeometry + geometry(space) -> DiffusionGeometry + +The `DiffusionGeometry` a tensor or space belongs to. Two tensors can only be combined +if this is the *same object* for both. + +# Examples +```jldoctest +julia> geometry(f) === dg +true + +julia> geometry(function_space(dg)) === dg +true +``` +""" geometry(t::AbstractTensor) = t.space.dg geometry(s::AbstractTensorSpace) = s.dg -"Batch (broadcast) dimensions preceding the coefficient axis." +""" + batch_shape(t) -> Tuple + +The batch (broadcast) dimensions preceding the coefficient axis — `()` for a single +tensor. + +# Examples +```jldoctest +julia> batch_shape(f) +() + +julia> batch_shape(dg_function(dg, ones(3, 60))) # three functions at once +(3,) +``` +""" batch_shape(t::AbstractTensor) = size(t.coeffs)[1:end-1] -"Differential form degree of a space (`nothing` unless it is a `FormSpace`)." +""" + space_degree(space) -> Int or nothing + +The differential-form degree of a space: `0` for a [`FunctionSpace`](@ref), `k` for a +degree-`k` [`FormSpace`](@ref), and `nothing` for the spaces that are not forms. + +# Examples +```jldoctest +julia> space_degree(form_space(dg, 2)), space_degree(function_space(dg)) +(2, 0) + +julia> space_degree(vector_field_space(dg)) === nothing +true +``` +""" space_degree(::AbstractTensorSpace) = nothing diff --git a/src/tensors/base_tensor/base_tensor.jl b/src/tensors/base_tensor/base_tensor.jl index 7ec53c2..020fca2 100644 --- a/src/tensors/base_tensor/base_tensor.jl +++ b/src/tensors/base_tensor/base_tensor.jl @@ -5,7 +5,27 @@ # methods can dispatch on `ScalarFunction`. Julia's multiple dispatch replaces the # Python `__array_ufunc__` / `isinstance` machinery entirely. -"""Pointwise data of a tensor at the sample points, shape `(batch..., n·C)`.""" +""" + to_pointwise_basis(t) -> Array + +The tensor's values at the sample points, shape `(batch..., n·C)` (`C = +component_dim`) — the inverse of the `dg_*` factories, up to the basis truncation. +This is how you get numbers back out of the coefficient world. + +# Examples +```jldoctest +julia> values = to_pointwise_basis(f); # f was built from cos θ + +julia> size(values) +(60,) + +julia> maximum(abs.(values .- cos.(θ))) < 1e-4 # the 8-mode projection error +true + +julia> size(to_pointwise_basis(grad(f))) # 60 points × 2 components +(120,) +``` +""" to_pointwise_basis(t::AbstractTensor) = _to_pointwise_basis(t.coeffs, t.space) # ── Display ──────────────────────────────────────────────────────────────────── @@ -16,6 +36,35 @@ Base.show(io::IO, t::AbstractTensor) = ", batch_shape=", batch_shape(t), ")") # ── Linear structure ─────────────────────────────────────────────────────────── +""" + +(a::AbstractTensor, b::AbstractTensor) -> AbstractTensor + -(a::AbstractTensor, b::AbstractTensor) -> AbstractTensor + *(s::Number, t::AbstractTensor) -> AbstractTensor + /(t::AbstractTensor, s::Number) -> AbstractTensor + +Tensors of the same space form a vector space: addition and scalar multiplication act +on the coefficients, and batch shapes broadcast (numpy right-alignment). Adding +tensors of *different* spaces is an error — use [`pack`](@ref) if you want them +side by side. + +A `Number` added to a `ScalarFunction` is the constant function of that value; no +other tensor absorbs a scalar that way. + +# Examples +```jldoctest +julia> (f + f).coeffs ≈ (2f).coeffs +true + +julia> (f - f).coeffs ≈ zeros(function_space(dg)).coeffs +true + +julia> round((f + 1).coeffs[1]; digits=4) # φ₀ ≡ 1, so the constant lands there +1.0 + +julia> f + grad(f) +ERROR: AssertionError: Operands must belong to the same space: FunctionSpace(dim=8) vs VectorFieldSpace(dim=16) +``` +""" Base.:-(t::AbstractTensor) = wrap(t.space, -t.coeffs) function Base.:+(a::AbstractTensor, b::AbstractTensor) @@ -42,6 +91,30 @@ function _scale_by_batch(t::AbstractTensor, w::AbstractArray{<:Number}, op) return wrap(t.space, op.(ct, wt)) end +""" + *(w::AbstractArray, t::AbstractTensor) -> AbstractTensor + /(t::AbstractTensor, w::AbstractArray) -> AbstractTensor + +A plain numeric array is a bag of *batch-wise scalars*, never a coefficient vector: +its shape broadcasts against the tensor's batch axes and the coefficient axis is left +alone. This is what makes a spectral filter work on the batched eigenbasis that +[`spectrum`](@ref) returns. + +# Examples +Heat-kernel damping of each eigenmode by `exp(-λ)`: + +```jldoctest +julia> evals, evecs = spectrum(laplacian(dg, 0)); + +julia> batch_shape(evecs) # one eigenvector per batch slot +(8,) + +julia> damped = exp.(-evals) .* evecs; # `.*` is a synonym for `*` here + +julia> batch_shape(damped) +(8,) +``` +""" Base.:*(w::AbstractArray{<:Number}, t::AbstractTensor) = _scale_by_batch(t, w, *) Base.:*(t::AbstractTensor, w::AbstractArray{<:Number}) = _scale_by_batch(t, w, *) Base.:/(t::AbstractTensor, w::AbstractArray{<:Number}) = _scale_by_batch(t, w, /) @@ -124,6 +197,32 @@ function _pointwise_divide(t::AbstractTensor, f) return wrap(t.space, _from_pointwise_basis(quotient, t.space)) end +""" + *(t::AbstractTensor, f::ScalarFunction) -> AbstractTensor + /(t::AbstractTensor, f::ScalarFunction) -> AbstractTensor + +Multiplying (or dividing) any tensor by a `ScalarFunction` is the *pointwise* product: +the values are multiplied at each sample point and the result is projected back onto +the basis. Two functions multiply pointwise as well. + +Because it round-trips through the sample points, `f * g` is not exactly bilinear in +the coefficients — it is the projection of the true product, and modes above the +truncation are lost. + +# Examples +```jldoctest +julia> squared = f * f; # cos²θ + +julia> maximum(abs.(to_pointwise_basis(squared) .- cos.(θ) .^ 2)) < 1e-3 +true + +julia> f * grad(f) # scaling a vector field by a function +VectorField(space=VectorFieldSpace(dim=16), shape=(16,), batch_shape=()) +``` + +Division floors the denominator at `1e-12`, so dividing by a function with zeros is +survivable but meaningless there. +""" # Function × Function is a pointwise product (most specific dispatch). Base.:*(f::ScalarFunction, g::ScalarFunction) = _pointwise_product(f, g) Base.:/(f::ScalarFunction, g::ScalarFunction) = _pointwise_divide(f, g) @@ -133,7 +232,20 @@ Base.:*(t::AbstractTensor, f::ScalarFunction) = _pointwise_product(t, f) Base.:*(f::ScalarFunction, t::AbstractTensor) = _pointwise_product(t, f) Base.:/(t::AbstractTensor, f::ScalarFunction) = _pointwise_divide(t, f) -# ── Pointwise exponentiation of a function ───────────────────────────────────── +""" + ^(f::ScalarFunction, p::Number) -> ScalarFunction + +Pointwise power: raise the function's *values* to `p` and project back onto the basis. + +# Examples +```jldoctest +julia> maximum(abs.(to_pointwise_basis(f^2) .- to_pointwise_basis(f * f))) < 1e-12 +true + +julia> to_pointwise_basis(f^0) ≈ ones(60) # every value to the zeroth power +true +``` +""" function Base.:^(f::ScalarFunction, p::Number) data = to_pointwise_basis(f) return dg_function(geometry(f), data .^ p) diff --git a/src/tensors/base_tensor/base_tensor_space.jl b/src/tensors/base_tensor/base_tensor_space.jl index fd4ef9d..e13a753 100644 --- a/src/tensors/base_tensor/base_tensor_space.jl +++ b/src/tensors/base_tensor/base_tensor_space.jl @@ -11,33 +11,153 @@ using LinearAlgebra: Symmetric, Diagonal, eigen # Coefficient functions truncated to the tensor rank n_coefficients. _u_coeffs(dg::DiffusionGeometry) = function_basis(dg)[:, 1:n_coefficients(dg)] -"""Number of pointwise components represented per sample point.""" +""" + component_dim(space) -> Int + +Number of pointwise components an element of `space` carries per sample point: 1 for a +function, `d` for a vector field or 1-form, `binomial(d, k)` for a k-form, `d²` for a +(0,2)-tensor, `d(d+1)/2` for a symmetric one. + +# Examples +```jldoctest +julia> component_dim(function_space(dg)), component_dim(vector_field_space(dg)) +(1, 2) + +julia> component_dim(tensor02_space(dg)), component_dim(tensor02sym_space(dg)) +(4, 3) +``` +""" component_dim(space::AbstractTensorSpace) = size(cdc_components(space), 2) -"""Number of coefficients needed to represent an element of `space`.""" +""" + space_dim(space) -> Int + +Length of the coefficient vector of an element of `space`: `n_coefficients × +component_dim` (functions use the full `n_function_basis` instead). This is the +number a [`wrap`](@ref) checks against, and the matrix size of every operator into or +out of the space. + +# Examples +```jldoctest +julia> space_dim(function_space(dg)) # 8 basis functions, 1 component +8 + +julia> space_dim(vector_field_space(dg)) # 8 × 2 +16 + +julia> space_dim(tensor02_space(dg)) # 8 × 2² +32 +``` +""" space_dim(space::AbstractTensorSpace) = n_coefficients(space.dg) * component_dim(space) # ── Riemannian metric ────────────────────────────────────────────────────────── -"""Metric tensor field `g(e_a, e_b)`, shape `(n, dim, dim)`.""" +""" + metric_tensor(space) -> Array + +The metric of the space written out in its own basis `{φ_i e_a}`: `g(φ_i e_a, φ_I e_b)` +at every point, shape `(n, space_dim, space_dim)`. Rarely needed directly — [`g`](@ref) +and [`metric_apply`](@ref) contract it against actual coefficients instead. + +# Examples +```jldoctest +julia> size(metric_tensor(vector_field_space(dg))) +(60, 16, 16) +``` +""" metric_tensor(space::AbstractTensorSpace) = metric(_u_coeffs(space.dg), cdc_components(space)) -"""Pointwise metric of two coefficient vectors, shape `(batch..., n)`.""" +""" + metric_apply(space, a_coeffs, b_coeffs) -> Array + +Pointwise metric `g(A, B)` of two *coefficient vectors* of `space`, shape +`(batch..., n)`. The implementation behind [`g`](@ref), which is the one to call — +this takes raw coefficients and skips the space check. + +# Examples +```jldoctest +julia> X = grad(f); + +julia> vals = metric_apply(vector_field_space(dg), X.coeffs, X.coeffs); + +julia> vals ≈ g(dg, X, X) +true +``` +""" metric_apply(space::AbstractTensorSpace, a_coeffs, b_coeffs) = _metric_apply(_u_coeffs(space.dg), regularise_fn(space.dg), a_coeffs, b_coeffs, cdc_components(space)) # ── Gram matrix and its spectral pseudo-inverse ──────────────────────────────── -"""Gram matrix `G_{ia,Ib} = ∫ φ_i φ_I g(e_a, e_b) dμ`, shape `(dim, dim)`.""" +""" + gram(space) -> Matrix + +Gram matrix `G_{ia,Ib} = ∫ φ_i φ_I g(e_a, e_b) dμ` of the space's basis, shape +`(space_dim, space_dim)`. It *is* the L² inner product: `⟨A, B⟩ = aᵀ G b`, which is why +it converts between an operator's weak and strong forms (see [`weak`](@ref) and +[`matrix`](@ref)). + +# Examples +The eigenfunction basis is L²-orthonormal, so the function space's Gram matrix is the +identity — but the vector-field basis `{φ_i ∇x_j}` is not orthonormal, and its Gram +matrix is where the metric of the data shows up: + +```jldoctest +julia> gram(function_space(dg)) ≈ I +true + +julia> gram(vector_field_space(dg)) ≈ I +false + +julia> size(gram(vector_field_space(dg))) +(16, 16) +``` +""" gram(space::AbstractTensorSpace) = gram(_u_coeffs(space.dg), cdc_components(space), measure(space.dg)) -"""`true` when the Gram matrix is (numerically) the identity (numpy allclose semantics).""" +""" + _is_orthonormal(space) -> Bool + +`true` when the space's [`gram`](@ref) matrix is numerically the identity (numpy +`allclose` semantics), in which case the basis conversions can skip inverting it. + +# Examples +```jldoctest +julia> DiffusionGeometryJ._is_orthonormal(function_space(dg)) +true + +julia> DiffusionGeometryJ._is_orthonormal(vector_field_space(dg)) +false +``` +""" function _is_orthonormal(space::AbstractTensorSpace) G = gram(space) E = Matrix{eltype(G)}(I, size(G)...) return all(abs.(G .- E) .<= 1e-10 .+ 1e-5 .* abs.(E)) end -"""Eigen-pairs of the Gram matrix retaining eigenvalues above `rcond`.""" +""" + _gram_spectrum(space) -> (values, vectors) + +Eigen-pairs of the [`gram`](@ref) matrix, keeping only eigenvalues above `dg.rcond`. +The truncation is what makes the geometry numerically stable: the discarded directions +are the ones the sampled data does not actually resolve. + +# Examples +Nothing falls below the default `rcond = 1e-5` here, so all 16 directions survive — +but the weakest are two orders of magnitude down on the rest. Those are the ambient +directions *normal* to the circle, which the metric barely sees: + +```jldoctest +julia> vals, vecs = DiffusionGeometryJ._gram_spectrum(vector_field_space(dg)); + +julia> length(vals) +16 + +julia> round(minimum(vals); sigdigits=2) +0.018 +``` +""" function _gram_spectrum(space::AbstractTensorSpace) G = gram(space) F = eigen(Symmetric(Matrix(G))) @@ -45,7 +165,29 @@ function _gram_spectrum(space::AbstractTensorSpace) return F.values[keep], F.vectors[:, keep] end -"""Moore–Penrose pseudo-inverse of the Gram matrix with spectral cutoff.""" +""" + gram_inv(space) -> Matrix + +Moore–Penrose pseudo-inverse of the [`gram`](@ref) matrix, with the `rcond` spectral +cutoff. Converts an operator's weak form back to its strong form. + +# Examples +`G⁺ G` is the projection onto the resolved subspace. Nothing is truncated in this +geometry, so here it is the identity — but it is a projection in general, and that is +why an operator round-tripped through its weak form loses whatever lived below `rcond`: + +```jldoctest +julia> G = gram(vector_field_space(dg)); Gi = gram_inv(vector_field_space(dg)); + +julia> P = Gi * G; + +julia> P * P ≈ P +true + +julia> P ≈ I # full rank here: nothing fell below rcond +true +``` +""" function gram_inv(space::AbstractTensorSpace) evals, evecs = _gram_spectrum(space) if length(evals) == 0 @@ -55,7 +197,24 @@ function gram_inv(space::AbstractTensorSpace) return evecs * Diagonal(1.0 ./ evals) * transpose(evecs) end -"""Orthonormal basis of the spectrally truncated subspace.""" +""" + orthonormal_basis(space) -> Matrix + +An L²-orthonormal basis of the subspace that survives the `rcond` cutoff, as columns of +a `(space_dim, k)` matrix. [`spectrum`](@ref) and [`inverse`](@ref) work in it, which +is how they avoid the unresolved directions of the space. + +# Examples +```jldoctest +julia> Φ = orthonormal_basis(vector_field_space(dg)); + +julia> size(Φ) # 16 coefficients, full rank here +(16, 16) + +julia> Φ' * gram(vector_field_space(dg)) * Φ ≈ I # orthonormal in the L² metric +true +``` +""" function orthonormal_basis(space::AbstractTensorSpace) evals, evecs = _gram_spectrum(space) length(evals) == 0 && return zeros(Float64, size(gram(space), 1), 0) @@ -63,11 +222,36 @@ function orthonormal_basis(space::AbstractTensorSpace) end # ── Elements ─────────────────────────────────────────────────────────────────── -"""Build an element of `space` from pointwise data (trailing shape `(n, C)`).""" +""" + from_pointwise(space, data) -> AbstractTensor + +Build an element of `space` from pointwise `data` (trailing shape `(n, C)`, with `C = +component_dim(space)`), by L²-projection onto the basis. The `dg_*` factories +([`dg_function`](@ref), [`dg_vector_field`](@ref), …) are the friendlier front doors. + +# Examples +```jldoctest +julia> from_pointwise(vector_field_space(dg), [-sin.(θ) cos.(θ)]) +VectorField(space=VectorFieldSpace(dim=16), shape=(16,), batch_shape=()) +``` +""" from_pointwise(space::AbstractTensorSpace, data::AbstractArray) = wrap(space, _from_pointwise_basis(data, space)) -"""Zero element with an optional batch shape.""" +""" + zeros(space, batch=()) -> AbstractTensor + +The zero element of `space`, optionally batched. + +# Examples +```jldoctest +julia> zeros(function_space(dg)) +ScalarFunction(space=FunctionSpace(dim=8), shape=(8,), batch_shape=()) + +julia> zeros(vector_field_space(dg), (3,)) +VectorField(space=VectorFieldSpace(dim=16), shape=(3, 16), batch_shape=(3,)) +``` +""" function Base.zeros(space::AbstractTensorSpace, batch::Tuple=()) return wrap(space, zeros(Float64, batch..., space_dim(space))) end diff --git a/src/tensors/base_tensor/metric_gram.jl b/src/tensors/base_tensor/metric_gram.jl index 3acdd9d..7911469 100644 --- a/src/tensors/base_tensor/metric_gram.jl +++ b/src/tensors/base_tensor/metric_gram.jl @@ -8,7 +8,17 @@ using OMEinsum: @ein_str, @optein_str Metric tensor field `g(e_a, e_b)`, shape `(n, n1 C, n1 C)`. `matrices` `(n, C, C)` are the local inner-product matrices (compound determinants -for k-forms, `gamma_02`/`gamma_02_sym` for (0,2)-tensors). `u_n1` is `(n, n1)`. +for k-forms, [`gamma_02`](@ref)/[`gamma_02_sym`](@ref) for (0,2)-tensors). `u_n1` is +`(n, n1)`. Called through [`metric_tensor`](@ref), which supplies the right matrices for +each space. + +# Examples +```jldoctest +julia> u = function_basis(dg)[:, 1:2]; + +julia> size(metric(u, gamma_coords(dg.cache))) # 2 basis functions × 2 components +(60, 4, 4) +``` """ function metric(u_n1::AbstractMatrix, matrices::AbstractArray{<:Any,3}) n, n1 = size(u_n1) @@ -35,7 +45,13 @@ end _metric_apply(u_n1, regularise_func, a_coeffs, b_coeffs, matrices) -> Array Pointwise metric `g(A, B)(p) = A^c(p) g_cd(p) B^d(p)` for two coefficient vectors, -returning `(batch..., n)`. `u_n1` is `(n, n1)`; `matrices` `(n, C, C)`. +returning `(batch..., n)`. `u_n1` is `(n, n1)`; `matrices` `(n, C, C)`. The kernel +behind [`metric_apply`](@ref) and hence [`g`](@ref) — the result is regularised before +it is returned. + +Takes the space apart into raw arrays, so it is called with `u_n1` and `matrices` rather +than a space; see [`metric_apply`](@ref) for a worked example through the public entry +point. """ function _metric_apply(u_n1::AbstractMatrix, regularise_func, a_coeffs, b_coeffs, matrices::AbstractArray{<:Any,3}) @@ -62,7 +78,16 @@ end gram(u_n1, matrices, measure) -> Matrix Gram matrix of a tensor basis, shape `(n1 C, n1 C)`. -`G_{ia, Ib} = ⟨φ_i e_a, φ_I e_b⟩ = ∫ φ_i φ_I g(e_a, e_b) dμ`. +`G_{ia, Ib} = ⟨φ_i e_a, φ_I e_b⟩ = ∫ φ_i φ_I g(e_a, e_b) dμ`. Called through +[`gram`](@ref), which supplies the right matrices for each space. + +# Examples +```jldoctest +julia> u = function_basis(dg)[:, 1:2]; + +julia> size(gram(u, gamma_coords(dg.cache), measure(dg))) +(4, 4) +``` """ function gram(u_n1::AbstractMatrix, matrices::AbstractArray{<:Any,3}, measure::AbstractVector) n1 = size(u_n1, 2) diff --git a/src/tensors/direct_sum/direct_sum.jl b/src/tensors/direct_sum/direct_sum.jl index 670c1e1..4acccaa 100644 --- a/src/tensors/direct_sum/direct_sum.jl +++ b/src/tensors/direct_sum/direct_sum.jl @@ -5,7 +5,28 @@ # matrix and its inverse are block-diagonal, and the metric is the sum of the # summands' metrics. -"""Direct sum of tensor spaces sharing a `DiffusionGeometry`.""" +""" + DirectSumSpace(dg, spaces) + +Direct sum of tensor spaces sharing a `DiffusionGeometry`, built with `+`. Coefficients +of the summands are concatenated; the Gram matrix is block-diagonal, and the metric is +the sum of the summands' metrics. Nested sums are flattened. + +This is what lets a *system* of tensors be one vector — a mixed-degree operator built +with [`block`](@ref) has a direct sum as its domain. + +# Examples +```jldoctest +julia> ds = function_space(dg) + vector_field_space(dg) +DirectSumSpace(spaces=[FunctionSpace, VectorFieldSpace], dim=24) + +julia> space_dim(ds) == space_dim(function_space(dg)) + space_dim(vector_field_space(dg)) +true + +julia> (function_space(dg) + vector_field_space(dg)) + form_space(dg, 1) # flattened +DirectSumSpace(spaces=[FunctionSpace, VectorFieldSpace, FormSpace], dim=40) +``` +""" struct DirectSumSpace <: AbstractTensorSpace dg::DiffusionGeometry spaces::Tuple{Vararg{AbstractTensorSpace}} @@ -24,7 +45,25 @@ struct DirectSumSpace <: AbstractTensorSpace end end -"""An element of a direct sum, storing concatenated coefficients.""" +""" + DirectSumElement + +An element of a [`DirectSumSpace`](@ref), storing the summands' coefficients end to end. +Build one with [`pack`](@ref) and take it apart with [`unpack`](@ref). + +# Examples +```jldoctest +julia> ds = function_space(dg) + vector_field_space(dg); + +julia> e = pack(ds, f, grad(f)) +DirectSumElement(space=DirectSumSpace(spaces=[FunctionSpace, VectorFieldSpace], dim=24), shape=(24,), batch_shape=()) + +julia> a, X = unpack(e); + +julia> a.coeffs == f.coeffs, X.coeffs == grad(f).coeffs +(true, true) +``` +""" struct DirectSumElement{A<:AbstractArray} <: AbstractTensor space::DirectSumSpace coeffs::A @@ -54,14 +93,61 @@ function wrap(ds::DirectSumSpace, coeffs::AbstractArray) return DirectSumElement(ds, coeffs) end +""" + split_coeffs(ds, coeffs) -> Tuple + +Split a direct-sum coefficient array into one array per summand (a tuple, in the +summands' order). The coefficient-level counterpart of [`unpack`](@ref). + +# Examples +```jldoctest +julia> ds = function_space(dg) + vector_field_space(dg); + +julia> parts = split_coeffs(ds, pack(ds, f, grad(f)).coeffs); + +julia> size.(parts) +((8,), (16,)) +``` +""" split_coeffs(ds::DirectSumSpace, coeffs::AbstractArray) = Tuple(copy(selectdim(coeffs, ndims(coeffs), r)) for r in _coeff_ranges(ds)) -"""Split a direct-sum element into its wrapped component tensors.""" +""" + unpack(e::DirectSumElement) -> Tuple + +Split a direct-sum element back into its component tensors, each wrapped in its own +space. The inverse of [`pack`](@ref). + +# Examples +```jldoctest +julia> ds = function_space(dg) + vector_field_space(dg); + +julia> unpack(pack(ds, f, grad(f))) +(ScalarFunction(space=FunctionSpace(dim=8), shape=(8,), batch_shape=()), VectorField(space=VectorFieldSpace(dim=16), shape=(16,), batch_shape=())) +``` +""" unpack(e::DirectSumElement) = Tuple(wrap(s, c) for (s, c) in zip(e.space.spaces, split_coeffs(e.space, e.coeffs))) -"""Combine one tensor per summand into a `DirectSumElement`.""" +""" + pack(ds::DirectSumSpace, tensors...) -> DirectSumElement + +Combine one tensor per summand (in order, sharing a batch shape) into a single element +of the direct sum. Inverted by [`unpack`](@ref). + +# Examples +```jldoctest +julia> ds = function_space(dg) + vector_field_space(dg); + +julia> e = pack(ds, f, grad(f)); + +julia> length(e.coeffs) # 8 + 16 +24 + +julia> pack(ds, f) # one tensor short +ERROR: AssertionError: Expected 2 tensors, received 1 +``` +""" function pack(ds::DirectSumSpace, tensors...) @assert length(tensors) == length(ds.spaces) "Expected $(length(ds.spaces)) tensors, received $(length(tensors))" bshape = batch_shape(tensors[1]) diff --git a/src/tensors/forms/forms.jl b/src/tensors/forms/forms.jl index c0a7742..d844782 100644 --- a/src/tensors/forms/forms.jl +++ b/src/tensors/forms/forms.jl @@ -6,7 +6,25 @@ # and live in `operators/tensor_actions.jl`. The pure algebra — the wedge product, # the tensor product of 1-forms, and the musical `sharp` — lives here. -"""Space of differential k-forms Ωᵏ(M); `binomial(dim, k)` components per point.""" +""" + FormSpace(dg, degree) + +Space of differential `degree`-forms `Ωᵏ(M)`, with `binomial(ambient_dim, degree)` +components per point. Its metric is the compound determinant of Γ (see +[`gamma_compound`](@ref)). Get it with [`form_space`](@ref), which also handles degree 0. + +# Examples +```jldoctest +julia> form_space(dg, 1) +FormSpace(degree=1, dim=16) + +julia> component_dim(form_space(dg, 2)) # binomial(2, 2) = 1 +1 + +julia> form_space(dg, 3) # beyond the ambient dimension +ERROR: AssertionError: Form degree k=3 must be between 1 and 2 +``` +""" struct FormSpace <: AbstractTensorSpace dg::DiffusionGeometry degree::Int @@ -16,7 +34,32 @@ struct FormSpace <: AbstractTensorSpace end end -"""A differential k-form `ω ∈ Ωᵏ(M)`, coefficients in the basis {φ_i dx_J}.""" +""" + Form + +A differential k-form `ω ∈ Ωᵏ(M)`, stored as coefficients in the basis `{φ_i dx_J}` +(the `J` in lexicographic order — see [`get_wedge_basis_indices`](@ref)). Build one with +[`d`](@ref), [`flat`](@ref) or [`dg_form`](@ref). + +The exterior algebra is here: [`wedge`](@ref) (also spelled `^`), [`d`](@ref), +[`codifferential`](@ref), the Laplacians, and [`hodge_decomposition`](@ref). A 1-form +can be applied to a vector field, giving the pointwise values of `ω(X)`. + +# Examples +```jldoctest +julia> ω = d(f) +Form(space=FormSpace(degree=1, dim=16), shape=(16,), batch_shape=()) + +julia> degree(ω) +1 + +julia> size(ω(grad(f))) # ω(X) = g(ω♯, X), a value per point +(60,) + +julia> d(ω) # dω is a 2-form +Form(space=FormSpace(degree=2, dim=8), shape=(8,), batch_shape=()) +``` +""" struct Form{A<:AbstractArray} <: AbstractTensor space::FormSpace coeffs::A @@ -52,6 +95,34 @@ zero function when the combined degree exceeds the ambient dimension. `^` is the same product, spelled as in Python. Beware that `*` on two `Form`s is the *tensor* product, not the wedge. + +# Examples +The wedge is antisymmetric on 1-forms, so `α ∧ α = 0` and `α ∧ β = -(β ∧ α)`: + +```jldoctest +julia> α = d(f); β = d(dg_function(dg, sin.(θ))); + +julia> wedge(α, β) +Form(space=FormSpace(degree=2, dim=8), shape=(8,), batch_shape=()) + +julia> (α ^ β).coeffs ≈ -(β ^ α).coeffs +true + +julia> maximum(abs.((α ^ α).coeffs)) +0.0 +``` + +Past the top degree there is nowhere to go, and the result is the zero function: + +```jldoctest +julia> α = d(f); γ = d(α); # a 1-form and a 2-form, in a 2-dim frame + +julia> wedge(α, γ) # degree 3 > 2, so it vanishes +ScalarFunction(space=FunctionSpace(dim=8), shape=(8,), batch_shape=()) + +julia> all(iszero, wedge(α, γ).coeffs) +true +``` """ function wedge(a::Form, b::Form) dg = geometry(a) @@ -92,6 +163,14 @@ A function is a degree-0 form, so the wedge with one degenerates to the pointwis product (as in Python). There is deliberately no `wedge` on two `ScalarFunction`s: `f^g` already reads as a pointwise power, and `Base.:^(::ScalarFunction, ::Number)` is that power — writing `f * g` for the product keeps `^` unambiguous. + +# Examples +```jldoctest +julia> ω = d(f); + +julia> wedge(f, ω).coeffs ≈ (f * ω).coeffs +true +``` """ wedge(f::ScalarFunction, ω::Form) = f * ω wedge(ω::Form, f::ScalarFunction) = ω * f @@ -110,8 +189,19 @@ Base.:^(ω::Form, f::ScalarFunction) = wedge(ω, f) # ── Tensor product of 1-forms → general (0,2)-tensor ─────────────────────────── """ -Tensor product of two 1-forms, `(α ⊗ β)_{ij} = α_i β_j`, giving a `Tensor02`. -This is *not* the wedge product — that is `wedge(α, β)`, or `α ^ β`. +Tensor product of two 1-forms, `(α ⊗ β)_{ij} = α_i β_j`, giving a [`Tensor02`](@ref). +This is *not* the wedge product — that is [`wedge`](@ref), or `α ^ β`. + +# Examples +```jldoctest +julia> α = d(f); β = d(dg_function(dg, sin.(θ))); + +julia> T = α * β +Tensor02(space=Tensor02Space(dim=32), shape=(32,), batch_shape=()) + +julia> transpose_tensor(T).coeffs ≈ (β * α).coeffs # (α ⊗ β)ᵀ = β ⊗ α +true +``` """ function Base.:*(a::Form, b::Form) dg = geometry(a) diff --git a/src/tensors/functions/functions.jl b/src/tensors/functions/functions.jl index 83b2adb..1ff1b49 100644 --- a/src/tensors/functions/functions.jl +++ b/src/tensors/functions/functions.jl @@ -5,12 +5,53 @@ # expanded in the full coefficient basis {φ_i} (n_function_basis columns), so they # use their own `basis_count` and `component_dim = 1`. -"""Space of scalar functions A ≅ L²(M, μ).""" +""" + FunctionSpace(dg) + +Space of scalar functions `A ≅ L²(M, μ)`. Unlike the other spaces it expands in the +*full* eigenfunction basis (`n_function_basis` coefficients, not `n_coefficients`), and +its metric is just the pointwise product — there is nothing to contract. + +Get it with [`function_space`](@ref) rather than constructing one. + +# Examples +```jldoctest +julia> function_space(dg) +FunctionSpace(dim=8) + +julia> component_dim(function_space(dg)), space_dim(function_space(dg)) +(1, 8) + +julia> gram(function_space(dg)) ≈ I # the eigenfunctions are L²-orthonormal +true +``` +""" struct FunctionSpace <: AbstractTensorSpace dg::DiffusionGeometry end -"""A scalar function `f ∈ A`, coefficients in the basis {φ_i}.""" +""" + ScalarFunction + +A scalar function `f ∈ A`, stored as coefficients in the eigenfunction basis `{φ_i}`. +Build one with [`dg_function`](@ref); get its values back with +[`to_pointwise_basis`](@ref). + +Functions are the scalars of this algebra: multiplying any tensor by one is the +pointwise product, and `+`/`-` against a `Number` adds a constant function. + +# Examples +```jldoctest +julia> f +ScalarFunction(space=FunctionSpace(dim=8), shape=(8,), batch_shape=()) + +julia> degree(f) # a function is a 0-form +0 + +julia> grad(f), d(f), laplacian(f) +(VectorField(space=VectorFieldSpace(dim=16), shape=(16,), batch_shape=()), Form(space=FormSpace(degree=1, dim=16), shape=(16,), batch_shape=()), ScalarFunction(space=FunctionSpace(dim=8), shape=(8,), batch_shape=())) +``` +""" struct ScalarFunction{A<:AbstractArray} <: AbstractTensor space::FunctionSpace coeffs::A diff --git a/src/tensors/tensor02/tensor02.jl b/src/tensors/tensor02/tensor02.jl index a2980b4..e12e3f1 100644 --- a/src/tensors/tensor02/tensor02.jl +++ b/src/tensors/tensor02/tensor02.jl @@ -5,12 +5,51 @@ # need the operator layer and are deferred to Phase 5. Symmetrisation and # transpose are pure algebra and live here. -"""Space of general (0,2)-tensors Ω⁰²(M); `dim²` components per point.""" +""" + Tensor02Space(dg) + +Space of general (0,2)-tensors, with `ambient_dim²` components per point. Its metric is +the Kronecker square of Γ (see [`gamma_02`](@ref)). Get it with +[`tensor02_space`](@ref). + +# Examples +```jldoctest +julia> tensor02_space(dg) +Tensor02Space(dim=32) + +julia> component_dim(tensor02_space(dg)) # d² = 4 +4 +``` +""" struct Tensor02Space <: AbstractTensorSpace dg::DiffusionGeometry end -"""A (0,2)-tensor `T ∈ Ω⁰²(M)`, coefficients in the basis {φ_i dx_j ⊗ dx_k}.""" +""" + Tensor02 + +A (0,2)-tensor, stored as coefficients in the basis `{φ_i dx_j ⊗ dx_k}` (row-major in +`(j, k)`). Produced by [`levi_civita`](@ref), by the tensor product of two 1-forms +(`α * β`), or by [`dg_tensor02`](@ref). + +It has two lives. Called on one vector field it acts as an *operator* `𝔛(M) → 𝔛(M)` +(see [`t02_operator`](@ref)); called on two it is the *bilinear form* `T(X, Y)`, +returning pointwise values. + +# Examples +```jldoctest +julia> T = d(f) * d(dg_function(dg, sin.(θ))); + +julia> T(grad(f)) # one argument → a vector field +VectorField(space=VectorFieldSpace(dim=16), shape=(16,), batch_shape=()) + +julia> size(T(grad(f), grad(f))) # two → a value per point +(60,) + +julia> symmetrise(T) +Tensor02Sym(space=Tensor02SymSpace(dim=24), shape=(24,), batch_shape=()) +``` +""" struct Tensor02{A<:AbstractArray} <: AbstractTensor space::Tensor02Space coeffs::A @@ -25,14 +64,45 @@ function wrap(space::Tensor02Space, coeffs::AbstractArray) return Tensor02(space, coeffs) end -"""Symmetric part of the tensor as a `Tensor02Sym`.""" +""" + symmetrise(T::Tensor02) -> Tensor02Sym + +The symmetric part `½(T + Tᵀ)`, re-expressed in the `d(d+1)/2` symmetric components. +Inverted by [`full_tensor`](@ref), which expands it back. + +# Examples +```jldoctest +julia> T = d(f) * d(dg_function(dg, sin.(θ))); + +julia> S = symmetrise(T) +Tensor02Sym(space=Tensor02SymSpace(dim=24), shape=(24,), batch_shape=()) + +julia> full_tensor(S).coeffs ≈ (0.5 * (T + transpose_tensor(T))).coeffs +true +``` +""" function symmetrise(T::Tensor02) dg = geometry(T) coeffs_sym = symmetrise_tensor_coeffs(T.coeffs, n_coefficients(dg), ambient_dim(dg)) return wrap(tensor02sym_space(dg), coeffs_sym) end -"""Transpose `(Tᵀ)_{ij} = T_{ji}`.""" +""" + transpose_tensor(T) -> Tensor02 + +Transpose `(Tᵀ)_{ij} = T_{ji}`. A symmetric tensor is its own transpose. + +# Examples +```jldoctest +julia> α = d(f); β = d(dg_function(dg, sin.(θ))); + +julia> transpose_tensor(α * β).coeffs ≈ (β * α).coeffs +true + +julia> transpose_tensor(transpose_tensor(α * β)).coeffs ≈ (α * β).coeffs +true +``` +""" function transpose_tensor(T::Tensor02) dg = geometry(T) n1, d = n_coefficients(dg), ambient_dim(dg) diff --git a/src/tensors/tensor02sym/tensor02sym.jl b/src/tensors/tensor02sym/tensor02sym.jl index fae7331..418aac7 100644 --- a/src/tensors/tensor02sym/tensor02sym.jl +++ b/src/tensors/tensor02sym/tensor02sym.jl @@ -4,12 +4,46 @@ # NOTE (phasing): the operator form and bilinear action delegate to the full # (0,2)-tensor and are deferred to Phase 5 with the operator layer. -"""Space of symmetric (0,2)-tensors Sym²(Ω¹(M)); `dim(dim+1)/2` components/point.""" +""" + Tensor02SymSpace(dg) + +Space of symmetric (0,2)-tensors, with `d(d+1)/2` components per point (see +[`get_symmetric_basis_indices`](@ref) for their order). Its metric carries the ×2 +weights that make the stored components account for the mirrored ones (see +[`gamma_02_sym`](@ref)). Get it with [`tensor02sym_space`](@ref). + +# Examples +```jldoctest +julia> tensor02sym_space(dg) +Tensor02SymSpace(dim=24) + +julia> component_dim(tensor02sym_space(dg)) # (T₁₁, T₁₂, T₂₂) +3 +``` +""" struct Tensor02SymSpace <: AbstractTensorSpace dg::DiffusionGeometry end -"""A symmetric (0,2)-tensor `S ∈ Sym²(Ω¹(M))`.""" +""" + Tensor02Sym + +A symmetric (0,2)-tensor. The codomain of [`hessian`](@ref), and what +[`symmetrise`](@ref) returns. It acts on vector fields exactly as a [`Tensor02`](@ref) +does, by expanding to one first. + +# Examples +```jldoctest +julia> H = hessian(f) +Tensor02Sym(space=Tensor02SymSpace(dim=24), shape=(24,), batch_shape=()) + +julia> size(H(grad(f), grad(f))) # Hess f (∇f, ∇f), a value per point +(60,) + +julia> transpose_tensor(H) === H # symmetric: its own transpose +true +``` +""" struct Tensor02Sym{A<:AbstractArray} <: AbstractTensor space::Tensor02SymSpace coeffs::A @@ -24,7 +58,26 @@ function wrap(space::Tensor02SymSpace, coeffs::AbstractArray) return Tensor02Sym(space, coeffs) end -"""Full (0,2)-tensor view, expanding the symmetric coefficients.""" +""" + full_tensor(S::Tensor02Sym) -> Tensor02 + +The same tensor written out in all `d²` components, mirroring each off-diagonal onto +its transpose. The inverse of [`symmetrise`](@ref) on symmetric tensors. + +# Examples +```jldoctest +julia> H = hessian(f); + +julia> T = full_tensor(H) +Tensor02(space=Tensor02Space(dim=32), shape=(32,), batch_shape=()) + +julia> symmetrise(T).coeffs ≈ H.coeffs # a round trip +true + +julia> transpose_tensor(T).coeffs ≈ T.coeffs # and it really is symmetric +true +``` +""" function full_tensor(S::Tensor02Sym) dg = geometry(S) coeffs_full = expand_symmetric_tensor_coeffs(S.coeffs, n_coefficients(dg), ambient_dim(dg)) diff --git a/src/tensors/vector_fields/vector_fields.jl b/src/tensors/vector_fields/vector_fields.jl index 71c6932..13bfe4d 100644 --- a/src/tensors/vector_fields/vector_fields.jl +++ b/src/tensors/vector_fields/vector_fields.jl @@ -8,12 +8,49 @@ using OMEinsum: @ein_str -"""Space of vector fields 𝔛(M); `dim` components per point.""" +""" + VectorFieldSpace(dg) + +Space of vector fields `𝔛(M)`, with `ambient_dim` components per point — the +coefficients of the frame `{φ_i ∇x_j}`. Get it with [`vector_field_space`](@ref). + +# Examples +```jldoctest +julia> vector_field_space(dg) +VectorFieldSpace(dim=16) + +julia> component_dim(vector_field_space(dg)) # one component per coordinate +2 +``` +""" struct VectorFieldSpace <: AbstractTensorSpace dg::DiffusionGeometry end -"""A vector field `X ∈ 𝔛(M)`, coefficients in the basis {φ_i ∇x_j}.""" +""" + VectorField + +A vector field `X ∈ 𝔛(M)`, stored as coefficients in the basis `{φ_i ∇x_j}`. Build one +with [`grad`](@ref) or [`dg_vector_field`](@ref). + +A vector field is also a *derivation*: call it on a function to get the directional +derivative `X(f) = g(X, ∇f)`. + +# Examples +```jldoctest +julia> X = grad(f) +VectorField(space=VectorFieldSpace(dim=16), shape=(16,), batch_shape=()) + +julia> X(f) # ∇f(f) = |∇f|², a function +ScalarFunction(space=FunctionSpace(dim=8), shape=(8,), batch_shape=()) + +julia> vf_operator(X)(f).coeffs ≈ X(f).coeffs # the call is that operator, applied +true + +julia> flat(X) # lower the index: the 1-form X♭ +Form(space=FormSpace(degree=1, dim=16), shape=(16,), batch_shape=()) +``` +""" struct VectorField{A<:AbstractArray} <: AbstractTensor space::VectorFieldSpace coeffs::A @@ -29,7 +66,27 @@ end degree(::VectorField) = 1 -"""Lower the index with the metric to get a 1-form (shares the coefficients).""" +""" + flat(X::VectorField) -> Form + +Lower the index with the metric (the musical isomorphism ♭): the 1-form `X♭ = g(X, ·)`. +Vector fields and 1-forms have the same coefficients in this basis, so this is a +reinterpretation, not a computation. Inverted by [`sharp`](@ref). + +# Examples +```jldoctest +julia> X = grad(f); + +julia> flat(X) +Form(space=FormSpace(degree=1, dim=16), shape=(16,), batch_shape=()) + +julia> sharp(flat(X)).coeffs == X.coeffs # ♯ ∘ ♭ = id +true + +julia> flat(X).coeffs ≈ d(f).coeffs # and (∇f)♭ = df +true +``` +""" flat(X::VectorField) = wrap(form_space(X.space.dg, 1), X.coeffs) # ── Least-squares reconstruction from ambient arrows ─────────────────────────── @@ -49,7 +106,15 @@ flat(X::VectorField) = wrap(form_space(X.space.dg, 1), X.coeffs) vector_field_to_quiver(dg) -> Array The linear map from vector-field coefficients to ambient arrows, as a -`(n, D, n_coefficients, dim)` tensor (`D` the ambient dimension). +`(n, D, n_coefficients, dim)` tensor (`D` the ambient dimension). The Jacobian of +[`to_ambient`](@ref) for a vector field, and the matrix +[`vector_field_from_reconstruction`](@ref) least-squares solves against. + +# Examples +```jldoctest +julia> size(vector_field_to_quiver(dg)) # (n, ambient D, n_coefficients, dim) +(60, 2, 8, 2) +``` """ function vector_field_to_quiver(dg::DiffusionGeometry) gamma = gamma_ambient(dg.cache) # (n, D, dim) @@ -62,7 +127,21 @@ end Fit a vector field to ambient arrows `data` of trailing shape `(n, D)` by least squares against [`vector_field_to_quiver`](@ref). The left inverse of -`to_ambient`; batched over any leading axes. +[`to_ambient`](@ref); batched over any leading axes. Reached through +`dg_vector_field(dg, data; mode=:reconstruct)`. + +Unlike the `:pullback` mode, this *round-trips*: the arrows you get back out of +`to_ambient` are the arrows you put in (to the extent they were representable). + +# Examples +```jldoctest +julia> arrows = to_ambient(grad(f)); # (60, 2) ambient quiver + +julia> X = vector_field_from_reconstruction(dg, arrows); + +julia> isapprox(to_ambient(X), arrows; atol=1e-8) +true +``` """ function vector_field_from_reconstruction(dg::DiffusionGeometry, data::AbstractArray) n, n1, d = npoints(dg), n_coefficients(dg), ambient_dim(dg) diff --git a/src/utils/basis_conversions.jl b/src/utils/basis_conversions.jl index e6df269..5119dd0 100644 --- a/src/utils/basis_conversions.jl +++ b/src/utils/basis_conversions.jl @@ -13,7 +13,35 @@ using OMEinsum: @ein_str, @optein_str _from_pointwise_basis(data, space; basis_count=nothing) -> Array Convert pointwise `data` (trailing shape `(n, C)`) to coefficients in the function -basis (trailing length `basis_count · C`, defaulting to `space_dim ÷ C`). +basis (trailing length `basis_count · C`, defaulting to `space_dim ÷ C`). The engine +behind [`from_pointwise`](@ref); inverted (up to the truncation) by +[`_to_pointwise_basis`](@ref). + +It solves the weak formulation `G c = ∫ φᵢ · data dμ`, so the result is the +*L²-projection* of the data onto the span of the basis, not an interpolation. + +# Examples +`cos θ` lives entirely in the λ = 1 eigenspace of the circle's Laplacian, so only +coefficients 2 and 3 are non-zero — and by Parseval they carry all of its L² norm, +`∫cos²θ = 1/2`: + +```jldoctest +julia> c = DiffusionGeometryJ._from_pointwise_basis(reshape(cos.(θ), 60, 1), function_space(dg)); + +julia> findall(>(1e-3), abs.(c)) # supported on the λ = 1 eigenspace alone +2-element Vector{Int64}: + 2 + 3 + +julia> round(sum(abs2, c); digits=4) +0.5 +``` + +!!! note "Individual coefficients are not reproducible" + λ = 1 has multiplicity two on the circle, so the eigensolver is free to return any + orthonormal pair spanning that plane. `c[2]` and `c[3]` therefore differ between + BLAS versions (only `c[2]² + c[3]²` is fixed). Never assert on a single coefficient + of a function that straddles a degenerate eigenspace. """ function _from_pointwise_basis(data::AbstractArray, space; basis_count=nothing) C = component_dim(space) @@ -39,7 +67,22 @@ end _to_pointwise_basis(coeffs, space; basis_count=nothing) -> Array Convert coefficients (trailing length `basis_count · C`) to pointwise data -(trailing length `n · C`). `basis_count` defaults to `n_coefficients`. +(trailing length `n · C`). `basis_count` defaults to `n_coefficients`. The engine +behind [`to_pointwise_basis`](@ref) — it just evaluates `Σᵢ cᵢ φᵢ` at the points. + +# Examples +Evaluating the coefficients of `f` back at the 60 sample points recovers `cos θ`, up +to the error of the 8-function truncation: + +```jldoctest +julia> values = DiffusionGeometryJ._to_pointwise_basis(f.coeffs, function_space(dg); basis_count=8); + +julia> size(values) +(60,) + +julia> maximum(abs.(values .- cos.(θ))) < 1e-4 +true +``` """ function _to_pointwise_basis(coeffs::AbstractArray, space; basis_count=nothing) dgv = space.dg diff --git a/src/utils/basis_utils.jl b/src/utils/basis_utils.jl index 75a0568..762b100 100644 --- a/src/utils/basis_utils.jl +++ b/src/utils/basis_utils.jl @@ -12,6 +12,21 @@ using Combinatorics: combinations, with_replacement_combinations Multi-indices for the symmetric (0,2)-tensor basis: `combinations_with_replacement(1:d, 2)`. Shape `(d*(d+1)÷2, 2)`, rows in lexicographic order. + +Row `a` names the basis element `dx_j ⊙ dx_k`, so it is the layout the symmetric +coefficients of a [`Tensor02Sym`](@ref) are stored in. + +# Examples +```jldoctest +julia> get_symmetric_basis_indices(3) +6×2 Matrix{Int64}: + 1 1 + 1 2 + 1 3 + 2 2 + 2 3 + 3 3 +``` """ function get_symmetric_basis_indices(d::Integer) rows = collect(with_replacement_combinations(1:d, 2)) @@ -23,6 +38,24 @@ end Multi-indices for the degree-`k` wedge basis: `combinations(1:d, k)`. Shape `(binomial(d, k), k)`, rows in lexicographic order. + +Row `a` names the basis k-form `dx_{J₁} ∧ … ∧ dx_{J_k}`, which is the component +order a [`Form`](@ref) of degree `k` stores its coefficients in. + +# Examples +```jldoctest +julia> get_wedge_basis_indices(3, 2) # dx₁∧dx₂, dx₁∧dx₃, dx₂∧dx₃ +3×2 Matrix{Int64}: + 1 2 + 1 3 + 2 3 + +julia> size(get_wedge_basis_indices(4, 2)) # binomial(4, 2) basis elements +(6, 2) + +julia> get_wedge_basis_indices(3, 0) # the (single) empty combination +1×0 Matrix{Int64} +``` """ function get_wedge_basis_indices(d::Integer, k::Integer) k == 0 && return Matrix{Int}(undef, 1, 0) @@ -39,6 +72,21 @@ end Lexicographic rank (1-based) of strictly-increasing `k`-combinations drawn from `1:n`. `idx` is either a length-`k` vector (returns an `Int`) or an `(m, k)` matrix of such combinations (returns a length-`m` vector). + +The inverse of [`get_wedge_basis_indices`](@ref): it maps a multi-index back to the +row of the wedge basis that holds it. + +# Examples +```jldoctest +julia> lex_rank([1, 3], 3) # dx₁∧dx₃ is the 2nd basis 2-form of ℝ³ +2 + +julia> lex_rank(get_wedge_basis_indices(3, 2), 3) # the basis ranks itself, in order +3-element Vector{Int64}: + 1 + 2 + 3 +``` """ function lex_rank(idx::AbstractVector{<:Integer}, n::Integer) k = length(idx) @@ -64,6 +112,34 @@ Indices and signs for the wedge product of a `k1`-form and a `k2`-form: `dx_I ∧ dx_J = sgn * dx_K`. `target`, `left`, `right` are 1-based ranks into the degree-`(k1+k2)`, `k1`, and `k2` wedge bases respectively; `signs` are `±1`. Returns empty vectors when `k1 + k2 > d`. + +This is the sparse contraction [`wedge`](@ref) runs: each entry says "component +`left` of α times component `right` of β, signed, accumulates into component +`target` of α ∧ β". + +# Examples +Wedging two 1-forms in ℝ³. Target 1 is `dx₁∧dx₂`, and it collects `α₁β₂ - α₂β₁` — +the antisymmetry is in the signs: + +```jldoctest +julia> target, left, right, signs = get_wedge_product_indices(3, 1, 1); + +julia> [target left right Int.(signs)] +6×4 Matrix{Int64}: + 1 1 2 1 + 1 2 1 -1 + 2 1 3 1 + 2 3 1 -1 + 3 2 3 1 + 3 3 2 -1 +``` + +Beyond the top degree there is nothing to build: + +```jldoctest +julia> get_wedge_product_indices(2, 1, 2) # a 3-form in ℝ² is zero +(Int64[], Int64[], Int64[], Int8[]) +``` """ function get_wedge_product_indices(d::Integer, k1::Integer, k2::Integer) ktot = k1 + k2 @@ -110,6 +186,31 @@ end Expand symmetric (0,2)-tensor coefficients (last axis `n1 · d_sym`, `d_sym = d(d+1)/2`) to the full `n1 · d²` basis, mirroring each off-diagonal component onto its transpose. Leading axes are treated as batch dims. + +The coefficient-level machinery behind [`full_tensor`](@ref); left-inverted by +[`symmetrise_tensor_coeffs`](@ref). + +# Examples +In `d = 2` the symmetric components `(T₁₁, T₁₂, T₂₂)` become the four components of +the full tensor, with `T₁₂` copied into the `T₂₁` slot: + +```jldoctest +julia> expand_symmetric_tensor_coeffs([1.0, 2.0, 3.0], 1, 2) +4-element Vector{Float64}: + 1.0 + 2.0 + 2.0 + 3.0 +``` + +Leading axes are batch dims, so a stack of tensors expands in one call: + +```jldoctest +julia> expand_symmetric_tensor_coeffs([1.0 2.0 3.0; 4.0 5.0 6.0], 1, 2) +2×4 Matrix{Float64}: + 1.0 2.0 2.0 3.0 + 4.0 5.0 5.0 6.0 +``` """ function expand_symmetric_tensor_coeffs(coeffs::AbstractArray, n_coefficients::Integer, d::Integer) n1 = n_coefficients @@ -139,6 +240,30 @@ end Project full (0,2)-tensor coefficients (last axis `n1 · d²`) onto the symmetric subspace, returning `n1 · d_sym` coefficients. Off-diagonal components average `(T_{jk} + T_{kj})/2`. Leading axes are batch dims. + +The coefficient-level machinery behind [`symmetrise`](@ref). + +# Examples +The full `d = 2` tensor `[1 2; 3 4]` (row-major) symmetrises to `(1, 2.5, 4)`: + +```jldoctest +julia> symmetrise_tensor_coeffs([1.0, 2.0, 3.0, 4.0], 1, 2) +3-element Vector{Float64}: + 1.0 + 2.5 + 4.0 +``` + +It is a left inverse of [`expand_symmetric_tensor_coeffs`](@ref) — expanding, then +symmetrising, is the identity on symmetric coefficients: + +```jldoctest +julia> symmetrise_tensor_coeffs(expand_symmetric_tensor_coeffs([1.0, 2.0, 3.0], 1, 2), 1, 2) +3-element Vector{Float64}: + 1.0 + 2.0 + 3.0 +``` """ function symmetrise_tensor_coeffs(coeffs::AbstractArray, n_coefficients::Integer, d::Integer) n1 = n_coefficients @@ -169,6 +294,34 @@ end For degree `k` (`1 ≤ k ≤ d-1`): the degree-`k` and degree-`(k+1)` wedge bases, plus `children[Jp, r]` = the 1-based rank in `idx_k` of `Jp` with its `r`-th entry removed, and the alternating Laplace-expansion `signs = (-1)^r`. + +This is the combinatorial skeleton of the exterior derivative: `(dω)_{Jp}` sums the +`k`-form components of the children of `Jp`, signed alternately. + +# Examples +The children of the basis 2-forms of ℝ³. Dropping an entry of `Jp = (1, 2)` leaves +`(2)` (rank 2) or `(1)` (rank 1), which enter with signs `+1` and `-1`: + +```jldoctest +julia> idx_k, idx_kp1, children, signs = kp1_children_and_signs(3, 1); + +julia> idx_kp1 # the 2-forms dx₁∧dx₂, dx₁∧dx₃, dx₂∧dx₃ +3×2 Matrix{Int64}: + 1 2 + 1 3 + 2 3 + +julia> children # the 1-form rank left after dropping entry r +3×2 Matrix{Int64}: + 2 1 + 3 1 + 3 2 + +julia> signs +2-element Vector{Int64}: + 1 + -1 +``` """ function kp1_children_and_signs(d::Integer, k::Integer) @assert 1 <= k <= d - 1 "Form degree k=$k must be between 1 and $(d-1)" @@ -193,7 +346,35 @@ end All `k!` permutations of `1:k`, as rows of a `(k!, k)` matrix in lexicographic order, together with their parities `(-1)^{inversions}`. Used to expand a k-form's -wedge-basis coefficients into a fully antisymmetric `(d,…,d)` tensor. +wedge-basis coefficients into a fully antisymmetric `(d,…,d)` tensor (see +[`to_ambient`](@ref)). + +The parity counts *inversions*, so the identity permutation always carries `+1`. +(The Python reference counts concordant pairs instead, which flips the sign of the +whole polyvector for `k ≡ 2, 3 (mod 4)`; see `docs/src/upstream-bugs.md`.) + +# Examples +```jldoctest +julia> perms, signs = permutations_with_signs(3); + +julia> perms +6×3 Matrix{Int64}: + 1 2 3 + 1 3 2 + 2 1 3 + 2 3 1 + 3 1 2 + 3 2 1 + +julia> signs # (1 2 3) is even; one transposition is odd +6-element Vector{Int64}: + 1 + -1 + -1 + 1 + 1 + -1 +``` """ function permutations_with_signs(k::Integer) k == 0 && return zeros(Int, 1, 0), Int[1] diff --git a/src/utils/batch_utils.jl b/src/utils/batch_utils.jl index 0f25a04..8104248 100644 --- a/src/utils/batch_utils.jl +++ b/src/utils/batch_utils.jl @@ -13,6 +13,21 @@ Broadcast batch shapes under numpy's rule: axes are matched from the *right*, an length-1 axis stretches to meet any other length. Julia's `Base.Broadcast` matches from the left, so `(3, 4)` with `(4,)` broadcasts here but not there — hence this function rather than a delegation. Throws `DimensionMismatch` if the shapes conflict. + +# Examples +```jldoctest +julia> DiffusionGeometryJ.broadcast_batch_shape((3, 4), (4,)) # right-aligned: (4,) fills the last axis +(3, 4) + +julia> DiffusionGeometryJ.broadcast_batch_shape((3, 1), (1, 4)) # length-1 axes stretch +(3, 4) + +julia> DiffusionGeometryJ.broadcast_batch_shape((), (5,)) # an unbatched tensor joins any batch +(5,) + +julia> DiffusionGeometryJ.broadcast_batch_shape((3,), (4,)) +ERROR: DimensionMismatch: Incompatible batch shapes: ((3,), (4,)) +``` """ function broadcast_batch_shape(shapes::Tuple...) rank = maximum(length, shapes; init=0) @@ -37,6 +52,17 @@ end compatible_batches(shape1, shape2) -> Bool `true` if two batch shapes are broadcast-compatible (numpy `broadcast_shapes`). +The predicate every binary tensor operation checks before it broadcasts; see +[`broadcast_batch_shape`](@ref) for the rule. + +# Examples +```jldoctest +julia> compatible_batches((3, 4), (4,)) +true + +julia> compatible_batches((3,), (4,)) +false +``` """ function compatible_batches(shape1::Tuple, shape2::Tuple) try @@ -54,6 +80,15 @@ Prepend singleton axes so the array's batch rank is `rank`, right-aligning the b axes it already has. Julia's own broadcasting then reproduces numpy's semantics: an `ntail=1` array with batch `(4,)` becomes `(1, 4, L)`, which broadcasts against `(3, 4, L)`. The array's `ntail` trailing (non-batch) axes are left alone. + +# Examples +```jldoctest +julia> size(DiffusionGeometryJ.pad_batch_dims(zeros(4, 6), 2)) # batch (4,), coeff axis 6 +(1, 4, 6) + +julia> size(DiffusionGeometryJ.pad_batch_dims(zeros(4, 6), 1)) # already rank 1: untouched +(4, 6) +``` """ function pad_batch_dims(array::AbstractArray, rank::Integer; ntail::Integer=1) npad = rank - (ndims(array) - ntail) @@ -68,6 +103,17 @@ end Right-align the batch axes of two arrays against each other by padding the shorter batch with singleton axes, so that `a .* b` (or any elementwise op) follows numpy's rule. Both arrays keep their `ntail` trailing (non-batch) axes. + +# Examples +```jldoctest +julia> a, b = DiffusionGeometryJ.align_batch_pair(zeros(3, 4, 6), zeros(4, 6)); + +julia> size(a), size(b) # b gained a leading singleton axis +((3, 4, 6), (1, 4, 6)) + +julia> size(a .+ b) # so Julia's broadcasting now agrees with numpy +(3, 4, 6) +``` """ function align_batch_pair(a::AbstractArray, b::AbstractArray; ntail::Integer=1) rank = max(ndims(a), ndims(b)) - ntail @@ -79,7 +125,18 @@ end Materialise `array` with its batch axes expanded to `target` (numpy right-alignment), keeping its `ntail` trailing axes. Use when the batch axes must be flattened for an -einsum; prefer [`pad_batch_dims`] when plain broadcasting will do. +einsum; prefer [`pad_batch_dims`](@ref) when plain broadcasting will do. + +# Examples +```jldoctest +julia> size(DiffusionGeometryJ.broadcast_batch_to(ones(4, 6), (3, 4))) # copied along the new axis +(3, 4, 6) + +julia> A = ones(4, 6); + +julia> DiffusionGeometryJ.broadcast_batch_to(A, (4,)) === A # already the target: no copy +true +``` """ function broadcast_batch_to(array::AbstractArray, target::Tuple; ntail::Integer=1) nb = ndims(array) - ntail @@ -95,6 +152,12 @@ end Expand the batch axes to `target` and collapse them (row-major) into a single leading axis, giving `(prod(target), tail...)` — the layout the einsum kernels want. + +# Examples +```jldoctest +julia> size(DiffusionGeometryJ.broadcast_flatten_batch(ones(4, 6), (3, 4))) # batch (3, 4) → 12 rows +(12, 6) +``` """ function broadcast_flatten_batch(array::AbstractArray, target::Tuple; ntail::Integer=1) expanded = broadcast_batch_to(array, target; ntail=ntail) @@ -106,7 +169,19 @@ end infer_batch_shape(array, expected_tail) -> (array, batch_shape) Validate the trailing dimensions of `array` against `expected_tail` (a tuple) and -return the array together with its leading batch shape. +return the array together with its leading batch shape. This is the check every +[`wrap`](@ref) runs: it is what rejects a coefficient vector of the wrong length. + +# Examples +```jldoctest +julia> _, batch = DiffusionGeometryJ.infer_batch_shape(zeros(3, 5, 16), (16,)); + +julia> batch +(3, 5) + +julia> DiffusionGeometryJ.infer_batch_shape(zeros(3, 15), (16,); name="VectorField") +ERROR: AssertionError: VectorField coefficients must have trailing shape (16,), got (3, 15) +``` """ function infer_batch_shape(array::AbstractArray, expected_tail::Tuple; name::AbstractString="Tensor") nt = length(expected_tail) @@ -122,6 +197,20 @@ end Collapse the batch axes (everything but the last) into a single leading axis, row-major, giving a `(prod(batch), L)` matrix plus the original batch shape. +Undone by [`restore_batch_dims`](@ref). + +# Examples +```jldoctest +julia> flat, batch = DiffusionGeometryJ.flatten_batch_dims(zeros(2, 3, 16)); + +julia> size(flat), batch +((6, 16), (2, 3)) + +julia> flat, batch = DiffusionGeometryJ.flatten_batch_dims(zeros(16)); # unbatched: one row, empty batch + +julia> size(flat), batch +((1, 16), ()) +``` """ function flatten_batch_dims(array::AbstractArray) @assert ndims(array) > 0 "Coefficient arrays must have at least one dimension" @@ -134,8 +223,21 @@ end """ restore_batch_dims(flat, batch_shape) -> Array -Inverse of [`flatten_batch_dims`]: reshape a `(prod(batch), L)` matrix back to +Inverse of [`flatten_batch_dims`](@ref): reshape a `(prod(batch), L)` matrix back to `batch_shape × (L,)` (row-major). An empty batch shape yields a length-`L` vector. + +# Examples +```jldoctest +julia> A = reshape(1.0:96.0, 2, 3, 16); + +julia> flat, batch = DiffusionGeometryJ.flatten_batch_dims(A); + +julia> DiffusionGeometryJ.restore_batch_dims(flat, batch) == A # a round trip +true + +julia> size(DiffusionGeometryJ.restore_batch_dims(zeros(1, 16), ())) # empty batch → a plain vector +(16,) +``` """ function restore_batch_dims(flat::AbstractMatrix, batch_shape::Tuple) L = size(flat, 2) diff --git a/src/utils/reshape_utils.jl b/src/utils/reshape_utils.jl index a05ca22..bd0a3d5 100644 --- a/src/utils/reshape_utils.jl +++ b/src/utils/reshape_utils.jl @@ -11,7 +11,34 @@ Reshape `A` with numpy C-order (row-major) semantics, so the result matches `numpy.reshape(A, dims)` element-for-element. Used to flatten the multi-index -weak-operator tensors into the matrices the parity fixtures store. +weak-operator tensors into the matrices the parity fixtures store, and to lay out +every coefficient vector in the package (`(i, a) ↦ i·C + a`, basis index slowest). + +# Examples +Reading order is by rows, not columns — this is what `Base.reshape` would *not* do: + +```jldoctest +julia> DiffusionGeometryJ.np_reshape([1 2 3; 4 5 6], 3, 2) +3×2 Matrix{Int64}: + 1 2 + 3 4 + 5 6 + +julia> reshape([1 2 3; 4 5 6], 3, 2) # Julia's column-major reshape, for contrast +3×2 Matrix{Int64}: + 1 5 + 4 3 + 2 6 +``` + +It round-trips, so a flattened coefficient vector can always be unpacked again: + +```jldoctest +julia> A = reshape(1:12, 2, 2, 3); + +julia> DiffusionGeometryJ.np_reshape(DiffusionGeometryJ.np_reshape(A, 2, 6), 2, 2, 3) == A +true +``` """ function np_reshape(A::AbstractArray, dims::Integer...) nin = ndims(A) diff --git a/src/visualisation/api.jl b/src/visualisation/api.jl index ec3557c..0a049c4 100644 --- a/src/visualisation/api.jl +++ b/src/visualisation/api.jl @@ -12,7 +12,7 @@ """Diverging blue → grey → red colour map, the package default for signed fields.""" const DIVERGING_COLORMAP = ["#166dde", "#d3d3d3", "#e32636"] -"""Cyclic variant of [`DIVERGING_COLORMAP`], for angle-valued fields on `[0, 2π]`.""" +"""Cyclic variant of [`DIVERGING_COLORMAP`](@ref), for angle-valued fields on `[0, 2π]`.""" const CYCLIC_COLORMAP = ["#166dde", "#d3d3d3", "#e32636", "#166dde"] """ @@ -23,11 +23,11 @@ Plot a tensor over the point cloud it lives on, choosing the visual from its typ | argument | visual | |:-----------------------------|:----------------------------------------------| -| [`ScalarFunction`] | scatter, coloured by value | -| [`VectorField`], 1-[`Form`] | quiver of ambient arrows | -| 2-[`Form`] | oriented discs (`d = 2` filled, `d = 3` pairs) | -| 3-[`Form`] | scatter sized and coloured by `ω₁₂₃` | -| [`Tensor02`], [`Tensor02Sym`]| ellipses (`d = 2`) or ellipsoids (`d = 3`) | +| [`ScalarFunction`](@ref) | scatter, coloured by value | +| [`VectorField`](@ref), 1-[`Form`](@ref) | quiver of ambient arrows | +| 2-[`Form`](@ref) | oriented discs (`d = 2` filled, `d = 3` pairs) | +| 3-[`Form`](@ref) | scatter sized and coloured by `ω₁₂₃` | +| [`Tensor02`](@ref), [`Tensor02Sym`](@ref)| ellipses (`d = 2`) or ellipsoids (`d = 3`) | The points come from `immersion_coords(geometry(tensor))` and the values from `to_ambient(tensor)`, so `dgplot(X)` is the whole of the Python idiom @@ -42,12 +42,29 @@ ratio — coordinates on an embedded manifold carry no meaning of their own. Pas into an axis you own and leaves it alone. Requires a Makie backend (`using GLMakie` or `using CairoMakie`). Other keyword arguments -are forwarded to the underlying recipe — see [`dgscatter`], [`dgquiver`], -[`dg2form`], [`dg3form`] and [`dgellipsoids`] for the attributes each accepts. +are forwarded to the underlying recipe — see [`dgscatter`](@ref), [`dgquiver`](@ref), +[`dg2form`](@ref), [`dg3form`](@ref) and [`dgellipsoids`](@ref) for the attributes each accepts. + +# Examples +```julia +using GLMakie, DiffusionGeometryJ + +θ = range(0, 2π; length=61)[1:60] +dg = from_point_cloud([cos.(θ) sin.(θ)]; knn_kernel=16, n_function_basis=8) +f = dg_function(dg, cos.(θ)) + +dgplot(f) # the function, as a coloured scatter +dgplot(grad(f); scale=0.3) # its gradient, as a quiver +dgplot(hessian(f)) # its Hessian, as ellipses + +fig, ax, _ = dgplot(f) # the Makie figure/axis/plot, as usual +dgplot!(ax, grad(f)) # draw the quiver over it +fig +``` """ function dgplot end -"""See [`dgplot`].""" +"""See [`dgplot`](@ref).""" function dgplot! end """ @@ -56,11 +73,11 @@ function dgplot! end Scatter `points` (an `(n, d)` matrix or vector of `Point`s) coloured by `values`. `colorrange` defaults to the symmetric range `(-a, a)` with `a = maximum(abs, values)`, so zero sits at the centre of the diverging colour map. Pass `cyclic = true` for -angle-valued data, which instead uses `(0, 2π)` and [`CYCLIC_COLORMAP`]. +angle-valued data, which instead uses `(0, 2π)` and [`CYCLIC_COLORMAP`](@ref). """ function dgscatter end -"""See [`dgscatter`].""" +"""See [`dgscatter`](@ref).""" function dgscatter! end """ @@ -77,7 +94,7 @@ arrow length. """ function dgquiver end -"""See [`dgquiver`].""" +"""See [`dgquiver`](@ref).""" function dgquiver! end """ @@ -85,14 +102,14 @@ function dgquiver! end Visualise a 2-form given by ambient skew-symmetric `matrices` of shape `(n, d, d)`. -In `d = 2` the form reduces (via [`hodge_star_2_form`]) to a scalar, drawn as filled +In `d = 2` the form reduces (via [`hodge_star_2_form`](@ref)) to a scalar, drawn as filled discs whose radius tracks magnitude and whose colour tracks sign. In `d = 3` it reduces to an axial vector, drawn as a pair of discs normal to it — red on the positive face, blue on the negative — so the orientation is visible from either side. """ function dg2form end -"""See [`dg2form`].""" +"""See [`dg2form`](@ref).""" function dg2form! end """ @@ -103,7 +120,7 @@ Visualise a 3-form on a 3-dimensional immersion, given by ambient `tensors` of s """ function dg3form end -"""See [`dg3form`].""" +"""See [`dg3form`](@ref).""" function dg3form! end """ @@ -119,7 +136,7 @@ Asymmetric input is symmetrised — an ellipsoid can only represent the symmetri """ function dgellipsoids end -"""See [`dgellipsoids`].""" +"""See [`dgellipsoids`](@ref).""" function dgellipsoids! end """ @@ -132,7 +149,7 @@ of `λ`. Used for Hessian eigen-directions. """ function dgeiglines end -"""See [`dgeiglines`].""" +"""See [`dgeiglines`](@ref).""" function dgeiglines! end """ @@ -144,19 +161,27 @@ spacing of `points` so the quads sit flush with the cloud. """ function dgtangentplanes end -"""See [`dgtangentplanes`].""" +"""See [`dgtangentplanes`](@ref).""" function dgtangentplanes! end """ dganimate(ft, filename; framerate=30, kwargs...) Record the time evolution of a batched tensor `ft` — as returned by -[`solve_differential_operator`] with a leading time axis — to `filename` +[`solve_differential_operator`](@ref) with a leading time axis — to `filename` (`.mp4` or `.gif`, chosen by extension). Returns `filename`. The colour range is held fixed across frames (spanning all of `ft` unless `colorrange` is given) so the animation shows the field decaying rather than being renormalised each frame. Replaces the Python `methods/pde.py::gif_from_functions`. + +# Examples +```julia +using GLMakie, DiffusionGeometryJ + +u = solve_differential_operator(-laplacian(dg, 0), f, range(0, 2; length=60)) +dganimate(u, "heat.mp4"; framerate=30) +``` """ function dganimate end diff --git a/src/visualisation/hodge_star.jl b/src/visualisation/hodge_star.jl index a5030f9..8990e23 100644 --- a/src/visualisation/hodge_star.jl +++ b/src/visualisation/hodge_star.jl @@ -20,6 +20,32 @@ left-handed (which negates the result). Only `d = 2` and `d = 3` are supported. The input is assumed to be in ambient (Euclidean) coordinates, e.g. `to_ambient(ω)` for a `Form` of degree 2. + +Note the sign convention: `(*ω) = -orientation · ε^{ij} ω_ij / 2`, matching the Python +reference (whose 2-form plots this reproduces). + +# Examples +```jldoctest +julia> ω = reshape([0.0 1.0; -1.0 0.0], 1, 2, 2); # one point, d = 2 + +julia> hodge_star_2_form(ω) +1-element Vector{Float64}: + -1.0 + +julia> hodge_star_2_form(ω; orientation=-1) # the other handedness +1-element Vector{Float64}: + 1.0 +``` + +In three dimensions the dual is the axial vector `(ω₂₃, ω₃₁, ω₁₂)`, signed: + +```jldoctest +julia> ω = reshape([0.0 3.0 -2.0; -3.0 0.0 1.0; 2.0 -1.0 0.0], 1, 3, 3); + +julia> hodge_star_2_form(ω) +1×3 Matrix{Float64}: + -1.0 -2.0 -3.0 +``` """ function hodge_star_2_form(omega::AbstractArray; orientation::Integer = 1) ndims(omega) == 3 || throw(ArgumentError("omega must have shape (n, d, d), got $(size(omega))")) @@ -46,11 +72,23 @@ end """ hodge_star_2_form(ω::Form; orientation=1) -> Array -Hodge dual of a degree-2 `Form`, taken in its ambient representation. +Hodge dual of a degree-2 [`Form`](@ref), taken in its ambient representation. Equivalent to `hodge_star_2_form(to_ambient(ω); orientation)`. For a vector field `X` on a 3-dimensional immersion, `hodge_star_2_form(d(flat(X)))` is the curl of `X`. + +# Examples +The curl of a vector field on the sphere, as an ambient axial vector per point: + +```jldoctest +julia> X = grad(dg_function(dg3, sphere[:, 3])); + +julia> curl = hodge_star_2_form(d(flat(X))); + +julia> size(curl) +(200, 3) +``` """ function hodge_star_2_form(ω::Form; orientation::Integer = 1) degree(ω) == 2 || throw(ArgumentError("hodge_star_2_form needs a 2-form, got degree $(degree(ω))")) diff --git a/test/doctests.jl b/test/doctests.jl new file mode 100644 index 0000000..26e89e7 --- /dev/null +++ b/test/doctests.jl @@ -0,0 +1,20 @@ +# Run every `jldoctest` in the package's docstrings as part of the test suite. +# +# The examples in the docstrings are the package's spec-by-example: they assert real +# values (Laplacian eigenvalues of the circle, adjointness of δ against d, the wedge's +# antisymmetry), so a change that breaks one is a change in behaviour, not just in prose. +# +# The preamble that defines `dg`, `dg3`, `f` and the point clouds lives in +# `docs/doctest_setup.jl`, shared with `docs/make.jl`. + +using Documenter + +include(joinpath(@__DIR__, "..", "docs", "doctest_setup.jl")) +DocMeta.setdocmeta!(DiffusionGeometryJ, :DocTestSetup, DOCTEST_SETUP; recursive=true) + +@testset "doctests" begin + # `manual=true` also runs the `jldoctest` blocks in the manual pages under + # `docs/src`, so the examples in the guide are held to the suite too, not just + # the ones in docstrings. + doctest(DiffusionGeometryJ; manual=true) +end diff --git a/test/runtests.jl b/test/runtests.jl index f30e32e..16fa3c8 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -21,6 +21,9 @@ include("pysuite.jl") include("test_notebooks.jl") include("test_makie_ext.jl") + # The docstring examples: they assert behaviour, so they are tests. + include("doctests.jl") + # Ports of the upstream Python suite (property tests, no fixtures). @testset "python tests/test_src" begin include("test_pysrc.jl") diff --git a/test/test_pysrc.jl b/test/test_pysrc.jl index ca26cbb..01cc3ca 100644 --- a/test/test_pysrc.jl +++ b/test/test_pysrc.jl @@ -176,7 +176,7 @@ end # Upstream asserts only the shape here — its value check is commented out # (`NOTE: Disabling strict value check`). The check is enabled below; the # expectation upstream compares against is the thing that is wrong, not the - # builder. See docs/upstream-bugs.md §5. + # builder. See docs/src/upstream-bugs.md §5. @testset "hessian_02_sym_weak matrix" begin H_weak = weak(hessian(dg)) @test size(H_weak) == (n1 * d_sym, n0)