Regularize basis set using idealized SNR measure.
(int_order, ext_order, S_decomp, mag_or_fine, extended_remove)
| 2353 | |
| 2354 | |
| 2355 | def _regularize_in(int_order, ext_order, S_decomp, mag_or_fine, extended_remove): |
| 2356 | """Regularize basis set using idealized SNR measure.""" |
| 2357 | n_in, n_out = _get_n_moments([int_order, ext_order]) |
| 2358 | |
| 2359 | # The "signal" terms depend only on the inner expansion order |
| 2360 | # (i.e., not sensor geometry or head position / expansion origin) |
| 2361 | a_lm_sq, rho_i = _compute_sphere_activation_in(np.arange(int_order + 1)) |
| 2362 | degrees, orders = _get_degrees_orders(int_order) |
| 2363 | a_lm_sq = a_lm_sq[degrees] |
| 2364 | |
| 2365 | I_tots = np.zeros(n_in) # we might not traverse all, so use np.zeros |
| 2366 | in_keepers = list(range(n_in)) |
| 2367 | out_removes = _regularize_out(int_order, ext_order, mag_or_fine, extended_remove) |
| 2368 | out_keepers = list(np.setdiff1d(np.arange(n_in, S_decomp.shape[1]), out_removes)) |
| 2369 | remove_order = [] |
| 2370 | S_decomp = S_decomp.copy() |
| 2371 | use_norm = np.sqrt(np.sum(S_decomp * S_decomp, axis=0)) |
| 2372 | S_decomp /= use_norm |
| 2373 | eigs = np.zeros((n_in, 2)) |
| 2374 | |
| 2375 | # plot = False # for debugging |
| 2376 | # if plot: |
| 2377 | # import matplotlib.pyplot as plt |
| 2378 | # fig, axs = plt.subplots(3, figsize=[6, 12]) |
| 2379 | # plot_ord = np.empty(n_in, int) |
| 2380 | # plot_ord.fill(-1) |
| 2381 | # count = 0 |
| 2382 | # # Reorder plot to match MF |
| 2383 | # for degree in range(1, int_order + 1): |
| 2384 | # for order in range(0, degree + 1): |
| 2385 | # assert plot_ord[count] == -1 |
| 2386 | # plot_ord[count] = _deg_ord_idx(degree, order) |
| 2387 | # count += 1 |
| 2388 | # if order > 0: |
| 2389 | # assert plot_ord[count] == -1 |
| 2390 | # plot_ord[count] = _deg_ord_idx(degree, -order) |
| 2391 | # count += 1 |
| 2392 | # assert count == n_in |
| 2393 | # assert (plot_ord >= 0).all() |
| 2394 | # assert len(np.unique(plot_ord)) == n_in |
| 2395 | noise_lev = 5e-13 # noise level in T/m |
| 2396 | noise_lev *= noise_lev # effectively what would happen by earlier multiply |
| 2397 | for ii in range(n_in): |
| 2398 | this_S = S_decomp.take(in_keepers + out_keepers, axis=1) |
| 2399 | u, s, v = _safe_svd(this_S, full_matrices=False, **check_disable) |
| 2400 | del this_S |
| 2401 | eigs[ii] = s[[0, -1]] |
| 2402 | v = v.T[: len(in_keepers)] |
| 2403 | v /= use_norm[in_keepers][:, np.newaxis] |
| 2404 | eta_lm_sq = np.dot(v * 1.0 / s, u.T) |
| 2405 | del u, s, v |
| 2406 | eta_lm_sq *= eta_lm_sq |
| 2407 | eta_lm_sq = eta_lm_sq.sum(axis=1) |
| 2408 | eta_lm_sq *= noise_lev |
| 2409 | |
| 2410 | # Mysterious scale factors to match MF, likely due to differences |
| 2411 | # in the basis normalizations... |
| 2412 | eta_lm_sq[orders[in_keepers] == 0] *= 2 |
no test coverage detected