diff --git a/espm/datasets/eds_spim.py b/espm/datasets/eds_spim.py index 2310f26..907a465 100644 --- a/espm/datasets/eds_spim.py +++ b/espm/datasets/eds_spim.py @@ -10,7 +10,6 @@ - set fixed W for the :class:`espm.estimators.NMFEstimator` decomposition """ -import json import warnings import hyperspy.api as hs @@ -31,11 +30,12 @@ from scipy.optimize import curve_fit from tqdm import tqdm -from espm.conf import NUMBER_PERIODIC_TABLE, SYMBOLS_PERIODIC_TABLE from espm.estimators import NMFEstimator, SmoothNMF from espm.models import EDXS from espm.utils import ( get_explained_intensity_W, + get_npt, + get_spt, num_to_symbol, number_to_symbol_list, quant_spectrum, @@ -43,9 +43,6 @@ symbol_to_number_list, ) -# TODO: use cache after merging optimisation -NPT = json.load(open(NUMBER_PERIODIC_TABLE)) - def _check_decomposition(func) -> bool: # @wraps(func) @@ -407,15 +404,12 @@ def carto_fixed_W(self, brstlg_comps: int = 1) -> np.ndarray: ) elements = self.metadata.EDS_model.elements if self.problem_type == "no_brstlg": - W = np.diag(-1 * np.ones((len(elements),))) + W = np.diag(np.full((len(elements),), -1.0)) elif self.problem_type == "bremsstrahlung": - W1 = np.diag(-1 * np.ones((len(elements),))) - W2 = np.zeros((2, len(elements))) - W_elts = np.vstack((W1, W2)) - W3 = np.zeros((len(elements), brstlg_comps)) - W4 = -1 * np.ones((2, brstlg_comps)) - W_brstlg = np.vstack((W3, W4)) - W = np.hstack((W_elts, W_brstlg)) + L = len(elements) + W = np.zeros((L + 2, L + brstlg_comps)) + np.fill_diagonal(W[:L, :L], -1.0) + W[L:, L:] = -1.0 return W @@ -449,31 +443,26 @@ def convert_to_symbols(elements=[]): conv_elts = convert_to_symbols(elements=elements) if self.problem_type == "no_brstlg": - W = -1 * np.ones((len(raw_elts), len(phases_dict.keys()))) + W = np.full((len(raw_elts), len(phases_dict.keys())), -1.0) elif self.problem_type == "bremsstrahlung": - W = -1 * np.ones((len(raw_elts) + 2, len(phases_dict.keys()))) + W = np.full((len(raw_elts) + 2, len(phases_dict.keys())), -1.0) else: raise ValueError( "problem type should be either no_brstlg or bremsstrahlung" ) for p, phase in enumerate(phases_dict): for key in phases_dict[phase]: - if key == "b0": - if self.problem_type == "bremsstrahlung": - W[-2, p] = phases_dict[phase][key] - else: - warnings.warn( - "The chosen EDXS modelling does not incorporate the bremsstrahlung. Input bremsstrahlung parameters will be ignored." - ) - if key == "b1": - if self.problem_type == "bremsstrahlung": - W[-1, p] = phases_dict[phase][key] - else: - warnings.warn( - "The chosen EDXS modelling does not incorporate the bremsstrahlung. Input bremsstrahlung parameters will be ignored." - ) - if key in conv_elts: - W[indices[conv_elts.index(key)], p] = phases_dict[phase][key] + match key: + case "b0" | "b1": + if self.problem_type == "bremsstrahlung": + idx = -2 if key == "b0" else -1 + W[idx, p] = phases_dict[phase][key] + else: + warnings.warn( + "The chosen EDXS modelling does not incorporate the bremsstrahlung. Input bremsstrahlung parameters will be ignored." + ) + case _ if key in conv_elts: + W[indices[conv_elts.index(key)], p] = phases_dict[phase][key] return W def print_concentration_report( @@ -591,11 +580,11 @@ def estimate_mass_thickness( W_init = estimator.W_ elts = list(self.model.get_elements(include_ignored=False)) elts_indices = self.model.NMF_simplex() - new_elts_dict = {elts[i]: W_init[elts_indices[i]] for i in range(len(elts))} + new_elts_dict = dict(zip(elts, W_init[elts_indices])) - _ = 0 curr_mt = self.metadata.Sample.thickness * self.metadata.Sample.density - while _ < 5: + + for _ in range(5): # first init of the model brstlg_model, mask = self.model.bremsstrahlung_only_tools( mass_thickness=curr_mt, elements_dict=new_elts_dict, ranges=self.ranges @@ -607,12 +596,12 @@ def estimate_mass_thickness( brstlg_estimator.fit(masked_X[:, np.newaxis]) W_brstlg = np.vstack( ( - -1 - * np.ones( + np.full( ( W_init.shape[0] - brstlg_estimator.W_.shape[0], brstlg_estimator.W_.shape[1], - ) + ), + -1.0, ), brstlg_estimator.W_, ) @@ -629,14 +618,10 @@ def estimate_mass_thickness( W_init = estimator.W_ H_init = estimator.H_ - elts = list(self.model.get_elements(include_ignored=False)) - elts_indices = self.model.NMF_simplex() - new_elts_dict = {elts[i]: W_init[elts_indices[i]] for i in range(len(elts))} + new_elts_dict = dict(zip(elts, W_init[elts_indices])) total_weight = self._elements_dict_to_weights(new_elts_dict) curr_mt = self._extract_mass_thickness(H_init.sum(), total_weight) - _ += 1 - print( f"The current estimated mass-thickness is {curr_mt} g.cm^-2", flush=True, @@ -677,7 +662,7 @@ def _elements_dict_to_weights(self, elements_dict): Total weight of the elements in grams. """ total_weight = sum( - quantity * NPT["table"][element]["atomic_mass"] * 1.66053906660e-24 + quantity * get_npt()[element]["atomic_mass"] * 1.66053906660e-24 for element, quantity in elements_dict.items() ) return total_weight @@ -736,7 +721,7 @@ def select_background_windows(self, num_windows=4, ranges=None): cm = self._register_ranges init_ranges = self._generate_ranges(num_windows) self.spans = [] - for i in range(num_windows): + for _ in range(num_windows): self.spans.append(Signal1DRangeSelector(self)) for j, span in enumerate(self.spans): @@ -766,7 +751,7 @@ def _register_ranges(self, signal, left, right): def _compute_bremsstrahlung(self): mt = self.metadata.Sample.density * self.metadata.Sample.thickness - elts_dict = {elt: 1.0 for elt in self.metadata.Sample.elements} + elts_dict = dict.fromkeys(self.metadata.Sample.elements, 1.0) brstlg_model, mask = self.model.bremsstrahlung_only_tools( mass_thickness=mt, elements_dict=elts_dict, ranges=self.ranges ) @@ -798,10 +783,8 @@ def _plot_background(self, model): def _generate_ranges(self, num): axis = self.axes_manager.signal_axes[0] - bounds = (axis.low_value, axis.high_value) - values = np.linspace(bounds[0], bounds[1], num=2 * num + 2) - ranges_list = [(values[2 * i - 1], values[2 * i]) for i in range(1, num + 1)] - return ranges_list + values = np.linspace(axis.low_value, axis.high_value, num=2 * num + 2) + return values[1:-1].reshape(num, 2).tolist() def decomposition( self, @@ -996,8 +979,7 @@ def concentration_report(self, selected_elts=[], W_input=None, fit_error=True): G = self.learning_results.decomposition_algorithm.G_ H = self.learning_results.decomposition_algorithm.H_ N = get_explained_intensity_W(G, W, H) - sqN = np.sqrt(N) - percentages = sqN / N * 100 + percentages = 100.0 / np.sqrt(N) else: W = W_input @@ -1056,44 +1038,39 @@ def estimate_best_binning(self, inspect=False): L = self.axes_manager[2].size K = self.axes_manager[0].size * self.axes_manager[1].size - facx = np.arange(1, self.axes_manager[0].size // 2 + 1) - facy = np.arange(1, self.axes_manager[1].size // 2 + 1) + size_x = self.axes_manager[0].size + size_y = self.axes_manager[1].size + + facx = np.arange(1, size_x // 2 + 1) + facy = np.arange(1, size_y // 2 + 1) - binx = [self.axes_manager[0].size / i for i in facx] - biny = [self.axes_manager[1].size / i for i in facy] - vars_est = np.array([]) - biases_est = np.array([]) - bprod = [] - for i in zip(binx, biny): - bprod.append((i[0], i[1])) + binx = size_x / facx + biny = size_y / facy - for i in tqdm(bprod): + vars_est = [] + biases_est = [] + + bprod = list(zip(binx, biny)) + data = self.data + mean_d = np.mean(data) + + for i0, i1 in tqdm(bprod): # Bin the measurement dataset and upsample to bring it back to its original dimensionality - B = i[0] * i[1] - binned = self.rebin(scale=(i[0], i[1], 1)) - upsampled = binned.rebin( - new_shape=( - self.axes_manager[0].size, - self.axes_manager[1].size, - self.axes_manager[2].size, - ) - ) + B = i0 * i1 + binned = self.rebin(scale=(i0, i1, 1)) + upsampled = binned.rebin(new_shape=(size_x, size_y, L)) upsampled_data = upsampled.data - data = self.data # Estimator of variance (Lemma 4.3) - \widehat{Var} (\hat{y}_i) = \alpha ^2 y_{i}+ (1-\alpha)^2 \sum_{k \in \mathcal{K}} (w_k^2 n_{i,k}) - vars_est = np.append(vars_est, np.mean(upsampled_data * 1 / B)) + vars_est.append(np.mean(upsampled_data) / B) # Estimator of squared bias (Lemma 4.4) - \widehat{Bias^2}(\hat{y}_i) = (1-\alpha)^2\left((y_{n_i} - y_{i})^2 - \sum_{k\in \mathcal{K}} w_k^2y_{i,k} - y_{i} \right) - biases_est = np.append( - biases_est, - np.mean( - (data - upsampled_data) ** 2 - - 1 / B * upsampled_data - - (1 - 2 / B) * data - ), - ) + mean = np.mean((data - upsampled_data) ** 2 - upsampled_data / B) + + biases_est.append(mean - (1 - 2 / B) * mean_d) + vars_est = np.array(vars_est) + biases_est = np.array(biases_est) mprimes_est = vars_est * K / L + biases_est estimated_binning = ( bprod[np.argmin(mprimes_est)][0], @@ -1159,8 +1136,8 @@ def generate_part_fixed_H_matrix( part_f_H : np.ndarray A fixed H matrix for one phase. """ - part_f_H = (-1) * np.ones( - shape=(self.data.shape[0], self.data.shape[1]), dtype=float + part_f_H = np.full( + shape=(self.data.shape[0], self.data.shape[1]), fill_value=-1.0, dtype=float ) if value < 0: @@ -1205,13 +1182,7 @@ def set_fixed_H(self, areas_dict: dict[str, np.ndarray]) -> np.ndarray: A fixed H matrix for the SmoothNMF decomposition algorithm. """ - H = (-1) * np.ones( - shape=(len(areas_dict), self.data.shape[0], self.data.shape[1]), dtype=float - ) - - for i, p in enumerate(areas_dict): - H[i, :, :] = areas_dict[p] - + H = np.array(list(areas_dict.values())) return H.reshape((len(areas_dict), self.data.shape[0] * self.data.shape[1])) @_check_decomposition @@ -1247,20 +1218,16 @@ def num_elts(elements=[]): _elts_indices = self.model.NMF_simplex() skipped_elts = sym_elts(elements=skipped_elements) - elts = [] - elts_indices = [] - for i, elt in enumerate(_elts): - if elt in skipped_elts: - pass - else: - elts.append(elt) - elts_indices.append(_elts_indices[i]) + elts = [elt for elt in _elts if elt not in skipped_elts] + elts_indices = [ + idx for idx, elt in zip(_elts_indices, _elts) if elt not in skipped_elts + ] W = _W[elts_indices, :] WH = W @ H WH = WH.reshape([W.shape[0]] + list(self.data.shape[:-1])) - WH /= WH.sum(0)[np.newaxis, ...] / 100 + WH /= WH.sum(axis=0, keepdims=True) / 100.0 if self.axes_manager.navigation_dimension == 2: Signal = hs.signals.Signal2D @@ -1282,7 +1249,8 @@ def num_elts(elements=[]): self.axes_manager[i], ["units", "scale", "name", "offset"] ) q.metadata.Signal.quantity = "Atomic %" - wh = Signal(WH) + + wh = Signal(WH) for i in range(self.axes_manager.navigation_dimension): wh.axes_manager[1 + i].update_from( @@ -1954,9 +1922,7 @@ def plot_table( hover_data = [] cmap = plt.get_cmap("tab10") - # TODO: use cache when optimisation is merged - with open(SYMBOLS_PERIODIC_TABLE, "r") as f: - SPT = json.load(f)["table"] + SPT = get_spt() for i, elt in enumerate(elements): if (anum := str(SPT[elt]["number"])) not in table: diff --git a/espm/estimators/base.py b/espm/estimators/base.py index 04fdb2c..598241b 100644 --- a/espm/estimators/base.py +++ b/espm/estimators/base.py @@ -261,13 +261,11 @@ def fit_transform(self, X, y=None, W=None, H=None): # Initialize the algorithm # ############################ - if self.hspy_comp: - self.X_ = validate_data(self, X.T, dtype=[np.float64, np.float32]) - - else: - self.X_ = validate_data(self, X, dtype=[np.float64, np.float32]) + self.X_ = validate_data( + self, X.T if self.hspy_comp else X, dtype=[np.float64, np.float32] + ) - if self.hspy_comp == False: + if not self.hspy_comp: try: import inspect @@ -519,11 +517,8 @@ def get_losses(self): For debug purposes : return the evolution of losses. """ if self.true_D is not None and self.true_H is not None: - mse_list = [] - angles_list = [] - for i in range(self.n_components): - angles_list.append(f"ang_p{i}") - mse_list.append(f"mse_p{i}") + mse_list = [f"mse_p{i}" for i in range(self.n_components)] + angles_list = [f"ang_p{i}" for i in range(self.n_components)] names = ( ["full_loss"] + self.loss_names_ @@ -533,39 +528,33 @@ def get_losses(self): + ["true_KL_loss"] ) - dt_list = [] - for elt in names: - dt_list.append((elt, "float64")) - dt = np.dtype(dt_list) - - tup_list = [] - for i in range(len(self.losses_)): - tup_list.append( - (self.losses_[i],) - + tuple(self.detailed_losses_[i]) - + tuple(self.rel_[i]) - + tuple(self.angles_[i]) - + tuple(self.mse_[i]) - + (self.true_losses_[i],) + dt = np.dtype([(elt, "float64") for elt in names]) + + tup_list = [ + (loss, *det_loss, *rel, *ang, *m, true_loss) + for loss, det_loss, rel, ang, m, true_loss in zip( + self.losses_, + self.detailed_losses_, + self.rel_, + self.angles_, + self.mse_, + self.true_losses_, ) + ] array = np.array(tup_list, dtype=dt) # TODO : Check this part for optimization else: names = ["full_loss"] + self.loss_names_ + ["rel_W", "rel_H"] - dt_list = [] - for elt in names: - dt_list.append((elt, "float64")) - dt = np.dtype(dt_list) - - tup_list = [] - for i in range(len(self.losses_)): - tup_list.append( - (self.losses_[i],) - + tuple(self.detailed_losses_[i]) - + tuple(self.rel_[i]) + dt = np.dtype([(elt, "float64") for elt in names]) + + tup_list = [ + (loss, *det_loss, *rel) + for loss, det_loss, rel in zip( + self.losses_, self.detailed_losses_, self.rel_ ) + ] array = np.array(tup_list, dtype=dt) @@ -576,8 +565,8 @@ def remove_zeros_lines(self, X, epsilon): new_X = X.copy() sum_cols = X.sum(axis=0) sum_rows = X.sum(axis=1) - new_X[:, np.where(sum_cols == 0)] = epsilon - new_X[np.where(sum_rows == 0), :] = epsilon + new_X[:, sum_cols == 0] = epsilon + new_X[sum_rows == 0, :] = epsilon return new_X else: raise ValueError("Negative values in data") diff --git a/espm/estimators/dicotomy.py b/espm/estimators/dicotomy.py index 6e27cc9..9077f92 100644 --- a/espm/estimators/dicotomy.py +++ b/espm/estimators/dicotomy.py @@ -4,7 +4,7 @@ def dichotomy_simplex( - num, denum, log_shift=log_shift, tol=dicotomy_tol, maxit=maxit_dichotomy + num, denum, log_shift=log_shift, tol=dicotomy_tol, maxit=maxit_dichotomy, safe=False ): """ Function to solve the num/(x+denum) -1 = 0 equation. Here, x is the Lagragian multiplier which is used to apply the simplex constraint. @@ -18,9 +18,10 @@ def dichotomy_simplex( # do some test - assert (num >= 0).all() - assert (denum >= 0).all() - assert (np.sum(num, axis=0) > 0).all() + if safe: + assert (num >= 0).all() + assert (denum >= 0).all() + assert (np.sum(num, axis=0) > 0).all() if log_shift > 0: # Check that a solution is possible if denum.shape[0] * log_shift >= 1: @@ -30,23 +31,9 @@ def dichotomy_simplex( denum_max = np.inf # Ideally we want to do this, but we have to exclude the case where num==0. # a = np.max(num/2 - denum, axis=0) - if denum.shape[1] > 1: - a = [] - for n, d in zip(num.T, denum.T): - m = n > 0 - # The divided by 2 is just a factor to help a bit. - a.append(np.max(n[m] / 2 - d[m])) - a = np.array(a) - else: - # This else is just to preserve the size and make it work in any case... - # There might be a possiblity to write this more elegantly - d = denum[:, 0] - - def max_masked(n): - m = n > 0 - return np.max(n[m] / 2 - d[m]) - - a = np.apply_along_axis(max_masked, 0, num) + val = num / 2 - denum + val[num <= 0] = -np.inf + a = np.max(val, axis=0) # r = np.sum(num/denum, axis=0) # b = np.zeros(r.shape) @@ -57,11 +44,17 @@ def func(x): new_x = x + denum return np.sum(np.maximum(num / new_x, log_shift), axis=0) - 1 - return dicotomy(a, b, func, maxit, tol) + return dicotomy(a, b, func, maxit, tol, safe=safe) def dichotomy_simplex_acc( - a, b, minus_c, log_shift=log_shift, tol=dicotomy_tol, maxit=maxit_dichotomy + a, + b, + minus_c, + log_shift=log_shift, + tol=dicotomy_tol, + maxit=maxit_dichotomy, + safe=False, ): """ Function to solve the dicotomy for the function: @@ -71,8 +64,9 @@ def dichotomy_simplex_acc( The second part applies the dichotomy algorithm to solve the equation. """ # do some test - assert a >= 0 - assert (minus_c >= 0).all() + if safe: + assert a >= 0 + assert (minus_c >= 0).all() if log_shift > 0: # Check that a solution is possible @@ -91,11 +85,11 @@ def func(x): axis=0, ) - return dicotomy(nu_max, nu_min, func, maxit, tol) + return dicotomy(nu_max, nu_min, func, maxit, tol, safe=safe) def dichotomy_simplex_projected_gradient( - a, log_shift=log_shift, tol=dicotomy_tol, maxit=maxit_dichotomy + a, log_shift=log_shift, tol=dicotomy_tol, maxit=maxit_dichotomy, safe=False ): r""" Function to solve the dicotomy for the function: @@ -119,19 +113,19 @@ def dichotomy_simplex_projected_gradient( def func(x): return np.sum(np.maximum(a + x, log_shift), axis=0) - 1 - return dicotomy(nu_max, nu_min, func, maxit, tol) + return dicotomy(nu_max, nu_min, func, maxit, tol, safe=safe) -def dicotomy(a, b, func, maxit, tol): +def dicotomy(a, b, func, maxit, tol, safe=False): """ Dicotomy algorithm searching for func(x)=0. Parameters ---------- - a : float or numpy array + a : numpy array Lower bound of the interval such that func(a) > 0 - b : float or numpy array + b : numpy array Upper bound of the interval such that func(b) < 0 func : function Function to solve @@ -152,10 +146,11 @@ def dicotomy(a, b, func, maxit, tol): func_max = func(a) func_min = func(b) - assert np.sum(func_min >= 0) == 0 - assert np.sum(func_max <= 0) == 0 - assert np.sum(np.isnan(func_max)) == 0 - assert np.sum(np.isnan(func_min)) == 0 + if safe: + assert not np.any(func_min >= 0) + assert not np.any(func_max <= 0) + assert not np.any(np.isnan(func_max)) + assert not np.any(np.isnan(func_min)) # Dichotomy algorithm to solve the equation it = 0 @@ -163,20 +158,19 @@ def dicotomy(a, b, func, maxit, tol): func_new = func(new) # print("A : {}, B: {}, new : {}, fA : {}, fB : {}, fnew : {}".format(np.min(np.abs(a)),np.min(np.abs(b)),np.min(np.abs(new)),np.min(np.abs(func(a))),np.min(np.abs(func(b))),np.min(np.abs(func(new))))) # print("A : {}, B: {}, new : {}, fA : {}, fB : {}, fnew : {}".format(np.max(a),np.max(b),np.max(new),np.max(func(a)),np.max(func(b)),np.max(func(new)))) + func_a = func_max + while np.max(np.abs(func_new)) > tol: it = it + 1 - func_a = func(a) - # func_b = func(b) - # if f(a)*f(new) <0 then f(new) < 0 --> store in b minus_bool = func_a * func_new <= 0 # if f(a)*f(new) > 0 then f(new) > 0 --> store in a # plus_bool = func_a * func_new > 0 - plus_bool = np.logical_not(minus_bool) + b = np.where(minus_bool, new, b) + a = np.where(minus_bool, a, new) + func_a = np.where(minus_bool, func_a, func_new) - b[minus_bool] = new[minus_bool] - a[plus_bool] = new[plus_bool] new = (a + b) / 2 func_new = func(new) if it >= maxit: diff --git a/espm/estimators/smooth_nmf.py b/espm/estimators/smooth_nmf.py index 15da68a..0dbc97f 100644 --- a/espm/estimators/smooth_nmf.py +++ b/espm/estimators/smooth_nmf.py @@ -325,185 +325,194 @@ def _iteration(self, W, H): if self.n_iter_ == 0: if self.gamma is None: - if self.algo in ["l2_surrogate", "log_surrogate", "bmd"]: - self.gamma_ = sigmaL - else: - gamma_W = estimate_Lipschitz_bound_w( - self.log_shift, self.X_, self.G_, k=self.n_components - ) - gamma_H = estimate_Lipschitz_bound_h( - self.log_shift, - self.X_, - self.G_, - k=self.n_components, - lambda_L=self.lambda_L, - mu=self.mu, - epsilon_reg=self.epsilon_reg, - ) - self.gamma_ = [gamma_H, gamma_W] + match self.algo: + case "l2_surrogate" | "log_surrogate" | "bmd": + self.gamma_ = sigmaL + case _: + gamma_W = estimate_Lipschitz_bound_w( + self.log_shift, self.X_, self.G_, k=self.n_components + ) + gamma_H = estimate_Lipschitz_bound_h( + self.log_shift, + self.X_, + self.G_, + k=self.n_components, + lambda_L=self.lambda_L, + mu=self.mu, + epsilon_reg=self.epsilon_reg, + ) + self.gamma_ = [gamma_H, gamma_W] else: self.gamma_ = deepcopy(self.gamma) # 1. Update for H if self.linesearch: Hold = H.copy() - if self.algo == "l2_surrogate": - H = multiplicative_step_hq( - self.X_, - self.G_, - W, - H, - simplex_H=self.simplex_H, - log_shift=self.log_shift, - safe=self.debug, - dicotomy_tol=self.dicotomy_tol, - lambda_L=self.lambda_L, - L=self.L_, - sigmaL=self.gamma_, - fixed_H=self.fixed_H, - ) - elif self.algo == "log_surrogate": - H = multiplicative_step_h( - self.X_, - self.G_, - W, - H, - simplex_H=self.simplex_H, - mu=self.mu, - log_shift=self.log_shift, - epsilon_reg=self.epsilon_reg, - safe=self.debug, - dicotomy_tol=self.dicotomy_tol, - lambda_L=self.lambda_L, - L=self.L_, - l2=self.l2, - fixed_H=self.fixed_H, - sigmaL=self.gamma_, - ) - elif self.algo == "projected_gradient": - H = proj_grad_step_h( - self.X_, - self.G_, - W, - H, - simplex_H=self.simplex_H, - mu=self.mu, - log_shift=self.log_shift, - epsilon_reg=self.epsilon_reg, - safe=self.debug, - dicotomy_tol=self.dicotomy_tol, - lambda_L=self.lambda_L, - L=self.L_, - l2=self.l2, - fixed_H=self.fixed_H, - gamma=self.gamma_[0], - ) - elif self.algo == "bmd": - H = multiplicative_step_h( - self.X_, - self.G_, - W, - H, - simplex_H=self.simplex_H, - mu=self.mu, - log_shift=self.log_shift, - epsilon_reg=self.epsilon_reg, - safe=self.debug, - dicotomy_tol=self.dicotomy_tol, - lambda_L=self.lambda_L, - L=self.L_, - l2=self.l2, - fixed_H=self.fixed_H, - sigmaL=self.gamma_, - use_bregman=True, - ) - else: - raise ValueError("Unknown algorithm") - - if self.linesearch: - if self.algo in ["l2_surrogate", "log_surrogate", "bmd"]: - d = diff_surrogate( - Hold, H, L=self.L_, sigmaL=self.gamma_, algo=self.algo + match self.algo: + case "l2_surrogate": + H = multiplicative_step_hq( + self.X_, + self.G_, + W, + H, + simplex_H=self.simplex_H, + log_shift=self.log_shift, + safe=self.debug, + dicotomy_tol=self.dicotomy_tol, + lambda_L=self.lambda_L, + L=self.L_, + sigmaL=self.gamma_, + fixed_H=self.fixed_H, ) - if d > 0: - self.gamma_ = self.gamma_ / 1.05 - else: - self.gamma_ = self.gamma_ * 1.5 - else: - gradf_xt = gradH( + case "log_surrogate": + H = multiplicative_step_h( self.X_, self.G_, W, - Hold, + H, + simplex_H=self.simplex_H, mu=self.mu, + log_shift=self.log_shift, + epsilon_reg=self.epsilon_reg, + safe=self.debug, + dicotomy_tol=self.dicotomy_tol, lambda_L=self.lambda_L, L=self.L_, + l2=self.l2, + fixed_H=self.fixed_H, + sigmaL=self.gamma_, + ) + case "projected_gradient": + H = proj_grad_step_h( + self.X_, + self.G_, + W, + H, + simplex_H=self.simplex_H, + mu=self.mu, + log_shift=self.log_shift, epsilon_reg=self.epsilon_reg, + safe=self.debug, + dicotomy_tol=self.dicotomy_tol, + lambda_L=self.lambda_L, + L=self.L_, + l2=self.l2, + fixed_H=self.fixed_H, + gamma=self.gamma_[0], + ) + case "bmd": + H = multiplicative_step_h( + self.X_, + self.G_, + W, + H, + simplex_H=self.simplex_H, + mu=self.mu, log_shift=self.log_shift, + epsilon_reg=self.epsilon_reg, safe=self.debug, + dicotomy_tol=self.dicotomy_tol, + lambda_L=self.lambda_L, + L=self.L_, + l2=self.l2, + fixed_H=self.fixed_H, + sigmaL=self.gamma_, + use_bregman=True, ) - f_xt = self.loss(W, Hold, X=self.X_, average=False) - f_x = self.loss(W, H, X=self.X_, average=False) - g_xxt = quadratic_surrogate(H, Hold, f_xt, gradf_xt, self.gamma_[0]) - d = g_xxt - f_x - if d > 0: - self.gamma_[0] = self.gamma_[0] / 1.05 - else: - self.gamma_[0] = self.gamma_[0] * 1.5 + case _: + raise ValueError("Unknown algorithm") + + if self.linesearch: + match self.algo: + case "l2_surrogate" | "log_surrogate" | "bmd": + d = diff_surrogate( + Hold, H, L=self.L_, sigmaL=self.gamma_, algo=self.algo + ) + if d > 0: + self.gamma_ = self.gamma_ / 1.05 + else: + self.gamma_ = self.gamma_ * 1.5 + case _: + gradf_xt = gradH( + self.X_, + self.G_, + W, + Hold, + mu=self.mu, + lambda_L=self.lambda_L, + L=self.L_, + epsilon_reg=self.epsilon_reg, + log_shift=self.log_shift, + safe=self.debug, + ) + f_xt = self.loss(W, Hold, X=self.X_, average=False) + f_x = self.loss(W, H, X=self.X_, average=False) + g_xxt = quadratic_surrogate(H, Hold, f_xt, gradf_xt, self.gamma_[0]) + d = g_xxt - f_x + if d > 0: + self.gamma_[0] = self.gamma_[0] / 1.05 + else: + self.gamma_[0] = self.gamma_[0] * 1.5 # 2. Update for W - if self.algo in ["l2_surrogate", "log_surrogate"]: - W = multiplicative_step_w( - self.X_, - self.G_, - W, - H, - log_shift=self.log_shift, - safe=self.debug, - l2=self.l2, - simplex_W=self.simplex_W, - fixed_W=self.fixed_W, - physics_model=self.physics_model_, - ) - elif self.algo == "bmd": - W = multiplicative_step_w( - self.X_, - self.G_, - W, - H, - log_shift=self.log_shift, - safe=self.debug, - l2=self.l2, - simplex_W=self.simplex_W, - fixed_W=self.fixed_W, - use_bregman=True, - physics_model=self.physics_model_, - ) - else: - if self.linesearch: - Wold = W.copy() - W = proj_grad_step_w( - self.X_, - self.G_, - W, - H, - log_shift=self.log_shift, - safe=self.debug, - gamma=self.gamma_[1], - simplex_W=self.simplex_W, - ) - if self.linesearch: - gradf_xt = gradW( - self.X_, self.G_, Wold, H, log_shift=self.log_shift, safe=self.debug + match self.algo: + case "l2_surrogate" | "log_surrogate": + W = multiplicative_step_w( + self.X_, + self.G_, + W, + H, + log_shift=self.log_shift, + safe=self.debug, + l2=self.l2, + simplex_W=self.simplex_W, + fixed_W=self.fixed_W, + physics_model=self.physics_model_, + ) + case "bmd": + W = multiplicative_step_w( + self.X_, + self.G_, + W, + H, + log_shift=self.log_shift, + safe=self.debug, + l2=self.l2, + simplex_W=self.simplex_W, + fixed_W=self.fixed_W, + use_bregman=True, + physics_model=self.physics_model_, ) - f_xt = self.loss(Wold, H, X=self.X_, average=False) - f_x = self.loss(W, H, X=self.X_, average=False) - g_xxt = quadratic_surrogate(W, Wold, f_xt, gradf_xt, self.gamma_[1]) - d = g_xxt - f_x - if d > 0: - self.gamma_[1] = self.gamma_[1] / 1.05 - else: - self.gamma_[1] = self.gamma_[1] * 1.5 + case _: + if self.linesearch: + Wold = W.copy() + W = proj_grad_step_w( + self.X_, + self.G_, + W, + H, + log_shift=self.log_shift, + safe=self.debug, + gamma=self.gamma_[1], + simplex_W=self.simplex_W, + ) + if self.linesearch: + gradf_xt = gradW( + self.X_, + self.G_, + Wold, + H, + log_shift=self.log_shift, + safe=self.debug, + ) + f_xt = self.loss(Wold, H, X=self.X_, average=False) + f_x = self.loss(W, H, X=self.X_, average=False) + g_xxt = quadratic_surrogate(W, Wold, f_xt, gradf_xt, self.gamma_[1]) + d = g_xxt - f_x + if d > 0: + self.gamma_[1] = self.gamma_[1] / 1.05 + else: + self.gamma_[1] = self.gamma_[1] * 1.5 # KL_surr = KL_loss_surrogate(self.X_, W, H, Hold, eps=0) # log_surr = log_surrogate(H, Hold, mu=self.mu, epsilon=self.epsilon_reg) diff --git a/espm/estimators/surrogates.py b/espm/estimators/surrogates.py index 303415a..a0ed689 100644 --- a/espm/estimators/surrogates.py +++ b/espm/estimators/surrogates.py @@ -42,13 +42,14 @@ def smooth_l2_surrogate(Ht, L, H=None, sigmaL=sigmaL, lambda_L=1): """ HtTL = Ht @ L - t1 = np.sum(HtTL * Ht) + t1 = np.vdot(HtTL, Ht) if H is None: t2 = t1 t3 = 0 else: - t2 = np.sum(HtTL * H) - t3 = np.sum((Ht - H) ** 2) + t2 = np.vdot(HtTL, H) + diff = Ht - H + t3 = np.vdot(diff, diff) return lambda_L / 2 * (2 * t2 - t1 + sigmaL * t3) @@ -101,7 +102,7 @@ def smooth_dgkl_surrogate(Ht, L, H=None, sigmaL=sigmaL, lambda_L=1): """ HtTL = Ht @ L - t1 = np.sum(HtTL * Ht) + t1 = np.vdot(HtTL, Ht) def dgkl(p, q): return p * np.log(p / q) - p + q @@ -110,7 +111,7 @@ def dgkl(p, q): t2 = t1 t3 = 0 else: - t2 = np.sum(HtTL * H) + t2 = np.vdot(HtTL, H) maxH = np.max(H, axis=1) t3 = np.sum(maxH * np.sum(dgkl(Ht, H), axis=1)) return lambda_L / 2 * (2 * t2 - t1 + sigmaL * t3) @@ -170,4 +171,5 @@ def quadratic_surrogate(x, xt, f_xt, gradf_xt, sigma): :returns: the answer """ - return f_xt + np.sum((x - xt) * gradf_xt) + sigma * np.sum((x - xt) ** 2) + diff = x - xt + return f_xt + np.vdot(diff, gradf_xt) + sigma * np.vdot(diff, diff) diff --git a/espm/estimators/updates.py b/espm/estimators/updates.py index e2bec65..bc9b1d4 100644 --- a/espm/estimators/updates.py +++ b/espm/estimators/updates.py @@ -27,9 +27,9 @@ def multiplicative_step_w( """ if safe: # Allow for very small negative values! - assert np.sum(H < -log_shift / 2) == 0 - assert np.sum(W < -log_shift / 2) == 0 - assert np.sum(G < -log_shift / 2) == 0 + assert not np.any(H < -log_shift / 2) + assert not np.any(W < -log_shift / 2) + assert not np.any(G < -log_shift / 2) H = np.maximum(H, log_shift) W = np.maximum(W, log_shift) @@ -37,7 +37,7 @@ def multiplicative_step_w( if l2: GG = G.T @ G HH = H @ H.T - GGWHH = GG @ W @ HH + GGWHH = (GG @ W) @ HH GXH = G.T @ (X @ H.T) @@ -47,15 +47,13 @@ def multiplicative_step_w( GWH = GW @ H if use_bregman: # check if G is the identity matrix - if np.allclose(G, np.eye(G.shape[0])): + if G.shape[0] == G.shape[1] and np.array_equal(G, np.eye(G.shape[0])): sigmaR = np.sum(X, axis=1, keepdims=True) else: sigmaR = np.sum(X) num = sigmaR * W - gradg = ( - -G.T @ (X / GWH) @ H.T - + np.sum(G, axis=0, keepdims=True).T - @ np.sum(H, axis=1, keepdims=True).T + gradg = -G.T @ ((X / GWH) @ H.T) + np.outer( + np.sum(G, axis=0), np.sum(H, axis=1) ) denum = gradg * W + sigmaR @@ -67,24 +65,22 @@ def multiplicative_step_w( GWH = np.maximum(GWH, log_shift) op1 = X / GWH - mult1 = G.T @ op1 - num = W * (mult1 @ H.T) - denum = ( - np.sum(G, axis=0, keepdims=True).T @ np.sum(H, axis=1, keepdims=True).T - ) + num = W * (G.T @ (op1 @ H.T)) + denum = np.outer(np.sum(G, axis=0), np.sum(H, axis=1)) if simplex_W: - if physics_model != None: + if physics_model is not None: indices = physics_model.NMF_simplex() nu = dichotomy_simplex( num[indices, :], denum[indices, :], log_shift=log_shift, tol=dicotomy_tol, + safe=safe, ) denum[indices, :] = denum[indices, :] + nu else: nu = dichotomy_simplex( - num, denum, log_shift=log_shift, tol=dicotomy_tol + num, denum, log_shift=log_shift, tol=dicotomy_tol, safe=safe ) denum = denum + nu @@ -126,7 +122,7 @@ def multiplicative_step_h( as a vector to regularize the different phase of A differently. To calculate the regularized step, we make a linear approximation of the log. """ - if not (lambda_L == 0): + if lambda_L != 0: if L is None: raise ValueError("Please provide the laplacian") HL = H @ L @@ -134,9 +130,9 @@ def multiplicative_step_h( if safe: # Allow for very small negative values! # TODO: update this - assert np.sum(H < -log_shift / 2) == 0 - assert np.sum(W < -log_shift / 2) == 0 - assert np.sum(G < -log_shift / 2) == 0 + assert not np.any(H < -log_shift / 2) + assert not np.any(W < -log_shift / 2) + assert not np.any(G < -log_shift / 2) H = np.maximum(H, log_shift) W = np.maximum(W, log_shift) @@ -147,7 +143,7 @@ def multiplicative_step_h( if np.isscalar(mu): assert mu == 0 else: - assert (mu == 0).all() + assert not np.any(mu != 0) WGGW = GW.T @ GW WGX = GW.T @ X num = WGX @@ -171,18 +167,20 @@ def multiplicative_step_h( if len(np.shape(mu)) == 1: mu = np.expand_dims(mu, axis=1) denum = denum + mu / (H + epsilon_reg) - if not (lambda_L == 0): + if lambda_L != 0: maxH = np.max(H, axis=1, keepdims=True) num = num + lambda_L * sigmaL * maxH denum = denum + lambda_L * sigmaL * maxH + lambda_L * HL num = H * num if simplex_H: - nu = dichotomy_simplex(num, denum, log_shift=log_shift, tol=dicotomy_tol) + nu = dichotomy_simplex( + num, denum, log_shift=log_shift, tol=dicotomy_tol, safe=safe + ) else: nu = 0 if safe: - assert np.sum(denum < 0) == 0 - assert np.sum(num < 0) == 0 + assert not np.any(denum < 0) + assert not np.any(num < 0) # Add the shift... new_H = np.maximum(num / (denum + nu), log_shift) @@ -210,7 +208,7 @@ def initialize_algorithms( if G is None: skip_second = True # G = sparse.diags(np.ones(X.shape[0]).astype(X.dtype)) - G = np.diag(np.ones(X.shape[0]).astype(X.dtype)) + G = np.eye(X.shape[0], dtype=X.dtype) # elif callable(G) : # assert not(model_params is None), "You need to input model_parameters" @@ -239,7 +237,7 @@ def initialize_algorithms( if skip_second: W = D else: - if physics_model != None: + if physics_model is not None: # [np.where(G[:,:-2].sum(axis=1)<(np.max(G[:,:-2].sum(axis=1))*0.001))[0],:] # Divide in two parts the initial fitting, otherwise the bremsstrahlung (which has a low intensity) tends to be poorly learned # First fit the caracteristic Xrays, then subtract that contribution to obtain a rough estimate of the bremsstralung parameters @@ -274,10 +272,9 @@ def initialize_algorithms( def update_q(D, H, log_shift=log_shift): """Perform a Q step.""" - Htmp = np.expand_dims(H.T, axis=0) - Dtmp = np.expand_dims(D, axis=1) - Ntmp = np.expand_dims(D @ H, axis=2) - return Htmp * (Dtmp / (Ntmp + log_shift)) + return H.T[np.newaxis, :, :] * ( + D[:, np.newaxis, :] / ((D @ H + log_shift)[:, :, np.newaxis]) + ) def multiplicative_step_wq( @@ -291,30 +288,31 @@ def multiplicative_step_wq( if safe: # Allow for very small negative values! - assert np.sum(H < -log_shift / 2) == 0 - assert np.sum(W < -log_shift / 2) == 0 - assert np.sum(G < -log_shift / 2) == 0 + assert not np.any(H < -log_shift / 2) + assert not np.any(W < -log_shift / 2) + assert not np.any(G < -log_shift / 2) GW = G @ W - Q = update_q(GW, H, log_shift=log_shift) - - XQ = np.sum(np.expand_dims(X, axis=2) * Q, axis=1) + XQ = GW * ((X / (GW @ H + log_shift)) @ H.T) term1 = G.T @ (XQ / (GW + log_shift)) - term2 = np.sum(G, axis=0, keepdims=True).T @ np.sum(H, axis=1, keepdims=True).T + term2 = np.outer(np.sum(G, axis=0), np.sum(H, axis=1)) if simplex_W: - if physics_model != None: + if physics_model is not None: indices = physics_model.NMF_simplex() nu = dichotomy_simplex( term1[indices, :], term2[indices, :], log_shift=log_shift, tol=dicotomy_tol, + safe=safe, ) term2[indices, :] = term2[indices, :] + nu else: - nu = dichotomy_simplex(term1, term2, log_shift=log_shift, tol=dicotomy_tol) + nu = dichotomy_simplex( + term1, term2, log_shift=log_shift, tol=dicotomy_tol, safe=safe + ) term2 = term2 + nu return W / term2 * term1 @@ -336,15 +334,15 @@ def multiplicative_step_hq( """ Multiplicative step in H. """ - if not lambda_L == 0: + if lambda_L != 0: if L is None: raise ValueError("Please provide the laplacian") if safe: # Allow for very small negative values! - assert np.sum(H < -log_shift / 2) == 0 - assert np.sum(W < -log_shift / 2) == 0 - assert np.sum(G < -log_shift / 2) == 0 + assert not np.any(H < -log_shift / 2) + assert not np.any(W < -log_shift / 2) + assert not np.any(G < -log_shift / 2) GW = G @ W # Also called D GWH = GW @ H @@ -352,18 +350,20 @@ def multiplicative_step_hq( minus_c = H * (GW.T @ (X / (GWH + log_shift))) b = np.sum(GW, axis=0, keepdims=True).T - if not lambda_L == 0: + if lambda_L != 0: b = b + lambda_L * H @ L - lambda_L * sigmaL * H a = lambda_L * sigmaL if simplex_H: nu = dichotomy_simplex_acc( - a, b, minus_c, log_shift=log_shift, tol=dicotomy_tol + a, b, minus_c, log_shift=log_shift, tol=dicotomy_tol, safe=safe ) b = b + nu new_H = (-b + np.sqrt(b**2 + 4 * a * minus_c)) / (2 * a) else: # We recover the classic case: multiplicative_step_a if simplex_H: - nu = dichotomy_simplex(minus_c, b, log_shift=log_shift, tol=dicotomy_tol) + nu = dichotomy_simplex( + minus_c, b, log_shift=log_shift, tol=dicotomy_tol, safe=safe + ) b = b + nu new_H = minus_c / b @@ -401,30 +401,28 @@ def gradH( safe=False, l2=False, ): - if not (lambda_L == 0): + if lambda_L != 0: if L is None: raise ValueError("Please provide the laplacian") - HL = H @ L if safe: H = np.maximum(H, log_shift) W = np.maximum(W, log_shift) + D = G @ W if l2: - D = G @ W grad = D.T @ (D @ H - X) else: - D = G @ W DH = D @ H grad = -D.T @ (X / DH) + np.sum(D, axis=0, keepdims=True).T if not (np.isscalar(mu) and mu == 0): if len(np.shape(mu)) == 1: mu = np.expand_dims(mu, axis=1) - grad += mu / (H + epsilon_reg) + grad = grad + mu / (H + epsilon_reg) - if not (lambda_L == 0): - grad += (lambda_L * L @ H.T).T + if lambda_L != 0: + grad = grad + (lambda_L * L @ H.T).T return grad @@ -506,7 +504,7 @@ def proj_grad_step_h( # Dichotomy if simplex_H: nu = dichotomy_simplex_projected_gradient( - new_H, log_shift=log_shift, tol=dicotomy_tol + new_H, log_shift=log_shift, tol=dicotomy_tol, safe=safe ) else: nu = 0 diff --git a/espm/measures.py b/espm/measures.py index 9e00892..f6cad35 100644 --- a/espm/measures.py +++ b/espm/measures.py @@ -506,11 +506,7 @@ def KLdiv_loss(X, W, H, log_shift=log_shift, average=False): 2.921251732961556 """ - W = np.maximum(W, log_shift) - H = np.maximum(H, log_shift) - X = np.maximum(X, log_shift) - - Y = W @ H + Y = np.maximum(W @ H, log_shift) if average: x_lin = np.mean(Y) x_log = np.mean(X * np.log(Y)) @@ -642,6 +638,6 @@ def squared_distance(x, y=None): xx = (x * x).sum(axis=1) yy = (y * y).sum(axis=1) xy = np.dot(x, y.T) - d = abs(np.kron(np.ones((ry, 1)), xx).T + np.kron(np.ones((rx, 1)), yy) - 2 * xy) + d = np.abs(xx[:, np.newaxis] + yy[np.newaxis, :] - 2 * xy) return d / cx diff --git a/espm/models/absorption_edxs.py b/espm/models/absorption_edxs.py index 74bbf18..355dc67 100644 --- a/espm/models/absorption_edxs.py +++ b/espm/models/absorption_edxs.py @@ -7,6 +7,7 @@ """ from pathlib import Path +from functools import lru_cache import numpy as np from scipy.interpolate import interp1d @@ -15,6 +16,13 @@ from espm.utils import approx_density, atomic_to_weight_dict, number_to_symbol_dict +@lru_cache(maxsize=128) +def _get_mac_interp_func(key): + x_db = HSPY_MAC[key]["energies (keV)"] + y_db = HSPY_MAC[key]["mass_absorption_coefficient (cm2/g)"] + return interp1d(x_db, y_db, kind="cubic") + + @number_to_symbol_dict def absorption_coefficient(x, atomic_fraction=False, *, elements_dict={"Si": 1.0}): r""" @@ -45,9 +53,7 @@ def absorption_coefficient(x, atomic_fraction=False, *, elements_dict={"Si": 1.0 sum_elts = sum(elements_dict.values()) for key in elements_dict.keys(): - x_db = HSPY_MAC[key]["energies (keV)"] - y_db = HSPY_MAC[key]["mass_absorption_coefficient (cm2/g)"] - interp_func = interp1d(x_db, y_db, kind="cubic") + interp_func = _get_mac_interp_func(key) mu += elements_dict[key] * interp_func(x) / sum_elts if len(elements_dict.keys()) == 0: @@ -149,6 +155,13 @@ def absorption_mass_thickness( return (1 - np.exp(-chi)) / chi +@lru_cache(maxsize=32) +def _get_det_efficiency_interp_func(filename, kind): + array = np.loadtxt(DB_PATH / Path(filename)) + x_curve, y_curve = array[:, 0], array[:, 1] + return interp1d(x_curve, y_curve, kind=kind) + + def det_efficiency_from_curve(x, filename, kind="cubic"): r""" Interpolate a detection efficiency curve that is stored in ~/espm/tables. @@ -167,9 +180,7 @@ def det_efficiency_from_curve(x, filename, kind="cubic"): Detection efficiency :np.array 1D: Interpolated detection efficiency. """ - array = np.loadtxt(DB_PATH / Path(filename)) - x_curve, y_curve = array[:, 0], array[:, 1] - interp_func = interp1d(x_curve, y_curve, kind=kind) + interp_func = _get_det_efficiency_interp_func(filename, kind) return interp_func(x) diff --git a/espm/utils.py b/espm/utils.py index 4682b5d..9386661 100644 --- a/espm/utils.py +++ b/espm/utils.py @@ -117,19 +117,29 @@ def bin_spim(data, n, m): # return a matrix of shape (n,m,k) bs = data.shape[0] // n, data.shape[1] // m # blocksize averaged over k = data.shape[2] - return np.reshape( - np.array( - [ - np.sum( - data[k1 * bs[0] : (k1 + 1) * bs[0], k2 * bs[1] : (k2 + 1) * bs[1]], - axis=(0, 1), - ) - for k1 in range(n) - for k2 in range(m) - ] - ), - (n, m, k), - ) + + cropped = data[: n * bs[0], : m * bs[1]] + return cropped.reshape(n, bs[0], m, bs[1], k).sum(axis=(1, 3)) + + +_NPT_CACHE = None +_SPT_CACHE = None + + +def get_npt(): + global _NPT_CACHE + if _NPT_CACHE is None: + with open(NUMBER_PERIODIC_TABLE, "r") as f: + _NPT_CACHE = json.load(f)["table"] + return _NPT_CACHE + + +def get_spt(): + global _SPT_CACHE + if _SPT_CACHE is None: + with open(SYMBOLS_PERIODIC_TABLE, "r") as f: + _SPT_CACHE = json.load(f)["table"] + return _SPT_CACHE def number_to_symbol_dict(func): @@ -143,8 +153,7 @@ def number_to_symbol_dict(func): def inner(*args, **kwargs): elts_dict = kwargs["elements_dict"] new_dict = {} - with open(NUMBER_PERIODIC_TABLE, "r") as f: - NPT = json.load(f)["table"] + NPT = get_npt() for key in elts_dict.keys(): if is_symbol(key): @@ -175,8 +184,7 @@ def symbol_to_number_dict(func): def inner(*args, **kwargs): elts_dict = kwargs["elements_dict"] new_dict = {} - with open(SYMBOLS_PERIODIC_TABLE, "r") as f: - SPT = json.load(f)["table"] + SPT = get_spt() for key in elts_dict.keys(): if is_number(key): new_dict[int(key)] = elts_dict[key] @@ -206,8 +214,7 @@ def symbol_to_number_list(func): def inner(*args, **kwargs): elts_list = kwargs["elements"] new_list = [] - with open(SYMBOLS_PERIODIC_TABLE, "r") as f: - SPT = json.load(f)["table"] + SPT = get_spt() for key in elts_list: if is_number(key): new_list.append(int(key)) @@ -235,8 +242,7 @@ def number_to_symbol_list(func): def inner(*args, **kwargs): elts_list = kwargs["elements"] new_list = [] - with open(NUMBER_PERIODIC_TABLE, "r") as f: - NPT = json.load(f)["table"] + NPT = get_npt() for key in elts_list: if is_number(key): new_list.append(NPT[str(key)]["symbol"])