Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Modified Levenberg Marquardt LOCO #62

Open
wants to merge 3 commits into
base: master
Choose a base branch
from
Open

Modified Levenberg Marquardt LOCO #62

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
b74bc03
Add Levenberg-Marquardt regularization to the SVD LOCO function.
Nov 5, 2024
672996c
Correct implementation.
Nov 5, 2024
1fb2e65
Add documentation
Nov 7, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
29 changes: 24 additions & 5 deletions pySC/correction/loco.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -89,15 +89,15 @@ def loco_correction_lm(initial_guess0, orm_model, orm_measured, Jn, lengths, inc
return result.x


def loco_correction_ng(initial_guess0, orm_model, orm_measured, J, lengths, including_fit_parameters, s_cut, weights=1):
def loco_correction_ng(initial_guess0, orm_model, orm_measured, J, lengths, including_fit_parameters, s_cut, weights=1, LM_lambda=0):
initial_guess = initial_guess0.copy()
mask = _get_parameters_mask(including_fit_parameters, lengths)
residuals = objective(initial_guess[mask], orm_measured - orm_model, J[mask, :, :], weights)
r = residuals.reshape(np.transpose(orm_model).shape)
t2 = np.zeros([len(initial_guess[mask]), 1])
for i in range(len(initial_guess[mask])):
t2[i] = np.sum(np.dot(np.dot(J[i].T, weights), r.T))
return get_inverse(J[mask, :, :], t2, s_cut, weights)
return get_inverse(J[mask, :, :], t2, s_cut, weights, LM_lambda=LM_lambda)


def objective(masked_params, orm_residuals, masked_jacobian, weights):
Expand Down Expand Up @@ -157,10 +157,29 @@ def select_equally_spaced_elements(total_elements, num_elements):
return total_elements[::step]


def get_inverse(jacobian, B, s_cut, weights, plot=False):
def get_inverse(jacobian, B, s_cut, weights, LM_lambda=0, plot=False):
"""
Calculates $(J^t J + \lambda diag(J^t J))^{-1} J^t r_0$ using a modified
version of the Levenberg-Marquardt method. More information at p.114 of
https://inis.iaea.org/collection/NCLCollectionStore/_Public/54/003/54003559.pdf

Args:
jacobian: jacobian matrix obtained from calculate_jacobian()
bpm_ords: array of element indices of registered BPMs for which to calculate readings
(for convenience, otherwise `SC.ORD.BPM` is used)
s_cut: number of kept singular values when inverting the matrix
weights: weights applied to the elements of the jacobian
LM_lambda: Levenberg-Marquardt lambda parameter. If LM_lambda=0, the
function is equivalent to the Gauss-Newton method
plot: boolean flag to plot the resulting inverse matrix

Returns:
Array of correction values of the same size as B.
"""
n_resp_mats = len(jacobian)
sum_corr = np.sum(jacobian, axis=2) # Sum over i and j for all planes
matrix = np.dot(np.dot(sum_corr, weights), sum_corr.T)
Jt = np.sum(jacobian, axis=2) # Sum over i and j for all planes
Jt_dot_J = np.dot(np.dot(Jt, weights), Jt.T)
matrix = Jt_dot_J + LM_lambda * np.diag(np.diag(Jt_dot_J))
inv_matrix = sc_tools.pinv(matrix, num_removed=n_resp_mats - min(n_resp_mats, s_cut), plot=plot)
results = np.ravel(np.dot(inv_matrix, B))
# e = np.ravel(np.dot(matrix, results)) - np.ravel(B)
Expand Down
Loading