diff --git a/lcls_tools/common/data/emittance.py b/lcls_tools/common/data/emittance.py index d27996cf..ae675982 100644 --- a/lcls_tools/common/data/emittance.py +++ b/lcls_tools/common/data/emittance.py @@ -11,6 +11,7 @@ def compute_emit_bmag( beamsize_squared: np.ndarray, rmat: np.ndarray, twiss_design: np.ndarray = None, + weighted_fit: bool = False, maxiter: int = None, ): """ @@ -40,6 +41,10 @@ def compute_emit_bmag( to each measurement in the respective batch for the calculation of Bmag (useful for quad scans). + weighted_fit : bool, default = False + Boolean specifying whether or not to weight the fitting loss function by the inverse beam + size. If True, smaller beam sizes will be prioritized in the non-linear curve fitting. + maxiter : int, optional Maximum number of iterations to perform in nonlinear fitting (minimization algorithm). @@ -98,8 +103,14 @@ def loss_torch(params): params = torch.reshape(params, [*beamsize_squared.shape[:-2], 3]) sig = torch.stack(beam_matrix_tuple(params), dim=-1).unsqueeze(-1) # sig should now be shape batchshape x 3 x 1 (column vectors) + if weighted_fit: + weights = 1.0 / beamsize_squared.sqrt() + else: + weights = 1.0 total_abs_error = ( - (torch.sqrt(amat @ sig) - torch.sqrt(beamsize_squared)).abs().nansum() + (weights * (torch.sqrt(amat @ sig) - torch.sqrt(beamsize_squared))) + .abs() + .nansum() ) return total_abs_error @@ -119,8 +130,12 @@ def loss(params): params = np.reshape(params, [*beamsize_squared.shape[:-2], 3]) sig = np.expand_dims(np.stack(beam_matrix_tuple(params), axis=-1), axis=-1) # sig should now be shape batchshape x 3 x 1 (column vectors) + if weighted_fit: + weights = 1.0 / np.sqrt(beamsize_squared) + else: + weights = 1.0 total_abs_error = np.nansum( - np.abs(np.sqrt(amat @ sig) - np.sqrt(beamsize_squared)) + np.abs(weights * (np.sqrt(amat @ sig) - np.sqrt(beamsize_squared))) ) return total_abs_error diff --git a/tests/unit_tests/lcls_tools/common/data/test_emittance.py b/tests/unit_tests/lcls_tools/common/data/test_emittance.py index 8e49db50..afd7c60a 100644 --- a/tests/unit_tests/lcls_tools/common/data/test_emittance.py +++ b/tests/unit_tests/lcls_tools/common/data/test_emittance.py @@ -54,6 +54,109 @@ def test_emittance_calc(self): result["bmag"], np.array([[1.0, 1.0, 1.0], [1.0, 1.0, 1.0]]), rtol=1e-2 ) + def test_emittance_calc_with_weighted_fit(self): + # set up test values + beamsize_squared = np.array( + [ + [[0.00299705], [0.12034662], [0.03792509]], + [[0.00259244], [0.23703413], [0.09598636]], + ] + ) + + rmat = np.array( + [ + [ + [[-0.24254679, 1.04600053], [-1.17140665, 0.92885986]], + [[-1.16646303, 0.71788384], [-0.8039517, -0.36251133]], + [[-0.55801204, -0.55330746], [0.61964408, -1.17765613]], + ], + [ + [[-0.24254679, 1.04600053], [-1.17140665, 0.92885986]], + [[-1.16646303, 0.71788384], [-0.8039517, -0.36251133]], + [[-0.55801204, -0.55330746], [0.61964408, -1.17765613]], + ], + ] + ) + + twiss_design = np.array( + [ + [ + [0.29970473, -1.58494282], + [12.03466156, -9.17146332], + [3.79250871, 3.01307483], + ], + [ + [0.1296221, -0.66578778], + [11.85170667, -9.88871364], + [4.79931781, 2.87869131], + ], + ] + ) + + # compute emittance & bmag + result = compute_emit_bmag( + beamsize_squared=beamsize_squared, + rmat=rmat, + twiss_design=twiss_design, + weighted_fit=True, + ) + + # compare results with ground-truth + assert np.allclose(result["emittance"], np.array([[0.01], [0.02]]), rtol=1e-2) + assert np.allclose( + result["bmag"], np.array([[1.0, 1.0, 1.0], [1.0, 1.0, 1.0]]), rtol=1e-2 + ) + + def test_emittance_calc_with_nans(self): + # set up test values + beamsize_squared = np.array( + [ + [[0.00299705], [0.12034662], [0.03792509]], + [[0.00259244], [0.23703413], [0.09598636]], + ] + ) + + rmat = np.array( + [ + [ + [[-0.24254679, 1.04600053], [-1.17140665, 0.92885986]], + [[-1.16646303, 0.71788384], [-0.8039517, -0.36251133]], + [[-0.55801204, -0.55330746], [0.61964408, -1.17765613]], + ], + [ + [[-0.24254679, 1.04600053], [-1.17140665, 0.92885986]], + [[-1.16646303, 0.71788384], [-0.8039517, -0.36251133]], + [[-0.55801204, -0.55330746], [0.61964408, -1.17765613]], + ], + ] + ) + + twiss_design = np.array( + [ + [ + [0.29970473, -1.58494282], + [12.03466156, -9.17146332], + [3.79250871, 3.01307483], + ], + [ + [0.1296221, -0.66578778], + [11.85170667, -9.88871364], + [4.79931781, 2.87869131], + ], + ] + ) + + # compute emittance & bmag + result = compute_emit_bmag( + beamsize_squared=beamsize_squared, rmat=rmat, twiss_design=twiss_design + ) + + # compare results with ground-truth + assert np.allclose(result["emittance"], np.array([[0.01], [0.02]]), rtol=1e-2) + assert np.allclose( + result["bmag"], np.array([[1.0, 1.0, 1.0], [1.0, 1.0, 1.0]]), rtol=1e-2 + ) + def test_emittance_calc_with_broadcast_twiss(self): # set up test values beamsize_squared = np.array(