surface

tit.stats.surface

fsaverage surface backend for cluster-based permutation stats.

The volumetric engine (:mod:tit.stats.engine) clusters voxels with scipy.ndimage.label on a 3-D grid -- which has no meaning for surface vertices. This module is the surface counterpart: it stacks the per-subject fsaverage field projections (.msh) written by :func:tit.source.project_fields_to_fsaverage into a (n_vertices, n_subjects) matrix and runs the same statistics with the same cluster conventions, swapping the grid clustering for graph connected-components over the fsaverage triangle adjacency.

It deliberately reuses the engine's space-agnostic kernels (:func:~tit.stats.engine.correlation, :func:~tit.stats.engine.ttest_ind, :func:~tit.stats.engine.ttest_rel, :func:~tit.stats.engine.pval_from_histogram) so the surface and volume paths give identical numbers on identical inputs; only the spatial layer differs. The volumetric engine is left untouched -- :mod:tit.stats.permutation dispatches here when config.space == "fsaverage".

Cluster conventions mirror the engine exactly: a cluster is a connected set of vertices with p < cluster_threshold (sign-restricted for one-sided tests); singletons (size 1) are ignored; cluster mass is the signed sum of t over its vertices; the null is the per-permutation max cluster stat; per-cluster p-values come from :func:~tit.stats.engine.pval_from_histogram.

load_group_surface_data

load_group_surface_data(subjects: list[tuple[str, str]], field: str, spacing: int) -> tuple[ndarray, list[str]]

Stack per-subject fsaverage projections into (n_vertices, n_subjects).

Each column is one node field of a *_space-fsaverage<N>_fields.msh (:func:tit.source.fsaverage.write_fsaverage_msh).

Parameters

subjects : list of (str, str) (subject_id, simulation_name) pairs (bare ids, no sub- prefix). field : str Which projected field to load -- one of :data:VALID_FSAVG_FIELDS. spacing : int fsaverage subdivision factor (5, 6, or 7).

Returns

(numpy.ndarray, list of str) The (n_vertices, n_subjects) data matrix and the subject ids in column order.

Source code in tit/stats/surface.py

def load_group_surface_data(
    subjects: list[tuple[str, str]], field: str, spacing: int
) -> tuple[np.ndarray, list[str]]:
    """Stack per-subject fsaverage projections into ``(n_vertices, n_subjects)``.

    Each column is one node field of a ``*_space-fsaverage<N>_fields.msh``
    (:func:`tit.source.fsaverage.write_fsaverage_msh`).

    Parameters
    ----------
    subjects : list of (str, str)
        ``(subject_id, simulation_name)`` pairs (bare ids, no ``sub-`` prefix).
    field : str
        Which projected field to load -- one of :data:`VALID_FSAVG_FIELDS`.
    spacing : int
        fsaverage subdivision factor (5, 6, or 7).

    Returns
    -------
    (numpy.ndarray, list of str)
        The ``(n_vertices, n_subjects)`` data matrix and the subject ids in
        column order.
    """
    if field not in VALID_FSAVG_FIELDS:
        raise ValueError(
            f"{field!r} is not on the fsaverage surface; valid: {VALID_FSAVG_FIELDS}. "
            + FSAVG_NOT_PROJECTED
        )
    pm = get_path_manager()
    expected = _FSAVG_NODES[spacing]
    columns: list[np.ndarray] = []
    ids: list[str] = []
    for sid, sim in subjects:
        msh_path = Path(pm.sim_fsaverage_fields(sid, sim, spacing))
        if not msh_path.exists():
            raise FileNotFoundError(
                f"No fsaverage projection for {sid}/{sim}: {msh_path}. "
                "Re-run its fsaverage projection first: simulate it again with the "
                "Simulator's 'Map fields to fsaverage' on, or map its existing outputs "
                "with tit.source.fsaverage.project_fields_to_fsaverage."
            )
        arr = read_fsaverage_field(msh_path, field)
        if arr.shape[0] != expected:
            raise ValueError(
                f"{msh_path.name}: expected {expected} fsaverage{spacing} vertices, "
                f"got {arr.shape[0]}"
            )
        columns.append(arr)
        ids.append(sid)
    return np.column_stack(columns), ids

build_fsaverage_adjacency

build_fsaverage_adjacency(spacing: int)

Block-diagonal lh+rh fsaverage vertex adjacency (cached per spacing).

No edges cross the hemisphere boundary, so a slow-wave cluster can never bridge the two hemispheres through a spurious midline edge -- matching the [lh; rh] node ordering the field projections are written in.

Mesh source, tried in order (see resources/fsaverage/README.md for why the bundled source is preferred on correctness grounds, not just to avoid a download): the vendored resources/fsaverage/{spacing}/ GIFTI surfaces (SimNIBS's own bundled fsaverage -- the same source :func:tit.source.fsaverage._compute_fields projects each subject's fields onto, via SimNIBS's cross_subject_map); then nilearn's fetch_surf_fsaverage (a runtime download, the previous sole behavior); a :class:RuntimeError naming both failed sources if neither is available.

Source code in tit/stats/surface.py

def build_fsaverage_adjacency(spacing: int):
    """Block-diagonal lh+rh fsaverage vertex adjacency (cached per spacing).

    No edges cross the hemisphere boundary, so a slow-wave cluster can never
    bridge the two hemispheres through a spurious midline edge -- matching the
    ``[lh; rh]`` node ordering the field projections are written in.

    Mesh source, tried in order (see ``resources/fsaverage/README.md`` for
    why the bundled source is preferred on correctness grounds, not just to
    avoid a download): the vendored ``resources/fsaverage/{spacing}/`` GIFTI
    surfaces (SimNIBS's own bundled fsaverage -- the same source
    :func:`tit.source.fsaverage._compute_fields` projects each subject's
    fields onto, via SimNIBS's ``cross_subject_map``); then nilearn's
    ``fetch_surf_fsaverage`` (a runtime download, the previous sole
    behavior); a :class:`RuntimeError` naming both failed sources if
    neither is available.
    """
    if spacing in _ADJ_CACHE:
        return _ADJ_CACHE[spacing]
    from scipy import sparse

    mesh = _load_fsaverage_mesh_bundled(spacing)
    source = "bundled resources/fsaverage"
    if mesh is None:
        try:
            mesh = _load_fsaverage_mesh_nilearn(spacing)
            source = "nilearn fetch_surf_fsaverage"
        except Exception as exc:
            raise RuntimeError(
                f"No fsaverage{spacing} mesh available: bundled "
                f"resources/fsaverage/{spacing}/ is missing, and the "
                f"nilearn fallback failed ({exc!r}). Vendor the resource "
                "(see resources/fsaverage/README.md) or ensure network "
                "access for nilearn's one-time download."
            ) from exc

    coords_l, faces_l, coords_r, faces_r = mesh
    adj = sparse.block_diag(
        [
            _faces_to_adjacency(faces_l, len(coords_l)),
            _faces_to_adjacency(faces_r, len(coords_r)),
        ]
    ).tocsr()
    if adj.shape[0] != _FSAVG_NODES[spacing]:
        raise ValueError(
            f"fsaverage{spacing} adjacency has {adj.shape[0]} nodes, "
            f"expected {_FSAVG_NODES[spacing]} (mesh source: {source})"
        )
    _ADJ_CACHE[spacing] = adj
    return adj

run_surface_correlation

run_surface_correlation(config, callback_handler=None, stop_callback=None) -> CorrelationResult

Vertexwise field-vs-response correlation on the fsaverage surface.

Source code in tit/stats/surface.py

def run_surface_correlation(
    config, callback_handler=None, stop_callback=None
) -> CorrelationResult:
    """Vertexwise field-vs-response correlation on the fsaverage surface."""
    import time

    t0 = time.time()
    output_dir, log, log_file = _output_dir_and_log(
        "correlation", config.analysis_name, callback_handler
    )
    log.info(
        "CLUSTER-BASED PERMUTATION TESTING - SURFACE CORRELATION (fsaverage%d)",
        config.fsaverage_spacing,
    )
    log.info(
        "field=%s  corr=%s  stat=%s  thr=%.3f  perms=%d  alpha=%.3f",
        config.fsaverage_field,
        config.correlation_type.value,
        config.cluster_stat.value,
        config.cluster_threshold,
        config.n_permutations,
        config.alpha,
    )

    pairs = [(s.subject_id, s.simulation_name) for s in config.subjects]
    effect = np.array([float(s.effect_size) for s in config.subjects], dtype=np.float64)
    weights = (
        np.array([float(s.weight) for s in config.subjects], dtype=np.float64)
        if config.use_weights
        else None
    )
    data, ids = load_group_surface_data(
        pairs, config.fsaverage_field, config.fsaverage_spacing
    )
    n_nodes = data.shape[0]
    adjacency = build_fsaverage_adjacency(config.fsaverage_spacing)

    # != 0 (not > 0): TI_normal is signed, so a consistently-inward vertex
    # must stay in the analysis; background is exactly 0 for all fields.
    valid_mask = np.any(data != 0, axis=1)
    valid_idx = np.flatnonzero(valid_mask)
    data_valid = data[valid_idx]
    log.info("Subjects=%d  valid vertices=%d/%d", len(ids), valid_idx.size, n_nodes)

    ctype = config.correlation_type.value
    # Pre-rank once for Spearman so each permutation skips re-ranking.
    if ctype == "spearman":
        from scipy.stats import rankdata

        data_valid = np.apply_along_axis(rankdata, 1, data_valid)
        preranked = True
    else:
        preranked = False

    def _maps(eff, w):
        r, t, p = correlation(
            data_valid,
            eff,
            correlation_type=ctype,
            weights=w,
            voxel_data_preranked=preranked,
        )
        r_full = np.zeros(n_nodes)
        t_full = np.zeros(n_nodes)
        p_full = np.ones(n_nodes)
        r_full[valid_idx], t_full[valid_idx], p_full[valid_idx] = r, t, p
        return r_full, t_full, p_full

    r_full, t_full, p_full = _maps(effect, weights)

    log.info("[null] %d permutations", config.n_permutations)
    null = np.empty(config.n_permutations)
    seed = resolve_seed(config.seed)  # recorded in the sidecar
    rng = np.random.default_rng(seed)
    for i in range(config.n_permutations):
        if stop_callback and stop_callback():
            raise KeyboardInterrupt("stopped")
        order = rng.permutation(len(ids))
        _, pt, pp = _maps(effect[order], None if weights is None else weights[order])
        labels, n = _label_graph_signed(
            (pp < config.cluster_threshold) & valid_mask, pt, adjacency, "two-sided"
        )
        null[i] = _max_cluster_stat(labels, n, pt, config.cluster_stat.value, tail=0)

    sig_mask, sig_clusters, observed = _identify_surface_clusters(
        t_full,
        p_full,
        valid_mask,
        adjacency,
        config.cluster_threshold,
        null,
        config.cluster_stat.value,
        config.alpha,
        "two-sided",
        r_full=r_full,
    )
    for obs in observed[:5]:
        log.info(
            "  cluster %d: %s=%.2f size=%d p=%.4f",
            obs["id"],
            config.cluster_stat.value,
            obs["stat_value"],
            obs["size"],
            obs["p_value"],
        )

    metadata = _run_metadata(
        config,
        "correlation",
        {
            "correlation_type": ctype,
            "n_subjects": len(ids),
            "subjects": [s.subject_id for s in config.subjects],
            "responses": effect.tolist(),
        },
        seed,
    )
    msh_path = _write_surface_outputs(
        output_dir,
        config,
        metadata,
        {
            "r": r_full,
            "t": t_full,
            "p": p_full,
            "sig_mask": sig_mask,
            "mean_field": data.mean(axis=1),
        },
        sig_clusters,
        null,
        data,
        [(s.subject_id, s.simulation_name, s.effect_size) for s in config.subjects],
        "response",
    )
    log.info("Saved surface maps -> %s", msh_path)

    return CorrelationResult(
        success=True,
        output_dir=output_dir,
        n_subjects=len(ids),
        n_significant_voxels=int(sig_mask.sum()),
        n_significant_clusters=len(sig_clusters),
        cluster_threshold=_null_threshold(null, config.alpha, config.n_permutations),
        analysis_time=time.time() - t0,
        clusters=sig_clusters,
        log_file=log_file,
    )

run_surface_group_comparison

run_surface_group_comparison(config, callback_handler=None, stop_callback=None) -> GroupComparisonResult

Responder-vs-non-responder t-test on the fsaverage surface.

Source code in tit/stats/surface.py

def run_surface_group_comparison(
    config, callback_handler=None, stop_callback=None
) -> GroupComparisonResult:
    """Responder-vs-non-responder t-test on the fsaverage surface."""
    import time

    t0 = time.time()
    output_dir, log, log_file = _output_dir_and_log(
        "group_comparison", config.analysis_name, callback_handler
    )
    log.info(
        "CLUSTER-BASED PERMUTATION TESTING - SURFACE GROUP (fsaverage%d)",
        config.fsaverage_spacing,
    )

    # Responders first, then non-responders, so column [:n_resp] are responders.
    resp = [s for s in config.subjects if s.response == 1]
    non = [s for s in config.subjects if s.response == 0]
    ordered = resp + non
    n_resp = len(resp)
    paired = config.test_type.value == "paired"
    if paired:
        check_paired_sizes(n_resp, len(non))
    pairs = [(s.subject_id, s.simulation_name) for s in ordered]
    data, ids = load_group_surface_data(
        pairs, config.fsaverage_field, config.fsaverage_spacing
    )
    n_nodes = data.shape[0]
    n_total = data.shape[1]
    adjacency = build_fsaverage_adjacency(config.fsaverage_spacing)

    # != 0 (not > 0): TI_normal is signed, so a consistently-inward vertex
    # must stay in the analysis; background is exactly 0 for all fields.
    valid_mask = np.any(data != 0, axis=1)
    valid_idx = np.flatnonzero(valid_mask)
    alt = config.alternative.value

    def _test(mat):
        if paired:
            return ttest_rel(mat, n_resp, alternative=alt)
        return ttest_ind(mat, n_resp, n_total - n_resp, alternative=alt)

    # Same rule as engine.ttest_voxelwise: a vertex with zero within-group variance
    # has t = nan (constant groups) or +/-inf (perfect separation); it leaves the
    # analysis, and a run with no testable vertex left is an error, not 0 clusters.
    _, _, degenerate = exclude_degenerate(*_test(data[valid_idx]), log)
    valid_mask[valid_idx[degenerate]] = False
    if not valid_mask.any():
        raise no_testable_voxel_error(valid_idx.size)
    valid_idx = np.flatnonzero(valid_mask)
    data_valid = data[valid_idx]
    log.info(
        "resp=%d  non=%d  valid vertices=%d  paired=%s",
        n_resp,
        len(non),
        valid_idx.size,
        paired,
    )

    def _maps(mat):
        # A relabelling can make a vertex degenerate; an inf mass would swamp the null.
        t, p = _neutralise_degenerate(*_test(mat))
        t_full = np.zeros(n_nodes)
        p_full = np.ones(n_nodes)
        t_full[valid_idx], p_full[valid_idx] = t, p
        return t_full, p_full

    t_full, p_full = _maps(data_valid)

    log.info("[null] %d permutations", config.n_permutations)
    null = np.empty(config.n_permutations)
    seed = resolve_seed(config.seed)  # recorded in the sidecar
    rng = np.random.default_rng(seed)
    for i in range(config.n_permutations):
        if stop_callback and stop_callback():
            raise KeyboardInterrupt("stopped")
        if paired:
            flips = rng.choice([-1.0, 1.0], size=n_resp)
            mean = (data_valid[:, :n_resp] + data_valid[:, n_resp:]) / 2
            diff = (data_valid[:, :n_resp] - data_valid[:, n_resp:]) / 2
            perm = np.empty_like(data_valid)
            perm[:, :n_resp] = mean + flips * diff
            perm[:, n_resp:] = mean - flips * diff
        else:
            perm = data_valid[:, rng.permutation(n_total)]
        pt, pp = _maps(perm)
        labels, n = _label_graph_signed(
            (pp < config.cluster_threshold) & valid_mask, pt, adjacency, alt
        )
        null[i] = _max_cluster_stat(
            labels, n, pt, config.cluster_stat.value, tail=tail_from_alternative(alt)
        )

    sig_mask, sig_clusters, observed = _identify_surface_clusters(
        t_full,
        p_full,
        valid_mask,
        adjacency,
        config.cluster_threshold,
        null,
        config.cluster_stat.value,
        config.alpha,
        alt,
    )
    for obs in observed[:5]:
        log.info(
            "  cluster %d: %s=%.2f size=%d p=%.4f",
            obs["id"],
            config.cluster_stat.value,
            obs["stat_value"],
            obs["size"],
            obs["p_value"],
        )

    metadata = _run_metadata(
        config,
        "group_comparison",
        {
            "test_type": config.test_type.value,
            "alternative": alt,
            "n_subjects": n_total,
            "groups": {
                config.group1_name: [s.subject_id for s in resp],
                config.group2_name: [s.subject_id for s in non],
            },
        },
        seed,
    )
    msh_path = _write_surface_outputs(
        output_dir,
        config,
        metadata,
        {
            "t": t_full,
            "p": p_full,
            "sig_mask": sig_mask,
            "mean_responders": data[:, :n_resp].mean(axis=1),
            "mean_non_responders": data[:, n_resp:].mean(axis=1),
        },
        sig_clusters,
        null,
        data,
        [(s.subject_id, s.simulation_name, config.group1_name) for s in resp]
        + [(s.subject_id, s.simulation_name, config.group2_name) for s in non],
        "group",
    )
    log.info("Saved surface maps -> %s", msh_path)

    return GroupComparisonResult(
        success=True,
        output_dir=output_dir,
        n_responders=n_resp,
        n_non_responders=len(non),
        n_significant_voxels=int(sig_mask.sum()),
        n_significant_clusters=len(sig_clusters),
        cluster_threshold=_null_threshold(null, config.alpha, config.n_permutations),
        analysis_time=time.time() - t0,
        clusters=sig_clusters,
        log_file=log_file,
    )