Skip to content

Separation metrics (pamica.metrics)

Quality metrics for a decomposition, as free functions over plain arrays. They are backend-agnostic and independent of any fitted model object, so they work on sources from any source.

  • mir — Mutual Information Reduction, in nats: how much mutual information a linear unmixing removes from the data. Higher is a better separation. It needs the full raw-data-to-sources transform (the unmixing composed with the sphering matrix), which must be square and invertible, since the estimate includes a log-Jacobian term.
  • pairwise_mi — the symmetric mutual-information matrix between sources. Its diagonal is each source's own entropy, not a mutual information.
  • block_diagonal_order — a permutation that clusters mutually dependent components near the diagonal, for reading structure out of a pairwise_mi matrix.

Most callers should prefer the model accessors, which compose the transform for you and refuse an unusable fit: AMICA.mir and AMICA.pmi. Use these free functions when you have sources or an unmixing from somewhere else. To plot a pairwise_mi matrix, see plot_pmi_heatmap, which applies block_diagonal_order itself and masks the entropy diagonal.

from pamica.metrics import mir, pairwise_mi, block_diagonal_order

Provenance

mir is a direct port of getMIR.m from bigdelys/pre_ICA_cleaning (Apache-2.0); the licence is vendored in THIRD_PARTY_NOTICES.md. It agrees with the original to 1.7e-15 relative on the bundled sample data.

pairwise_mi and block_diagonal_order are a clean-room reimplementation. The comparable MATLAB code (minfojp.m in postAmicaUtility) is GPL-2.0-or-later and pamica is BSD-3-Clause, so that source was never read; the implementation works from the published description in Delorme et al. (2012), "Independent EEG sources are dipolar", PLoS ONE. It agrees with that reference at r=0.9887 on identical signals.

pamica.metrics.mir

Mutual Information Reduction (MIR) ICA separation-quality metric.

MIR measures how much a linear unmixing transform reduces the total mutual information among channels: it is the (signed) difference between the sum of per-channel differential entropies before and after unmixing, corrected by the log-Jacobian of the transform. A well-separated decomposition drives the unmixed channels toward independence, so MIR should be large and positive. The identity transform (or any scalar multiple of it) gives MIR exactly zero, since the bin-partition estimator below is scale-invariant; a generic non-identity rotation of dependent, non-Gaussian data need not (that directional sensitivity is the entire mechanism ICA exploits).

Ported from getMIR.m / its nested getent4 in bigdelys/pre_ICA_cleaning <https://github.com/bigdelys/pre_ICA_cleaning/blob/b6034f03889a6a418968ee119123f3df55251957/getMIR.m>_ (Apache License 2.0, full text reproduced in THIRD_PARTY_NOTICES.md per that license's Section 4(a)). This is a permitted derivative under Apache-2.0; the port keeps the original's binned-entropy estimator (including its bin-count-1 correction and eigenvalue-based log-Jacobian) rather than substituting a different entropy estimator.

A 2023 MATLAB adaptation of getMIR (the nested mir() function in sccn/NEMAR-pipeline's eeg_nemar_dataqual.m, currently dead/commented-out code there) added an explicit sphering step ahead of this computation, but its mir(data, linT) references eig(W) and /N, neither of which are its own parameters nor ever assigned anywhere in the enclosing function -- it would raise an undefined-variable error if uncommented, which is why it's dead code. It also depends on a robust_sphering_matrix routine that itself calls GPL-licensed helpers (block_geometric_median, hlp_memfree), out of scope for this BSD-3-Clause project. This port instead takes an explicit unmixing matrix: callers that want the sphere-then-unmixing composition (e.g. wiring MIR to a fitted AMICA model) compose it themselves before calling mir.

MIR follows the separation-quality-metric approach introduced by Delorme, Palmer, Onton, Oostenveld, and Makeig (2012), "Independent EEG sources are dipolar", PLoS ONE (this repo's paper.bib key delorme2012independent).

mir(unmixing, data, nbins=None)

Mutual Information Reduction of unmixing applied to data, in nats.

Parameters:

Name Type Description Default
unmixing ndarray

Square (n, n) linear transform applied to data (e.g. an ICA unmixing matrix, optionally composed with a sphering matrix).

 required
data ndarray

(n, N) data matrix unmixing is applied to.

 required
nbins int

Histogram bin count for the marginal-entropy estimator; see _marginal_entropies. Defaults to round(3 * log2(1 + N/10)).

 None

Returns:

Name Type Description
mir_nats float

Mutual information (in nats) removed by unmixing from data.
variance float

Variance of the MIR estimate.

Raises:

Type Description
ValueError

If unmixing is singular or near-singular (the log-Jacobian term is undefined for a non-invertible transform), or if data (or the transformed unmixing @ data) contains non-finite values or a constant channel -- a non-finite unmixing is also caught this way, indirectly, via the transformed data.
Source code in pamica/metrics/mir.py 
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
def mir(
    unmixing: np.ndarray, data: np.ndarray, nbins: int | None = None
) -> tuple[float, float]:
    """Mutual Information Reduction of `unmixing` applied to `data`, in nats.

    Parameters
    ----------
    unmixing : np.ndarray
        Square (n, n) linear transform applied to `data` (e.g. an ICA
        unmixing matrix, optionally composed with a sphering matrix).
    data : np.ndarray
        (n, N) data matrix `unmixing` is applied to.
    nbins : int, optional
        Histogram bin count for the marginal-entropy estimator; see
        `_marginal_entropies`. Defaults to ``round(3 * log2(1 + N/10))``.

    Returns
    -------
    mir_nats : float
        Mutual information (in nats) removed by `unmixing` from `data`.
    variance : float
        Variance of the MIR estimate.

    Raises
    ------
    ValueError
        If `unmixing` is singular or near-singular (the log-Jacobian term is
        undefined for a non-invertible transform), or if `data` (or the
        transformed `unmixing @ data`) contains non-finite values or a
        constant channel -- a non-finite `unmixing` is also caught this way,
        indirectly, via the transformed data.
    """
    eigvals = np.linalg.eigvals(unmixing)
    min_abs_eig = np.min(np.abs(eigvals))
    if min_abs_eig < _MIN_ABS_EIGENVALUE:
        raise ValueError(
            f"mir(): unmixing matrix is singular or near-singular (smallest "
            f"|eigenvalue| = {min_abs_eig:.3e}); the log-Jacobian term is "
            "undefined for a non-invertible transform. Verify `unmixing` is "
            "full rank (e.g. check for duplicated/collinear components from "
            "an unconverged ICA fit)."
        )
    hx, vx = _marginal_entropies(data, nbins)
    y = unmixing @ data
    hy, vy = _marginal_entropies(y, nbins)
    n_samples = data.shape[1]
    # eigvals (not slogdet), matching getMIR.m's own switch from det to eig.
    mir_nats = float(np.sum(np.log(np.abs(eigvals))) + np.sum(hx) - np.sum(hy))
    variance = float((np.sum(vx) + np.sum(vy)) / n_samples)
    return mir_nats, variance

pamica.metrics.pairwise_mi(sources, nbins=None)

Symmetric (n, n) pairwise mutual-information matrix, in nats.

Parameters:

Name Type Description Default
sources ndarray

(n_components, n_samples) signal matrix (e.g. ICA source activations or raw channels).

 required
nbins int

Histogram bin count per axis of the joint 2-D histogram. Defaults to round(3 * log2(1 + N/10)); see pamica.metrics._common.resolve_nbins.

 None

Returns:

Name Type Description
mi_matrix ndarray

(n, n) symmetric matrix; mi_matrix[i, j] is the mutual information (nats) between components i and j. The diagonal holds each component's own entropy (mi_matrix[i, i] ~= H(X_i)).

Raises:

Type Description
ValueError

If sources contains non-finite values, a constant channel, or nbins is too large for n_samples (too many joint-histogram cells to be estimated from that many samples).
Source code in pamica/metrics/pmi.py 
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
def pairwise_mi(sources: np.ndarray, nbins: int | None = None) -> np.ndarray:
    """Symmetric (n, n) pairwise mutual-information matrix, in nats.

    Parameters
    ----------
    sources : np.ndarray
        (n_components, n_samples) signal matrix (e.g. ICA source activations
        or raw channels).
    nbins : int, optional
        Histogram bin count per axis of the joint 2-D histogram. Defaults to
        ``round(3 * log2(1 + N/10))``; see `pamica.metrics._common.resolve_nbins`.

    Returns
    -------
    mi_matrix : np.ndarray
        (n, n) symmetric matrix; `mi_matrix[i, j]` is the mutual information
        (nats) between components `i` and `j`. The diagonal holds each
        component's own entropy (`mi_matrix[i, i] ~= H(X_i)`).

    Raises
    ------
    ValueError
        If `sources` contains non-finite values, a constant channel, or
        `nbins` is too large for `n_samples` (too many joint-histogram cells
        to be estimated from that many samples).
    """
    validate_signal_matrix(sources)
    n, n_samples = sources.shape
    nbins = resolve_nbins(n_samples, nbins)
    if nbins**2 > n_samples:
        raise ValueError(
            f"pairwise_mi: nbins={nbins} gives a {nbins}x{nbins} joint "
            f"histogram ({nbins**2} cells) for only {n_samples} samples; "
            "most cells would be empty or singleton, making the MI estimate "
            "meaningless. Reduce nbins or supply more samples."
        )

    mi_matrix = np.zeros((n, n))
    for i in range(n):
        for j in range(i, n):
            joint_counts, _, _ = np.histogram2d(sources[i], sources[j], bins=nbins)
            # Marginals are summed from this same joint histogram rather than
            # independently rebinned, so H(X), H(Y), H(X,Y) share one
            # consistent partition (and the diagonal i==j reduces exactly to
            # H(X_i), since the joint histogram of a channel with itself is
            # diagonal-only).
            marginal_x = joint_counts.sum(axis=1)
            marginal_y = joint_counts.sum(axis=0)
            h_x = _binned_entropy_from_counts(marginal_x, n_samples)
            h_y = _binned_entropy_from_counts(marginal_y, n_samples)
            h_xy = _binned_entropy_from_counts(joint_counts.ravel(), n_samples)
            mi = h_x + h_y - h_xy
            mi_matrix[i, j] = mi
            mi_matrix[j, i] = mi
    return mi_matrix

pamica.metrics.block_diagonal_order(mi_matrix)

Greedy nearest-neighbor-chain permutation clustering high-MI pairs near the diagonal.

Parameters:

Name Type Description Default
mi_matrix ndarray

(n, n) symmetric mutual-information matrix, e.g. from pairwise_mi.

 required

Returns:

Name Type Description
order ndarray

Length-n permutation of 0..n-1. Reordering both axes of mi_matrix by order tends to cluster high-MI pairs adjacent to the diagonal.

Raises:

Type Description
ValueError

If mi_matrix isn't square, contains non-finite values, or isn't symmetric.
Source code in pamica/metrics/pmi.py 
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
def block_diagonal_order(mi_matrix: np.ndarray) -> np.ndarray:
    """Greedy nearest-neighbor-chain permutation clustering high-MI pairs near the diagonal.

    Parameters
    ----------
    mi_matrix : np.ndarray
        (n, n) symmetric mutual-information matrix, e.g. from `pairwise_mi`.

    Returns
    -------
    order : np.ndarray
        Length-`n` permutation of `0..n-1`. Reordering both axes of
        `mi_matrix` by `order` tends to cluster high-MI pairs adjacent to the
        diagonal.

    Raises
    ------
    ValueError
        If `mi_matrix` isn't square, contains non-finite values, or isn't
        symmetric.
    """
    validate_mi_matrix(mi_matrix)
    n = mi_matrix.shape[0]
    if n <= 1:
        return np.arange(n)

    off_diag = mi_matrix.copy()
    np.fill_diagonal(off_diag, -np.inf)
    i, j = np.unravel_index(np.argmax(off_diag), off_diag.shape)
    chain = [int(i), int(j)]
    remaining = set(range(n)) - {int(i), int(j)}

    while remaining:
        best_value = -np.inf
        best_component = -1
        best_at_start = False
        for c in remaining:
            value_start = mi_matrix[chain[0], c]
            if value_start > best_value:
                best_value = value_start
                best_component = c
                best_at_start = True
            value_end = mi_matrix[chain[-1], c]
            if value_end > best_value:
                best_value = value_end
                best_component = c
                best_at_start = False
        if best_at_start:
            chain.insert(0, best_component)
        else:
            chain.append(best_component)
        remaining.remove(best_component)

    return np.array(chain)