calc

tit.calc

Temporal interference field calculation utilities.

Vectorised NumPy implementations of the TI/mTI modulation-amplitude envelope from Grossman et al. (2017), extended to an arbitrary even number of electrode pairs (mTI).

Public API

get_TI_vectors TI modulation-amplitude vectors for a single electrode pair (K=1). get_mTI_vectors Modulation-amplitude vectors for K >= 1 electrode pairs; the verified N>2 replacement for :func:get_nTI_vectors. get_TI_avg Direction-averaged modulation depth for K >= 1 electrode pairs. get_magnitude_am Direction-free magnitude-envelope AM, K >= 1 electrode pairs. get_nTI_vectors Deprecated. Recursive binary-tree N-field TI; not physically valid for N > 2. Delegates to :func:get_mTI_vectors.

Attribution

The K >= 2 envelope and the _fibonacci_sphere / _validate_field_list helpers originate from collaborator Larissa Albantakis's branch alba/mTI_testing and are ported here with attribution, not reimplemented. get_magnitude_am is likewise ported from that branch's _botzanowski_magnitude_am_components.

get_TI_vectors

get_TI_vectors(E1_org, E2_org)

Compute the TI modulation-amplitude vectors for two electric fields.

Sign-agnostic closed form (Hirata et al. 2024), equivalent to the original Grossman et al. (2017) preprocess-then-branch formulation but needs no magnitude-ordering swap or acute-angle sign flip of the inputs: TI = 2*min(|E1|,|E2|) (as a vector, sign-corrected) when min(|E1|,|E2|) <= sqrt(|E1.E2|), else TI = 2 * that same sign-corrected vector's component perpendicular to h, whichever of E1-E2/E1+E2 has the smaller norm.

Parameters

E1_org, E2_org : np.ndarray, shape (N, 3) Electric field vectors [V/m] from the two electrode pairs.

Returns

np.ndarray, shape (N, 3) TI vectors [V/m]: direction/magnitude of max envelope modulation.

References

Grossman, N. et al. (2017). Cell, 169(6), 1029-1041. Hirata, A. et al. (2024). Computers in Biology and Medicine, 178, 108697.

Source code in tit/calc.py

def get_TI_vectors(E1_org, E2_org):
    """Compute the TI modulation-amplitude vectors for two electric fields.

    Sign-agnostic closed form (Hirata et al. 2024), equivalent to the
    original Grossman et al. (2017) preprocess-then-branch formulation but
    needs no magnitude-ordering swap or acute-angle sign flip of the
    inputs: ``TI = 2*min(|E1|,|E2|)`` (as a vector, sign-corrected) when
    ``min(|E1|,|E2|) <= sqrt(|E1.E2|)``, else ``TI = 2 *`` that same
    sign-corrected vector's component perpendicular to ``h``, whichever of
    ``E1-E2``/``E1+E2`` has the smaller norm.

    Parameters
    ----------
    E1_org, E2_org : np.ndarray, shape (N, 3)
        Electric field vectors [V/m] from the two electrode pairs.

    Returns
    -------
    np.ndarray, shape (N, 3)
        TI vectors [V/m]: direction/magnitude of max envelope modulation.

    References
    ----------
    Grossman, N. et al. (2017). Cell, 169(6), 1029-1041.
    Hirata, A. et al. (2024). Computers in Biology and Medicine, 178, 108697.
    """
    assert E1_org.shape == E2_org.shape, "E1 and E2 must have same shape"
    assert E1_org.shape[1] == 3, "Vectors must be 3D"

    E1 = E1_org
    E2 = E2_org

    normE1 = np.linalg.norm(E1, axis=1)
    normE2 = np.linalg.norm(E2, axis=1)
    dot = np.sum(E1 * E2, axis=1)

    min_norm = np.minimum(normE1, normE2)
    regime1_mask = min_norm <= np.sqrt(np.abs(dot))

    # The smaller-magnitude field, sign-corrected for an acute angle; ties
    # go to E2, matching the strict '>' swap convention of the original
    # preprocess-then-branch form.
    use_E1_as_small = normE1 < normE2
    Es_raw = np.where(use_E1_as_small[:, None], E1, E2)
    sign = np.where(dot < 0, -1.0, 1.0)
    Es = sign[:, None] * Es_raw

    TI_vectors = np.zeros_like(E1)
    TI_vectors[regime1_mask] = 2.0 * Es[regime1_mask]

    # Regime 2 (oblique): TI = 2 * (Es perpendicular to h), h = whichever
    # of (E1-E2), (E1+E2) has the smaller norm.
    regime2_mask = ~regime1_mask
    if np.any(regime2_mask):
        a2 = E1[regime2_mask]
        b2 = E2[regime2_mask]
        dot2 = dot[regime2_mask]
        h_minus = a2 - b2
        h_plus = a2 + b2
        # Decide from sign(dot) directly rather than comparing the two
        # norms: |E1-E2|^2 - |E1+E2|^2 == -4*dot, so the sign is exact,
        # while for near-orthogonal fields (dot tiny relative to |E1|,
        # |E2|) the two norms round to the same float64 and the
        # comparison loses the sign to cancellation.
        use_minus = dot2 >= 0
        h = np.where(use_minus[:, None], h_minus, h_plus)
        h_norm = np.linalg.norm(h, axis=1)
        h_norm_safe = np.where(h_norm == 0, 1.0, h_norm)
        e_h = h / h_norm_safe[:, None]

        Es_r2 = Es[regime2_mask]
        Es_parallel_component = np.sum(Es_r2 * e_h, axis=1)[:, None] * e_h
        Es_perp = Es_r2 - Es_parallel_component
        TI_vectors[regime2_mask] = 2.0 * Es_perp

    return TI_vectors

get_mTI_vectors

get_mTI_vectors(fields, channels=None, psi=None)

Compute mTI modulation-amplitude vectors for K >= 1 electrode pairs.

fields is [E_1a, E_1b, ..., E_Ka, E_Kb], 2K arrays of shape (N, 3), paired positionally into K channels by default. Pass channels to override the grouping -- e.g. many electrode pairs sharing just two carriers (Lee et al. 2022) becomes one channel of summed fields. K=1 dispatches exactly to :func:get_TI_vectors; K>=2 returns best_direction * md from the verified :func:_mti_modulation_depth envelope. hf_peak/hf_sar (:mod:tit.fields) are unaffected by channels: they always sum over every carrier field.

Parameters

fields : list of np.ndarray, each shape (N, 3) Carrier field vectors, referenced by index from channels. channels : sequence of (group_a, group_b), or None Per-channel index groups into fields (see :func:_resolve_channels); None is consecutive pairing, identical to today's behaviour. psi : array-like, shape (K,), or None Per-pair envelope phase offset (radians); None means phase-aligned pairs (psi_k=0), the standard case. Ignored at K=1 (phase-invariant there).

Returns

np.ndarray, shape (N, 3) Modulation-amplitude vectors [V/m]; norm is the modulation depth.

Raises

ValueError Invalid fields, channels, or psi; see :func:_resolve_channels and :func:_validate_psi.

References

Grossman, N. et al. (2017). Cell, 169(6), 1029-1041 (K=1 closed form).

Source code in tit/calc.py

def get_mTI_vectors(fields, channels=None, psi=None):
    """Compute mTI modulation-amplitude vectors for K >= 1 electrode pairs.

    ``fields`` is ``[E_1a, E_1b, ..., E_Ka, E_Kb]``, 2K arrays of shape
    ``(N, 3)``, paired positionally into K channels by default. Pass
    ``channels`` to override the grouping -- e.g. many electrode pairs
    sharing just two carriers (Lee et al. 2022) becomes one channel of
    summed fields. K=1 dispatches exactly to :func:`get_TI_vectors`; K>=2
    returns ``best_direction * md`` from the verified
    :func:`_mti_modulation_depth` envelope. ``hf_peak``/``hf_sar``
    (:mod:`tit.fields`) are unaffected by ``channels``: they always sum
    over every carrier field.

    Parameters
    ----------
    fields : list of np.ndarray, each shape (N, 3)
        Carrier field vectors, referenced by index from ``channels``.
    channels : sequence of (group_a, group_b), or None
        Per-channel index groups into ``fields`` (see
        :func:`_resolve_channels`); ``None`` is consecutive pairing,
        identical to today's behaviour.
    psi : array-like, shape (K,), or None
        Per-pair envelope phase offset (radians); ``None`` means
        phase-aligned pairs (``psi_k=0``), the standard case. Ignored
        at K=1 (phase-invariant there).

    Returns
    -------
    np.ndarray, shape (N, 3)
        Modulation-amplitude vectors [V/m]; norm is the modulation depth.

    Raises
    ------
    ValueError
        Invalid ``fields``, ``channels``, or ``psi``; see
        :func:`_resolve_channels` and :func:`_validate_psi`.

    References
    ----------
    Grossman, N. et al. (2017). Cell, 169(6), 1029-1041 (K=1 closed form).
    """
    arrs = _resolve_channels(fields, channels)
    n_pairs = len(arrs) // 2
    _validate_psi(psi, n_pairs)

    if n_pairs == 1:
        return get_TI_vectors(arrs[0], arrs[1])

    result = _mti_modulation_depth(arrs, psi=psi)
    return result["best_direction"] * result["md"][:, None]

get_TI_avg

get_TI_avg(fields, channels=None, psi=None)

Direction-averaged modulation depth for K >= 1 electrode pairs.

TI_max (:func:get_mTI_vectors) maximizes the envelope over direction -- a best case for a neuron aligned with the optimal axis. TI_avg instead averages the same coarse Fibonacci-sphere sweep over all sampled directions, giving what a randomly-oriented neuron sees on average. Local refinement (accurate for a single best direction only) is skipped as irrelevant to an average.

Parameters

fields : list of np.ndarray, each shape (N, 3) Carrier field vectors, referenced by index from channels. channels : sequence of (group_a, group_b), or None Per-channel index groups into fields; see :func:get_mTI_vectors and :func:_resolve_channels. psi : array-like, shape (K,), or None Per-pair envelope phase offset (radians); see :func:get_mTI_vectors.

Returns

np.ndarray, shape (N,) Modulation depth [V/m], averaged over sampled directions.

Source code in tit/calc.py

def get_TI_avg(fields, channels=None, psi=None):
    """Direction-averaged modulation depth for K >= 1 electrode pairs.

    ``TI_max`` (:func:`get_mTI_vectors`) maximizes the envelope over
    direction -- a best case for a neuron aligned with the optimal axis.
    ``TI_avg`` instead averages the same coarse Fibonacci-sphere sweep
    over all sampled directions, giving what a randomly-oriented neuron
    sees on average. Local refinement (accurate for a single best
    direction only) is skipped as irrelevant to an average.

    Parameters
    ----------
    fields : list of np.ndarray, each shape (N, 3)
        Carrier field vectors, referenced by index from ``channels``.
    channels : sequence of (group_a, group_b), or None
        Per-channel index groups into ``fields``; see
        :func:`get_mTI_vectors` and :func:`_resolve_channels`.
    psi : array-like, shape (K,), or None
        Per-pair envelope phase offset (radians); see
        :func:`get_mTI_vectors`.

    Returns
    -------
    np.ndarray, shape (N,)
        Modulation depth [V/m], averaged over sampled directions.
    """
    arrs = _resolve_channels(fields, channels)
    n_pairs = len(arrs) // 2
    psi_arr = _validate_psi(psi, n_pairs)
    return _mti_modulation_depth_avg(arrs, psi_arr)

get_magnitude_am

get_magnitude_am(fields)

Direction-free amplitude-modulation envelope of ||E(t)||.

Not the direction-maximized modulation depth from :func:get_mTI_vectors -- the AM envelope of the field magnitude itself, no direction search, computed on the full 3-vectors: P = 0.5*sum_i ||E_i||^2, Q = |sum_k E_ka . E_kb| (3D dot products per pair), result = sqrt(2*(P+Q)) - sqrt(2*max(P-Q, 0)). At K=1 this reduces to abs(|E1+E2| - |E1-E2|).

Parameters

fields : list of np.ndarray, each shape (N, 3) Field vectors for 2K sub-channels (K electrode pairs), K >= 1.

Returns

np.ndarray, shape (N,) Magnitude-AM envelope [V/m].

See Also

get_mTI_vectors : Direction-maximized modulation-amplitude vectors -- a different quantity from this magnitude envelope.

References

Botzanowski, B. et al. (2025). Bioelectronic Medicine, 11(1), 7.

Source code in tit/calc.py

def get_magnitude_am(fields):
    """Direction-free amplitude-modulation envelope of ``||E(t)||``.

    Not the direction-maximized modulation depth from
    :func:`get_mTI_vectors` -- the AM envelope of the field *magnitude*
    itself, no direction search, computed on the full 3-vectors:
    ``P = 0.5*sum_i ||E_i||^2``, ``Q = |sum_k E_ka . E_kb|`` (3D dot
    products per pair), result ``= sqrt(2*(P+Q)) - sqrt(2*max(P-Q, 0))``.
    At K=1 this reduces to ``abs(|E1+E2| - |E1-E2|)``.

    Parameters
    ----------
    fields : list of np.ndarray, each shape (N, 3)
        Field vectors for 2K sub-channels (K electrode pairs), K >= 1.

    Returns
    -------
    np.ndarray, shape (N,)
        Magnitude-AM envelope [V/m].

    See Also
    --------
    get_mTI_vectors : Direction-maximized modulation-amplitude vectors --
        a different quantity from this magnitude envelope.

    References
    ----------
    Botzanowski, B. et al. (2025). Bioelectronic Medicine, 11(1), 7.
    """
    arrs = _validate_field_list(fields)
    n_pairs = len(arrs) // 2

    P = np.zeros(arrs[0].shape[0], dtype=np.float64)
    for e in arrs:
        P += np.sum(e * e, axis=1)
    P *= 0.5

    dot_sum = np.zeros(arrs[0].shape[0], dtype=np.float64)
    for k in range(n_pairs):
        dot_sum += np.sum(arrs[2 * k] * arrs[2 * k + 1], axis=1)
    Q = np.abs(dot_sum)

    env_max = np.sqrt(2.0 * np.maximum(P + Q, 0.0))
    env_min = np.sqrt(2.0 * np.maximum(P - Q, 0.0))
    return env_max - env_min

get_nTI_vectors

get_nTI_vectors(fields)

Deprecated: recursive binary-tree N-field TI. Use :func:get_mTI_vectors.

This paired fields via TI(TI(E1,E2), TI(E3,E4), ...), feeding already-modulated envelope vectors back into :func:get_TI_vectors -- a formula derived only for two carrier fields. Measured against the verified :func:_mti_modulation_depth envelope on random fields: signed mean error +38.6% (range -90% to +416%) at N=4, +103% at N=8.

Parameters

fields : list of np.ndarray, each shape (N, 3)

Returns

np.ndarray, shape (N, 3)

Source code in tit/calc.py

def get_nTI_vectors(fields):
    """Deprecated: recursive binary-tree N-field TI. Use :func:`get_mTI_vectors`.

    This paired fields via ``TI(TI(E1,E2), TI(E3,E4), ...)``, feeding
    already-modulated envelope vectors back into :func:`get_TI_vectors` --
    a formula derived only for two carrier fields. Measured against the
    verified :func:`_mti_modulation_depth` envelope on random fields: signed
    mean error +38.6% (range -90% to +416%) at N=4, +103% at N=8.

    Parameters
    ----------
    fields : list of np.ndarray, each shape (N, 3)

    Returns
    -------
    np.ndarray, shape (N, 3)
    """
    warnings.warn(
        "get_nTI_vectors is deprecated and physically invalid for N>2 "
        "(measured +38.6% mean error at N=4 vs. the verified envelope); "
        "use get_mTI_vectors instead.",
        DeprecationWarning,
        stacklevel=2,
    )
    return get_mTI_vectors(fields)