From c3bbd8461fec40ecb98a70c86484c67a4e635e24 Mon Sep 17 00:00:00 2001 From: Jeffrey Ward Date: Sat, 11 Jul 2026 22:14:35 -0400 Subject: [PATCH 1/3] docs: add doctests to every function and a Documenter site MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Every documented function now carries runnable `# Examples` (223 jldoctest blocks), and the 28 exports that had no docstring at all (degree, cdc_components, wrap, the accessors, the dg_* factories) have one. The examples assert behaviour rather than illustrate it: the circle's Laplacian spectrum, adjointness of δ against d, antisymmetry of the wedge, ♯∘♭ = id, and each weak-form builder reproducing exactly the matrix its operator carries. They run under `Pkg.test()` (test/doctests.jl) and in the docs build, sharing the preamble in docs/doctest_setup.jl. Writing them turned up three things worth recording: - A function's individual coefficients are not reproducible when it straddles a degenerate eigenspace: λ=1 has multiplicity two on the circle, so the eigensolver returns an arbitrary orthonormal pair spanning that plane. Only the Parseval norm is fixed. Noted on _from_pointwise_basis. - betti_spectrum's harmonic eigenvalue *magnitude* is machine-dependent — the w=1e10 penalty multiplies the float noise with it — while the gap that gives the Betti number is not. BettiSpectrum now warns to count the gap, never a threshold. - laplacian(dg, 0) and codifferential ∘ d are not the same matrix: composing passes through the Gram pseudo-inverse, so the two agree on the resolved modes and part company on the truncation tail. Documented in Conventions. DiffusionGeometry, MarkovTriple, ImmersedMarkovTriple and GammaCache also gain compact `show` methods — printing a geometry previously dumped 123 KB of eigenbasis to the terminal. docs/plotting.md and docs/upstream-bugs.md move under docs/src so they are part of the built site; the three references to them are updated. CI gains a docs job that builds the site (and so runs the doctests) with checkdocs=:exports, which fails the build on any undocumented export. --- .github/workflows/ci.yml | 18 + .gitignore | 1 + PORTING_PLAN.md | 2 +- Project.toml | 3 +- README.md | 18 +- docs/Project.toml | 6 + docs/doctest_setup.jl | 40 ++ docs/make.jl | 34 ++ docs/src/api/diffusion.md | 50 ++ docs/src/api/geometry.md | 91 ++++ docs/src/api/methods.md | 38 ++ docs/src/api/operators.md | 109 ++++ docs/src/api/tensors.md | 84 ++++ docs/src/api/utils.md | 51 ++ docs/src/conventions.md | 86 ++++ docs/src/index.md | 92 ++++ docs/{ => src}/plotting.md | 31 ++ docs/{ => src}/upstream-bugs.md | 0 src/core/diffusion/carre_du_champ.jl | 105 ++++ src/core/diffusion/diffusion_process.jl | 145 ++++++ src/core/diffusion/markov_triples.jl | 126 ++++- src/core/diffusion/regularise.jl | 35 ++ src/core/geometry/cache.jl | 127 ++++- src/core/geometry/diffusion_geometry.jl | 475 +++++++++++++++++- src/methods/geodesics.jl | 26 + src/methods/pde.jl | 24 + src/methods/topology.jl | 96 ++++ .../differential_operators/derivative.jl | 16 + .../geometry_operators.jl | 255 +++++++++- .../differential_operators/hessian.jl | 48 +- .../differential_operators/laplacian.jl | 17 + .../differential_operators/levi_civita.jl | 16 + .../differential_operators/lie_bracket.jl | 16 + src/operators/tensor_actions.jl | 319 +++++++++++- src/operators/types/bilinear.jl | 115 ++++- src/operators/types/direct_sum.jl | 46 +- src/operators/types/linear.jl | 251 ++++++++- src/tensors/base_tensor/abstract.jl | 199 +++++++- src/tensors/base_tensor/base_tensor.jl | 116 ++++- src/tensors/base_tensor/base_tensor_space.jl | 206 +++++++- src/tensors/base_tensor/metric_gram.jl | 37 +- src/tensors/direct_sum/direct_sum.jl | 94 +++- src/tensors/forms/forms.jl | 98 +++- src/tensors/functions/functions.jl | 45 +- src/tensors/tensor02/tensor02.jl | 78 ++- src/tensors/tensor02sym/tensor02sym.jl | 59 ++- src/tensors/vector_fields/vector_fields.jl | 89 +++- src/utils/basis_conversions.jl | 47 +- src/utils/basis_utils.jl | 183 ++++++- src/utils/batch_utils.jl | 108 +++- src/utils/reshape_utils.jl | 29 +- src/visualisation/api.jl | 63 ++- src/visualisation/hodge_star.jl | 40 +- test/doctests.jl | 17 + test/runtests.jl | 3 + test/test_pysrc.jl | 2 +- 56 files changed, 4382 insertions(+), 143 deletions(-) create mode 100644 docs/Project.toml create mode 100644 docs/doctest_setup.jl create mode 100644 docs/make.jl create mode 100644 docs/src/api/diffusion.md create mode 100644 docs/src/api/geometry.md create mode 100644 docs/src/api/methods.md create mode 100644 docs/src/api/operators.md create mode 100644 docs/src/api/tensors.md create mode 100644 docs/src/api/utils.md create mode 100644 docs/src/conventions.md create mode 100644 docs/src/index.md rename docs/{ => src}/plotting.md (93%) rename docs/{ => src}/upstream-bugs.md (100%) create mode 100644 test/doctests.jl diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 24fa896..7a8d4c5 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -23,3 +23,21 @@ 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 + 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..c58eebf 100644 --- a/README.md +++ b/README.md @@ -257,7 +257,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 @@ -302,6 +302,22 @@ 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, 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`. + ## 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..7cae892 --- /dev/null +++ b/docs/doctest_setup.jl @@ -0,0 +1,40 @@ +# 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 + + # The unexported helpers documented on the "Utilities" page. Their examples are + # written as if they were exported, so bring them into scope here. + using DiffusionGeometryJ: np_reshape, broadcast_batch_shape, pad_batch_dims, + align_batch_pair, broadcast_batch_to, broadcast_flatten_batch, + infer_batch_shape, flatten_batch_dims, restore_batch_dims, + _from_pointwise_basis, _to_pointwise_basis, _is_orthonormal, _gram_spectrum, + gamma_coords_regularised, regularise_fn, + coeffs, space, geometry, batch_shape, space_degree + + # 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. + 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] + global sphere = [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..d6acf11 --- /dev/null +++ b/docs/make.jl @@ -0,0 +1,34 @@ +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 + warnonly=[:missing_docs], + 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..9d1dc8a --- /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> 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/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..26fe1b4 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> _is_orthonormal(function_space(dg)) +true + +julia> _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 = _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..277fa0c 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,19 @@ 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. + +# Examples +```jldoctest +julia> X = grad(f); + +julia> vals = metric_apply(vector_field_space(dg), X.coeffs, X.coeffs); + +julia> size(vals) +(60,) +``` """ function _metric_apply(u_n1::AbstractMatrix, regularise_func, a_coeffs, b_coeffs, matrices::AbstractArray{<:Any,3}) @@ -62,7 +84,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..faf3f10 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 = _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 = _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..41981b9 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> broadcast_batch_shape((3, 4), (4,)) # right-aligned: (4,) fills the last axis +(3, 4) + +julia> broadcast_batch_shape((3, 1), (1, 4)) # length-1 axes stretch +(3, 4) + +julia> broadcast_batch_shape((), (5,)) # an unbatched tensor joins any batch +(5,) + +julia> 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(pad_batch_dims(zeros(4, 6), 2)) # batch (4,), coeff axis 6 +(1, 4, 6) + +julia> size(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 = 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(broadcast_batch_to(ones(4, 6), (3, 4))) # copied along the new axis +(3, 4, 6) + +julia> A = ones(4, 6); + +julia> 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(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 = infer_batch_shape(zeros(3, 5, 16), (16,)); + +julia> batch +(3, 5) + +julia> 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 = flatten_batch_dims(zeros(2, 3, 16)); + +julia> size(flat), batch +((6, 16), (2, 3)) + +julia> flat, batch = 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 = flatten_batch_dims(A); + +julia> restore_batch_dims(flat, batch) == A # a round trip +true + +julia> size(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..949d5f8 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> 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> np_reshape(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..69c033d --- /dev/null +++ b/test/doctests.jl @@ -0,0 +1,17 @@ +# 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 + doctest(DiffusionGeometryJ; manual=false) +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) From 12ff60a1768bd8d86a950929ae7e20f5c21254d4 Mon Sep 17 00:00:00 2001 From: Jeffrey Ward Date: Sat, 11 Jul 2026 22:56:26 -0400 Subject: [PATCH 2/3] docs: make the docs checks real, and the examples honest about scope MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Follow-up to c3bbd84, from a review of it. Four of the guarantees that commit advertised were not actually enforced, and the examples quietly leaned on names a reader cannot call. - `warnonly=[:missing_docs]` in docs/make.jl downgraded `checkdocs=:exports` to a warning, so the claim that the build fails on an undocumented export was vacuous. Dropped it; the strict build passes, so all 179 exports were already documented and nothing was hiding behind the warning. - The CI docs job passed GITHUB_TOKEN but granted no `contents: write` and had no tag trigger, so `deploydocs` could never push to gh-pages and versioned docs would never build. Added the permissions and a `v*` tag trigger. - test/doctests.jl ran with `manual=false`, so the examples in index.md, conventions.md and the API pages were exercised only by the docs job, not by `Pkg.test()` — while README claimed they "run as part of the test suite, so they cannot drift". Now `manual=true`, which makes the README true as written. - The doctest preamble imported twenty unexported names so the examples could call them bare. Copying `_is_orthonormal(function_space(dg))` out of the rendered docs got you an UndefVarError. That last one split in two. The five accessors — `coeffs`, `space`, `geometry`, `batch_shape`, `space_degree` — are the read side of every tensor a user holds and were unexported by oversight (the c3bbd84 message already called them exports), so they are exported now. The genuinely internal helpers (`np_reshape`, the batch utilities, `_from_pointwise_basis`) stay private and their examples spell them `DiffusionGeometryJ.np_reshape(...)`, as a reader would have to. The preamble now imports nothing unexported, so an example that reaches for a private name fails the doctests rather than misleading. Also: `_metric_apply`'s example did not call `_metric_apply` — it re-ran the public `metric_apply` example verbatim. Replaced with a cross-reference, since the kernel takes raw arrays rather than a space. --- .github/workflows/ci.yml | 7 ++++ docs/doctest_setup.jl | 16 ++++----- docs/make.jl | 3 +- docs/src/conventions.md | 2 +- src/DiffusionGeometryJ.jl | 2 ++ src/tensors/base_tensor/base_tensor_space.jl | 6 ++-- src/tensors/base_tensor/metric_gram.jl | 12 ++----- src/utils/basis_conversions.jl | 4 +-- src/utils/batch_utils.jl | 34 ++++++++++---------- src/utils/reshape_utils.jl | 4 +-- test/doctests.jl | 5 ++- 11 files changed, 49 insertions(+), 46 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 7a8d4c5..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: @@ -26,6 +27,12 @@ jobs: 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 diff --git a/docs/doctest_setup.jl b/docs/doctest_setup.jl index 7cae892..ae4a53b 100644 --- a/docs/doctest_setup.jl +++ b/docs/doctest_setup.jl @@ -10,14 +10,10 @@ const DOCTEST_SETUP = quote using DiffusionGeometryJ using LinearAlgebra - # The unexported helpers documented on the "Utilities" page. Their examples are - # written as if they were exported, so bring them into scope here. - using DiffusionGeometryJ: np_reshape, broadcast_batch_shape, pad_batch_dims, - align_batch_pair, broadcast_batch_to, broadcast_flatten_batch, - infer_batch_shape, flatten_batch_dims, restore_batch_dims, - _from_pointwise_basis, _to_pointwise_basis, _is_orthonormal, _gram_spectrum, - gamma_coords_regularised, regularise_fn, - coeffs, space, geometry, batch_shape, space_degree + # 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. @@ -30,11 +26,11 @@ const DOCTEST_SETUP = quote # 200 points on the unit 2-sphere in ℝ³ (Fibonacci lattice), for the examples # that need a surface: 2-forms, curvature, the Hodge star. - let n = 200, ϕ = (1 + sqrt(5)) / 2 + 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] - global sphere = [r .* cos.(ψ) r .* sin.(ψ) z] + [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 index d6acf11..1c0a361 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -10,7 +10,8 @@ makedocs(; modules=[DiffusionGeometryJ], doctest=true, checkdocs=:exports, # every exported name must appear in an @docs block - warnonly=[:missing_docs], + # 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", diff --git a/docs/src/conventions.md b/docs/src/conventions.md index 9d1dc8a..15119e8 100644 --- a/docs/src/conventions.md +++ b/docs/src/conventions.md @@ -41,7 +41,7 @@ carries its own [`broadcast_batch_shape`](@ref). julia> batch_shape(dg_function(dg, ones(3, 60))) # three functions in one tensor (3,) -julia> broadcast_batch_shape((3, 4), (4,)) # right-aligned +julia> DiffusionGeometryJ.broadcast_batch_shape((3, 4), (4,)) # right-aligned (3, 4) ``` 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/tensors/base_tensor/base_tensor_space.jl b/src/tensors/base_tensor/base_tensor_space.jl index 26fe1b4..e13a753 100644 --- a/src/tensors/base_tensor/base_tensor_space.jl +++ b/src/tensors/base_tensor/base_tensor_space.jl @@ -123,10 +123,10 @@ gram(space::AbstractTensorSpace) = gram(_u_coeffs(space.dg), cdc_components(spac # Examples ```jldoctest -julia> _is_orthonormal(function_space(dg)) +julia> DiffusionGeometryJ._is_orthonormal(function_space(dg)) true -julia> _is_orthonormal(vector_field_space(dg)) +julia> DiffusionGeometryJ._is_orthonormal(vector_field_space(dg)) false ``` """ @@ -149,7 +149,7 @@ but the weakest are two orders of magnitude down on the rest. Those are the ambi directions *normal* to the circle, which the metric barely sees: ```jldoctest -julia> vals, vecs = _gram_spectrum(vector_field_space(dg)); +julia> vals, vecs = DiffusionGeometryJ._gram_spectrum(vector_field_space(dg)); julia> length(vals) 16 diff --git a/src/tensors/base_tensor/metric_gram.jl b/src/tensors/base_tensor/metric_gram.jl index 277fa0c..7911469 100644 --- a/src/tensors/base_tensor/metric_gram.jl +++ b/src/tensors/base_tensor/metric_gram.jl @@ -49,15 +49,9 @@ returning `(batch..., n)`. `u_n1` is `(n, n1)`; `matrices` `(n, C, C)`. The kern behind [`metric_apply`](@ref) and hence [`g`](@ref) — the result is regularised before it is returned. -# Examples -```jldoctest -julia> X = grad(f); - -julia> vals = metric_apply(vector_field_space(dg), X.coeffs, X.coeffs); - -julia> size(vals) -(60,) -``` +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}) diff --git a/src/utils/basis_conversions.jl b/src/utils/basis_conversions.jl index faf3f10..5119dd0 100644 --- a/src/utils/basis_conversions.jl +++ b/src/utils/basis_conversions.jl @@ -26,7 +26,7 @@ coefficients 2 and 3 are non-zero — and by Parseval they carry all of its L² `∫cos²θ = 1/2`: ```jldoctest -julia> c = _from_pointwise_basis(reshape(cos.(θ), 60, 1), function_space(dg)); +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}: @@ -75,7 +75,7 @@ Evaluating the coefficients of `f` back at the 60 sample points recovers `cos θ to the error of the 8-function truncation: ```jldoctest -julia> values = _to_pointwise_basis(f.coeffs, function_space(dg); basis_count=8); +julia> values = DiffusionGeometryJ._to_pointwise_basis(f.coeffs, function_space(dg); basis_count=8); julia> size(values) (60,) diff --git a/src/utils/batch_utils.jl b/src/utils/batch_utils.jl index 41981b9..8104248 100644 --- a/src/utils/batch_utils.jl +++ b/src/utils/batch_utils.jl @@ -16,16 +16,16 @@ function rather than a delegation. Throws `DimensionMismatch` if the shapes conf # Examples ```jldoctest -julia> broadcast_batch_shape((3, 4), (4,)) # right-aligned: (4,) fills the last axis +julia> DiffusionGeometryJ.broadcast_batch_shape((3, 4), (4,)) # right-aligned: (4,) fills the last axis (3, 4) -julia> broadcast_batch_shape((3, 1), (1, 4)) # length-1 axes stretch +julia> DiffusionGeometryJ.broadcast_batch_shape((3, 1), (1, 4)) # length-1 axes stretch (3, 4) -julia> broadcast_batch_shape((), (5,)) # an unbatched tensor joins any batch +julia> DiffusionGeometryJ.broadcast_batch_shape((), (5,)) # an unbatched tensor joins any batch (5,) -julia> broadcast_batch_shape((3,), (4,)) +julia> DiffusionGeometryJ.broadcast_batch_shape((3,), (4,)) ERROR: DimensionMismatch: Incompatible batch shapes: ((3,), (4,)) ``` """ @@ -83,10 +83,10 @@ axes it already has. Julia's own broadcasting then reproduces numpy's semantics: # Examples ```jldoctest -julia> size(pad_batch_dims(zeros(4, 6), 2)) # batch (4,), coeff axis 6 +julia> size(DiffusionGeometryJ.pad_batch_dims(zeros(4, 6), 2)) # batch (4,), coeff axis 6 (1, 4, 6) -julia> size(pad_batch_dims(zeros(4, 6), 1)) # already rank 1: untouched +julia> size(DiffusionGeometryJ.pad_batch_dims(zeros(4, 6), 1)) # already rank 1: untouched (4, 6) ``` """ @@ -106,7 +106,7 @@ rule. Both arrays keep their `ntail` trailing (non-batch) axes. # Examples ```jldoctest -julia> a, b = align_batch_pair(zeros(3, 4, 6), zeros(4, 6)); +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)) @@ -129,12 +129,12 @@ einsum; prefer [`pad_batch_dims`](@ref) when plain broadcasting will do. # Examples ```jldoctest -julia> size(broadcast_batch_to(ones(4, 6), (3, 4))) # copied along the new axis +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> broadcast_batch_to(A, (4,)) === A # already the target: no copy +julia> DiffusionGeometryJ.broadcast_batch_to(A, (4,)) === A # already the target: no copy true ``` """ @@ -155,7 +155,7 @@ axis, giving `(prod(target), tail...)` — the layout the einsum kernels want. # Examples ```jldoctest -julia> size(broadcast_flatten_batch(ones(4, 6), (3, 4))) # batch (3, 4) → 12 rows +julia> size(DiffusionGeometryJ.broadcast_flatten_batch(ones(4, 6), (3, 4))) # batch (3, 4) → 12 rows (12, 6) ``` """ @@ -174,12 +174,12 @@ return the array together with its leading batch shape. This is the check every # Examples ```jldoctest -julia> _, batch = infer_batch_shape(zeros(3, 5, 16), (16,)); +julia> _, batch = DiffusionGeometryJ.infer_batch_shape(zeros(3, 5, 16), (16,)); julia> batch (3, 5) -julia> infer_batch_shape(zeros(3, 15), (16,); name="VectorField") +julia> DiffusionGeometryJ.infer_batch_shape(zeros(3, 15), (16,); name="VectorField") ERROR: AssertionError: VectorField coefficients must have trailing shape (16,), got (3, 15) ``` """ @@ -201,12 +201,12 @@ Undone by [`restore_batch_dims`](@ref). # Examples ```jldoctest -julia> flat, batch = flatten_batch_dims(zeros(2, 3, 16)); +julia> flat, batch = DiffusionGeometryJ.flatten_batch_dims(zeros(2, 3, 16)); julia> size(flat), batch ((6, 16), (2, 3)) -julia> flat, batch = flatten_batch_dims(zeros(16)); # unbatched: one row, empty batch +julia> flat, batch = DiffusionGeometryJ.flatten_batch_dims(zeros(16)); # unbatched: one row, empty batch julia> size(flat), batch ((1, 16), ()) @@ -230,12 +230,12 @@ Inverse of [`flatten_batch_dims`](@ref): reshape a `(prod(batch), L)` matrix bac ```jldoctest julia> A = reshape(1.0:96.0, 2, 3, 16); -julia> flat, batch = flatten_batch_dims(A); +julia> flat, batch = DiffusionGeometryJ.flatten_batch_dims(A); -julia> restore_batch_dims(flat, batch) == A # a round trip +julia> DiffusionGeometryJ.restore_batch_dims(flat, batch) == A # a round trip true -julia> size(restore_batch_dims(zeros(1, 16), ())) # empty batch → a plain vector +julia> size(DiffusionGeometryJ.restore_batch_dims(zeros(1, 16), ())) # empty batch → a plain vector (16,) ``` """ diff --git a/src/utils/reshape_utils.jl b/src/utils/reshape_utils.jl index 949d5f8..bd0a3d5 100644 --- a/src/utils/reshape_utils.jl +++ b/src/utils/reshape_utils.jl @@ -18,7 +18,7 @@ every coefficient vector in the package (`(i, a) ↦ i·C + a`, basis index slow Reading order is by rows, not columns — this is what `Base.reshape` would *not* do: ```jldoctest -julia> np_reshape([1 2 3; 4 5 6], 3, 2) +julia> DiffusionGeometryJ.np_reshape([1 2 3; 4 5 6], 3, 2) 3×2 Matrix{Int64}: 1 2 3 4 @@ -36,7 +36,7 @@ It round-trips, so a flattened coefficient vector can always be unpacked again: ```jldoctest julia> A = reshape(1:12, 2, 2, 3); -julia> np_reshape(np_reshape(A, 2, 6), 2, 2, 3) == A +julia> DiffusionGeometryJ.np_reshape(DiffusionGeometryJ.np_reshape(A, 2, 6), 2, 2, 3) == A true ``` """ diff --git a/test/doctests.jl b/test/doctests.jl index 69c033d..26e89e7 100644 --- a/test/doctests.jl +++ b/test/doctests.jl @@ -13,5 +13,8 @@ include(joinpath(@__DIR__, "..", "docs", "doctest_setup.jl")) DocMeta.setdocmeta!(DiffusionGeometryJ, :DocTestSetup, DOCTEST_SETUP; recursive=true) @testset "doctests" begin - doctest(DiffusionGeometryJ; manual=false) + # `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 From a9e4ae7a434e31677aed329aa0da6352d4e52616 Mon Sep 17 00:00:00 2001 From: Jeffrey Ward Date: Sat, 11 Jul 2026 23:01:41 -0400 Subject: [PATCH 3/3] docs: bring the README back in line with the source MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Audited every API name and code block in the README against the package by executing them, rather than reading. The prose was accurate but had fallen behind three merged features and carried one broken link. - The install URL pointed at `wardjm/DiffusionGeometryJ.jl`; the repository is `wardjm/DiffusionGeometryJ` (no `.jl`), so `Pkg.add` as documented would have failed. - Batching went unmentioned entirely — no `batch_shape`, and no note that batch axes broadcast numpy-style from the right rather than by Julia's rule. That is a real trap when porting, and it was only recorded in the Conventions page. - `.*` and `./` (#24) and the accessors were undocumented. A numeric array scaling a tensor means batch-wise scalars, which is what makes `exp.(-evals) .* evecs` a spectral filter; that now appears where a reader will look for it. - The `show` methods added in c3bbd84 are shown, and the coefficients-not-values distinction that motivates `to_pointwise_basis` is stated in both the tensor section and Conventions. The Documentation section now says the guide pages' examples run in the suite too (true as of the previous commit) and links the Conventions and upstream-bugs pages. Everything else — the constructor table, the operator and Betti sections, the degree conventions, `hodge_decomposition`'s two return shapes — was checked against a live geometry and is correct as written. --- README.md | 57 +++++++++++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 53 insertions(+), 4 deletions(-) diff --git a/README.md b/README.md index c58eebf..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 @@ -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 @@ -314,9 +357,15 @@ julia --project=docs docs/make.jl # writes docs/build, and runs every doc 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, 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`. +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