[parashardhapola]

Summary

Two figures that users look for first: the samples on their principal components, as
DESeq2's plotPCA computes and draws them, and the dot plot of genes by cell type
(the fraction of cells that express each gene and their mean expression), as scanpy
draws it. principal_components returns the coordinates, shares of variance and
loadings that plot_pca draws, so that a guide can use them; plot_dotplot draws from
the sums a cube keeps, without the cells, and gives scanpy's values on the cells the cube
was pooled from. Separately, the incomplete beta function, which hung or ran 10^7 steps
of its continued fraction at arguments near the edges of the double range, now takes
well under a millisecond at every argument, with values checked against mpmath; on the
existing validation grid every result is unchanged to the bit but one, which moved to
the exact value.

Changes

• cytepool.principal_components.principal_components(data, *, group=None, top=500, min_cells=None):
  the samples on the principal components of their top genes of largest variance, for
  a samples x genes matrix, a VoomResult or one group of a Cube (its samples with at
  least min_cells cells and some counts, the genes filter_by_expression keeps,
  log-CPM with TMM factors and an average count of 2, as plot_mds prepares a cube).
• cytepool.principal_components.PrincipalComponents: coordinates (samples x
  components), variance_explained (DESeq2's percentVar), loadings (genes used x
  components), genes (positions of the genes used, by decreasing variance), samples
  (a cube's sample table rows, or ids from 1), settings, n_components; read-only.
• cytepool.plotting.plot_pca(data, *, group, color_by, marker_by, labels, top, components=(1, 2), min_cells, ax):
  plotPCA's figure, from a PrincipalComponents or from data: samples coloured and
  marked by sample variables as plot_mds draws them, axes labelled "PC1: 45%
  variance", equal scales.
• cytepool.plotting.plot_dotplot(cube, genes, *, groups, samples, standard_scale, mean_only_expressed, ax):
  scanpy's dot plot from a cube's detected, log_sum and n_cells, pooled over all
  samples or some (ids or a mask), for genes by feature id or symbol, or in named groups
  drawn under brackets; standard_scale="gene" or "group"; mean_only_expressed.
  The colour bar's label is written above the bar, as scanpy writes it, so that it
  never runs into the key of the dots' areas below the bar; a test checks that they
  do not overlap.
• The incomplete beta function (cytepool/stats/_scalar_special_functions.py): Temme's
  uniform expansion near the center of two shapes from 1e10; forms of the continued
  fraction and its prefactor that stay in range for shapes from 2^500, for a + b beyond
  the double range, for shapes below 1e-280, and for shapes more than 1e17 apart; exact
  products of shapes from 1e290; the small-shape series of 1 - I divided by the shape
  where its exponent would be absorbed. deviance_term halves arguments whose sum
  overflows, which also stops the incomplete gamma function from hanging there. Normal
  arguments take the old code path.
• Validation: validation/pca_and_dot_plots.py (new); validation/distributions.py
  gains the incomplete beta function at the edges of the double range and near the
  center of large shapes against mpmath, with the time per value.
• Tests: tests/test_pca_and_dot_plots.py (new) with fixtures in
  tests/data/pca_and_dot_plots/ (49 KB); tests/test_special.py gains the extreme
  arguments (fixture tests/data/distributions/betainc_extreme.npz, every point), a
  time bound, closed forms at huge and tiny shapes, and the incomplete gamma function
  where shape and argument overflow; the masked-input inventory gains the three
  functions (for arrays and cubes) and the record.
• Documentation: an API page for cytepool.principal_components in the API index; two
  capability matrix rows; a section "Samples and cells at a glance" in "The cube and its
  additive layers"; Pearson (1901) and Temme (1982) in references.bib (DOIs checked
  against Crossref) and on the Bibliography page; plot_pca and plot_dotplot under
  [planned] in the coverage exemptions, with the guides that will show them; the
  Validation page (a section and summary row for the new script, the distributions
  section, totals, dates, timings, the command, and a note that scanpy's dot plot
  conventions were read from its BSD-licensed source); a README status row.

Design notes

Principal components, as DESeq2's plotPCA

The help page of plotPCA (DESeq2 1.46.0) says: the ntop genes "selected by highest
row variance", returnData for "the data.frame of PC1 and PC2 with intgroup
covariates", and pcsToUse. What it leaves open was found by running it as a black
box on DESeqTransform objects of chosen values:

• The coordinates are those of R's prcomp on the centred top genes (U S of the
  singular value decomposition), and percentVar is each squared singular value over
  their sum. prcomp's help page says the signs are arbitrary and may differ between
  builds of R. CytePool makes each component's largest absolute coordinate positive (of
  near ties, within a relative 1e-9, the first), as plot_mds orients its dimensions;
  the validation compares up to sign (2 of 10 components differ from R's).
• Genes of equal variance keep their order: of four tied genes at ranks 19 to 22 with
  ntop = 20, the first two are kept. ntop beyond the number of genes takes them all.
• The axes read "PC1: 55% variance" (the share rounded to a whole percent), the figure
  has equal scales (coord_fixed, ratio 1), and two variables of intgroup colour the
  samples by their interaction ("a:x").
• Two samples: DESeq2 returns a second component of rounding noise (share 9.5e-31) and
  draws it. CytePool returns min(n_samples - 1, n_genes) components, the dimensions the
  centred values have, and refuses to draw a second one (a deliberate difference). It
  also refuses values that are not finite, where DESeq2 ranks a gene with NA last, and
  values that vary only by rounding.

A cube's group is prepared as plot_mds prepares it, as the edgeR User's Guide
prepares counts for plotMDS; the validation has R prepare the same samples with
edgeR (filterByExpr, normLibSizes, cpm(log = TRUE, prior.count = 2)) before
plotPCA. plot_pca draws its points with a private _draw_samples, which repeats
plot_mds's drawing so that this change touches only new functions of
plotting.py; plot_mds can move onto it later.

Dot plots, as scanpy's dotplot

scanpy's documentation (1.12.4) gives the fraction of cells whose value exceeds
expression_cutoff (0) and the mean in each group, standard_scale ("subtract the
minimum and divide each by its maximum", per variable or per group),
mean_only_expressed, size_exponent 1.5, largest_dot 200 and var_names as a dict
for groups of genes. scanpy is BSD-licensed, and three conventions its documentation
leaves open were read from its source: the largest fraction shown is rounded up to a
tenth, the key of the areas steps by 0.2, or by 0.1 and 0.05 when the largest fraction
is at most 0.6 and 0.3, and a scaled row of equal means is 0. No code was copied; a
notice in THIRD_PARTY_NOTICES.md may still be wanted (see below).

From a cube: a cell counts as expressing a gene when its count is positive (detected),
which is scanpy's default cutoff on log-normalized values, since log(1 + x) > 0 exactly
when the count is; the fraction and mean over the cells of all samples pooled are sums
over samples divided by summed cell numbers, which equals scanpy's on the cells. The mean
over expressing cells, log_sum / detected, needs a transform that maps a zero count to
zero, so it is refused unless the cube's normalization records a pseudocount of 1, as
from_cells records it. Options that need the cells (another expression_cutoff, a
dendrogram of the cells) are not offered; styling options are left to the Axes, and the
colours follow CytePool's palette (blue from light to dark) rather than scanpy's reds.
The figure exposes its values through its artists: the dots' sizes are the areas and
their array the values the colours show.

The incomplete beta function

On a grid of 0, 5e-324, 1e-300, 1e-8, 0.5, 1 - 2^-53, 1, 2, 30, 1e8, 1e300 and the
largest double for a, b and x (1,728 points), the old code did not return at 32 points
within 5 s, took over 1 ms at 144 (up to 3.4 s), and gave NaN at 158 points of the
domain. The causes:

• Loader's deviance term looped forever when x + m overflowed (a + b, or a + x for the
  incomplete gamma function, beyond 1.8e308): the series' step became NaN and its stop
  test, an update that changes nothing, never held. It now halves its arguments
  (bd0 is homogeneous of degree one) where their sum reaches 2^1023.
• The continued fraction overflowed for shapes beyond about 1e154 (a (1 + lam), lam a,
  (a + b) x and k (b - k) x), and with denominators beyond 2^1022 the reciprocal that
  Lentz's method keeps was subnormal and inexact, so the ratio of steps never reached 1
  within the tolerance: 10^7 steps (0.15 s, longer on subnormals), then NaN. For shapes
  from 2^500 it now runs on products ordered to stay in range and keeps the denominators
  themselves instead of their reciprocals; normal shapes run the old routine.
• Near the center of two large shapes the fraction needs about the cube root of the
  shapes in steps (counted: 2,500 at 1e8, 11,600 at 1e10, 55,000 at 1e12, 1.1 million at
  1e16), beyond the 10^7 cap from about 1e20. Within one standard deviation of the
  center (r eta^2 / 2 <= 1/2) of two shapes from 1e10, Temme's (1982) uniform expansion
  replaces it: I = erfc(-eta sqrt(r/2)) / 2 - G exp(-r eta^2 / 2) c0(eta) / sqrt(2 pi r),
  with c0's Taylor series (p - q) / (3 s) + (1 - s^2) / (12 s^2) eta + (p - q)(2 + s^2) /
  (135 s^3) eta^2, derived by series reversion and checked against mpmath; the omitted
  terms are of relative order min(a, b)^(-3/2). For such shapes outside the band the
  fraction takes at most about 1,300 steps, and for shapes below 1e10 at most about
  11,600, at the center.
• At tiny shapes the values were wrong rather than slow: the exponent s of the series of
  1 - I (about a log x) was absorbed by Stirling corrections of order 0.003 that the sum
  adds and subtracts, and subnormal shapes rounded the fraction's first elements. Below
  1e-3, where |s| < 1e-4, s is now formed divided by a; below 1e-280 the fraction is
  evaluated divided by a.
• Two inaccuracies came to light on the same grid: with a + b beyond the double range
  the deviance term of (a + b)(1 - x) = b + lam overflowed and was taken as 0 (1.7e-10
  relative in log C at a = 1.8e308, b = 1e300), and for shapes more than 1e17 apart with
  a tiny x the direct deviance form subtracted logarithms near 690 (2.1e-12 at a = 30,
  b = 1e300, x = 1e-300). Both are now formed at half scale and from the rounded
  product (a + b) x.

Now every one of the 1,728 values takes at most 0.33 ms per call from Python (0.21 ms,
best of three, in the validation), the grid evaluates as one array in 0.39 ms, no value
in the domain is NaN, and all 882 points of the extreme family agree with mpmath to
2.0e-13. Results on the 19,493 points of the existing validation grids were compared
bit for bit between the old and the new package: one point changed, F(x = 1; 10^12,
10^12), whose cdf and sf were 0.5000000000000071 and 0.49999999999999295 and are now
0.5, which mpmath (and the symmetry of I_1/2(a, a)) gives exactly. A first version of the
small-shape change moved the t quantile at p = 0.26, df = 0.001 from 9.7e-14 to 2.1e-13
of mpmath's value; the change was narrowed to where the old sum loses s, and that
quantile is again unchanged to the bit. The incomplete gamma function, which hung at
gammainc(1.8e308, 1.8e308) for the same reason, gives 1/2 there, and
log Q(1e300, 1.8e308) = -1.7976929347904657e308, as mpmath gives.

The reference at the edges of the double range needed care: where one tail is within
1e-30 of 1, the other is its complement, so the continued fraction is evaluated in
2,400 bits (enough for a + b with a subnormal and the largest double) to 400 digits;
near the center of large shapes, where the fraction converges too slowly, the beta
density is integrated numerically.

Validation

validation/pca_and_dot_plots.py (R 4.4.3, DESeq2 1.46.0, edgeR 4.4.0, ggplot2 4.0.3;
scanpy 1.12.4, anndata 0.13.4, matplotlib 3.11.2):

Check family	Reference	Largest difference	Tolerance	Outcome
Principal components: coordinates up to sign and shares of variance, with the top 100 genes (components 1 and 2, 3 and 4) and with 500; genes of equal variance; a cube's cell type prepared by edgeR	plotPCA(returnData = TRUE)	1.1e-14	1e-10	3 pass
The figure: points up to sign, axis labels, equal scales, samples coloured by condition and batch	plotPCA	6.9e-15	1e-10	1 pass
Two samples: the first component; DESeq2's second is rounding noise	plotPCA(returnData = TRUE)	8.9e-16	n/a	1 deliberate difference
Components whose sign differs from R's	plotPCA(returnData = TRUE)	2 of 10	n/a	1 information
Dot plots: areas and the fractions they encode, colour values plain, scaled by gene or group and over expressing cells, for all samples and for 3 samples and 3 cell types in a chosen order, genes in groups; names of rows and columns, the key of the areas	scanpy.pl.dotplot	1.0e-15	1e-10	6 pass
Brackets of genes in groups: columns and names	scanpy.pl.dotplot	identical	0	1 pass

validation/distributions.py (mpmath 1.3.0), the rows that changed (the other 46
checks print the values they printed before):

Check family	Reference	Largest difference	Tolerance	Outcome
The incomplete beta function, its complement and their logarithms: every combination of 11 shapes and 7 arguments from 5e-324 to 1.8e308, and 35 points near the center of shapes from 1e10 to 1e20	mpmath, 50 digits	2.04e-13	1e-12	4 pass
Time per value at the 1,728 combinations of 12 values as shapes and argument	timed	2.1e-4 s (median 2.5e-5 s)	1e-2 s	1 pass

REF=/home/parashar/.cache/cytepool-reference
RSCRIPT=$REF/r/bin/Rscript SCANPY_PYTHON=$REF/python-scanpy/bin/python PYTHONPATH=. \
    $REF/python/bin/python validation/pca_and_dot_plots.py
PYTHONPATH=. $REF/python/bin/python validation/distributions.py

Both exit without a failure (19.8 s, 671 MiB; 12 min 28 s, 174 MiB, of which mpmath
took 128 s at the edges of the double range).

Tests

• tests/test_pca_and_dot_plots.py, 37 tests in 1.0 s warm: DESeq2's coordinates,
  shares, ties, cube and figure from the fixtures; scanpy's areas, colour values,
  labels, key and brackets for the six cases; by hand, three samples on two genes whose
  components are the genes, two samples, one gene, and a gene expressed in half the
  cells of two samples (areas and means); invariants: the order of samples and genes,
  negated, shifted and scaled values, orthonormal loadings, a merged group and pooled
  samples against the same cells given one label, scaled means from 0 to 1; errors and
  their messages.
• tests/test_special.py, 9 new: the extreme fixture against mpmath (all 882 points,
  9.8 KB), 3,456 values at the grid within 2 s (under 1 ms here), I_1/2(a, a) = 1/2 up to the largest double, the
  reflection I_x(a, b) = 1 - I_(1-x)(b, a) inside and outside the band, 1 - x^a at
  shapes down to 5e-324, the incomplete gamma function where shape and argument overflow.
• The masked-input inventory: 33 new parametrized checks.

Checks

• The full suite on master with the "Compare conditions" guides: 4,901 passed, 18
  skipped; the same with bounds checks, and with the rounding changed
  (OPENBLAS_CORETYPE=Sandybridge, 4 threads, NUMBA_CPU_NAME=generic, a fresh
  cache); Python 3.14: 4,922 passed; Python 3.15: 4,816 passed, 107 skipped.
• validation/pca_and_dot_plots.py and validation/distributions.py run again on
  that tree: every check passes.
• ruff check and ruff format --check on cytepool tests validation docs: clean;
  mypy cytepool: no issues in 72 files.
• make -C docs api and make -C docs html (Node.js 24): built without warnings;
  tests/test_documentation.py: 74 passed.