Metadata-Version: 2.4
Name: vitali
Version: 0.1.0
Classifier: License :: OSI Approved :: Apache Software License
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Rust
Classifier: Topic :: Scientific/Engineering :: Mathematics
Classifier: Intended Audience :: Science/Research
Requires-Dist: numpy>=1.24
License-File: LICENSE
Summary: Multiscale geometric analysis of point clouds: geometric multi-resolution analysis and diffusion wavelets (pure Rust, Python bindings)
Keywords: multiscale,manifold-learning,wavelets,diffusion-wavelets,gmra,dimension-estimation
Author-email: Alejandro Soto Franco <sotofranco.eng@gmail.com>
License: Apache-2.0
Requires-Python: >=3.10
Description-Content-Type: text/markdown; charset=UTF-8; variant=GFM
Project-URL: Documentation, https://github.com/alejandro-soto-franco/vitali#readme
Project-URL: Repository, https://github.com/alejandro-soto-franco/vitali

# vitali

Multiscale geometric analysis of point clouds, in Rust, with Python bindings and a
MATLAB bridge.

Two constructions live here, both from Mauro Maggioni and coauthors, and neither
with a maintained implementation in any package ecosystem before this crate.

**Geometric multi-resolution analysis** builds a multiscale, piecewise-affine
approximation of a point cloud: a tree partition of the data, a local PCA scaling
basis and centre at every node, and a geometric wavelet basis relating each node to
its parent. Allard, Chen and Maggioni, *Multi-scale geometric methods for data sets
II: geometric multi-resolution analysis*, Applied and Computational Harmonic
Analysis, 2012; Chen and Maggioni, *Multiscale geometric wavelets for the analysis
of point clouds*, CISS, 2010.

**Diffusion wavelets** build an orthonormal multiresolution directly from a
diffusion operator, rather than from translates and dilates of a fixed mother
wavelet: rank-revealing Gram-Schmidt on the operator's columns gives a scaling basis,
the operator is compressed onto it, the compressed operator is squared to represent
the next dyadic power, and the process repeats. Coifman and Maggioni, *Diffusion
wavelets*, Applied and Computational Harmonic Analysis, 2006.

The crate is named for the Vitali covering lemma. The greedy net that builds every
scale of the GMRA partition is a Vitali covering argument.

## Prior art in August 2026

Geometric multi-resolution analysis has three GitHub repositories, none above one
star, one a 2015 Python toolkit and one an R package last pushed in February 2021.
Nothing on PyPI, in the Julia General registry, or on crates.io.

Diffusion wavelets have nothing in any ecosystem. `Wavelets.jl`,
`ContinuousWavelets.jl` and PyWavelets are classical wavelets on regular grids. The
construction that builds an orthonormal multiresolution from a diffusion operator
exists in the 2007 MATLAB sources and nowhere else.

## Install

```bash
cargo add vitali     # Rust
pip install vitali   # Python
```

## Rust

```rust
use vitali::gmra;

let points = gmra::datasets::d_sphere(4000, 4, 100, 0);
let tree = gmra::build(&points, &gmra::Options::default());

println!("{} nodes, depth {}", tree.nodes.len(), tree.depth());

let scales = gmra::reconstruct_all_scales(&tree, &points);
for (j, approx) in scales.iter().enumerate() {
    let err = (approx - &points).norm() / points.norm();
    println!("scale {}: relative error {:.3e}", j + 1, err);
}
```

Diffusion wavelets take a sparse operator:

```rust
use vitali::diffusion_wavelets::{Options, build_tree, circle_diffusion};

let t = circle_diffusion(256);
let tree = build_tree(&t, &Options::default());

for level in &tree.scaling {
    println!("level {}: dimension {}", level.level, level.ext_basis.ncols());
}
```

Point clouds are `D x N` in Rust, one column per point, following the reference
MATLAB distribution.

## Python

Point clouds are `(n_points, n_features)` in Python, one row per point, following
scikit-learn. The boundary transposes.

```python
import numpy as np
import vitali

X = vitali.datasets.d_sphere(4000, d=4, n_features=100, seed=0)
tree = vitali.gmra(X, min_leaf_size=50)

tree.median_scaling_dimension()   # 4.0, the manifold's true dimension
tree.depth
tree.num_nodes

scales = tree.reconstruct_all_scales(X)
for j, approx in enumerate(scales, start=1):
    err = np.linalg.norm(approx - X) / np.linalg.norm(X)
    print(f"scale {j}: relative error {err:.3e}")
```

A node exposes its own geometry:

```python
node = tree.node(tree.leaves()[0])
node.centre          # (n_features,)
node.scaling_basis   # (local_dim, n_features), one direction per row
node.sigmas          # every singular value, before truncation
node.wavelet_basis   # the basis relating this node to its parent
node.radius
```

Diffusion wavelets take a dense array or any scipy sparse matrix, which is converted
to CSR without ever being densified:

```python
import scipy.sparse as sp
import vitali

T = sp.csr_matrix(vitali.circle_diffusion(256))
tree = vitali.diffusion_wavelets(T, precision=1e-10)

tree.scaling_dimensions()          # dimension at every level
tree.scaling()[0].ext_basis        # Phi_1 in the original delta basis
tree.wavelet()[1].ext_basis        # Psi_2 in the original delta basis
```

## MATLAB

MATLAB reaches vitali through its Python interface, so the Rust runs underneath
without a MEX build. Point `pyenv` at an interpreter your MATLAB release supports,
once per machine:

```matlab
pyenv(Version = "/path/to/python", ExecutionMode = "OutOfProcess");
```

Out-of-process execution keeps the extension's libraries away from the ones MATLAB
loads itself, and is the sensible default.

Two conversions do all the work, and both are in `matlab/` in this repository.
`tonumpy` respects MATLAB's column-major order, which is where a hand-written
conversion usually goes wrong:

```matlab
addpath("matlab");

X = rand(4000, 100);              % n_points by n_features
tree = py.vitali.gmra(tonumpy(X), pyargs("min_leaf_size", int64(50)));

k = double(tree.median_scaling_dimension())
d = double(tree.depth)

scales = cell(tree.reconstruct_all_scales(tonumpy(X)));
approx = fromnumpy(scales{end});   % back to a MATLAB double matrix
err = norm(approx - X, "fro") / norm(X, "fro")
```

Three points decide whether a call works:

- **Keyword arguments go through `pyargs`.** Every option is keyword-only.
- **Integers need an explicit cast.** MATLAB numeric literals are doubles, and
  `min_leaf_size`, `max_dim`, `levels`, `stop_fcns`, `manifold_dimension` and every
  `seed` are integers. Pass `int64(50)`, never `50`.
- **Orientation is transposed against the MATLAB toolboxes.** `GMRA.m` and
  `DWPTree.m` take `D x N`, one column per point. vitali's Python surface takes
  `(n_points, n_features)`. Data prepared by `GenerateData_and_SetParameters.m`
  passes as `tonumpy(X.')`.

A diffusion operator goes across as a dense matrix:

```matlab
T = fromnumpy(py.vitali.circle_diffusion(int64(256)));   % or your own operator
tree = py.vitali.diffusion_wavelets(tonumpy(T), pyargs("precision", 1e-10));

levels = cell(tree.scaling());
Phi1 = fromnumpy(levels{1}.ext_basis);
```

`matlab/vitali_demo.m` runs the whole sequence end to end and prints the recovered
dimension and the per-scale reconstruction error.

## Differences from the reference distribution

**Partitioning.** The reference distribution partitions with METIS-based nested
dissection or a cover tree, both reached only through precompiled MEX binaries with
no source in the distribution: `nesdis.mexa64` and `covertree.mexa64` do not load
under GNU Octave (`'covertree' undefined near line 4`, checked directly). This crate
recursively splits each node's points at the median of their projection onto the
node's top local PCA direction, which is a principal-direction divisive partition
tree. Node
for node, reconstruction error and dimension counts here are not directly comparable
against a run using either reference partitioner. The comparisons below are against
independently generated data of known dimension, and against the reference
distribution's own reported figures as a target range.

**Reconstruction.** The reference distribution's inverse transform accumulates
wavelet detail along the root-to-leaf path, corrected at every step by a
tangential-correction mechanism (`addTangentialCorrections`, its own default) that
keeps the accumulated sum aligned with each node's own direct local PCA
reconstruction. Without that correction the accumulated sum drifts, because two
nodes' independently computed local PCA planes are not nested in general, unlike
diffusion wavelets' scaling spaces: an early build of the wavelet-sum reconstruction
here measured 26x worse error at the finest scale than the same node's own direct
reconstruction (0.66 against 0.025 relative error) on curved test data, worsening
with depth. Reproducing that bookkeeping (`ScalBasisChange`, `ScalWavBasisChange`,
`ScalCenterConsts`) is out of scope. `reconstruct_all_scales` instead reconstructs
each scale directly from that scale's own node, `centre + Phi' Phi (x - centre)`,
which is the value the tangential correction exists to make the accumulated sum
equal. Wavelet bases and translation constants are still constructed and available
on every node, for compression-oriented use.

**Precision.** This crate is `f64` throughout. See the D-Cube result below.

**`op_threshold`.** `DWPTree.m`'s docstring states that `OpThreshold` defaults to 0;
its shipped code defaults it to `Precision` when left unset. `Options::op_threshold`
follows the shipped code as an explicit field with no implicit default, so a caller
wanting the docstring's behaviour passes `0.0` and a caller wanting the shipped
behaviour passes `precision` itself.

## Validation

35 tests. Diffusion wavelets are checked against the four structural identities the
construction is stated by, and against MATLAB-exported fixtures; GMRA is checked
against manifolds of known intrinsic dimension.

### Structural identities of the diffusion wavelets

On the circle at `N = 256`:

| Identity | Residual |
|---|---|
| Orthonormality, `Phi_j' Phi_j = I` | 2e-15 |
| Nesting, `span Phi_j` inside `span Phi_{j-1}` | 3e-15 |
| Wavelet-scaling complement | 6e-16 |
| Level-`j` operator equals the compression of `T^(2^(j-1))` | 5e-15 |

### Dimension recovery under GMRA

Median scaling dimension across the whole tree, against each manifold's true
dimension, fully adaptive (no fixed `manifold_dimension`, `internal_threshold = 0.5`,
`leaf_threshold = 1e-2`, both reference defaults), `N = 4000`:

| Manifold | True dimension | Median recovered |
|---|---|---|
| D-Sphere | 4 | 4 |
| D-Cube | 6 | 5 |
| SwissRoll | 2 | 2 |
| S-Manifold | 2 | 2 |

### Reconstruction error per scale under GMRA

The reference distribution's own double-precision run on D-Cube reconstructs to
machine precision at every scale with no rise anywhere. Its single-precision run
shows a non-monotone error rising to about 0.04 near the deepest scale. With
`manifold_dimension` fixed at the true dimension (matching the reference
distribution's own driver for these synthetic cases) and a leaf size comfortably
above the dimension being estimated locally (`min_leaf_size = 50` for 6-dimensional
data), D-Cube at `N = 4000`, `D = 50` gives:

```
[3.55e-16, 3.22e-16, 2.93e-16, 2.99e-16, 2.73e-16, 2.45e-16, 2.26e-16, 2.01e-16]
```

Machine precision at every one of the eight scales, and no rise between consecutive
scales above 1e-12, which is the threshold the test asserts. Variation between
neighbouring scales sits at the 1e-18 level, which is rounding rather than a trend.
That reproduces the reference
distribution's double-precision finding independently and through a different
partitioner, and supports reading the reported single-precision spike as an artefact
of that floating-point path rather than a property of the construction.

A smaller leaf size (`min_leaf_size = 20`, twenty points for a 6-dimensional local
PCA) reproduces a related, smaller artefact: one leaf under-samples and keeps rank 5
rather than 6, leaving every dimension-fixed internal scale at machine precision and
the finest scale at 0.37% relative error. Both runs are kept as tests. The effect is
small, explained, and specific to under-sampled local PCA.

Curved data reconstructs approximately, with error decreasing by scale and settling
above machine precision at the finest one. `N = 2000`, `D = 50`, `min_leaf_size = 20`,
`manifold_dimension` fixed at the true dimension, seed 1:

```
D-Sphere (d = 4): [0.431, 0.246, 0.213, 0.172, 0.126, 0.074, 0.060, 0.025]
SwissRoll (d = 2): [0.390, 0.240, 0.179, 0.151, 0.102, 0.067, 0.035, 0.004]
```

`cargo test -p vitali --release --test reconstruction -- --nocapture` prints all of
these, and `crates/vitali/tests/` states the assertion behind each one.

## Layout

```
crates/vitali/      the library, published to crates.io
crates/vitali-py/   PyO3 bindings, published to PyPI as vitali
matlab/             tonumpy, fromnumpy and a demo script
```

## Licence

Apache-2.0.

