Skip to content

calc

tit.calc

Shared TI Field Calculation Utilities Used by optimization tools

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.

Source code in tit/calc.py 
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 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
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
127
128
129
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.
    """
    # 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 mTI vectors for N E-fields using recursive binary-tree pairwise TI.

N must be even. Fields are paired sequentially: (E1,E2), (E3,E4), etc. Then intermediate TI results are paired recursively until one 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.

Source code in tit/calc.py 
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
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
def get_nTI_vectors(fields):
    """Compute mTI vectors for N E-fields using recursive binary-tree pairwise TI.

    N must be even. Fields are paired sequentially: (E1,E2), (E3,E4), etc.
    Then intermediate TI results are paired recursively until one 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.
    """
    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 channel E-fields.

This 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

Source code in tit/calc.py 
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
def get_mTI_vectors(E1_org, E2_org, E3_org, E4_org):
    """
    Calculate multi-temporal interference (mTI) vectors from four channel E-fields.

    This 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
    """
    # 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