Metadata-Version: 2.4
Name: vsh-tools
Version: 1.0.0
Summary: Vector and scalar spherical harmonic analysis of astrometric catalogues.
Author-email: Neo <niu.liu@nju.edu.cn>
License-Expression: MIT
Project-URL: Homepage, https://git.nju.edu.cn/neo/vsh
Project-URL: Repository, https://git.nju.edu.cn/neo/vsh
Requires-Python: >=3.9
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: numpy>=1.22
Requires-Dist: scipy>=1.9
Requires-Dist: astropy>=5.0
Provides-Extra: plot
Requires-Dist: matplotlib>=3.7; extra == "plot"
Provides-Extra: dev
Requires-Dist: pytest>=7.4; extra == "dev"
Provides-Extra: release
Requires-Dist: build>=1.0; extra == "release"
Requires-Dist: twine>=5.0; extra == "release"
Dynamic: license-file

# VSH Tools: Vector Spherical Harmonic Analysis

Real vector spherical harmonics (VSH) and scalar spherical harmonics (SSH)
for astrometric catalogue analysis. Dense fits support full observation
covariance, including cross-source correlations. Chunked fits use QR/SVD
with independent sources and optional within-source component correlations,
avoiding a full design matrix or dense observation covariance.

Reference: Mignard & Klioner (2012), *Astronomy & Astrophysics* **547**, A59,
[doi:10.1051/0004-6361/201219927](https://doi.org/10.1051/0004-6361/201219927).
Source repository: [git.nju.edu.cn/neo/vsh](https://git.nju.edu.cn/neo/vsh).

## Local installation

From the repository root, activate your virtual environment and install:

```bash
source ~/.venvs/astro/bin/activate
bash build_code.sh local_install
```

The script performs an editable installation. The equivalent command is
`python -m pip install --editable .`. To include the test dependencies:

```bash
python -m pip install --editable '.[dev]'
```

These examples describe the refactored code in this checkout, not necessarily
the version currently published on PyPI.

Building source and wheel distributions requires the release tools:

```bash
python -m pip install --editable '.[dev,release]'
bash build_code.sh build
```

This command builds and checks both archives without uploading them. See
[the release procedure](docs/releasing.md) for versioning, tagging and publishing.

## Numerical conventions

- Positions `ra`, `dec` are in radians. Convert degree coordinates explicitly
  with `np.deg2rad` before calling the new interfaces.
- `dra` is already the tangent component `delta_RA * cos(dec)`. Its error must
  describe that same component. Do not multiply either by `cos(dec)` again.
- Coordinates and observations must be finite. Filter missing data explicitly,
  applying the same selection to coordinates, observations, and covariance.
- All observation components have consistent units, for example microarcseconds
  or microarcseconds per year. Coefficients retain those units; covariance has
  their squared units. NumPy float64 arrays and SciPy numerical routines are used.
- A VSH observation vector is ordered as **all RA components, then all Dec
  components**. Its covariance has shape `(2*N, 2*N)` in that same order. SSH
  uses `(N,)` observations and `(N,N)` covariance.
- Full covariance, including cross-source correlations, is supported by the
  dense method. Both methods also accept standard-error columns. Omitting
  both covariance and errors assumes known unit variance; it does not estimate
  a noise scale.
- The real expansion uses the `m=0` column and, for `m>0`, `2*Re` and `-2*Im`
  columns. VSH names begin `t10, s10, t11r, s11r, t11i, s11i` in a full fit.
  SSH begins `y00, y10, y11r, y11i`. Use the naming functions to obtain the
  complete order. VSH starts at degree 1; SSH includes degree 0.
- Coefficient covariance is **not** multiplied by reduced chi-square. A
  rank-deficient fit raises an exception instead of returning a minimum-norm
  solution with misleading uncertainties.

Dense covariance alone occupies about 0.8 GB at 5,000 sources and 3.2 GB at
10,000 sources for VSH. Factorization and other work arrays require additional
memory; the entire design matrix and thin SVD also reside in memory.

## A VSH fit

The following example generates a synthetic vector field. It is runnable in a
Python session or notebook; no notebook file is required.

```python
import numpy as np
from vsh_tools import (
    covariance_from_errors, fit_vsh, predict_vsh, vsh_coefficient_names,
)
from vsh_tools.statistics import fit_statistics
from vsh_tools.parameters import physical_from_fit, vector_apex
from vsh_tools.power import power_from_fit, significance_from_fit

rng = np.random.default_rng(42)
n = 120
ra = rng.uniform(0, 2*np.pi, n)
dec = np.arcsin(rng.uniform(-1, 1, n))
l_max = 2
names = vsh_coefficient_names(l_max)
true_coefficients = rng.normal(0, 20, len(names))

sigma_ra = np.full(n, 10.0)
sigma_dec = np.full(n, 12.0)
covariance = covariance_from_errors(sigma_ra, sigma_dec, rho=0.2)
noise = np.linalg.cholesky(covariance) @ rng.standard_normal(2*n)
observed = predict_vsh(ra, dec, true_coefficients, l_max) + noise.reshape(2, n)
dra, ddec = observed

fit = fit_vsh(ra, dec, dra, ddec, l_max, covariance=covariance)
sigma = np.sqrt(np.diag(covariance)).reshape(2, n)
stats = fit_statistics(fit, sigma=sigma)
physical = physical_from_fit(fit)
power = power_from_fit(fit)
significance = significance_from_fit(fit)

print(dict(zip(fit.coefficient_names, fit.coefficients)))
print(stats.chi2, stats.dof, stats.reduced_chi2)
print(dict(zip(physical.names, physical.values)))
print(power.degrees, power.total)
print(significance.components["TS"].p_value)

# For a full VSH fit, the first three physical parameters are G1, G2, G3.
glide_apex = vector_apex(physical.values[:3], physical.covariance[:3, :3], degrees=True)
print(glide_apex.amplitude, glide_apex.ra, glide_apex.dec)
```

`covariance_from_errors` constructs a dense matrix with independent sources
and optional within-source RA-Dec correlations. If a catalogue supplies a
full cross-source covariance, pass that matrix directly to the fit.

`fit.residuals` has shape `(2,N)` and contains observed minus computed values.
`fit.covariance` is the coefficient covariance, not the observation covariance.
For a pure T or S model, pass `fit_type="T"` or `"S"`; physical conversion
returns only the corresponding rotation or glide family and quadrupole terms.

## Chunked fitting

For independent sources, supply the errors directly. The following fit is
equivalent to the dense example above, including its within-source correlation:

```python
chunked = fit_vsh(ra, dec, dra, ddec, l_max, method="chunked",
                  sigma_ra=sigma_ra, sigma_dec=sigma_dec, rho=0.2, block_size=50)
```

`block_size` counts sources. A final partial block is supported. Full
cross-source covariance requires `method="dense"`; it cannot be replaced by
independent blocks without changing the statistical model. For SSH, pass
`sigma` instead of the two vector error arrays. With `store_residuals=False`,
coefficients, covariance and chi-square remain available, but residual arrays
and descriptive residual statistics are omitted.

Chunked fitting needs approximately `O(block_size*p + p**2)` working memory,
in addition to input columns and any stored residuals, where `p` is the number
of coefficients. Runtime still grows as `O(N*p**2)`; blocking alone does not
make very high-degree or Gaia-sized analyses inexpensive.

## Iterative clipping and selected covariance

Each iteration fits the current selection, computes residuals for **all original
sources**, and classifies all sources again. Excluded sources may re-enter.
Convergence requires identical masks, not merely the same number of sources.

Continuing from the arrays above:

```python
from vsh_tools import fit_vsh_clipped

clipped = fit_vsh_clipped(
    ra, dec, dra, ddec, l_max,
    covariance=covariance, threshold=4.5, mode="absolute", max_iter=20,
)
if not clipped.converged:
    raise RuntimeError("Clipping did not converge; inspect clipped.history.")

fit = clipped.fit
mask = clipped.mask
selected = np.flatnonzero(mask)
rows = np.r_[selected, n + selected]
selected_covariance = covariance[np.ix_(rows, rows)]
selected_sigma = np.sqrt(np.diag(selected_covariance)).reshape(2, selected.size)
stats = fit_statistics(fit, sigma=selected_sigma)
```

The threshold here is illustrative, not an automatically selected cutoff.
`mode="median"` uses `threshold * median(scores)` over all original sources
at each iteration. VSH scores are two-component marginal Mahalanobis distances;
SSH scores are absolute residuals divided by their marginal standard errors.

`clipped.residuals` covers all original sources; `clipped.fit.residuals` covers
only the selected sources. At the iteration limit the returned fit still
matches its returned mask, and `converged` is false.

For SSH selection, use `covariance[np.ix_(selected, selected)]` and a one-dimensional
sigma array. When using error columns, subset those columns and `rho` with the
same source mask. Unit-weight fits need no covariance submatrix; use unit
sigmas if normalized summaries are wanted.

Clipping also accepts `method="chunked"` and the source-error arguments of
`fit_vsh`/`fit_ssh`. Its median cutoff is computed over all original sources,
not separately within each block. `store_scores=False` omits the returned
score array; `store_residuals=False` omits both full and selected residuals.

## Bootstrap uncertainties

Using the selected data from the clipping example:

```python
from vsh_tools.bootstrap import bootstrap_vsh, physical_from_bootstrap

boot = bootstrap_vsh(
    ra[mask], dec[mask], dra[mask], ddec[mask], l_max,
    covariance=selected_covariance,
    method="parametric", n_resamples=1000, rng=42,
)
summary = boot.summary(confidence_level=0.95)
print(summary.errors)
print(summary.interval)  # Shape (2,p): lower and upper percentile bounds.
physical_bootstrap = physical_from_bootstrap(boot)
```

- `pairs` (default) resamples whole source records. It requires independent
  sources, retains within-source covariance, and rejects any cross-source
  covariance. A rank-deficient resample raises with its draw number; failed
  draws are not silently dropped or replaced.
- `parametric` holds positions fixed and uses the specified Gaussian covariance,
  including cross-source correlations. For a linear fit, direct draws from
  `N(fit.coefficients, fit.covariance)` have exactly the same coefficient
  distribution as simulating and refitting catalogues. This shortcut avoids
  repeated dense factorizations. It reproduces the assumed noise model rather
  than providing an independent calibration of that model.

`boot.samples` stores all coefficient draws. Bootstrap covariance uses denominator
`B-1`; `summary.rms` measures RMS about the original estimate with denominator
`B`. Original formal errors are not overwritten. Percentile intervals are
marginal, not simultaneous or BCa intervals. Clipping is not repeated within
bootstrap: using a selected catalogue conditions on that selection.

## SSH fitting

```python
from vsh_tools import fit_ssh, predict_ssh
from vsh_tools.bootstrap import bootstrap_ssh

# Reuse the sky positions and RNG from the first example.
values = 5.0 + rng.normal(0, 1, n)
scalar_fit = fit_ssh(ra, dec, values, l_max=0)
constant = scalar_fit.coefficients[0] / np.sqrt(4*np.pi)
scalar_model = predict_ssh(ra, dec, scalar_fit.coefficients, l_max=0)
scalar_power = power_from_fit(scalar_fit)
scalar_significance = significance_from_fit(scalar_fit)
scalar_boot = bootstrap_ssh(ra, dec, values, 0, n_resamples=1000, rng=42)
print(constant, scalar_power.components["Y"])
```

`fit_ssh_clipped` provides the corresponding iterative scalar workflow.
SSH uses `(l_max+1)**2` coefficients and one observation per source.

## Catalogue input and unified analysis

`vsh_catalog_from_table` and `ssh_catalog_from_table` in `vsh_tools.catalog`
adapt Astropy `Table` and `QTable` objects. Column unit metadata takes priority;
positions are converted to radians and observations/errors to a common unit.
The vector RA component and its error must already include `cos(dec)`.

Error columns are retained in `catalog.errors` with optional
`catalog.correlation`; `catalog.covariance` remains `None`. An explicitly
supplied full covariance is retained instead. `catalog.select(mask)` subsets
the weights and preserves `row_indices` and `source_ids`. Invalid input raises
by default; `invalid="drop"` explicitly removes invalid rows.

`analyze_vsh` and `analyze_ssh` accept a table, a table file path or a prepared
catalogue and assemble fitting, clipping and post-fit results:

```python
from astropy import units as u
from astropy.table import Table
from vsh_tools import analyze_vsh

table = Table()
table["ra"], table["dec"] = ra * u.rad, dec * u.rad
table["dra"], table["ddec"] = dra * u.uas, ddec * u.uas
table["sx"], table["sy"] = sigma_ra * u.uas, sigma_dec * u.uas
analysis = analyze_vsh(table, l_max, method="chunked", block_size=50,
                       catalog_options={"errors": ("sx", "sy"), "correlation": 0.2})
print(analysis.physical.values, analysis.reduced_chi2)
```

For optional clipping, pass `clipping=True` and an explicit `threshold`.
`analysis.mask` indexes the prepared valid catalogue; `analysis.original_mask`
indexes the original input table. `analysis.write("result.npz")` saves arrays
and JSON metadata, refuses to overwrite existing files by default, and can be
read using `np.load(..., allow_pickle=False)`. Input observations are not
duplicated in the archive. With `store_residuals=False`, `analysis.statistics`
is `None`; global chi-square, power and significance remain available.

## Power and significance

Power is the integral of the squared field over the sphere, without a `1/(4*pi)`
normalization. Each `m=0` coefficient has weight 1 and every positive-m real or
imaginary coefficient has weight 2. `sqrt(power/(4*pi))` is the spherical RMS.
Power includes noise bias; it is not automatically debiased.

`significance_from_fit` defaults to the full-covariance Wald statistic.
Single-family VSH tests have `2*l+1` degrees of freedom; joint `TS` tests have
`2*(2*l+1)`. SSH includes its one-dof monopole. The paper's approximations are
available as `method="zonal"`, `"diagonal"`, and `"mean"`. Prefer returned
`p_value` to the approximate Wilson-Hilferty `z` score.

Chi-square reference distributions require a fixed, identifiable linear model
with Gaussian noise and known covariance. They do not account for clipping,
adaptive degree selection, or multiple testing. Estimating an overall noise
scale from residuals changes the appropriate reference test.

`sampled_field_power` estimates field power from sample values, with optional
sky-area weights. An unweighted nonuniform catalogue need not reproduce the
continuous harmonic integral. Evaluate actual O-C values to obtain sampled
residual power; subtracting harmonic power from sampled data power is not a
general residual-power identity.

## Refactored modules

| Module | Responsibility |
| --- | --- |
| `basis.py` | Real SSH/VSH basis arrays and coefficient names |
| `covariance.py` | Dense covariance and direct source-pair whitening |
| `solver.py` | Dense SVD and chunked QR/SVD least squares |
| `model.py` | Predictions and O-C residuals |
| `fitting.py` | Fixed-sample SSH/VSH fits and `HarmonicFitResult` |
| `clipping.py` | Iterative source classification with re-entry |
| `statistics.py` | Residual summaries and fit diagnostics |
| `parameters.py` | Physical VSH parameters, covariance propagation, vector apexes |
| `power.py` | Degree-wise power and significance |
| `bootstrap.py` | Pairs/parametric resampling and uncertainty summaries |
| `catalog.py` | Table/QTable units, weights, filtering and source identities |
| `analysis.py` | Unified catalogue analysis and optional NPZ output |
| `legacy/` | Historical interfaces and algorithms |

Common basis, covariance, prediction, and fitting interfaces are also available
directly from `vsh_tools`; the unified analysis exports are loaded lazily.
Import post-fit analyses from their respective modules, as above. Legacy
modules live in `vsh_tools.legacy`; the new numerical path
does not call the old normal-equation solver or the old `vsh_fit` workflow.

## Examples and validation

Run the complete synthetic VLBI example, including cross-source covariance,
outlier clipping, statistics, physical parameters, power, and bootstrap:

```bash
python examples/vlbi_workflow.py
```

Run all module, legacy and integration tests from the repository root after
installing the development dependencies:

```bash
python -m pytest -q
```
