diff --git a/espm/datasets/eds_spim.py b/espm/datasets/eds_spim.py index e6a955e..2310f26 100644 --- a/espm/datasets/eds_spim.py +++ b/espm/datasets/eds_spim.py @@ -26,20 +26,24 @@ from hyperspy.roi import BaseROI, RectangularROI from hyperspy.signal_tools import Signal1DRangeSelector from hyperspy.ui_registry import get_gui +from matplotlib.collections import LineCollection from prettytable import PrettyTable from scipy.optimize import curve_fit from tqdm import tqdm -from espm.conf import NUMBER_PERIODIC_TABLE +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, num_to_symbol, number_to_symbol_list, + quant_spectrum, + symbol_to_number_dict, symbol_to_number_list, ) +# TODO: use cache after merging optimisation NPT = json.load(open(NUMBER_PERIODIC_TABLE)) @@ -69,6 +73,10 @@ def __init__(self, *args, **kwargs): self.model_ = None self.custom_init_ = None self.ranges = None + self.average_spectrum_ = None + self.energy_axis_ = None + self.elements_ = None + self._set_default_analysis_params() ############## @@ -214,10 +222,39 @@ def model(self) -> EDXS: self.model_ = EDXS(**mod_pars, custom_init=self.custom_init_) return self.model_ + @property + def energy_axis(self): + if self.energy_axis_ is None: + self.energy_axis_ = self.axes_manager.signal_axes[0].axis + return self.energy_axis_ + + @property + def average_spectrum(self): + if self.average_spectrum_ is None: + nav_axes = self.axes_manager.navigation_axes + if len(nav_axes) > 0: + self.average_spectrum_ = self.mean(axis=nav_axes).data + else: + self.average_spectrum_ = self.data + return self.average_spectrum_ + + @property + def elements(self): + if self.elements_ is None: + + @symbol_to_number_list + def convert_to_numbers(elements): + return elements + + elements = convert_to_numbers(elements=self.metadata.Sample.elements) + self.elements_ = [str(el) for el in elements] + return self.elements_ + def build_G( self, problem_type: str = "bremsstrahlung", ignored_elements: list[str] = ["Cu"], + use_calibration: bool = False, *, elements_dict: dict[str, float] = {}, ) -> None: @@ -231,6 +268,10 @@ def build_G( - "bremsstrahlung" : the G matrix is a callable with both characteristic X-rays and a bremsstrahlung model. - "no_brstlg" : the G matrix is a matrix with only characteristic X-rays. - "identity" : the G matrix is None which is equivalent to an identity matrix for espm functions. + ignored_elements : list, optional + List of chemical elements to ignore when building the G matrix. + use_calibration : bool, optional + If True, build G using the polynomially calibrated table from `auto_calibrate`. elements_dict : dict, optional Dictionary containing atomic numbers and a corresponding cut-off energies. It is used to separate the characteristic X-rays of the given elements into two energies ranges and assign them each a column in the G matrix instead of having one column per element. For example elements_dict = {"26",3.0} will separate the characteristic X-rays of the element Fe into two energies ranges and assign them each a column in the G matrix. This is useful to circumvent issues with the absorption. @@ -239,13 +280,21 @@ def build_G( None """ self._check_metadata_G() + + if use_calibration: + assert self.model.calibrated_db_dict is not None, ( + "Auto calibration has not been performed. Please use `auto_calibrate` before running `build_G` with `use_calibration=True`." + ) + self.problem_type = problem_type self.separated_lines = elements_dict + g_pars = { "g_type": problem_type, "ignored_elements": ignored_elements, "elements": self.metadata.Sample.elements, "elements_dict": elements_dict, + "use_calibration": use_calibration, } self.model.generate_g_matr(**g_pars) @@ -598,7 +647,7 @@ def estimate_mass_thickness( f"Estimated mass-thickness : {curr_mt} g.cm^-2" ) - axis = self.axes_manager.signal_axes[0].axis + axis = self.energy_axis self._plot.signal_plot.ax.plot( axis, estimator.G_ @ estimator.W_ @ estimator.H_, "b-", label="Full model" ) @@ -744,7 +793,7 @@ def _compute_bremsstrahlung(self): def _plot_background(self, model): # The full range from self.compute_background misses both ends # The axis needs to be trimmed accordingly - axis = self.axes_manager.signal_axes[0].axis[1:-1] + axis = self.energy_axis[1:-1] self._plot.signal_plot.ax.plot(axis, model) def _generate_ranges(self, num): @@ -925,7 +974,7 @@ def convert_elts(elements=[]): if indices: component = sum([G[:, idx] @ W[idx, :] @ H for idx in indices]) spectrum_1D._plot.signal_plot.ax.plot( - self.axes_manager.signal_axes[0].axis, + self.energy_axis, component, label=f"{conv_elts_dict[elt]}", linestyle=line_styles[_ % len(line_styles)], @@ -961,34 +1010,29 @@ def convert_elts(elements=[]): elts = self.model.get_elements(False) elts_indices = self.model.NMF_simplex() - if selected_elts: - conv_elts = convert_elts(elements=elts) - conv_elts_dict = {conv_elts[i]: num for i, num in enumerate(elts_indices)} - new_elts_indices = [] - for elt in selected_elts: - if elt in conv_elts_dict: - new_elts_indices.append(conv_elts_dict[elt]) - - W = W[new_elts_indices, :] * 100 / W[new_elts_indices, :].sum(axis=0) - if fit_error: - errors = percentages[new_elts_indices, :] - errors[errors > 10000] = np.inf - else: - errors = np.zeros_like(W) + conv_elts = convert_elts(elements=elts) - return selected_elts, W, errors + if selected_elts: + conv_elts_dict = dict(zip(conv_elts, elts_indices)) + indices = [ + conv_elts_dict[elt] for elt in selected_elts if elt in conv_elts_dict + ] + returned_elts = [elt for elt in selected_elts if elt in conv_elts_dict] + W = W[indices, :] * 100 / W[indices, :].sum(axis=0) else: - conv_elts = convert_elts(elements=elts) + indices = elts_indices + returned_elts = conv_elts - W = W[elts_indices, :] * 100 # /W[indices,:].sum(axis = 0) - if fit_error: - errors = percentages[elts_indices, :] - errors[errors > 10000] = np.inf - else: - errors = np.zeros_like(W) + W = W[indices, :] * 100 - return conv_elts, W, errors + if fit_error: + errors = percentages[indices, :] + errors[errors > 10000] = np.inf + else: + errors = np.zeros_like(W) + + return returned_elts, W, errors def estimate_best_binning(self, inspect=False): r""" @@ -1638,6 +1682,387 @@ def get_full_el_list(self): els_names = [num_to_symbol(el) for el in els] return els_names + #################### + # Auto Calibration # + #################### + + def fit_single_peak(self, window, energy): + energy_axis = self.energy_axis + + mask = (energy_axis >= energy - window) & (energy_axis <= energy + window) + if not np.any(mask): + return None + + xdata = energy_axis[mask] + ydata = self.average_spectrum[mask] + + x0_guess = xdata[np.argmax(ydata)] + A_guess = max(np.max(ydata) - np.min(ydata), 1e-6) + sigma_guess = ( + self.model.width_slope * energy + self.model.width_intercept + ) / 2.3548 + C_guess = np.min(ydata) + + p0 = (A_guess, x0_guess, sigma_guess, C_guess) + bounds = ( + (0.0, energy - window, 0.9 * sigma_guess, 0.0), + (np.inf, energy + window, 1.1 * sigma_guess, np.inf), + ) + + try: + popt, _ = curve_fit( + gaussian_with_bg, + xdata, + ydata, + p0, + # sigma=np.sqrt(np.maximum(ydata, 1.0)), + # absolute_sigma=True, + bounds=bounds, + ) + except Exception: + return None + + return popt + + def fit_table(self, window, filter_cs, filter_conc): + table = self.model.db_dict + + calibrated_table = {} + + concentrations, _ = quant_spectrum( + self.mean(axis=nav_axes) + if len(nav_axes := self.axes_manager.navigation_axes) > 0 + else self + ) + + @symbol_to_number_dict + def symbol_to_number(elements_dict): + return {str(k): v for k, v in elements_dict.items()} + + concentrations = symbol_to_number(elements_dict=concentrations) + + for element in self.elements: + if ( + element not in concentrations + or (conc := concentrations[element]) < filter_conc + ): + continue + + lines = table[element] + + max_cs = max([line["cs"] for line in lines.values()]) + + calibrated_table[element] = {} + + for line_name, line_data in lines.items(): + energy = line_data["energy"] + cs = line_data["cs"] + + if cs < filter_cs * max_cs: + continue + + fit = self.fit_single_peak(window, energy) + if fit is None: + continue + + A, x0, sigma, C = fit + + calibrated_table[element][line_name] = { + "energy": x0, + "cs": cs, + "theoretical": energy, + "sigma": sigma, + "amplitude": A, + "bg": C, + "conc": conc, + "relative_cs": cs / max_cs, + } + + return calibrated_table + + def auto_calibrate( + self, + window, + degree=2, + filter_cs=0.0, + filter_conc=0.0, + weighted=lambda conc, cs: conc * cs, + ): + r""" + Automatically calibrate the dataset energy scale linearly, and further (non-)linearly correct the table. + + This method performs the following steps: + 1. Fits Gaussian peaks to available X-ray emission lines on the average spectrum of the dataset. + 2. Applies an initial affine (degree 1) linear correction directly to the dataset's energy axis. + 3. Fits a higher-degree polynomial correction mapping theoretical peak energies to empirical line positions. + 4. Updates `self.model` for downstream matrix generation. + + `build_G` by default only uses the linear correction done in step 2. In order to benefit from further corrections, + `use_calibration=True` as to be passed to `build_G` + + Parameters + ---------- + window : float + Half-width of the energy window (in keV) around each theoretical peak position used for Gaussian fitting. + This is the only user input necessary. + degree : int, optional + Degree of the polynomial fit for non-linear energy calibration refinement applied to the table, not + the energy axis (default is 2). + filter_cs : float, optional + Relative cross-section filter threshold in [0.0, 1.0]. Lines with cross-sections smaller than + `filter_cs * max_cs` for a given element are excluded from calibration (default is 0.0, i.e. no filtering). + filter_conc : float, optional + Elemental concentration filter threshold. Elements with estimated concentration below + `filter_conc` by percentage are excluded from calibration (default is 0.0, i.e. no filtering). + weighted : callable or None, optional + Weighting function `f(conc, cs)` taking elemental concentration and relative cross-section + to compute line weights for polynomial fitting (default is `lambda conc, cs: conc * cs`). + If set to `None` or `False`, unweighted fitting is performed. + + Returns + ------- + calibrated_db : dict + Dictionary mapping element symbols to their calibrated emission lines and fitted energy parameters. + energy_poly : np.ndarray + 1D array of polynomial coefficients mapping theoretical energies to calibrated line positions. + sigma_poly : np.ndarray + 1D array of polynomial coefficients mapping theoretical energies to calibrated width of bell curves. + """ + + calibrated = self.fit_table(window, filter_cs, filter_conc) + + theoretical = [] + empirical = [] + sigma = [] + weight = [] + + for lines in calibrated.values(): + for line in lines.values(): + theoretical.append(line["theoretical"]) + empirical.append(line["energy"]) + sigma.append(line["sigma"]) + if weighted: + weight.append(weighted(line["conc"], line["relative_cs"])) + + theoretical = np.asarray(theoretical) + empirical = np.asarray(empirical) + sigma = np.asarray(sigma) + + w = weight if weighted else None + + affine = np.polyfit(empirical, theoretical, 1, w=w) + + a, b = affine + + axis = self.axes_manager[-1] + axis.scale *= a + axis.offset = a * axis.offset + b + self.energy_axis_ = None + self.average_spectrum_ = None + self.model_ = None + + empirical = np.polyval(affine, empirical) + + energy_poly = np.polyfit(theoretical, empirical, degree, w=w) + sigma_poly = np.polyfit(theoretical, sigma, degree, w=w) + + calibrated_db = { + element: { + line: { + **( + calibrated[element][line] + if element in calibrated and line in calibrated[element] + else {} + ), + "energy": np.polyval(energy_poly, data["energy"]), + "cs": data["cs"], + # "sigma": np.polyval(sigma_poly, data["energy"]), + "theoretical": data["energy"], + } + for line, data in v.items() + } + for element, v in self.model.db_dict.items() + if element in self.elements + } + + self.model.calibrated_db_dict = calibrated_db + self.model.energy_calibration_poly = energy_poly + self.model.sigma_calibration_poly = sigma_poly + + return (calibrated_db, energy_poly, sigma_poly) + + def plot_table( + self, + table=None, + elements=None, + linestyle="--", + bell=False, + legend=True, + ax=None, + ): + """Plot X-ray emission line energies and optional line profiles for elements. + + This method plots vertical lines corresponding to theoretical or calibrated X-ray + emission line energies for a specified list of elements. Optionally, fitted Gaussian + bell curves can be displayed for each line if Gaussian parameters (amplitude, sigma, background) + are present in the provided emission database `table`. Interactive hover tooltips display the + element symbol, line name, and energy value upon mouse movement over the plotted lines. + + Parameters + ---------- + table : dict or None, optional + Dictionary mapping atomic numbers (as strings, e.g. "29") to emission line data dictionaries + (containing `"energy"`, and optionally `"amplitude"`, `"sigma"`, `"bg"`). + If `None`, defaults to `self.model.db_dict`. + elements : list of str or None, optional + List of element symbols (e.g. `["Cu", "Fe"]`) for which to plot lines. + If `None`, defaults to `self.metadata.Sample.elements`. + linestyle : str, optional + Line style for the vertical emission line markers (default is `"--"`). + bell : bool, optional + If `True`, plots Gaussian bell curves around emission line energies using Gaussian line + shape parameters stored in `table` (default is `False`). + legend : bool, optional + If `True`, displays a legend mapping element symbols to line colors (default is `True`). + ax : matplotlib.axes.Axes or None, optional + Matplotlib Axes instance on which to render the plot. If `None`, a new figure and axes are + created, the average spectrum of the dataset is plotted in black, and axis labels are set + (default is `None`). + + Returns + ------- + ax : matplotlib.axes.Axes + The Matplotlib Axes object containing the plotted X-ray emission lines and spectrum. + """ + ax_was_none = ax is None + + if table is None: + table = self.model.db_dict + if elements is None: + elements = self.metadata.Sample.elements + if ax_was_none: + fig, ax = plt.subplots() + + xs = [] + colours = [] + + legend_entries = {} + + bell_segments = [] + bell_colours = [] + + 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"] + + for i, elt in enumerate(elements): + if (anum := str(SPT[elt]["number"])) not in table: + continue + lines = table[anum] + + colour = cmap(i) + legend_entries[elt] = colour + + for name, data in lines.items(): + energy = data["energy"] + + xs.append(energy) + colours.append(colour) + + hover_data.append( + {"x": energy, "name": f"{elt} {name}", "color": colour} + ) + + if bell and "amplitude" in data: + sigma = data["sigma"] + A = data["amplitude"] + C = data["bg"] + # m = data["background_slope"] + + mask = (self.energy_axis > energy - 3 * sigma) & ( + self.energy_axis < energy + 3 * sigma + ) + x_axis = self.energy_axis[mask] + + y_gauss = gaussian_with_bg(x_axis, A, energy, sigma, C) + + points = np.column_stack([x_axis, y_gauss]) + bell_segments.append(points) + bell_colours.append(colour) + + vline_collection = ax.vlines( + x=xs, + ymin=0, + ymax=1, + colors=colours, + linestyles=linestyle, + linewidths=1, + pickradius=5, + transform=ax.get_xaxis_transform(), + ) + + if bell: + bell_collection = LineCollection( + bell_segments, colors=bell_colours, linewidths=1.5 + ) + ax.add_collection(bell_collection) + + if legend: + for elt, col in legend_entries.items(): + ax.plot([], [], color=col, label=elt, lw=1) + ax.legend() + + if ax_was_none: + ax.plot( + self.energy_axis, + self.average_spectrum, + linewidth=2, + color="k", + label="Dataset", + ) + + annot = ax.annotate( + "", + xy=(0, 0), + xytext=(10, 10), + textcoords="offset points", + bbox={"boxstyle": "round", "fc": "w", "alpha": 0.9, "ec": "gray"}, + ) + annot.set_visible(False) + + def on_hover(event): + if event.inaxes == ax: + contained, info = vline_collection.contains(event) + + if contained: + line_idx = info["ind"][0] + item = hover_data[line_idx] + + annot.xy = (item["x"], event.ydata) + annot.set_text(f"{item['name']}: {item['x']:.3f}") + annot.get_bbox_patch().set_edgecolor(item["color"]) + + if not annot.get_visible(): + annot.set_visible(True) + fig.canvas.draw_idle() + return + + if annot.get_visible(): + annot.set_visible(False) + fig.canvas.draw_idle() + + fig = ax.figure + fig.canvas.mpl_connect("motion_notify_event", on_hover) + + if ax_was_none: + ax.set_xlabel("Energy (keV)") + ax.set_ylabel("Intensity") + + return ax + ####################### # Auxiliary functions # @@ -1697,3 +2122,7 @@ def build_G(model, g_params): def Gauss(x, a, x0, sigma): return a * np.exp(-((x - x0) ** 2) / (2 * sigma**2)) + + +def gaussian_with_bg(x, A, x0, sigma, C): + return A * np.exp(-((x - x0) ** 2) / (2 * sigma**2)) + C # + m * (x - x0) diff --git a/espm/models/edxs.py b/espm/models/edxs.py index 4d09287..4fc8a2d 100644 --- a/espm/models/edxs.py +++ b/espm/models/edxs.py @@ -61,81 +61,63 @@ def __init__( # Tranfer the ranges from eds_espm to the physical model self.ranges = None - def __add_elts_G(self, reference_elt={}, *, elements=[]): + self.calibrated_db_dict = None + self.energy_calibration_poly = None + self.sigma_calibration_poly = None + + def __add_elts_G(self, use_calibration, reference_elt={}, *, elements=[]): for elt in elements: - if self.lines: - energies, cs = read_lines_db(elt, self.db_dict) - else: - energies, cs = read_compact_db(elt, self.db_dict) - if elt in reference_elt: - peaks_low = np.zeros((self.x.shape[0], 1)) - peaks_high = np.zeros((self.x.shape[0], 1)) - for i, energy in enumerate(energies): - if (energy > np.min(self.x)) and (energy < np.max(self.x)): - if type(self.params_dict["Det"]) == str: - D = det_efficiency_from_curve( - energy, self.params_dict["Det"] - ) - else: - D = det_efficiency(energy, self.params_dict["Det"]) - - A = absorption_correction( - energy, **self.params_dict["Abs"], elements_dict={elt: 1.0} - ) + energies, cs = ( + read_lines_db(elt, self.calibrated_db_dict) + if use_calibration + else read_lines_db(elt, self.db_dict) + if self.lines + else read_compact_db(elt, self.db_dict) + ) + lines = [ + ( + energy, + np.polyval(self.sigma_calibration_poly, energy) + if use_calibration + else (self.width_slope * energy + self.width_intercept) / 2.3548, + c, + ) + for energy, c in zip(energies, cs) + ] - width = self.width_slope * energy + self.width_intercept - if energy < reference_elt[elt]: - peaks_low += ( - ( - cs[i] - * gaussian(self.x, energy, width / 2.3548)[ - np.newaxis - ].T - ) - * D - * A - ) - else: - peaks_high += ( - ( - cs[i] - * gaussian(self.x, energy, width / 2.3548)[ - np.newaxis - ].T - ) - * D - * A - ) - peaks = np.hstack((peaks_low, peaks_high)) + shape = (self.x.shape[0], 1) + if elt in reference_elt: + peaks_low = np.zeros(shape) + peaks_high = np.zeros(shape) else: - peaks = np.zeros((self.x.shape[0], 1)) - for i, energy in enumerate(energies): - # The actual detected width is calculated at each energy - if (energy > np.min(self.x)) and (energy < np.max(self.x)): - if type(self.params_dict["Det"]) == str: - D = det_efficiency_from_curve( - energy, self.params_dict["Det"] - ) - else: - D = det_efficiency(energy, self.params_dict["Det"]) - - A = absorption_correction( - energy, **self.params_dict["Abs"], elements_dict={elt: 1.0} - ) + peaks = np.zeros(shape) - width = self.width_slope * energy + self.width_intercept + for energy, sigma, cs in lines: + # The actual detected width is calculated at each energy + if not np.min(self.x) < energy < np.max(self.x): + continue - peaks += ( - ( - cs[i] - * gaussian(self.x, energy, width / 2.3548)[np.newaxis].T - ) - * D - * A - ) + if type(self.params_dict["Det"]) == str: + D = det_efficiency_from_curve(energy, self.params_dict["Det"]) + else: + D = det_efficiency(energy, self.params_dict["Det"]) + + A = absorption_correction( + energy, **self.params_dict["Abs"], elements_dict={elt: 1.0} + ) + + delta = cs * gaussian(self.x, energy, sigma)[np.newaxis].T * D * A + if elt in reference_elt: + if energy < reference_elt[elt]: + peaks_low += delta + else: + peaks_high += delta + else: + peaks += delta + + if elt in reference_elt: + peaks = np.hstack((peaks_low, peaks_high)) - # print(np.max(peaks, axis = 0)) - # print(str(elt)) if np.all((np.max(peaks, axis=0)) > 0.0): self.G = np.concatenate((self.G, peaks), axis=1) if elt in reference_elt: @@ -149,32 +131,43 @@ def __add_elts_G(self, reference_elt={}, *, elements=[]): ) raise ValueError("Empty G column") - def _add_ignored_elts(self, elements=[]): + def _add_ignored_elts(self, use_calibration, elements=[]): for elt in elements: - if self.lines: - energies, cs = read_lines_db(elt, self.db_dict) - else: - energies, cs = read_compact_db(elt, self.db_dict) + energies, cs = ( + read_lines_db(elt, self.calibrated_db_dict) + if use_calibration + else read_lines_db(elt, self.db_dict) + if self.lines + else read_compact_db(elt, self.db_dict) + ) + lines = [ + ( + energy, + np.polyval(self.sigma_calibration_poly, energy) + if use_calibration + else (self.width_slope * energy + self.width_intercept) / 2.3548, + c, + ) + for energy, c in zip(energies, cs) + ] peaks_list = [] - for i, energy in enumerate(energies): - # The actual detected width is calculated at each energy - if (energy > np.min(self.x)) and (energy < np.max(self.x)): - if type(self.params_dict["Det"]) == str: - D = det_efficiency_from_curve(energy, self.params_dict["Det"]) - else: - D = det_efficiency(energy, self.params_dict["Det"]) + for energy, sigma, cs in lines: + if not np.min(self.x) < energy < np.max(self.x): + continue + if type(self.params_dict["Det"]) == str: + D = det_efficiency_from_curve(energy, self.params_dict["Det"]) + else: + D = det_efficiency(energy, self.params_dict["Det"]) - A = absorption_correction( - energy, **self.params_dict["Abs"], elements_dict={elt: 1.0} - ) + A = absorption_correction( + energy, **self.params_dict["Abs"], elements_dict={elt: 1.0} + ) - width = self.width_slope * energy + self.width_intercept + peaks_list.append((cs * gaussian(self.x, energy, sigma)) * D * A) - peaks_list.append( - (cs[i] * gaussian(self.x, energy, width / 2.3548)) * D * A - ) peaks = np.array(peaks_list).T + if np.all((np.max(peaks, axis=0)) > 0.0): self.G = np.concatenate((self.G, peaks), axis=1) for i in range(peaks.shape[1]): @@ -194,6 +187,7 @@ def generate_g_matr( *, elements=[], elements_dict={}, + use_calibration=False, **kwargs, ): r""" @@ -245,24 +239,27 @@ def convert_elts(elements=ignored_elements): conv_ignored_elts = convert_elts(elements=ignored_elements) - valid_elts = self.__check_elts_in_G(elements) - valid_ignored = self.__check_elts_in_G(conv_ignored_elts) + valid_elts = self.__check_elts_in_G(elements, use_calibration) + valid_ignored = self.__check_elts_in_G(conv_ignored_elts, use_calibration) - if g_type == "bremsstrahlung": - self.bkgd_in_G = True - else: - self.bkgd_in_G = False + self.bkgd_in_G = g_type == "bremsstrahlung" # None is a default value for the G matrix and thus G will be considered to be the identity matrix in most of espm functions. if len(valid_elts) == 0 or g_type == "identity": self.G = None # model based on elements_list - elif (g_type == "bremsstrahlung") or (g_type == "no_brstlg"): + elif self.bkgd_in_G or (g_type == "no_brstlg"): # The number of shells depend on the element, it is then not straightforward to pre-determine the size of g_matr self.G = np.zeros((self.x.shape[0], 0)) # For each element we unpack all shells and then unpack all lines of each shell. - self.__add_elts_G(reference_elt=elements_dict, elements=valid_elts) - self._add_ignored_elts(elements=valid_ignored) + self.__add_elts_G( + use_calibration=use_calibration, + reference_elt=elements_dict, + elements=valid_elts, + ) + self._add_ignored_elts( + use_calibration=use_calibration, elements=valid_ignored + ) # Appends a pure continuum spectrum is needed if self.bkgd_in_G: @@ -279,10 +276,12 @@ def convert_elts(elements=ignored_elements): self.bkgd_in_G = False norms = np.sum(self.G, axis=0, keepdims=True) - if g_type == "bremsstrahlung": + + if self.bkgd_in_G: norms[0][:-2] = np.mean(norms[0][:-2]) else: norms[0] = np.mean(norms[0]) + self.norm = norms self.G /= self.norm else: @@ -290,23 +289,26 @@ def convert_elts(elements=ignored_elements): 'g_type has to be one of those : "bremsstrahlung", "no_brstlg" or "identity". G will be None, corresponding to "identity". ' ) - def __check_elts_in_G(self, elements): + def __check_elts_in_G(self, elements, use_calibration): """ Check if the elements of the metadata are in the range of the energy axis. """ valid_elts = [] for elt in elements: - if self.lines: - energies, cs = read_lines_db(elt, self.db_dict) - else: - energies, cs = read_compact_db(elt, self.db_dict) + energies, _ = ( + read_lines_db(elt, self.calibrated_db_dict) + if use_calibration + else read_lines_db(elt, self.db_dict) + if self.lines + else read_compact_db(elt, self.db_dict) + ) energy_range = [np.min(self.x), np.max(self.x)] - found = 0 + found = False for energy in energies: - if (energy > energy_range[0]) and (energy < energy_range[1]): + if energy_range[0] < energy < energy_range[1]: valid_elts.append(elt) - found = 1 + found = True break if not found: print(f"No peak is present in the energy range for element : {elt}") diff --git a/espm/utils.py b/espm/utils.py index e8eea0b..4682b5d 100644 --- a/espm/utils.py +++ b/espm/utils.py @@ -1,5 +1,7 @@ r"""Utils for the ESPM package""" +import contextlib +import io import json from functools import wraps @@ -10,7 +12,6 @@ import seaborn import skimage as ski from exspy.material import atomic_to_weight, density_of_mixture -from IPython.utils import io from scipy.optimize import nnls from scipy.sparse import block_diag, lil_matrix from sklearn.linear_model import LinearRegression as LR @@ -525,19 +526,13 @@ def quant_spectrum(s1, skip_elements=[]): ] s.build_G() - est = espm.estimators.SmoothNMF(n_components=1, G=s.G(), verbose=0) - with io.capture_output() as captured: + est = espm.estimators.SmoothNMF(n_components=1, G=s.G, verbose=0) + with contextlib.redirect_stdout(io.StringIO()): est.fit_transform(X=s1.data[:, np.newaxis], H=np.array([1.0])[:, np.newaxis]) s.learning_results.decomposition_algorithm = est - with io.capture_output() as captured: - s.print_concentration_report(selected_elts=selected_elements) - # print(captured) - return dict( - [ - [i.split(":")[0][:-1], float(i.split(":")[1])] - for i in captured.stdout.splitlines()[2:] - ] - ), s + conv_elts, W, _ = s.concentration_report(selected_elts=selected_elements) + quant_dict = {el: float(W[i, 0]) for i, el in enumerate(conv_elts)} + return quant_dict, s def cluster_analysis_concentration_report(s, cluster_source=None, print_std=False): diff --git a/notebooks/calibration.ipynb b/notebooks/calibration.ipynb new file mode 100644 index 0000000..c4cc8da --- /dev/null +++ b/notebooks/calibration.ipynb @@ -0,0 +1,943 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "5f404f75", + "metadata": {}, + "source": [ + "### TL;DR\n", + "\n", + "In order to use auto calibration, do the following:\n", + "\n", + "1. Call `EDSespm::plot_table` to see how shifted the dataset is.\n", + "2. Call `EDSespm::auto_calibrate` with `window` you obtained in step 1.\n", + " This does the following:\n", + " - Afflinely correct the energy axis of the dataset\n", + " - Correct the line table with higher degree polynomial\n", + "3. Call `EDSespm::build_G` with `use_calibration=True` if you want to use the higher degree correction." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f937247d", + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib qt\n", + "# %matplotlib widget\n", + "\n", + "import json\n", + "\n", + "import hyperspy.api as hs\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "from matplotlib.collections import LineCollection\n", + "\n", + "from espm.conf import SYMBOLS_PERIODIC_TABLE\n", + "from espm.datasets.eds_spim import gaussian_with_bg" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "90f42db7", + "metadata": {}, + "outputs": [], + "source": [ + "FILENAME = \"../playground/X3-13MAY22_MAP06.bcf\" # Change this to your dataset\n", + "BIN = 32 # Scale for rebinning\n", + "\n", + "WINDOW = 0.05\n", + "FILTER_CS = 0.1 # Keep lines with cs >= FILTER_CS * CS_max\n", + "FILTER_CONC = 1.0 # Keep elements with concentration >= FILTER_CONC %\n", + "DEGREE = 2 # The degree of polynomial used to calibrate the energy axis\n", + "\n", + "MAX_DEGREE = 10" + ] + }, + { + "cell_type": "markdown", + "id": "8be084a1", + "metadata": {}, + "source": [ + "Matplotlib setup, please ignore." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "755008d2", + "metadata": {}, + "outputs": [], + "source": [ + "with open(SYMBOLS_PERIODIC_TABLE, \"r\") as f:\n", + " SPT = json.load(f)[\"table\"]\n", + "\n", + "cmap = plt.get_cmap(\"tab10\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "677097f0", + "metadata": {}, + "outputs": [], + "source": [ + "def plot_table(ax, table, elements, energy_axis, l1=\"--\", l2=\"-.\", bell=False):\n", + " xs = []\n", + " colours = []\n", + " linestyles = []\n", + "\n", + " legend = {}\n", + "\n", + " bell_segments = []\n", + " bell_colours = []\n", + "\n", + " hover_data = []\n", + "\n", + " idx = set()\n", + "\n", + " for i, elt in enumerate(elements):\n", + " elt_num_str = str(SPT[elt][\"number\"])\n", + " if elt_num_str not in table:\n", + " continue\n", + " db_entries = table[elt_num_str]\n", + "\n", + " colour = cmap(i)\n", + " legend[elt] = colour\n", + "\n", + " for name, data in db_entries.items():\n", + " energy = data[\"energy\"]\n", + "\n", + " xs.append(energy)\n", + " colours.append(colour)\n", + " linestyles.append(l1 if energy not in idx else l2)\n", + "\n", + " hover_data.append({\"x\": energy, \"name\": f\"{elt} {name}\", \"color\": colour})\n", + "\n", + " if bell and \"amplitude\" in data:\n", + " sigma = data[\"sigma\"]\n", + " A = data[\"amplitude\"]\n", + " C = data[\"bg\"]\n", + " # m = data[\"background_slope\"]\n", + "\n", + " mask = (energy_axis > energy - 3 * sigma) & (\n", + " energy_axis < energy + 3 * sigma\n", + " )\n", + " x_axis = energy_axis[mask]\n", + "\n", + " y_gauss = gaussian_with_bg(x_axis, A, energy, sigma, C)\n", + "\n", + " points = np.column_stack([x_axis, y_gauss])\n", + " bell_segments.append(points)\n", + " bell_colours.append(colour)\n", + "\n", + " idx.add(energy)\n", + "\n", + " vline_collection = ax.vlines(\n", + " x=xs,\n", + " ymin=-0.1,\n", + " ymax=2 * BIN * BIN,\n", + " colors=colours,\n", + " linestyles=linestyles,\n", + " linewidths=1,\n", + " pickradius=5,\n", + " )\n", + "\n", + " if bell:\n", + " bell_collection = LineCollection(\n", + " bell_segments, colors=bell_colours, linewidths=1.5\n", + " )\n", + " ax.add_collection(bell_collection)\n", + "\n", + " for elt, col in legend.items():\n", + " ax.plot([], [], color=col, label=elt, lw=1)\n", + " ax.legend()\n", + "\n", + " annot = ax.annotate(\n", + " \"\",\n", + " xy=(0, 0),\n", + " xytext=(10, 10),\n", + " textcoords=\"offset points\",\n", + " bbox={\"boxstyle\": \"round\", \"fc\": \"w\", \"alpha\": 0.9, \"ec\": \"gray\"},\n", + " )\n", + " annot.set_visible(False)\n", + "\n", + " def on_hover(event):\n", + " if event.inaxes == ax:\n", + " contained, info = vline_collection.contains(event)\n", + "\n", + " if contained:\n", + " line_idx = info[\"ind\"][0]\n", + " item = hover_data[line_idx]\n", + "\n", + " annot.xy = (item[\"x\"], event.ydata)\n", + " annot.set_text(f\"{item['name']}: {item['x']:.3f}\")\n", + " annot.get_bbox_patch().set_edgecolor(item[\"color\"])\n", + "\n", + " if not annot.get_visible():\n", + " annot.set_visible(True)\n", + " fig.canvas.draw_idle()\n", + " return\n", + "\n", + " if annot.get_visible():\n", + " annot.set_visible(False)\n", + " fig.canvas.draw_idle()\n", + "\n", + " fig = ax.figure\n", + " fig.canvas.mpl_connect(\"motion_notify_event\", on_hover)\n", + "\n", + "\n", + "def plot(funcs, title):\n", + " fig, ax = plt.subplots()\n", + "\n", + " for func in funcs:\n", + " func(ax)\n", + "\n", + " ax.set_title(title)\n", + " ax.set_xlabel(\"Energy (keV)\")\n", + " ax.set_ylabel(\"Intensity (counts)\")\n", + "\n", + " ax.legend()\n", + "\n", + " def on_press(event):\n", + " if event.key == \" \":\n", + " toolbar = event.canvas.toolbar\n", + " if toolbar.mode != \"pan/zoom\":\n", + " toolbar.pan()\n", + "\n", + " def on_release(event):\n", + " if event.key == \" \":\n", + " toolbar = event.canvas.toolbar\n", + " if toolbar.mode == \"pan/zoom\":\n", + " toolbar.pan()\n", + "\n", + " fig.canvas.mpl_connect(\"key_press_event\", on_press)\n", + " fig.canvas.mpl_connect(\"key_release_event\", on_release)\n", + "\n", + " plt.show()\n", + "\n", + "\n", + "def plot_avg(ax, avg_spectrum, energy_axis, **kwargs):\n", + " ax.plot(energy_axis, avg_spectrum, **kwargs)" + ] + }, + { + "cell_type": "markdown", + "id": "83fabfdb", + "metadata": {}, + "source": [ + "### Load Dataset" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0c17427c", + "metadata": {}, + "outputs": [], + "source": [ + "signals = hs.load(FILENAME)\n", + "signal = signals[1].rebin(scale=(BIN, BIN, 1)).isig[0.2:]\n", + "signal.set_signal_type(\"EDS_espm\")\n", + "\n", + "signal.set_analysis_parameters(\n", + " thickness=10e-5,\n", + " density=4.1,\n", + " detector_type=\"SDD_efficiency.txt\",\n", + " width_slope=0.01,\n", + " width_intercept=0.065,\n", + " geom_eff=None,\n", + " xray_db=\"200keV_xrays.json\",\n", + ")\n", + "signal.change_dtype(\"float64\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1fe3bc61", + "metadata": {}, + "outputs": [], + "source": [ + "average_spectrum = signal.average_spectrum\n", + "energy_axis = signal.energy_axis\n", + "theoretical_table = signal.model.db_dict\n", + "\n", + "elements = signal.metadata.Sample.elements\n", + "# elements.remove(\"Re\")\n", + "\n", + "elements" + ] + }, + { + "cell_type": "markdown", + "id": "ce0c46b2", + "metadata": {}, + "source": [ + "### Result You should Expect :)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "98c5641b", + "metadata": {}, + "outputs": [], + "source": [ + "s = signal.deepcopy()\n", + "s.auto_calibrate(WINDOW)\n", + "s_filter = signal.deepcopy()\n", + "s_filter.auto_calibrate(WINDOW, filter_cs=FILTER_CS, filter_conc=FILTER_CONC)\n", + "\n", + "plot(\n", + " [\n", + " lambda ax: plot_avg(\n", + " ax, average_spectrum, energy_axis, linewidth=3, label=\"Input\"\n", + " ),\n", + " lambda ax: plot_table(ax, theoretical_table, elements, energy_axis),\n", + " ]\n", + " + [\n", + " lambda ax: plot_avg(\n", + " ax,\n", + " s.average_spectrum,\n", + " s.energy_axis,\n", + " linewidth=3,\n", + " label=\"Calibrated\",\n", + " ),\n", + " lambda ax: plot_table(\n", + " ax,\n", + " s.model.calibrated_db_dict,\n", + " elements,\n", + " s.energy_axis,\n", + " l1=\":\",\n", + " l2=\"-\",\n", + " ),\n", + " ]\n", + " + [\n", + " lambda ax: plot_avg(\n", + " ax,\n", + " s_filter.average_spectrum,\n", + " s_filter.energy_axis,\n", + " linewidth=3,\n", + " label=\"Calibrated (Filtered)\",\n", + " ),\n", + " lambda ax: plot_table(\n", + " ax,\n", + " s_filter.model.calibrated_db_dict,\n", + " elements,\n", + " s_filter.energy_axis,\n", + " l1=\":\",\n", + " l2=\"-\",\n", + " ),\n", + " ],\n", + " \"Original vs Calibrated\",\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "9841729f", + "metadata": {}, + "source": [ + "### Average Spectrum vs Theoretical Lines" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "74bbcead", + "metadata": {}, + "outputs": [], + "source": [ + "signal.plot_table()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "2198f0f3", + "metadata": {}, + "source": [ + "### Independently and Naively Fitted Lines\n", + "\n", + "Dashed lines (- - -) are theoretical lines, and dotted lines (. . . .) are fitted lines" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "275fe830", + "metadata": {}, + "outputs": [], + "source": [ + "plot(\n", + " [\n", + " lambda ax: plot_avg(\n", + " ax, average_spectrum, energy_axis, color=\"k\", linewidth=3, label=\"Input\"\n", + " ),\n", + " lambda ax: plot_table(ax, theoretical_table, elements, energy_axis),\n", + " lambda ax: plot_table(\n", + " ax,\n", + " signal.fit_table(WINDOW, filter_cs=0.0, filter_conc=0.0),\n", + " elements,\n", + " energy_axis,\n", + " l1=\":\",\n", + " l2=\"-\",\n", + " bell=True,\n", + " ),\n", + " ],\n", + " \"Independent Naive Line Fit\",\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "738f7e58", + "metadata": {}, + "source": [ + "### Independently and Fitted Lines with Cross Section Filter\n", + "\n", + "Dashed lines (- - -) are theoretical lines, and dotted lines (. . . .) are fitted lines\n", + "\n", + "> Only lines with cross section $\\geq$ `FILTER_CS` $\\cdot$ (max cs of the element) are kept" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cc38b6dd", + "metadata": {}, + "outputs": [], + "source": [ + "plot(\n", + " [\n", + " lambda ax: plot_avg(\n", + " ax, average_spectrum, energy_axis, color=\"k\", linewidth=3, label=\"Input\"\n", + " ),\n", + " lambda ax: plot_table(ax, theoretical_table, elements, energy_axis),\n", + " lambda ax: plot_table(\n", + " ax,\n", + " signal.fit_table(WINDOW, filter_cs=FILTER_CS, filter_conc=0.0),\n", + " elements,\n", + " energy_axis,\n", + " l1=\":\",\n", + " l2=\"-\",\n", + " bell=True,\n", + " ),\n", + " ],\n", + " f\"Independent Line Fit with CS Filter {FILTER_CS}\",\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "0df5c9e4", + "metadata": {}, + "source": [ + "### Independently and Fitted Lines with Cross Section and Concentration Filter\n", + "\n", + "Dashed lines (- - -) are theoretical lines, and dotted lines (. . . .) are fitted lines\n", + "\n", + "> Only lines with cross section $\\geq$ `FILTER_CS` $\\cdot$ (max cs of the element) are kept\n", + "\n", + "> Only elements with concentration $\\geq$ `FILTER_CONC` % are kept" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "13f1e26b", + "metadata": {}, + "outputs": [], + "source": [ + "from espm.utils import quant_spectrum\n", + "\n", + "quant_spectrum(signal.mean(axis=(0, 1)))[0]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e44b6171", + "metadata": {}, + "outputs": [], + "source": [ + "plot(\n", + " [\n", + " lambda ax: plot_avg(\n", + " ax, average_spectrum, energy_axis, color=\"k\", linewidth=3, label=\"Input\"\n", + " ),\n", + " lambda ax: plot_table(ax, theoretical_table, elements, energy_axis),\n", + " lambda ax: plot_table(\n", + " ax,\n", + " signal.fit_table(window=0.05, filter_cs=FILTER_CS, filter_conc=FILTER_CONC),\n", + " elements,\n", + " energy_axis,\n", + " l1=\":\",\n", + " l2=\"-\",\n", + " bell=True,\n", + " ),\n", + " ],\n", + " f\"Independent Line Fit with CS Filter {FILTER_CS} and Concentration Filter {FILTER_CONC}\",\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "e4a1c989", + "metadata": {}, + "source": [ + "### Lines Calibrated using (Unweighted) Polynomial Estimation [Bad]\n", + "\n", + "Dashed lines (- - -) are theoretical lines, and dotted lines (. . . .) are fitted lines\n", + "\n", + "> Every line is independently fitted and a polynomial is fitted with the original and calibrated position." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "41085c0e", + "metadata": {}, + "outputs": [], + "source": [ + "s = signal.deepcopy()\n", + "s.auto_calibrate(\n", + " WINDOW,\n", + " degree=DEGREE,\n", + " weighted=False,\n", + ")\n", + "\n", + "plot(\n", + " [\n", + " lambda ax: plot_avg(\n", + " ax, average_spectrum, energy_axis, color=\"k\", linewidth=3, label=\"Input\"\n", + " ),\n", + " lambda ax: plot_avg(\n", + " ax,\n", + " s.average_spectrum,\n", + " s.energy_axis,\n", + " linewidth=3,\n", + " label=\"Calibrated\",\n", + " color=\"r\",\n", + " ),\n", + " lambda ax: plot_table(ax, theoretical_table, elements, energy_axis),\n", + " lambda ax: plot_table(\n", + " ax,\n", + " s.model.calibrated_db_dict,\n", + " elements,\n", + " s.energy_axis,\n", + " l1=\":\",\n", + " l2=\"-\",\n", + " bell=True,\n", + " ),\n", + " ],\n", + " f\"Unweighted Polynomial Fitting (degree {DEGREE})\",\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "9447e678", + "metadata": {}, + "source": [ + "### Lines Calibrated using (Weighted) Polynomial Estimation [Good]\n", + "\n", + "Dashed lines (- - -) are theoretical lines, and dotted lines (. . . .) are fitted lines\n", + "\n", + "> Every line is independently fitted and a polynomial is fitted with the original and calibrated position, with amplitude of each peak as weight." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c3ed48a4", + "metadata": {}, + "outputs": [], + "source": [ + "s = signal.deepcopy()\n", + "s.auto_calibrate(\n", + " WINDOW,\n", + " degree=DEGREE,\n", + ")\n", + "\n", + "plot(\n", + " [\n", + " lambda ax: plot_avg(\n", + " ax, average_spectrum, energy_axis, color=\"k\", linewidth=3, label=\"Input\"\n", + " ),\n", + " lambda ax: plot_avg(\n", + " ax,\n", + " s.average_spectrum,\n", + " s.energy_axis,\n", + " linewidth=3,\n", + " label=\"Calibrated\",\n", + " color=\"r\",\n", + " ),\n", + " lambda ax: plot_table(ax, theoretical_table, elements, energy_axis),\n", + " lambda ax: plot_table(\n", + " ax,\n", + " s.model.calibrated_db_dict,\n", + " elements,\n", + " s.energy_axis,\n", + " l1=\":\",\n", + " l2=\"-\",\n", + " bell=True,\n", + " ),\n", + " ],\n", + " f\"Weighted Polynomial Fitting (degree {DEGREE})\",\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "b3647594", + "metadata": {}, + "source": [ + "### Comparing Polynomial of Different Degrees\n", + "\n", + "We see that most low degree polynomial agree that the calibration should be approximately linear, and degree 2 polynomial behaves the best. (For a polynomial with too high degree, overfitting will happen, which is not good)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6370eec5", + "metadata": {}, + "outputs": [], + "source": [ + "def plot_polys(filter_cs, filter_conc):\n", + " _, ax = plt.subplots()\n", + "\n", + " lines = []\n", + "\n", + " for d in range(1, MAX_DEGREE + 1):\n", + " s = signal.deepcopy()\n", + " _, y_poly, _ = s.auto_calibrate(\n", + " degree=d, window=0.05, filter_cs=filter_cs, filter_conc=filter_conc\n", + " )\n", + " ax.plot(energy_axis, np.polyval(y_poly, energy_axis), label=f\"degree {d}\")\n", + "\n", + " ax.plot(\n", + " energy_axis, energy_axis, label=\"identity\", color=\"k\", linewidth=3, zorder=50\n", + " )\n", + "\n", + " x = []\n", + " y = []\n", + " c = []\n", + "\n", + " calibrated = signal.fit_table(\n", + " window=0.05, filter_cs=filter_cs, filter_conc=filter_conc\n", + " )\n", + "\n", + " for i, el in enumerate(calibrated):\n", + " lines = calibrated[el]\n", + " for _, line in lines.items():\n", + " x.append(line[\"theoretical\"])\n", + " y.append(line[\"energy\"])\n", + " c.append(cmap(i))\n", + "\n", + " ax.scatter(x, y, c=c, zorder=100)\n", + "\n", + " ax.set_aspect(\"equal\")\n", + " ax.set_box_aspect(1)\n", + " ax.set_xlim(0, 20)\n", + " ax.set_ylim(0, 20)\n", + " ax.legend(loc=\"upper right\")\n", + "\n", + " ax.set_title(f\"Theoretical vs Calibrated Energy ({filter_cs}, {filter_conc})\")\n", + " ax.set_xlabel(\"Theoretical Energy (keV)\")\n", + " ax.set_ylabel(\"Calibrated Energy (keV)\")\n", + "\n", + " plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "8a2a7182", + "metadata": {}, + "source": [ + "#### Polynomial Fittings Without Filtering by Cross Section" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f98a1c53", + "metadata": {}, + "outputs": [], + "source": [ + "plot_polys(0.0, 0.0)" + ] + }, + { + "cell_type": "markdown", + "id": "b514f249", + "metadata": {}, + "source": [ + "#### Polynomial Fittings With Filtering by Cross Section" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3bef849a", + "metadata": {}, + "outputs": [], + "source": [ + "plot_polys(FILTER_CS, FILTER_CONC)" + ] + }, + { + "cell_type": "markdown", + "id": "deccd2ae", + "metadata": {}, + "source": [ + "### Visualising the `G` Built by Different Calibration Methods\n", + "\n", + "At the end, we will plot a graph with the MSE of each method. We will see that weighted poly fit with degree 2 is generally very good, and naive method is overfitted." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "22e842a4", + "metadata": {}, + "outputs": [], + "source": [ + "from espm.estimators.smooth_nmf import SmoothNMF\n", + "\n", + "\n", + "def run_decomp(name, build_g_func):\n", + " s = signal.deepcopy()\n", + " build_g_func(s)\n", + "\n", + " estimator = SmoothNMF(\n", + " n_components=3,\n", + " G=s.G,\n", + " max_iter=500,\n", + " tol=1e-5,\n", + " init=\"nndsvdar\",\n", + " random_state=42,\n", + " hspy_comp=True,\n", + " )\n", + "\n", + " s.decomposition(algorithm=estimator)\n", + "\n", + " G_est = estimator.G_\n", + " W_est = estimator.W_\n", + " H_est = estimator.H_\n", + "\n", + " reconstructed = G_est @ W_est @ H_est\n", + " reconstructed_mean = reconstructed.mean(1)\n", + "\n", + " mse = np.mean((reconstructed_mean - s.average_spectrum) ** 2)\n", + " mae = np.mean(np.abs(reconstructed_mean - s.average_spectrum))\n", + "\n", + " print(f\"{name} MSE: {mse:.6f}, MAE: {mae:.6f}\")\n", + " return reconstructed_mean, mse, mae, G_est" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d2dbc164", + "metadata": {}, + "outputs": [], + "source": [ + "rec_uncal, mse_uncal, mae_uncal, G_uncal = run_decomp(\n", + " \"Uncalibrated\", lambda s: s.build_G()\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b8d87e8d", + "metadata": {}, + "outputs": [], + "source": [ + "def build_calibrated_naive(s):\n", + " s.model.db_dict = s.fit_table(window=WINDOW, filter_cs=0.0, filter_conc=0.0)\n", + " s.build_G()\n", + "\n", + "\n", + "rec_naive, mse_naive, mae_naive, G_naive = run_decomp(\n", + " \"peak fit\", build_calibrated_naive\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ef9b79f0", + "metadata": {}, + "outputs": [], + "source": [ + "rec_poly_unweighted = {}\n", + "mse_poly_unweighted = {}\n", + "mae_poly_unweighted = {}\n", + "G_poly_unweighted = {}\n", + "\n", + "\n", + "def build_poly_unweighted(s, d):\n", + " s.auto_calibrate(\n", + " WINDOW,\n", + " degree=d,\n", + " weighted=False,\n", + " )\n", + " s.build_G(use_calibration=True)\n", + "\n", + "\n", + "for d in range(1, MAX_DEGREE + 1):\n", + " rec, mse, mae, G = run_decomp(\n", + " f\"Unweighted Poly Fit degree {d}\", lambda s, d=d: build_poly_unweighted(s, d)\n", + " )\n", + " rec_poly_unweighted[d] = rec\n", + " mse_poly_unweighted[d] = mse\n", + " mae_poly_unweighted[d] = mae\n", + " G_poly_unweighted[d] = G" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3c178aa2", + "metadata": {}, + "outputs": [], + "source": [ + "rec_poly, mse_poly, mae_poly, G_poly = {}, {}, {}, {}\n", + "\n", + "\n", + "def build_poly(s, d):\n", + " s.auto_calibrate(WINDOW, degree=d)\n", + " s.build_G(use_calibration=True)\n", + "\n", + "\n", + "for d in range(1, MAX_DEGREE + 1):\n", + " rec, mse, mae, G = run_decomp(\n", + " f\"Weighted Poly Fit degree {d}\", lambda s, d=d: build_poly(s, d)\n", + " )\n", + " rec_poly[d] = rec\n", + " mse_poly[d] = mse\n", + " mae_poly[d] = mae\n", + " G_poly[d] = G" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "59c2a9a4", + "metadata": {}, + "outputs": [], + "source": [ + "s = signal.deepcopy()\n", + "s.auto_calibrate(WINDOW)\n", + "s_unweighted = signal.deepcopy()\n", + "s_unweighted.auto_calibrate(WINDOW, weighted=False)\n", + "\n", + "plot(\n", + " [\n", + " lambda ax: plot_avg(\n", + " ax, average_spectrum, energy_axis, linewidth=3, label=\"Input\"\n", + " ),\n", + " lambda ax: plot_avg(\n", + " ax,\n", + " s.average_spectrum,\n", + " s.energy_axis,\n", + " linewidth=3,\n", + " label=\"Calibrated\",\n", + " ),\n", + " lambda ax: plot_avg(\n", + " ax,\n", + " rec_uncal,\n", + " energy_axis,\n", + " linewidth=3,\n", + " label=f\"Uncalibrated {mse_uncal:.3f}\",\n", + " linestyle=\"-.\",\n", + " ),\n", + " lambda ax: plot_avg(\n", + " ax,\n", + " rec_naive,\n", + " energy_axis,\n", + " linewidth=3,\n", + " label=f\"Naive {mse_naive:.3f}\",\n", + " linestyle=\"-.\",\n", + " ),\n", + " ]\n", + " + [\n", + " lambda ax, d=d: plot_avg(\n", + " ax,\n", + " rec_poly[d],\n", + " s.energy_axis,\n", + " label=f\"Weighted Poly deg {d} {mse_poly[d]:.3f}\",\n", + " linewidth=3,\n", + " )\n", + " for d in range(2, 2 + 1)\n", + " ]\n", + " # + [\n", + " # lambda ax, d=d: plot_avg(\n", + " # ax,\n", + " # rec_poly_unweighted[d],\n", + " # s.energy_axis,\n", + " # label=f\"Unweighted Poly deg {d} {mse_poly_unweighted[d]:.3f}\",\n", + " # linestyle=\"--\",\n", + " # )\n", + " # for d in range(1, MAX_DEGREE + 1)\n", + " # ]\n", + " + [\n", + " lambda ax: plot_table(\n", + " ax, theoretical_table, elements, energy_axis, l1=\"-.\", l2=\"-.\"\n", + " )\n", + " ]\n", + " + [\n", + " lambda ax: plot_table(\n", + " ax,\n", + " s.model.calibrated_db_dict,\n", + " elements,\n", + " s.energy_axis,\n", + " l1=\":\",\n", + " l2=\"-\",\n", + " # bell=True,\n", + " ),\n", + " ],\n", + " # + [\n", + " # lambda ax: plot_table(\n", + " # ax,\n", + " # s_unweighted.model.calibrated_db_dict,\n", + " # elements,\n", + " # s_unweighted.energy_axis,\n", + " # l1=\"--\",\n", + " # l2=\"-.\",\n", + " # # bell=True,\n", + " # ),\n", + " # ],\n", + " \"G of Different Methods\",\n", + ")" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "espm (3.12.x)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.13" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/pyproject.toml b/pyproject.toml index 3ddfc86..0ea09cd 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -46,7 +46,6 @@ dependencies = [ "exspy>=0.3.2", "hyperspy>=2.4.0", "intervaltree>=3.2.1", - "ipython>=8.39.0", "numpy>=2.2.6", "scikit-image>=0.25.2", "scikit-learn>=1.7.2",