pynibs.optimization package¶

Routines to compute optimal sets of electric fields, either by optimizing

single-channel TMS coil placements, or
input currents for multichannel TMS coil arrays.

Used in [1]

[1]

Numssen, O., Martin, C. W., Worbs, T., Thielscher, A., Weise, K., & Knösche, T. R. (2025). Optimizing and assessing multichannel TMS focality. Cold Spring Harbor Laboratory. 10.1101/2025.09.19.677136

Submodules¶

pynibs.optimization.coil_opt module¶

The coil_opt.py module provides functions for optimizing coil placements in Transcranial Magnetic Stimulation (TMS) based on given electric field matrices and fMRI statistics. It includes functions for identifying optimal coil placement regions, calculating the gain map between a reference and an optimized sequence of electric fields, and performing virtual online optimization to determine the congruence factor.

The module includes the following functions:

rowvec_diff(): This function returns the coil configuration out of all available configurations exhibiting the highest minimum difference to the already selected configurations.
get_optimal_coil_positions(): This function determines a set of optimal coil positions for TMS regression analysis.
online_optimization(): This function performs virtual online optimization to determine the congruence factor. After an initial set of coil positions, the algorithm iteratively optimizes the next coil position based on the virtually measured MEP data.
calc_opt_gain_map(): This function calculates the gain map between a reference e_matrix (e.g. from random sampling) and an optimized sequence of electric fields for mapping.
optimal_coilplacement_region(): This function identifies the optimal coil placement regions based on given electric field (E-field) matrices and fMRI statistics.

Each function in this module is documented with docstrings providing more detailed information about its purpose, parameters, and return values.

This module is primarily used for handling and optimizing coil placements in TMS studies.

pynibs.optimization.coil_opt.calc_opt_gain_map(e_matrix_ref, e_matrix_opt, points, con, fn_out=None, threshold=0.75)¶

Calculates the gain map between a reference e_matrix (e.g. from random sampling) and an optimized sequence of electric fields for mapping.

Parameters:

e_matrix_ref (np.ndarray of float [n_ref, n_ele_roi]) – Electric field matrix of reference simulations. E-fields in ROI are in the rows. (n_ref does not have to match n_opt)
e_matrix_opt (np.ndarray of float [n_opt, n_ele_roi]) – Electric field matrix of optimal simulations. E-fields in ROI are in the rows. (n_ref does not have to match n_opt)
points (np.ndarray of float [n_points_roi, 3]) – Node coordinates of the ROI.
con (np.ndarray of float [n_ele_roi, 3]) – Connectivity matrix of ROI surface.
fn_out (str, optional, default: None) – Filename of .hdf5 and .xdmf file for plots with paraview. (Without file extension)
threshold (float, optional, default: 0.75) – Threshold of correlation the focality is quantified by area.

Returns:

focality_ref (np.ndarray of float [n_ele_roi, 3]) – Focality measure (area) of PSF in each element > threshold for reference case.
focality_opt (np.ndarray of float [n_ele_roi, 3]) – Focality measure (area) of PSF in each element > threshold for optimal case.
focality_dif (np.ndarray of float [n_ele_roi, 3]) – Difference between focality_opt and focality_ref quantifying absolute gain in mm^2. Values < 0 : Optimal solution has smaller PSF than reference if hotspot would be in this element. Values > 0 : Optimal solution has larger PSF than reference if hotspot would be in this element.
<file> (.hdf5 and .xmdf) – Geo and data files for visualization in paraview.

pynibs.optimization.coil_opt.get_optimal_coil_positions(e_matrix, criterion, n_stim, ele_idx_1=None, ele_idx_2=None, fn_out_hdf5=None, n_cpu=4, zap_idx_opt=None, regression_cmap=None, regression_fit_parameters=None, metrics_weights=None, overwrite=True, verbose=True, fn_coilpos_hdf5=None, start_zap_idx=-1, fim_fit_fun=None, fim_p2p_amps=None, fim_didt_list=None, fim_rmt_mso=None, fim_mso_didt_conversion_factor=1.43, fim_visited_positions_e_mat=None, fim_regression_n_refit=10, fim_debug_screenshot_dir_fn=None, fim_roi_pts=None, fim_roi_tris=None, fim_use_gpu=False)¶

Determine set of optimal coil positions for TMS regression analysis.

Parameters:

e_matrix (np.ndarray of float) – (n_stim, n_ele) Matrix containing the electric field values in the ROI.
criterion (str) – Optimization criterion: * “mc_cols”: Minimization of mutual coherence between columns * “mc_rows”: Minimization of mutual coherence between rows * “svd”: Minimization of condition number * “dist”: Equal distant sampling * “dist_svd”: Minimization of condition number and equidistant sampling * “dist_mc_cols”: Minimization of mutual coherence between columns and equidistant sampling * “dist_mc_rows”: Minimization of mutual coherence between rows and equidistant sampling * “coverage”: Maximizes the electric field coverage * “variability”: Maximizes variability between elements
n_stim (int) – Maximum number of stimulations.
ele_idx_1 (np.ndarray of int, optional) – Element indices the first optimization goal is performed for, If None, all elements are considered.
ele_idx_2 (np.ndarray of int, optional) – Element indices the first optimization goal is performed for. If None, all elements are considered.
n_cpu (int) – Number of threads.
fn_out_hdf5 (str, optional) –
Returns the list of optimal zap indices if fn_out_hdf5 is None, otherwise, save the results in .hdf5 file. Filename of output .hdf5 file where the zap index lists are saved in subfolder “zap_index_lists”
- ”zap_index_lists/0”: [213]
- ”zap_index_lists/1”: [213, 5]
- etc
zap_idx_opt (list of int, optional) – List of already selected optimal coil positions (those are ignored in the optimization and will not be picked again).
fim_fit_fun (function object) – Function object defined in interval [0, 1] (only needed for fim optimization).
regression_fit_parameters (dict [n_ele], optional, optional) – The parameter estimates that should be used for the FIM optimization (whole ROI). The keys are the parameter names of fun (only needed for fim and dist optimization).
regression_cmap (np.ndarray of float [n_ele], optional, optional) – Congruence factor in each ROI element. Used to weight fim and dist optimization (only needed for fim and dist optimization).
metrics_weights (list of float [2], default: [0.5, 0.5]) – Weights of optimization criteria in case of multiple goal functions (e_all_coil_pos.g. fim_svd). Higher weight means higher importance for the respective criteria. By default, both optimization criteria are weighted equally [0.5, 0.5].
overwrite (bool, default: True) – Overwrite existing solutions or read existing hdf5 file and continue optimization.
verbose (bool, default: True) – Print output messages.
fn_coilpos_hdf5 (str) – File containing the corresponding coil positions and orientations (centers, m0, m1, m2).
start_zap_idx (int, default: 0) – First zap index to start greedy search.
fim_didt_list (np.ndarray[float], (len(zap_idx_op))) – List of realized dI/dt of each of the already stimulated coil configurations in ‘zap_idx_opt’. Not required for any other metric than FIM.
fim_rmt_mso (int) – Resting motor threshold used as the lower boundary of the FIM optimal e-field scaling. Unit in %MSO Not required for any other metric than FIM.
fim_mso_didt_conversion_factor (float, default: 1.43) – Factor to convert between realized current (dI/dt) and percentage of maximum stimulator output (%MSO). Defaults to 1.43 describing the factor of a MagVenture Pro with an MCF-B65 coil. Not required for any other metric than FIM.
fim_visited_positions_e_mat (np.ndarray[float], (len(zap_idx_opt], n_ele), optional) – The efield matrix computed using the actually approached coil configurations. This is intended to make the FIM method even more precise by taking into account slight deviations in the approached coil configuration instead of relying on the static pre-computed coil configurations. Not required for any other metric than FIM.
fim_p2p_amps (np.ndarray[float], (len(zap_idx_opt))) – EMG peak to peak amplitudes associated with the already collected (optimal) coil positions. Not required for any other metric than FIM.
fim_didt_list – Realized current (didt) in TMS coil (as returned by the stimulator). of the already collected (optimal) coil positions. Not required for any other metric than FIM.
fim_regression_n_refit (float) – Number of refits used in the mag(E)<>p2p regression during FIM optimization. Not required for any other metric than FIM.
fim_debug_screenshot_dir_fn (str) – String representation of the fully qualified path to a directory where a 3D rendering of the FIM optimal coil positions, ie coil positions that can reach the FIM optimal e-field strength at the current target hotspot given the MSO bounds [fim_rmt_mso, 100], should be saved. Not required for any other metric than FIM.
fim_roi_pts (np.ndarray, (n_points x 3)) – Points (vertices) of ROI surface mesh (where the congruence scores should be computed on). Not required for any other metric than FIM.
fim_roi_tris (np.ndarray, (n_tris x 3)) – Connectivity list of the ‘fim_roi_points’. Not required for any other metric than FIM.
fim_use_gpu (bool) – True: Use cupy and CUDA acceleration for the computation of the correlation matrix. False: Use Python multiprocessing for the computation of the corrleation matrix. Not required for any other metric than FIM.

Returns:

zap_idx_e_opt (list of int) – (n_stim) Optimal zap indices.
<File> .hdf5 file – Output file containing the zap index lists.

pynibs.optimization.coil_opt.online_optimization(fn_subject_hdf5, fn_roi_ss_indices_hdf5, fn_out_hdf5, fn_stimsites_hdf5, e_matrix, mep, mesh_idx, roi_idx, n_zaps_init=3, criterion_init='mc_rows', criterion='coverage', n_cpu=4, threshold=0.8, weights=None, eps0=0.01, eps0_dist=1, exponent=5, perc=99, n_refit=0, fun=None, verbose=True)¶

Performs virtual online optimization to determine the congruence factor. After an initial set of coil positions, the algorithm iteratively optimizes the next coil position based on the virtually measured MEP data.

Parameters:

fn_subject_hdf5 (str) – Filename of subject .hdf5 file.
fn_roi_ss_indices_hdf5 (str) – Filename of .hdf5 file containing the element indices of the subsampled ROI in f["roi_indices"].
e_matrix (np.ndarray of float) – (n_zaps, n_ele) Electric field matrix.
mep (np.ndarray of float) – (n_zaps) Motor evoked potentials for every stimulation.
fn_out_hdf5 (str) – Filename of .hdf5 output file containing the coil positions and the congruence factor maps for every iteration.
fn_stimsites_hdf5 (str) – Filename of the .hdf5 file containing the stimulation sites in “centers”, “m0”, “m1”, “m2”.
mesh_idx (int) – Mesh index.
roi_idx (int) – ROI index.
n_zaps_init (int, default: 3) – Number of initial samples optimized using optimization criterion specified in criterion_init.
criterion_init (str, default: "mc_rows") – Optimization criterion for which the initial samples are optimized (e.g. “mc_rows”, “svd”, …).
criterion (str, default: "coverage") – Optimization criterion for which the online optimization is performed (e.g. “coverage”, “mc_rows”, “svd”, …).
n_cpu (int, optional, dfault: 4) – Number of CPU cores to use.
threshold (float, default: 0.1) – Threshold between [0 … 1] of the maximal congruence factor. Elements where c > threshold * max(c) are included in the online optimization to select the next optimal coil position.
weights (list of float [2], default: [0.5, 0.5]) – Weights of optimization criteria in case of multiple goal functions (e.g. fim_svd). Higher weight means higher importance for the respective criteria. By default, both optimization criteria are weighted equally [0.5, 0.5].
eps0 (float, default: 0.01) – First error threshold to terminate the online optimization. The normalized root-mean-square deviation is calculated between the current and the previous solution. If the error is lower than eps0 for 3 times in a row, the online optimization terminates and returns the results.
eps0_dist (float, default: 1) – Second error threshold to terminate the online optimization. The geodesic distance in mm of the hotspot is calculated between the current and the previous solution. If the error is lower than eps0_dist for 3 times in a row, the online optimization terminates and returns the results.
exponent (float, default: 5) – Exponent the congruence factor map is scaled c**exponent.
perc (float, default: 99) – Percentile the congruence factor map is normalized (between 0 and 100).
n_refit (int, default: 0) – Number of refit iterations. No refit is applied if n_refit=0.
fun (function object, default: pynibs.linear) – Function to use to determine the congruence factor (e.g. pynibs.linear, pynibs.expio.sigmoid, …).
verbose (bool, default: True) – Plot output messages.

Returns:

Results output file containing the coil positions and the congruence factor maps for every iteration.

Return type:

<file> .hdf5 file

pynibs.optimization.coil_opt.optimal_coilplacement_region(e, fmri_stats, best_n=1, metric='dot', non_negative=False, pos_matrices=None, res_folder=None)¶

Identify the optimal coil placement regions based on given electric field (E-field) matrices and fMRI statistics.

Parameters:

e (np.ndarray) – The E-field matrix containing field strengths for various positions.
fmri_stats (ndarray) – fMRI statistics used to identify optimal stimulation regions. The values are normalized to the maximum value.
best_n (int, default: 1) – Number of top positions to identify.
metric (str, default: 'dot') –
The evaluation metric to use.
- ’dot’ calculates the dot product,
- ’cor’ computes the correlation coefficient.
non_negative (bool, optional) – If True, only considers non-negative fMRI statistics.
pos_matrices (ndarray, optional) – Position matrices for each coil position.
res_folder (str, optional) – Directory to save the results. If not specified, no files are saved.

Returns:

best_coil_id – N best coil ids, accoring to best_n

Return type:

int

Raises:

ValueError – If the specified metric is not recognized.

Notes

The function normalizes the input E-field matrix and fMRI statistics for evaluation. If res_folder is specified, visualizations of the results are saved to the directory. The function supports saving and writing TMS navigator instrument marker files using simnibs.localite.

pynibs.optimization.coil_opt.rowvec_diff(candidate_coil_idcs, selected_coil_idcs, efields_diff_mat)¶

Given a difference matrix (e.g. of row vectors/coil configurations) this function returns the coil configuration out of all available configurations exhibiting the highest minimum difference to the already selected configurations.

Parameters:

candidate_coil_idcs (np.ndarry[int]) – List of indices of coil configurations that are still available to pick for the optiized sequence.
selected_coil_idcs (np.ndarray[int]) – List of indices of coil configurations that have already been selected for the optimized sequence.
efields_diff_mat (np.ndarray[float], [n_coil,n_coil]) – Difference matrix, where each cell denotes the magnitude of the difference vector between two coil configurations (determined by row_idx,col_idx).

Returns:

coil_idx – index of coil configuration with maximal minimal difference to the set of already selected coil configurations.

Return type:

int

pynibs.optimization.multichannel module¶

Functions to optimize single channel currents for multichannel TMS arrays.

pynibs.optimization.multichannel.calc_e(l, j)¶

Calculate E-field. Multiplication of lead-field and coil currents.

Parameters:

l (np.ndarray) – [n_elm*3, n_coils] Flattened lead-field matrix. x,y,z direction in each point as consecutive rows.
j (np.ndarray) – [n_coils]

Returns:

E – [n_elm, 3] Resulting e-field for all elements in ROI.

Return type:

np.ndarray of dimenstion

pynibs.optimization.multichannel.calc_ll(l)¶

Calculate matrix multiplication of Lead-filed * transpose(Lead-field).

Parameters:: l (np.ndarray) – [n_elm*3, n_coils] Flattened lead-field matrix. x,y,z direction in each point as consecutive rows.
Returns:: LL – [n_coils, n_coils] Resulting square matrix LL.
Return type:: np.ndarray

pynibs.optimization.multichannel.calculate_strain(currents)¶

Computes planar mTMS strain from single channel currents.

Written by Mikael Laine, 2025.

Parameters:: currents (list) –
Returns:: pulse_strain
Return type:: float

pynibs.optimization.multichannel.calculate_triangle_area(p1, p2, p3)¶

Calculate the area of a triangle given its three vertices. v

Parameters:

p1 (np.ndarray) – (n, 3) Coordinates of the triangle’s vertices.
p2 (np.ndarray) – (n, 3) Coordinates of the triangle’s vertices.
p3 (np.ndarray) – (n, 3) Coordinates of the triangle’s vertices.

Returns:

area – (n, ) Area of the triangles.

Return type:

np.ndarray

pynibs.optimization.multichannel.constraint_min_current(x, n_stim, channels=(2, 3), min_current=0.7, frac=0.95, factor=100.0, verbose=False)¶

This function constrains the currents to be above a certain threshold. Should be used with the ‘SLSQP’ method of scipy.optimize.minimize as an inequality constraint.

Parameters:

x (np.ndarray of float) – (n_channel * n_stim_opt, ) Vector to scale channels for each.
n_stim (int) – Number of stimulations to compute score for.
channels (tuple of int, default: (2, 3)) – Channels to constrain.
min_current (float, default: 0.7) – Minimum current.
frac (float, default: 0.95) – Fraction of channels that must be within the bounds.
factor (float, default: 100.) – Factor to scale the constraint.
verbose (bool, default: False) – Print additional information

Returns:

res – >=0 if the constraint is met, <0 otherwise

Return type:

float

pynibs.optimization.multichannel.foc_opt(l, ll, k, target_stim=0.1, current_constraint=None, verbose=False)¶

Focality optimization routine. Minimizing goal function value via gradient descend of scipy.minimize(). The goal function is the squared local lead-field matrix at target position k. The optimization is done in the space of the hyperellipsoid cylinder.

Initially, an unconstrained optimization is performed.
In case current_constraints are set and the resulting currents exceed the constraints,
a constrained optimization is performed. If the constrained optimization fails, the function returns None for e, eps, j, and phi_opt.

l is np.vstack([E_chan1.flatten(), E_chan2.flatten(),...]).T with E_chan1.shape = (n_elms, 3)

ll is L.T @ L

Parameters:

l (np.ndarray of float) – [n_elm*3, n_coils] Flattened lead-field matrix. x,y,z direction in each point as consecutive rows.
ll (np.ndarray of float) – [n_coils, n_coils] Result of matrix multiplication of lead-field matrix * transpose(lead-field matrix)
k (int) – index of optimization target element in l.
target_stim (int, default = 80) – desired target stimulation. scaling value for unit current.
verbose (bool, default = False) – If True, displays optimization progress and debug information. This parameter controls the verbosity of the optimization process. If set to False, it suppresses the output.
current_constraint (int, optional) – If set, a current constraint is applied for the optimization, ensuring the resulting currents don’t exceed this value.

Returns:

e (np.ndarray) – [n_elm, 3] The resulting electric field from the optimization process.
eps (float) – The result of the goal function for the given phi.
j (np.ndarray) – [n_coils] The optimal currents for the given phi.
phi_opt (np.ndarray of dimension) – [n_coils] The resulting optimal angles of the hyperellipsoid.
constraint_opt (str) – The type of constraint optimization used. ‘foc_opt’ for unconstrained, ‘foc_opt_constrained’ for constrained.
eps_raw (float) – The result of the goal function for the given phi.

pynibs.optimization.multichannel.foc_opt_worker_fun(idx, l, ll, points_wb, connections_wb, l_wb, current_constraint, target_stim)¶: Worker fun to be used with multiprocessing.

pynibs.optimization.multichannel.foc_targets(r, oa, thr_oa, thr_r, result_path)¶

Calculate the prevalence of the target field and transform it into a measure in cm^2.

Parameters:

r (np.ndarray of dimension) – [n_nodes,] Target-to-max Ratio focality measure evaluated at all nodes.
oa (np.ndarray of dimension) – [n_nodes,] Overstimulated Area focality measure evaluated all nodes.
thr_oa (int) – threshold for focality mask: OA <= thr_OA.
thr_r (float) – threshold for focality mask: R >= thr_R.
result_path (string) – path to results folder.

Returns:

valid_array (np.ndarray of dimension) – [n_nodes] Array where valid_indices = 1, invalid_indices = 0, non-evaluated = nan.
valid_indices (np.ndarray) – Indices for which OA <= thr_OA and R >= thr_R is true.

pynibs.optimization.multichannel.get_current_init(n_channel, n_stim, bounds=None, only_pos=None, seed=None)¶: Get initial guess for the optimization problem.

pynibs.optimization.multichannel.get_focality_score(j, ll)¶: Calculate the focality score.

pynibs.optimization.multichannel.get_maximality_currents(e, current_constraint=100)¶

Compute the current set that maximizes |E|

Parameters:

e (np.ndarray of float) – (n_chans, 3) Leadfield in target element.
current_constraint (float, default=100) – Maximum current constraint.

Returns:

j_max (list of float) – [n_chans] Current-set to maximize E.
e_max (Maximum |E| for given current_constraint.)

pynibs.optimization.multichannel.get_score(x, e, n_stim, n_ele, n_channel, x_opt=None, opt_target='elms', fun=None, fun_args=None)¶

Normalize the score by the number of elements.

Parameters:

x (np.ndarray of float) – (n_channel * n_stim_opt, ) Vector to scale channels for each.
e (np.ndarray of float) – (n_ele*3, n_channels) E-field.
n_stim (int) – Number of stimulations to compute score for.
n_ele (int) – Number of elements.
n_channel (int) – Number of channels.
x_opt (np.ndarray of float, optinonal) – (n_pre_opt,) Previously optimized channel currents.

Returns:

score – The normalized score.

Return type:

float

pynibs.optimization.multichannel.get_score_raw(x, e, n_stim, n_ele, n_channel, x_opt=None, opt_target='elms', fun=None, fun_args=None, res_scale=0.001, metric='mean', verbose=False)¶

Compute score for e-efield cross correlations. Non-normalized score is returned, so you need to do sth like

score = 1 / ((n_ele ** 2 - n_ele) / 2) * get_score()

Parameters:

x (np.ndarray of float) – (n_channel * n_stim_opt, ) Vector to scale channels for each.
e (np.ndarray of float) – (n_ele*3, n_channels) E-field.
n_stim (int) – Number of stimulations to compute score for.
n_ele (int) – Number of elements.
n_channel (int) – Number of channels.
x_opt (np.ndarray of float, optinonal) – (n_pre_opt,) Previously optimized channel currents.
opt_target (str, default: 'elms') – Optimization target. ‘elms’ for optimizing decorrelations of elements, ‘stims’ for stimulations.
fun (function, optional) – Function to scale the e-field.
fun_args (tuple, optional) – Arguments for the function.
res_scale (float, default: 1) – Scaling factor for the result.
verbose (book, default: False) – Print out progress.

Returns:

score – non-normalized score np.nansum(np.abs(np.triu(np.corrcoef(e_mag), k=1)))

Return type:

float

pynibs.optimization.multichannel.get_score_raw_single_channel(x, e, x_opt=None)¶

Compute score for e-efield cross correlations. Non-normalized score is returned, so you need to do sth like

score = 1 / ((n_ele ** 2 - n_ele) / 2) * get_score()

Parameters:

x (np.ndarray of float) – (n_placements, ) selection of coil placements.
e (np.ndarray of float) – (n_ele*3, n_channels) E-field.
x_opt (np.ndarray of float, optinonal) – (n_pre_opt,) Previously optimized channel currents.

Returns:

score – non-normalized score np.nansum(np.abs(np.triu(np.corrcoef(e_mag), k=1)))

Return type:

float

pynibs.optimization.multichannel.goal(phi, l, ll, target_stim=80, u=0, s=0)¶

Goal Function for focality optimization.

Parameters:

phi (np.ndarray) – [n_coils] Angles of hyperellipsoid cylinder.
l (np.ndarray of float) – [n_elm*3, n_coils] Flattened lead-field matrix. x,y,z direction in each point as consecutive rows.
ll (np.ndarray of float) – [n_coils, n_coils] Result of matrix multiplication of lead-field matrix * transpose(lead-field matrix)
target_stim (int, default: 80) – Desired target stimulation. scaling value for unit current.
u (np.ndarray of float) – [n_coils, n_coils] Left singular vectors of lead-field matrix.
s (np.ndarray of float) – [n_coils] Singular values of lead-field matrix.

Returns:

eps (int) – Result value of goal function for given phi.
J (np.ndarray of dimension) – [n_coils] Resulting currents for given phi.

pynibs.optimization.multichannel.maximality_focality_optimization(j_foc, flat_e_data, target_stim=80, current_constraint=100, steps=100)¶

Find balanced current solution between maximality and focality that yields a target stimulation while keeping within a current constraint.

Parameters:

j_foc (np.ndarray of float) – Focality current set.
flat_e_data (np.ndarray of float) – (n_chans, 3) Flattened leadfield matrix for the element.
target_stim (float, default=80) – Desired target stimulation.
current_constraint (float, default=100) – Maximum current constraint.
steps (int, default=100) – How many steps to consider for the balancing. Basically the resolution of the optimization.

Returns:

j_inter (np.ndarray of float) – Current set that balances maximality and focality.
dist_foc_int (float) – Distance between focality current and the intersection current.
dist_max_int (float) – Distance between maximality current and the intersection current.

pynibs.optimization.multichannel.optimization_all_k(l, points, connections, radius, center, ll, current_constraint=None, target_stim=80)¶

Run optimization on all target points inside sphere with certain radius. Save results as hdf5.

Parameters:

l (np.ndarray of dimension) – [n_elm*3, n_coils] Flattened lead-field matrix. x,y,z direction in each point as consecutive rows.
points (np.ndarray of dimension) – [n_points, 3] Point elements for describing triangle surface. x,y,z coordinates.
connections (np.ndarray of dimension) – [n_elm x 3] Triangle surface elements. defined by three point element indices.
radius (int) – Sphere radius.
center (list of float) – [3] Center [x, y, z] of target sphere.
ll (np.ndarray) – [n_coils, n_coils] Result of matrix multiplication of lead-field matrix * transpose(lead-field matrix).
current_constraint (int, optional) – If set, applies a current constraint during the optimization, ensuring the resulting currents don’t exceed a certain threshold (100 in this case).
target_stim (int) – desired target stimulation value.

Returns:

ee (np.ndarray of dimension [n_targets, n_elm, 3]) – Electric fields computed for each target in the optimization process. Each row represents the electric field values across the ROI for a specific target.
eps (list of int) – Result values of the goal function for each target. This represents how well the optimization achieved the target stimulation for each point.
jj (np.ndarray of dimension [n_points x n_coils]) – Currents for each point in the sphere. Each row corresponds to the resulting currents for a specific point.
phi (list of np.ndarray) – Optimal angles of the hyperellipsoid for each target point, resulting from the optimization.
curr_cons (list of float) – The maximum stimulation for each target after applying current constraints (if any).
idx (list of int) – Indices of the points inside the target sphere, selected based on the radius and center.

pynibs.optimization.multichannel.optimization_roi(l, points_wb, connections_wb, l_wb, opt_indices=None, current_constraint=None, target_stim=100, hdf5_file_path=None, n_cpus=1)¶

Run optimization on all target points inside sphere with certain radius. Save results as hdf5.

Parameters:

l (np.ndarray of dimension) – [n_elm*3, n_coils] Flattened lead-field matrix. x,y,z direction in each point as consecutive rows.
points_wb (np.ndarray of float) – [n_points, 3] Point elements for describing triangle surface. x,y,z coordinates.
connections_wb (np.ndarray of float) – [n_elm x 3] Triangle surface elements. defined by three point element indices.
l_wb (np.ndarray of float) – [n_elm*3, n_coils] Flattened lead-field matrix. x,y,z direction in each point as consecutive rows.
opt_indices (list of int, optional) – [n_points] Indices of the elements to optimize.
current_constraint (int, optional) – If set, applies a current constraint during the optimization, ensuring the resulting currents don’t exceed a certain threshold (100 in this case).
target_stim (int) – desired target stimulation value
hdf5_file_path (str, optional) – Path to save the results in HDF5 format.
n_cpus (int, default: 1) – Number of CPUs to use for parallel processing.

Returns:

eps (list of int) – Result values of the goal function for each target. This represents how well the optimization achieved the target stimulation for each point.
j (np.ndarray of dimension [n_points x n_coils]) – Currents for each point in the sphere. Each row corresponds to the resulting currents for a specific point.
phi (list of np.ndarray) – Optimal angles of the hyperellipsoid for each target point, resulting from the optimization.
max_e (list of float) – The maximum stimulation for each target after applying current constraints (if any).
constraint_types (list of string) – The type of constraint optimization used. ‘foc_opt’ for unconstrained, ‘foc_opt_constrained’ for constrained.
eps_raw (list of float) – The result of the goal function for the given phi.
eval_res (dict) – Dictionary containing various evaluation metrics for the optimization process, including overstimulated area and target-to-max ratio.

pynibs.optimization.multichannel.optimize_currents(e, n_stim, currents_prev=None, seed=None, maxiter=200, basiniter=50, method='SLSQP', opt_target='elms', fun=None, fun_args=None, bounds=None, x0=None, constraints=None, verbose=False)¶

Optimize the currents for a multichannel TMS array by minimizing e-fields cross-correlation.

Parameters:

e (np.ndarray of float) – (n_elms * 3, n_channel) or (n_elms, 3, n_channel). E in ROI for currents = 1.
n_stim (int) – Number of stimulations.
currents_prev (np.ndarray of float, optional) – (n_channels, n_stims_prev) Previous currents to append to.
seed (int, optional) – Seed for random number generator.
maxiter (int, default=200) – Max iterations of the optimization.
method (str, default: 'SLSQP') – Optimization method.
opt_target (str, default: 'elms') – Optimization target. ‘elms’ for optimizing decorrelations of elements, ‘stims’ for stimulations.
fun (function, optional) – Function to scale the e-field.
fun_args (tuple, optional) – Arguments for the fun function.
bounds (tuple, optional) – Bounds of the optimization problem.
x0 (np.ndarray of float, optional) – Initial guess for the optimization problem.
constraints (list of dict, optional) – Constraints for the optimization problem.
verbose (bool, default: False) – Print additional information.

Returns:

currents (np.ndarray) – (n_channels, n_stims) The optimized currents to drive the multichannel array.
score (float) – Final score of the solution.
res

pynibs.optimization.multichannel.optimize_currents_single_channel(e, n_stim, currents_prev=None, seed=None, maxiter=200, method='SLSQP', verbose=False)¶

Optimize the coil placement selection for a single channel e-field set minimizing e-fields cross-correlation.

Parameters:

e (np.ndarray of float) – (n_elms3, n_placements).
n_stim (int) – Number of stimulations.
currents_prev (np.ndarray of float, optional) – (n_channels, n_stims_prev) Previous currents to append to.
seed (int, optional) – Seed for random number generator.
maxiter (int, default=200) – Max iterations of the optimization.
method (str, default: 'SLSQP') – Optimization method.
verbose (bool, default: False) – Print additional information.

Returns:

currents (np.ndarray) – (n_channels, n_stims) The optimized currents to drive the multichannel array.
score (float) – Final score of the solution.

pynibs.optimization.multichannel.overstimulated_area(e, k, connections, points, percentile=1, e_target=None, e_in_tris=True)¶

Calculate the prevalence of the target field and transform it into a measure in cm^2.

Parameters:

e (np.ndarray) – [n_nodes,] E-Field magnitudes (in triangles or nodes) .
k (int) – Index of the optimized target.
connections (np.ndarray) – [n_connections, 3] Triangle surface elements defined by three point element indices.
points (np.ndarray of dimension) – [n_points, 3] Coordinates of the mesh points.
percentile (int, default = 1) – percentile*e_target >= E, range 0-1.
e_target (float, optional) – E-field at the target position. If not provided, it is calculated from e.
e_in_tris (bool, default: True) – Is E provided in tris (True) or nodes (False)

Returns:

eva_rel (float) – Relative area for which |E|_target >= |E|_offtarget.
eva_abs (int) – Absolute count of the indices where |E|_target >= |E|_offtarget.
total_area_cm2 (float) – Total area in cm² for which |E|_target >= |E|_offtarget.

pynibs.optimization.multichannel.scd(connections, points, idx, mesh_skin, center_skin=None)¶

Calculate skin-cortex distance between midlayer target points and skin surface

Parameters:

connections (np.ndarray of dimension [n_ROI/3 x 3]) – Triangle surface elements, defined by three point element indices for each triangle.
points (np.ndarray of dimension [n_points x 3]) – Point elements describing the triangle surface. Each point has x, y, z coordinates.
idx (list or int) – Indices of the target points or a single index. If scalar, it is converted to a list of one element.
mesh_skin (simnibs.mesh_io.Mesh) – A SimNIBS skin surface mesh.
center_skin (np.ndarray of dimension, default: None) – (n_elements, 3) The centers of the skin mesh elements. If not provided, it is calculated using the elements_baricenters() method of the mesh_skin.

Returns:

target_idx (np.ndarray of dimension [n_targetsphere]) – indices of targets in sphere
scd (np.ndarray of dimension [n_targetsphere]) – skin cortex distance, or rather skin gm midlayer distance for each triangle mesh in target sphere

pynibs.optimization.multichannel.sphere(points, connections, radius, center)¶

Define target-sphere around M1-center with given radius. All triangle surface elements included that contain minimum one point inside sphere.

Parameters:

points (np.ndarray) – [n_points x 3] Point elements for describing triangle surface. x,y,z coordinates.
connections (np.ndarray) – [n_elm x 3] Triangle connectivity, indices of points.
radius (int) – Sphere radius.
center (list of int) – [3] Center [x, y, z] of sphere.

Returns:

valid_p (np.ndarray) – All points included in the sphere.
idx_c (np.ndarray) – Triangle surface indices included in the sphere.

pynibs.optimization.multichannel.strain_constr(x, max_strain, n_channel, n_stim)¶: Compute mechanical strain on 5 channel planar mTMS coil

pynibs.optimization.multichannel.strain_constr_min(x, max_strain, n_channel, n_stim)¶: Compute mechanical strain on 5 channel planar mTMS coil

pynibs.optimization.multichannel.target_to_max_ratio(e, k, e_target=None)¶

Evaluation-metric: E_target/E_max. Quantifies the difference between target |E| and max(|E|).

Parameters:

e (np.ndarray) – [n_elm, 1] or [n_elm, 3] E-Field.
k (int) – Index of optimized target.
e_target (float, optional) – E-field at the target position. If not provided, it is calculated from e.

Returns:

t2m – e_target/e_max Metric if max_stim = 1: e_target == e_max if max_stim < 1: e_target << e_max -> values < 1 indicate worse focality

Return type:

int

pynibs.optimization.multichannel.write_magventure_xls(opt_currents_fn, max_mso, xls_fn=None, rep=1)¶: Writes .xls file for MagVenture mTMS systems.

pynibs.optimization.opt_mep module¶

Functions to select optimal stimulation conditions based on the estimated y (i.e., MEPs in MEP studies).

pynibs.optimization.opt_mep.get_det_fim(x, fun, p, fim_matrix)¶

Updates the Fisher Information Matrix and returns the negative determinant based on the sample x. It is a score how much information the additional sample yields.

Parameters:

fun (callable) – Callable defined in interval [0, 1].
x (float) – Single sample location (interval [0, 1]).
p (dict) – Dictionary containing the parameter estimates. The keys are the parameter names of fun.
fim_matrix (np.ndarray of float) – (n_params, n_params) Fisher Information Matrix.

Returns:

det – Determinant of the Fisher Information Matrix after adding sample x.

Return type:

float

pynibs.optimization.opt_mep.get_fim_sample(fun, x, p)¶

Get Fisher Information Matrix of one single sample.

Parameters:

fun (callable) – Callable the fisher information matrix is calculated for. The sample is passed as the first argument.
x (float) – Sample passed to function.
p (dict) – Dictionary containing the parameter estimates. The keys are the parameter names of fun.

Returns:

fim_matrix – (n_params, n_params) Fisher information matrix.

Return type:

np.ndarray of float

pynibs.optimization.opt_mep.get_optimal_sample_fim(fun, p, x=None)¶

Determines optimal location of next sample by maximizing the determinant of the Fisher Information Matrix.

Parameters:

fun (callable) – Callable (interval [0, 1]).
x (np.ndarray of float, optional) – Previous sample locations (interval [0, 1]).
p (dict) – Dictionary containing the parameter estimates. The keys are the parameter names of fun.

Returns:

x_opt – Optimal location of next sample (interval [0, 1]).

Return type:

float

pynibs.optimization.opt_mep.init_fim_matrix(fun, x, p)¶

Initializes the Fisher Information Matrix based on the samples given in x.

Parameters:

fun (callable) – Callable defined in interval [0, 1].
x (np.ndarray of float) – Initial sample locations (interval [0, 1]).
p (dict) – Dictionary containing the parameter estimates. The keys are the parameter names of fun.

Returns:

fim_matrix – Fisher Information Matrix.

Return type:

np.ndarray of float [n_params x n_params]

pynibs.optimization.workhorses module¶

pynibs.optimization.workhorses.coil_wise_corr(idx_list, array, ele_idx_1, **kwargs)¶

pynibs.optimization.workhorses.coverage(idx_list, array, x, y, ele_idx_1, **kwargs)¶

Determine coverage score (likelihood) for given zap indices in idx_list

Parameters:

idx_list (list of lists [n_combs][n_zaps]) – Index lists of zaps containing different possible combinations. Usually only the last index changes.
array (np.ndarray of float [n_zaps x n_ele]) – Electric field for different coil positions and elements.
x (np.ndarray of float [200 x n_ele]) – x-values of coverage distributions, defined in interval [0, 1] (element wise normalized electric field).
y (np.ndarray of float [200 x n_ele]) – y-values of coverage distributions (element wise probability of already included e-fields).
ele_idx_1 (np.ndarray of float [n_roi]) – Element indices for which the coverage optimization is performed for.

Returns:

res – Coverage score (likelihood) for given electric field combinations. Lower values indicate that the new zap fills a gap which was not covered before.

Return type:

np.ndarray of float [n_combs]

pynibs.optimization.workhorses.coverage_prepare(idx_list, array, zap_idx, **kwargs)¶

Prepares coverage calculation. Determines coverage distributions for elements in idx_list given the zaps in zap_idx

Parameters:

idx_list (list [n_ele]) – Index lists of elements.
array (ndarray of float [n_zaps x n_ele]) – Electric field for different coil positions and elements.
zap_idx (ndarray of int) – Included zaps in coverage distribution.

Returns:

x (ndarray of float [200 x n_ele]) – x-values of coverage distributions, defined in interval [0, 1] (element wise normalized electric field).
y (ndarray of float [200 x n_ele]) – y-values of coverage distributions (element wise probability of already included e-fields).

pynibs.optimization.workhorses.dist(idx_list, array, ele_idx_1, **kwargs)¶

Determines distance score for given zap indices in idx_list.

Parameters:

idx_list (list of lists [n_combs][n_zaps]) – Index lists of zaps containing different possible combinations. Usually only the last index changes.
array (np.ndarray of float [n_zaps x n_ele]) – Electric field for different coil positions and elements.
ele_idx_1 (np.ndarray of float [n_ele]) – Element indices for which the optimization is performed.

Returns:

res – Distance based score. Lower values indicate more equidistant sampling (better)

Return type:

np.ndarray of float [n_combs]

pynibs.optimization.workhorses.dist_mc(idx_list, array, ele_idx_1, ele_idx_2, mode='cols', **kwargs)¶

Determines distance score and mutual coherence for given zap indices in idx_list. If c_max_idx is given, the distance based score is calculated only for this element. The condition number however is optimized for all elements in array.

Parameters:

idx_list (list of lists [n_combs][n_zaps]) – Index lists of zaps containing different possible combinations. Usually only the last index changes.
array (np.ndarray of float [n_zaps x n_ele]) – Electric field for different coil positions and elements.
mode (str, default: "cols") – Set if the mutual coherence is calculated w.r.t. columns or rows (“cols”, “rows”).
ele_idx_1 (np.ndarray of float [n_ele]) – Element indices for which the dist optimization is performed for.
ele_idx_2 (np.ndarray of float [n_ele]) – Element indices for which the mc optimization is performed for.

Returns:

res_dist (np.ndarray of float [n_combs]) – Distance based score. Lower values indicate more equidistant sampling (better).
res_mc (np.ndarray of float [n_combs]) – Mutual coherence. Lower values indicate more orthogonal e-field combinations (better).

pynibs.optimization.workhorses.dist_svd(idx_list, array, ele_idx_1, ele_idx_2, **kwargs)¶

Determines distance score and condition number for given zap indices in idx_list. If c_max_idx is given, the distance based score is calculated only for this element. The condition number however is optimized for all elements in array

Parameters:

idx_list (list of lists [n_combs][n_zaps]) – Index lists of zaps containing different possible combinations. Usually only the last index changes.
array (np.ndarray of float [n_zaps x n_ele]) – Electric field for different coil positions and elements.
ele_idx_1 (np.ndarray of float [n_ele]) – Element indices for which the dist optimization is performed for.
ele_idx_2 (np.ndarray of float [n_ele]) – Element indices for which the svd optimization is performed for.

Returns:

res_dist (np.ndarray of float [n_combs]) – Distance based score. Lower values indicate more equidistant sampling (better).
res_svd (np.ndarray of float [n_combs]) – Condition number. Lower values indicate more orthogonal e-field combinations (better).

pynibs.optimization.workhorses.fim(idx_list, array, ele_idx_1, e_opt, c=None, **kwargs)¶

Determine difference between e-fields and optimal e-field determined using the Fisher Information Matrix.

Parameters:

idx_list (list of lists [n_combs][n_zaps]) – Index lists of zaps containing different possible combinations. Usually only the last index changes.
array (np.ndarray of float [n_zaps x n_ele]) – Electric field for different coil positions and elements.
ele_idx_1 (np.ndarray of float [n_roi]) – Element indices for which the fim optimization is performed for.
e_opt (np.ndarray of float [n_roi]) – Optimal electric field value(s) (target) determined by FIM method.
c (np.ndarray of float [n_ele], optional) – Congruence factor map normalized to 1 (whole ROI) used to weight the difference between the optimal e-field and the candidate e-field. If None, no weighting is applied.

Returns:

res – Difference between e-fields and optimal e-field.

Return type:

np.ndarray of float [n_combs]

pynibs.optimization.workhorses.fim_mc(idx_list, array, ele_idx_1, ele_idx_2, e_opt, c=None, mode='rows', **kwargs)¶

Determine difference between e-fields and optimal e-field determined using the Fisher Information Matrix and mutual coherence.

Parameters:

idx_list (list of lists [n_combs][n_zaps]) – Index lists of zaps containing different possible combinations. Usually only the last index changes.
array (np.ndarray of float [n_zaps x n_ele]) – Electric field for different coil positions and elements.
ele_idx_1 (np.ndarray of float [n_roi_1]) – Element indices for which the fim optimization is performed for.
ele_idx_2 (np.ndarray of float [n_roi_2]) – Element indices for which the mc optimization is performed for.
e_opt (float) – Optimal electric field value (target) determined by FIM method.
c (np.ndarray of float [n_ele], optional) – Congruence factor map normalized to 1 (whole ROI) used to weight the difference between the optimal e-field and the candidate e-field. If None, no weighting is applied.
mode (str) –

Returns:

res_fim (np.ndarray of float [n_combs]) – Difference between e-fields and optimal e-field.
res_mc (np.ndarray of float [n_combs]) – Mutual coherence. Lower values indicate more orthogonal e-field combinations (better)

pynibs.optimization.workhorses.fim_svd(idx_list, array, ele_idx_1, ele_idx_2, e_opt, c=None, **kwargs)¶

Determine difference between e-fields and optimal e-field determined using the Fisher Information Matrix and condition number.

Parameters:

idx_list (list of lists [n_combs][n_zaps]) – Index lists of zaps containing different possible combinations. Usually only the last index changes.
array (np.ndarray of float [n_zaps x n_ele]) – Electric field for different coil positions and elements.
ele_idx_1 (np.ndarray of float [n_roi_1]) – Element indices for which the fim optimization is performed for.
ele_idx_2 (np.ndarray of float [n_roi_2]) – Element indices for which the svd optimization is performed for.
e_opt (float) – Optimal electric field value (target) determined by FIM method.
c (np.ndarray of float [n_ele], optional) – Congruence factor map normalized to 1 (whole ROI) used to weight the difference between the optimal e-field and the candidate e-field. If None, no weighting is applied.

Returns:

res_fim (np.ndarray of float [n_combs]) – Difference between e-fields and optimal e-field.
res_svd (np.ndarray of float [n_combs]) – Condition number. Lower values indicate more orthogonal e-field combinations (better)

pynibs.optimization.workhorses.mc(idx_list, array, ele_idx_1, mode='cols', **kwargs)¶

Determines mutual coherence for given zap indices in idx_list.

Parameters:

idx_list (list of lists [n_combs][n_zaps]) – Index lists of zaps containing different possible combinations. Usually only the last index changes.
array (np.ndarray of float [n_zaps x n_ele]) – Electric field for different coil positions and elements.
ele_idx_1 (np.ndarray of float [n_ele]) – Element indices for which the optimization is performed.
mode (str, default: "cols") – Set if the mutual coherence is calculated w.r.t. columns or rows (“cols”, “rows”).

Returns:

res – Mutual coherence. Lower values indicate more orthogonal e-field combinations (better).

Return type:

np.ndarray of float [n_combs]

pynibs.optimization.workhorses.roi_elmt_wise_corr(idx_list, e_mat, ele_idx_1, decorrelate_hotspot_only=False, backend=<module 'numpy' from '/home/docs/checkouts/readthedocs.org/user_builds/pynibs/envs/master/lib/python3.11/site-packages/numpy/__init__.py'>, **kwargs)¶

Compute element wise correlation for sets of e-field.

Parameters:

idx_list (list of int) – List of efield indicies in e_mat.
e_mat (np.ndarray) – All e-fields.
ele_idx_1 (list of int) – All element indices to compute corcoeff for.
decorrelate_hotspot_only (bool, default: False) – If true, ele_idx_1[-.1] is used for decorrelation.
backend (module) – Package to use to compute correlation - probably either numpy or cuda.

Returns:

res – (e_mat.shape[0], ) with correlations for all e-field sets in idx_list.

Return type:

np.ndarray

pynibs.optimization.workhorses.rowvec_diff(candidate_coil_idcs, selected_coil_idcs, efields_diff_mat)¶

Parameters:

candidate_coil_idcs (np.ndarry[int]) – List of indices of coil configurations that are still available to pick for the optiized sequence.
selected_coil_idcs (np.ndarray[int]) – List of indices of coil configurations that have already been selected for the optimized sequence.
efields_diff_mat (np.ndarray[float], [n_coil,n_coil]) – Difference matrix, where each cell denotes the magnitude of the difference vector between two coil configurations (determined by row_idx,col_idx).

Returns:

coil_idx – Index of coil configuration with maximal minimal difference to the set of already selected coil configurations.

Return type:

int

pynibs.optimization.workhorses.rowvec_diff_prepare(idx_list, array, ele_idx_1, **kwargs)¶

Computes the part of the difference matrix of row-vectors specified by the row indices in idx_list. Assumption: ‘idx_list’ must be sorted and valid within ‘array’.

Parameters:

idx_list – typing.List[int] List of row indices whose difference should be determined.
array – numpy.typing.ArrayLike [n_coil x n_ele] E-field matrix of all possible coil positions.
ele_idx_1 – numpy.typing.ArrayLike [n_ele] Indices of the ROI elements that should be considered for optimization.

Returns:

numpy.typing.ArrayLike [array.shape[0] x array.shape[0]] = [n_rows x n_rows] The difference matrix with the lenght of the difference vectors between pairs of row vectors specified by idx_list. All other (not calculated) paris of row vectors have a score of 0 in this matrix.

pynibs.optimization.workhorses.smooth(idx_list, array, ele_idx_1, **kwargs)¶

pynibs.optimization.workhorses.svd(idx_list, array, ele_idx_1, **kwargs)¶

Determines condition number for given zap indices in idx_list.

Parameters:

idx_list (list of lists [n_combs][n_zaps]) – Index lists of zaps containing different possible combinations. Usually only the last index changes.
array (np.ndarray of float [n_zaps x n_ele]) – Electric field for different coil positions and elements.
ele_idx_1 (np.ndarray of float [n_ele]) – Element indices for which the optimization is performed.

Returns:

res – Condition number. Lower values indicate more orthogonal e-field combinations (better).

Return type:

np.ndarray of float [n_combs]

pynibs.optimization.workhorses.var(idx_list, array, ele_idx_1, **kwargs)¶

pynibs.optimization.workhorses.variability(idx_list, array, ele_idx_1, **kwargs)¶

Determines variability score for given zap indices in idx_list.

Parameters:

idx_list (list of lists [n_combs][n_zaps]) – Index lists of zaps containing different possible combinations. Usually only the last index changes.
array (np.ndarray of float [n_zaps x n_ele]) – Electric field for different coil positions and elements.
ele_idx_1 (np.ndarray of float [n_ele]) – Element indices for which the optimization is performed.

Returns:

res – Condition number. Lower values indicate more orthogonal e-field combinations (better).

Return type:

np.ndarray of float [n_combs]

pynibs.optimization package

pynibs.optimization package¶

Submodules¶

pynibs.optimization.coil_opt module¶

pynibs.optimization.multichannel module¶

pynibs.optimization.opt_mep module¶

pynibs.optimization.workhorses module¶

Table of Contents

Previous topic

Next topic

This Page