Adds align profile, quick Bayesian exit, fixes bug in data background and BlittingSupport (#212)

Author: StephenNneji

cpp/RAT

@@ -17,7 +17,7 @@ dependencies = [
     "numpy>=1.20",
     "prettytable>=3.9.0",
     "pydantic>=2.7.2",
-    "scipy>=1.13.1",
+    "scipy>=1.13.1, <=1.17.1",
     "strenum>=0.4.15 ; python_full_version < '3.11'",
     "tqdm>=4.66.5",
 ]

pyproject.toml

@@ -243,7 +243,7 @@ def make_problem(project: ratapi.Project) -> ProblemDefinition:
         contrast_background_param = []
 
         if background.type == TypeOptions.Data:
-            contrast_background_param.append(project.data.index(background.source, True))
+            # Background data is appended to contrast data so empty index
             if background.value_1 != "":
                 contrast_background_param.append(project.background_parameters.index(background.value_1))
             # If we are using a data background, we add the background data to the contrast data
@@ -492,7 +492,7 @@ def check_indices(problem: ProblemDefinition) -> None:
 
     source_param_lists = {
         "constant": "backgroundParams",
-        "data": "data",
+        "data": "backgroundParams",
         "function": "customFiles",
     }
 
@@ -501,11 +501,13 @@ def check_indices(problem: ProblemDefinition) -> None:
 
         # check source param is in range of the relevant parameter list
         param_list = getattr(problem, source_param_lists[background_type])
-        source_index = background_data[0]
-        if not 0 < source_index <= len(param_list):
+        # For data, background data is appended to the contrast data so first index could
+        # either be empty or be the index of an offset.
+        if background_data and not 0 < background_data[0] <= len(param_list):
+            first_entry_name = "data offset" if background_type == "data" else "source"
             raise IndexError(
                 f'Entry {i} of contrastBackgroundParams has type "{background_type}" '
-                f"and source index {source_index}, "
+                f"and {first_entry_name} index {background_data[0]}, "
                 f'which is outside the range of "{source_param_lists[background_type]}".'
             )
 

ratapi/inputs.py

@@ -66,6 +66,33 @@ def zip_if_several(*params) -> tuple | list[tuple]:
             return zip(*params, strict=False)
         return [params]
 
+    def __opaque_to_string(opaque):
+        """Extract java string from MATLABOpaque object.
+
+        Parameters
+        ----------
+        opaque : MATLABOpaque
+            object containing java string.
+
+        Returns
+        -------
+        str
+            string extracted from object.
+
+        Raises
+        ------
+        ValueError
+            Raised if there is no java string array in the object.
+        """
+        entries = opaque[0]
+        for entry in entries:
+            if isinstance(entry, ndarray):
+                break
+        else:
+            raise ValueError("No array in MatlabOpaque")
+
+        return bytes(entry[7:]).decode("ascii")
+
     def read_param(names, constrs, values, fits):
         """Read in a parameter list from the relevant keys, and fix constraints for non-fit parameters.
 
@@ -217,11 +244,11 @@ def fix_invalid_constraints(name: str, constrs: tuple[float, float], value: floa
     # which is given as the byte data of a Java string; this consists of 7 metadata bytes (ignored)
     # and then the actual string characters (index [7:]) in ascii format (.decode("ascii"))
     if len(mat_project["contrastNames"]) == 1 and isinstance(mat_project["contrastNames"], MatlabOpaque):
-        mat_project["contrastNames"] = bytes(mat_project["contrastNames"][0][3][7:]).decode("ascii")
+        mat_project["contrastNames"] = __opaque_to_string(mat_project["contrastNames"])
     else:
         for i, contrast_name in enumerate(mat_project["contrastNames"]):
             if isinstance(contrast_name, MatlabOpaque):
-                mat_project["contrastNames"][i] = bytes(contrast_name[0][3][7:]).decode("ascii")
+                mat_project["contrastNames"][i] = __opaque_to_string(contrast_name)
 
     # if just one contrast, resolNames is a string; fix that here
     if isinstance(mat_project["resolNames"], str):

ratapi/utils/convert.py

@@ -60,6 +60,7 @@ def _extract_plot_data(event_data: PlotEventData, q4: bool, show_error_bar: bool
         # Plot the reflectivity on plot (1,1)
         results["ref"].append([r[:, 0], r[:, 1] * mult])
 
+        results["error"].append([])
         if event_data.dataPresent[i]:
             sd_x = data[:, 0]
             sd_y, sd_e = map(lambda x: x * mult, (data[:, 1], data[:, 2]))
@@ -73,7 +74,7 @@ def _extract_plot_data(event_data: PlotEventData, q4: bool, show_error_bar: bool
                 valid = np.ones(len(sd_e)).astype(bool)
                 sd_e = errors
 
-            results["error"].append([sd_x[valid], sd_y[valid], sd_e[valid]])
+            results["error"][-1].extend([sd_x[valid], sd_y[valid], sd_e[valid]])
 
         results["sld"].append([])
         for j in range(len(sld)):
@@ -111,6 +112,7 @@ def plot_ref_sld_helper(
     show_legend: bool = True,
     shift_value: float = 100,
     animated=False,
+    align_profile=False,
 ):
     """Clear the previous plots and updates the ref and SLD plots.
 
@@ -160,6 +162,8 @@ def plot_ref_sld_helper(
     ref_plot.cla()
     sld_plot.cla()
 
+    if align_profile:
+        _align_profiles(data, confidence_intervals)
     plot_data = _extract_plot_data(data, q4, show_error_bar, shift_value)
     for i, name in enumerate(data.contrastNames):
         ref_plot.plot(plot_data["ref"][i][0], plot_data["ref"][i][1], label=name, linewidth=1, animated=animated)
@@ -234,6 +238,54 @@ def plot_ref_sld_helper(
         plt.pause(0.005)
 
 
+def _align_profiles(data: PlotEventData, confidence_intervals: dict | None = None):
+    """Align SLD profiles and resampled layers.
+
+    Aligns the A/L SLD profiles so that the substrates line up by padding the
+    start of any shorter than the longest profile.
+
+    Parameters
+    ----------
+    data : PlotEventData
+        The plot event data that contains all the information
+        to generate the ref and sld plots
+    confidence_intervals : dict or None, default None
+        The Bayesian confidence intervals for reflectivity and SLD.
+        Only relevant if the procedure used is Bayesian (NS or DREAM)
+    """
+    slds = data.sldProfiles
+    size = (len(slds), len(slds[0]))
+
+    # Find the length of the longest profile.
+    lengths = [[sld.shape[0] for sld in sld_row] for sld_row in slds]
+    max_value = np.max(lengths)
+    max_index = np.unravel_index(np.argmax(lengths), shape=size)
+
+    max_x = slds[max_index[0]][max_index[1]][:, 0]
+    max_x_value = max_x[-1]
+
+    for i in range(size[0]):
+        for j in range(size[1]):
+            cur_sld = slds[i][j]
+            diff = max_value - cur_sld.shape[0]
+            if diff:
+                pad = np.zeros(diff)
+                max_y = np.concatenate((pad, cur_sld[:, 1]))
+                slds[i][j] = np.column_stack((max_x, max_y))
+
+                cur_resample_layer = data.resampledLayers[i][j]
+                if not np.all(cur_resample_layer):
+                    total_length = sum(cur_resample_layer[:, 0])
+                    offset = max_x_value - total_length
+                    data.resampledLayers[i][j] = np.vstack(([offset, 0, 0], cur_resample_layer))
+
+                if confidence_intervals is not None:
+                    cur_ci = confidence_intervals["sld"][i][j]
+                    inter_a = np.concatenate((pad, cur_ci[0]))
+                    inter_b = np.concatenate((pad, cur_ci[1]))
+                    confidence_intervals["sld"][i][j] = (inter_a, inter_b)
+
+
 def plot_ref_sld(
     project: ratapi.Project,
     results: ratapi.outputs.Results | ratapi.outputs.BayesResults,
@@ -294,7 +346,7 @@ def plot_ref_sld(
     data.reflectivity = copy.deepcopy(results.reflectivity)
     data.shiftedData = results.shiftedData
     data.sldProfiles = copy.deepcopy(results.sldProfiles)
-    data.resampledLayers = results.resampledLayers
+    data.resampledLayers = copy.deepcopy(results.resampledLayers)
     data.dataPresent = ratapi.inputs.make_data_present(project)
     data.subRoughs = results.contrastParams.subRoughs
     data.resample = ratapi.inputs.make_resample(project)
@@ -355,6 +407,7 @@ def plot_ref_sld(
         show_grid=show_grid,
         show_legend=show_legend,
         shift_value=shift_value,
+        align_profile=(project.geometry == "air/substrate" and project.model != "custom xy"),
     )
 
     if return_fig:
@@ -533,13 +586,12 @@ def update_foreground(self, data):
         self.figure.canvas.restore_region(self.bg)
         plot_data = _extract_plot_data(data, self.q4, self.show_error_bar, self.shift_value)
 
-        offset = 2
-        for i in range(
-            0,
-            len(self.figure.axes[0].lines),
-        ):
-            self.figure.axes[0].lines[i].set_data(plot_data["ref"][i // offset][0], plot_data["ref"][i // offset][1])
-            self.figure.axes[0].draw_artist(self.figure.axes[0].lines[i])
+        offset = 0
+        for i in range(len(data.contrastNames)):
+            for _ in range(int(data.dataPresent[i]) + 1):
+                self.figure.axes[0].lines[offset].set_data(plot_data["ref"][i][0], plot_data["ref"][i][1])
+                self.figure.axes[0].draw_artist(self.figure.axes[0].lines[offset])
+                offset += 1
 
         i = 0
         for j in range(len(plot_data["sld"])):
@@ -553,12 +605,14 @@ def update_foreground(self, data):
                 self.figure.axes[1].draw_artist(self.figure.axes[1].lines[i])
                 i += 1
 
-        for i, container in enumerate(self.figure.axes[0].containers):
-            self.adjust_error_bar(
-                container, plot_data["error"][i][0], plot_data["error"][i][1], plot_data["error"][i][2]
-            )
-            self.figure.axes[0].draw_artist(container[2][0])
-            self.figure.axes[0].draw_artist(container[0])
+        i = 0
+        for error in plot_data["error"]:
+            if error:
+                container = self.figure.axes[0].containers[i]
+                self.adjust_error_bar(container, error[0], error[1], error[2])
+                self.figure.axes[0].draw_artist(container[2][0])
+                self.figure.axes[0].draw_artist(container[0])
+                i += 1
 
         self.figure.canvas.blit(self.figure.bbox)
         self.figure.canvas.flush_events()

ratapi/utils/plotting.py

@@ -534,17 +534,18 @@ def test_background_params_source_indices(self, test_problem, background_type, b
         test_problem = request.getfixturevalue(test_problem)
         test_problem.contrastBackgroundParams = bad_value
         test_problem.contrastBackgroundTypes = [background_type]
-
         source_param_lists = {
             "constant": "backgroundParams",
-            "data": "data",
+            "data": "backgroundParams",
             "function": "customFiles",
         }
 
+        first_entry_name = "data offset" if background_type == "data" else "source"
+
         with pytest.raises(
             IndexError,
             match=f'Entry 0 of contrastBackgroundParams has type "{background_type}" '
-            f"and source index {bad_value[0][0]}, "
+            f"and {first_entry_name} index {bad_value[0][0]}, "
             f'which is outside the range of "{source_param_lists[background_type]}".',
         ):
             check_indices(test_problem)

tests/test_inputs.py

@@ -204,9 +204,14 @@ def test_plot_ref_sld(mock: MagicMock, input_project, reflectivity_calculation_r
         for sld, sld_results in zip(sldProfile, result_sld_profile, strict=False):
             assert (sld == sld_results).all()
 
+    for resampledLayers, reflectivity_results in zip(
+        data.resampledLayers, reflectivity_calculation_results.resampledLayers, strict=False
+    ):
+        for resam, resam_results in zip(resampledLayers, reflectivity_results, strict=False):
+            assert (resam == resam_results).all()
+
     assert data.modelType == input_project.model
     assert data.shiftedData == reflectivity_calculation_results.shiftedData
-    assert data.resampledLayers == reflectivity_calculation_results.resampledLayers
     assert data.dataPresent.size == 0
     assert (data.subRoughs == reflectivity_calculation_results.contrastParams.subRoughs).all()
     assert data.resample.size == 0