calc

tit.calc

Temporal interference field calculation utilities.

Vectorised NumPy implementations of the TI modulation-amplitude algorithm from Grossman et al. (2017), extended to multi-channel (mTI) configurations. Used by the simulation engine and optimisation tools.

Public API

get_TI_vectors Compute TI modulation-amplitude vectors for a single electrode pair. get_nTI_vectors Generalised N-channel TI via recursive binary-tree pairing. get_mTI_vectors 4-channel mTI (convenience wrapper around :func:get_TI_vectors).

get_TI_vectors

get_TI_vectors(E1_org, E2_org)

Calculate the temporal interference (TI) modulation amplitude vectors.

This function implements the Grossman et al. 2017 algorithm for computing TI vectors that represent both the direction and magnitude of maximum modulation amplitude when two sinusoidal electric fields interfere.

PHYSICAL INTERPRETATION: When two electric fields E1(t) = E1cos(2πf1t) and E2(t) = E2cos(2πf2t) with slightly different frequencies are applied simultaneously, they create a beating pattern. The TI vector indicates: - DIRECTION: Spatial direction of maximum envelope modulation - MAGNITUDE: Maximum envelope amplitude = 2 * effective_amplitude

ALGORITHM (Grossman et al. 2017): 1. Preprocessing: Ensure |E1| ≥ |E2| and acute angle α < π/2 2. Regime selection based on geometric relationship: - Regime 1 (parallel): |E2| ≤ |E1|cos(α) → TI = 2E2 - Regime 2 (oblique): |E2| > |E1|cos(α) → TI = 2E2_perpendicular_to_h where h = E1 - E2

Parameters

E1_org : np.ndarray, shape (N, 3) Electric field vectors from electrode pair 1 [V/m] E2_org : np.ndarray, shape (N, 3) Electric field vectors from electrode pair 2 [V/m]

Returns

TI_vectors : np.ndarray, shape (N, 3) TI modulation amplitude vectors [V/m] Direction: Maximum modulation direction Magnitude: Maximum envelope amplitude

References

Grossman, N. et al. (2017). Noninvasive Deep Brain Stimulation via Temporally Interfering Electric Fields. Cell, 169(6), 1029-1041.

See Also

get_mTI_vectors : 4-channel mTI (two pairs of pairs). get_nTI_vectors : Generalised N-channel recursive TI.

Source code in tit/calc.py

def get_TI_vectors(E1_org, E2_org):
    """
    Calculate the temporal interference (TI) modulation amplitude vectors.

    This function implements the Grossman et al. 2017 algorithm for computing
    TI vectors that represent both the direction and magnitude of maximum
    modulation amplitude when two sinusoidal electric fields interfere.

    PHYSICAL INTERPRETATION:
    When two electric fields E1(t) = E1*cos(2πf1*t) and E2(t) = E2*cos(2πf2*t)
    with slightly different frequencies are applied simultaneously, they create
    a beating pattern. The TI vector indicates:
    - DIRECTION: Spatial direction of maximum envelope modulation
    - MAGNITUDE: Maximum envelope amplitude = 2 * effective_amplitude

    ALGORITHM (Grossman et al. 2017):
    1. Preprocessing: Ensure |E1| ≥ |E2| and acute angle α < π/2
    2. Regime selection based on geometric relationship:
       - Regime 1 (parallel): |E2| ≤ |E1|cos(α) → TI = 2*E2
       - Regime 2 (oblique): |E2| > |E1|cos(α) → TI = 2*E2_perpendicular_to_h
       where h = E1 - E2

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

    Returns
    -------
    TI_vectors : np.ndarray, shape (N, 3)
        TI modulation amplitude vectors [V/m]
        Direction: Maximum modulation direction
        Magnitude: Maximum envelope amplitude

    References
    ----------
    Grossman, N. et al. (2017). Noninvasive Deep Brain Stimulation via
    Temporally Interfering Electric Fields. Cell, 169(6), 1029-1041.

    See Also
    --------
    get_mTI_vectors : 4-channel mTI (two pairs of pairs).
    get_nTI_vectors : Generalised N-channel recursive TI.
    """
    # Input validation
    assert E1_org.shape == E2_org.shape, "E1 and E2 must have same shape"
    assert E1_org.shape[1] == 3, "Vectors must be 3D"

    # Work with copies to avoid modifying input arrays
    E1 = E1_org.copy()
    E2 = E2_org.copy()

    # =================================================================
    # PREPROCESSING STEP 1: Magnitude ordering |E1| ≥ |E2|
    # =================================================================
    # Ensures consistency by always treating E1 as the "stronger" field
    # This simplifies the subsequent regime analysis
    idx_swap = np.linalg.norm(E2, axis=1) > np.linalg.norm(E1, axis=1)
    E1[idx_swap], E2[idx_swap] = E2[idx_swap], E1_org[idx_swap]

    # =================================================================
    # PREPROCESSING STEP 2: Acute angle constraint α < π/2
    # =================================================================
    # Ensures constructive interference by flipping E2 if dot product < 0
    # This avoids destructive interference scenarios
    idx_flip = np.sum(E1 * E2, axis=1) < 0
    E2[idx_flip] = -E2[idx_flip]

    # =================================================================
    # GEOMETRIC PARAMETERS CALCULATION
    # =================================================================
    # Calculate field magnitudes and angle between vectors
    normE1 = np.linalg.norm(E1, axis=1)
    normE2 = np.linalg.norm(E2, axis=1)

    # Safe cosine calculation to avoid division by zero and numerical errors
    denom = normE1 * normE2
    denom[denom == 0] = 1.0  # Prevent division by zero
    cosalpha = np.clip(np.sum(E1 * E2, axis=1) / denom, -1.0, 1.0)

    # =================================================================
    # REGIME SELECTION CRITERION
    # =================================================================
    # Critical condition from Grossman 2017: |E2| ≤ |E1| * cos(α)
    # This determines whether E2 is "small" relative to E1's projection
    regime1_mask = normE2 <= normE1 * cosalpha

    # Initialize output array
    TI_vectors = np.zeros_like(E1)

    # =================================================================
    # REGIME 1: PARALLEL ALIGNMENT (|E2| ≤ |E1| cos(α))
    # =================================================================
    # Physical interpretation: E2 is effectively "contained" within E1's projection
    # The TI amplitude is determined entirely by E2's magnitude and direction
    # Formula: TI = 2 * E2
    TI_vectors[regime1_mask] = 2.0 * E2[regime1_mask]

    # =================================================================
    # REGIME 2: OBLIQUE CONFIGURATION (|E2| > |E1| cos(α))
    # =================================================================
    # Physical interpretation: E2 has significant perpendicular component to E1
    # The TI is determined by the component of E2 perpendicular to h = E1 - E2
    # Formula: TI = 2 * E2_perpendicular_to_h
    regime2_mask = ~regime1_mask
    if np.any(regime2_mask):
        # Calculate difference vector h = E1 - E2
        h = E1[regime2_mask] - E2[regime2_mask]
        h_norm = np.linalg.norm(h, axis=1)

        # Handle degenerate case (h = 0) by setting unit norm
        h_norm[h_norm == 0] = 1.0
        e_h = h / h_norm[:, None]  # Unit vector along h

        # Project E2 onto h, then subtract to get perpendicular component
        # E2_perp = E2 - proj_h(E2) = E2 - (E2·ĥ)ĥ
        E2_parallel_component = np.sum(E2[regime2_mask] * e_h, axis=1)[:, None] * e_h
        E2_perp = E2[regime2_mask] - E2_parallel_component

        # The TI vector in regime 2 is twice the perpendicular component
        TI_vectors[regime2_mask] = 2.0 * E2_perp

    return TI_vectors

get_nTI_vectors

get_nTI_vectors(fields)

Compute TI vectors for N E-fields using recursive binary-tree pairing.

N must be even. Fields are paired sequentially — (E1, E2), (E3, E4), etc. — then intermediate TI results are paired recursively until a single result remains.

For 2 fields: TI(E1, E2) For 4 fields: TI(TI(E1,E2), TI(E3,E4)) For 6 fields: TI(TI(TI(E1,E2), TI(E3,E4)), TI(E5,E6)) For 8 fields: TI(TI(TI(E1,E2), TI(E3,E4)), TI(TI(E5,E6), TI(E7,E8)))

Parameters

fields : list of np.ndarray, each shape (N, 3) Electric field vectors, one per electrode pair.

Returns

result : np.ndarray, shape (N, 3) Combined TI modulation amplitude vectors.

Raises

ValueError If number of fields is not even or less than 2.

See Also

get_TI_vectors : Core 2-field TI calculation. get_mTI_vectors : 4-channel convenience wrapper.

Source code in tit/calc.py

def get_nTI_vectors(fields):
    """Compute TI vectors for *N* E-fields using recursive binary-tree pairing.

    *N* must be even.  Fields are paired sequentially — ``(E1, E2)``,
    ``(E3, E4)``, etc. — then intermediate TI results are paired recursively
    until a single result remains.

    For 2 fields: ``TI(E1, E2)``
    For 4 fields: ``TI(TI(E1,E2), TI(E3,E4))``
    For 6 fields: ``TI(TI(TI(E1,E2), TI(E3,E4)), TI(E5,E6))``
    For 8 fields: ``TI(TI(TI(E1,E2), TI(E3,E4)), TI(TI(E5,E6), TI(E7,E8)))``

    Parameters
    ----------
    fields : list of np.ndarray, each shape (N, 3)
        Electric field vectors, one per electrode pair.

    Returns
    -------
    result : np.ndarray, shape (N, 3)
        Combined TI modulation amplitude vectors.

    Raises
    ------
    ValueError
        If number of fields is not even or less than 2.

    See Also
    --------
    get_TI_vectors : Core 2-field TI calculation.
    get_mTI_vectors : 4-channel convenience wrapper.
    """
    n = len(fields)
    if n < 2 or n % 2 != 0:
        raise ValueError(
            f"get_nTI_vectors requires an even number of fields >= 2, got {n}"
        )

    # First round: pair adjacent fields
    current = []
    for i in range(0, n, 2):
        current.append(get_TI_vectors(fields[i], fields[i + 1]))

    # Recursive rounds: pair results until one remains
    while len(current) > 1:
        next_round = []
        i = 0
        while i + 1 < len(current):
            next_round.append(get_TI_vectors(current[i], current[i + 1]))
            i += 2
        # Odd element carries forward
        if i < len(current):
            next_round.append(current[i])
        current = next_round

    return current[0]

get_mTI_vectors

get_mTI_vectors(E1_org, E2_org, E3_org, E4_org)

Calculate multi-temporal interference (mTI) vectors from four E-fields.

Computes TI between channels 1 and 2 to get TI_A, TI between channels 3 and 4 to get TI_B, and finally TI between TI_A and TI_B to produce the mTI vector field.

Parameters

E1_org : np.ndarray, shape (N, 3) Electric field vectors for channel 1. E2_org : np.ndarray, shape (N, 3) Electric field vectors for channel 2. E3_org : np.ndarray, shape (N, 3) Electric field vectors for channel 3. E4_org : np.ndarray, shape (N, 3) Electric field vectors for channel 4.

Returns

mTI_vectors : np.ndarray, shape (N, 3) Multi-TI modulation amplitude vectors [V/m].

Raises

ValueError If any input array does not have shape (N, 3) or if the shapes are not identical.

See Also

get_TI_vectors : Core 2-field TI calculation. get_nTI_vectors : Generalised N-channel recursive TI.

Source code in tit/calc.py

def get_mTI_vectors(E1_org, E2_org, E3_org, E4_org):
    """Calculate multi-temporal interference (mTI) vectors from four E-fields.

    Computes TI between channels 1 and 2 to get ``TI_A``, TI between
    channels 3 and 4 to get ``TI_B``, and finally TI between ``TI_A`` and
    ``TI_B`` to produce the mTI vector field.

    Parameters
    ----------
    E1_org : np.ndarray, shape (N, 3)
        Electric field vectors for channel 1.
    E2_org : np.ndarray, shape (N, 3)
        Electric field vectors for channel 2.
    E3_org : np.ndarray, shape (N, 3)
        Electric field vectors for channel 3.
    E4_org : np.ndarray, shape (N, 3)
        Electric field vectors for channel 4.

    Returns
    -------
    mTI_vectors : np.ndarray, shape (N, 3)
        Multi-TI modulation amplitude vectors [V/m].

    Raises
    ------
    ValueError
        If any input array does not have shape ``(N, 3)`` or if the shapes
        are not identical.

    See Also
    --------
    get_TI_vectors : Core 2-field TI calculation.
    get_nTI_vectors : Generalised N-channel recursive TI.
    """
    # Validate shapes
    for i, arr in enumerate([E1_org, E2_org, E3_org, E4_org], start=1):
        if arr.ndim != 2 or arr.shape[1] != 3:
            raise ValueError(f"E{i}_org must have shape (N, 3), got {arr.shape}")

    if not (E1_org.shape == E2_org.shape == E3_org.shape == E4_org.shape):
        raise ValueError(
            "All input arrays must have identical shapes. "
            f"Got: {[E1_org.shape, E2_org.shape, E3_org.shape, E4_org.shape]}"
        )

    # Step 1: TI between (E1, E2)
    TI_A = get_TI_vectors(E1_org, E2_org)

    # Step 2: TI between (E3, E4)
    TI_B = get_TI_vectors(E3_org, E4_org)

    # Step 3: TI between (TI_A, TI_B) → mTI
    mTI_vectors = get_TI_vectors(TI_A, TI_B)

    return mTI_vectors