Offline Chart Analysis

These functions operate on a completed survey result. The elevation map is read from the measured log-likelihood; z_fidelity_map flags projection aliasing.

elevation_map

elevation_map(surface, n_grid: int = DEFAULT_MOVIE_N_GRID, quantile: float = DEFAULT_ELEVATION_QUANTILE, datum=None, threshold=None, extent=None) -> np.ma.MaskedArray

Render map elevation from the measured per-sample log-likelihood.

The elevation of a cell is a high quantile of the log-likelihood over the samples in that cell. The value comes from the alignment likelihood, not from the bias, so it does not depend on delta_T.

Parameters:

Name	Type	Description	Default
surface	chart result	A completed survey result (duck-typed). The function reads only lambda1, lambda2, logL, and h.	required
n_grid	int	Grid resolution along each axis.	DEFAULT_MOVIE_N_GRID
quantile	float	Per-cell log-likelihood quantile in (0, 1]. A value of 1.0 gives the strict per-cell maximum.	DEFAULT_ELEVATION_QUANTILE
datum	float or None	An offset. The function subtracts it from every elevation. Pass surface.logL.max() to reference the elevation to the per-gene maximum.	None
threshold	float or None	The function masks a cell whose elevation is below this value.	None
extent	tuple or None	The bounds (lam1_min, lam1_max, lam2_min, lam2_max) for the grid. When None, the bounds come from the sample positions with a margin of surface.h.	None

Returns:

Type	Description
(MaskedArray, shape(n_grid, n_grid), float64)	Elevation at each cell. Element [i, j] is lambda1 bin i and lambda2 bin j. The function masks a cell that has no samples.

Source code in src/hifuku/analysis2d.py

def elevation_map(
    surface,
    n_grid: int = DEFAULT_MOVIE_N_GRID,
    quantile: float = DEFAULT_ELEVATION_QUANTILE,
    datum=None,
    threshold=None,
    extent=None,
) -> np.ma.MaskedArray:
    """
    Render map elevation from the measured per-sample log-likelihood.

    The elevation of a cell is a high quantile of the log-likelihood over the
    samples in that cell.  The value comes from the alignment likelihood, not
    from the bias, so it does not depend on delta_T.

    Parameters
    ----------
    surface : chart result
        A completed survey result (duck-typed).  The function reads only ``lambda1``, ``lambda2``,
        ``logL``, and ``h``.
    n_grid : int
        Grid resolution along each axis.
    quantile : float
        Per-cell log-likelihood quantile in (0, 1].  A value of 1.0 gives the
        strict per-cell maximum.
    datum : float or None
        An offset.  The function subtracts it from every elevation.  Pass
        ``surface.logL.max()`` to reference the elevation to the per-gene
        maximum.
    threshold : float or None
        The function masks a cell whose elevation is below this value.
    extent : tuple or None
        The bounds ``(lam1_min, lam1_max, lam2_min, lam2_max)`` for the grid.
        When None, the bounds come from the sample positions with a margin of
        ``surface.h``.

    Returns
    -------
    np.ma.MaskedArray, shape (n_grid, n_grid), float64
        Elevation at each cell.  Element ``[i, j]`` is lambda1 bin ``i`` and
        lambda2 bin ``j``.  The function masks a cell that has no samples.
    """
    from scipy.stats import binned_statistic_2d

    lam1 = np.asarray(surface.lambda1, dtype=np.float64)
    lam2 = np.asarray(surface.lambda2, dtype=np.float64)
    logL = np.asarray(surface.logL, dtype=np.float64)

    if extent is None:
        h = float(surface.h)
        pos = np.stack([lam1, lam2], axis=1)
        lam1_min, lam1_max, lam2_min, lam2_max = _data_bounds(pos, margin=h)
    else:
        lam1_min, lam1_max, lam2_min, lam2_max = extent

    edges1 = np.linspace(lam1_min, lam1_max, n_grid + 1)
    edges2 = np.linspace(lam2_min, lam2_max, n_grid + 1)

    def _cell_quantile(v):
        return np.quantile(v, quantile) if v.size else np.nan

    stat, _, _, _ = binned_statistic_2d(
        lam1, lam2, logL, statistic=_cell_quantile, bins=[edges1, edges2]
    )
    elev = stat.astype(np.float64)
    if datum is not None:
        elev = elev - float(datum)

    mask = np.isnan(elev)
    if threshold is not None:
        mask = mask | (elev < float(threshold))

    return np.ma.array(elev, mask=mask)

z_fidelity_map

z_fidelity_map(surface, n_grid: int = DEFAULT_MOVIE_N_GRID, extent=None, min_samples: int = 2) -> np.ma.MaskedArray

Per-cell aliasing indicator from the variance of the residual z.

The chart projection is many-to-one, so distinct trees can share a chart point. Within a cell, this function splits the variance of z into a within-walker part and a between-walker part. The within-walker part is genuine out-of-plane roughness. The between-walker part is aliasing: several smooth sheets stacked at one chart point. The two parts sum to the total variance (law of total variance).

Parameters:

Name	Type	Description	Default
surface	chart result	A completed survey result (duck-typed). The function reads only lambda1, lambda2, z_residual, chain_id, n_chains, and h.	required
n_grid	int	Grid resolution along each axis.	DEFAULT_MOVIE_N_GRID
extent	tuple or None	The bounds (lam1_min, lam1_max, lam2_min, lam2_max) for the grid. When None, the bounds come from the sample positions with a margin of surface.h. Pass the same value as elevation_map to align the two.	None
min_samples	int	The function masks a cell with fewer than this many samples.	2

Returns:

Type	Description
(MaskedArray, shape(n_grid, n_grid), float64)	The aliasing indicator between_var / total_var in [0, 1]. A value near 1 marks aliasing. A value near 0 marks genuine roughness. The function masks a cell with too few samples or with zero z variance.

Notes

This is a measure-only diagnostic. z is a scalar residual, a cheap and partial probe of a high-dimensional orthogonal complement.

Source code in src/hifuku/analysis2d.py

def z_fidelity_map(
    surface,
    n_grid: int = DEFAULT_MOVIE_N_GRID,
    extent=None,
    min_samples: int = 2,
) -> np.ma.MaskedArray:
    """
    Per-cell aliasing indicator from the variance of the residual z.

    The chart projection is many-to-one, so distinct trees can share a chart
    point.  Within a cell, this function splits the variance of z into a
    within-walker part and a between-walker part.  The within-walker part is
    genuine out-of-plane roughness.  The between-walker part is aliasing: several
    smooth sheets stacked at one chart point.  The two parts sum to the total
    variance (law of total variance).

    Parameters
    ----------
    surface : chart result
        A completed survey result (duck-typed).  The function reads only ``lambda1``, ``lambda2``,
        ``z_residual``, ``chain_id``, ``n_chains``, and ``h``.
    n_grid : int
        Grid resolution along each axis.
    extent : tuple or None
        The bounds ``(lam1_min, lam1_max, lam2_min, lam2_max)`` for the grid.
        When None, the bounds come from the sample positions with a margin of
        ``surface.h``.  Pass the same value as ``elevation_map`` to align the two.
    min_samples : int
        The function masks a cell with fewer than this many samples.

    Returns
    -------
    np.ma.MaskedArray, shape (n_grid, n_grid), float64
        The aliasing indicator ``between_var / total_var`` in [0, 1].  A value
        near 1 marks aliasing.  A value near 0 marks genuine roughness.  The
        function masks a cell with too few samples or with zero z variance.

    Notes
    -----
    This is a measure-only diagnostic.  z is a scalar residual, a cheap and
    partial probe of a high-dimensional orthogonal complement.
    """
    lam1 = np.asarray(surface.lambda1, dtype=np.float64)
    lam2 = np.asarray(surface.lambda2, dtype=np.float64)
    z = np.asarray(surface.z_residual, dtype=np.float64)
    cid = np.asarray(surface.chain_id, dtype=np.int64)
    n_chains = int(surface.n_chains)

    valid = cid >= 0
    lam1, lam2, z, cid = lam1[valid], lam2[valid], z[valid], cid[valid]

    if extent is None:
        h = float(surface.h)
        pos = np.stack([lam1, lam2], axis=1)
        lam1_min, lam1_max, lam2_min, lam2_max = _data_bounds(pos, margin=h)
    else:
        lam1_min, lam1_max, lam2_min, lam2_max = extent

    edges1 = np.linspace(lam1_min, lam1_max, n_grid + 1)
    edges2 = np.linspace(lam2_min, lam2_max, n_grid + 1)
    ix = np.clip(np.digitize(lam1, edges1) - 1, 0, n_grid - 1)
    iy = np.clip(np.digitize(lam2, edges2) - 1, 0, n_grid - 1)
    cell = ix * n_grid + iy
    n_cells = n_grid * n_grid

    # One group per (cell, walker).
    gkey = cell * n_chains + cid
    uniq, inv = np.unique(gkey, return_inverse=True)
    n_g = np.bincount(inv).astype(np.float64)
    sum_z = np.bincount(inv, weights=z)
    sum_z2 = np.bincount(inv, weights=z * z)
    mu_g = sum_z / n_g
    var_g = np.maximum(sum_z2 / n_g - mu_g * mu_g, 0.0)
    cell_of_g = uniq // n_chains

    N_c = np.bincount(cell_of_g, weights=n_g, minlength=n_cells)
    within_num = np.bincount(cell_of_g, weights=n_g * var_g, minlength=n_cells)
    sum_nmu = np.bincount(cell_of_g, weights=n_g * mu_g, minlength=n_cells)
    sum_nmu2 = np.bincount(cell_of_g, weights=n_g * mu_g * mu_g, minlength=n_cells)

    with np.errstate(invalid="ignore", divide="ignore"):
        within_var = within_num / N_c
        mu_c = sum_nmu / N_c
        between_var = np.maximum(sum_nmu2 / N_c - mu_c * mu_c, 0.0)
        total_var = within_var + between_var
        aliasing = between_var / total_var

    aliasing = aliasing.reshape(n_grid, n_grid)
    total_var = total_var.reshape(n_grid, n_grid)
    N_c = N_c.reshape(n_grid, n_grid)

    eps = 1e-12
    mask = (N_c < min_samples) | ~(total_var > eps) | np.isnan(aliasing)
    return np.ma.array(aliasing, mask=mask)

grid_extent

grid_extent(surface) -> tuple[float, float, float, float]

Return (lam1_min, lam1_max, lam2_min, lam2_max) for a completed surface.

Bounds are computed from all particle positions across all stored snapshots plus all MCMC sample positions, with a margin of surface.h on each side. Use the returned tuple as the extent argument to imshow and to set axis limits:

ext = grid_extent(surf)
ax.imshow(data, extent=[ext[0], ext[1], ext[2], ext[3]], ...)

Source code in src/hifuku/analysis2d.py

def grid_extent(surface) -> tuple[float, float, float, float]:
    """
    Return (lam1_min, lam1_max, lam2_min, lam2_max) for a completed surface.

    Bounds are computed from all particle positions across all stored snapshots
    plus all MCMC sample positions, with a margin of surface.h on each side.
    Use the returned tuple as the ``extent`` argument to ``imshow`` and to set
    axis limits:

        ext = grid_extent(surf)
        ax.imshow(data, extent=[ext[0], ext[1], ext[2], ext[3]], ...)
    """
    h = float(surface.h)
    all_pos = surface.bias_round_positions.reshape(-1, 2)
    sample_pos = np.stack(
        [np.asarray(surface.lambda1), np.asarray(surface.lambda2)], axis=1
    )
    combined = np.concatenate([all_pos, sample_pos], axis=0)
    return _data_bounds(combined, margin=h)