loop roll update

ramav87 · ramav87 · commit 0fdf94883790 · 2026-04-13T07:11:49.000-04:00
diff --git a/examples/mcp/fit_beps_dataset_workflow.md b/examples/mcp/fit_beps_dataset_workflow.md
@@ -6,10 +6,13 @@ This is the small-model-friendly workflow for the common BEPS fitting task.
 
 1. Read the file with the SciFiReaders MCP server.
 2. Pass the original file path to the sidpy workflow tool.
-3. Fit SHO over the full DC sweep for one cycle.
+3. Fit SHO over the full DC sweep for one selected cycle only.
 4. Derive the loop input by projecting the fitted SHO amplitude and phase
    with BGlib `projectLoop`.
-5. Save the SHO and BEPS fit-parameter datasets.
+5. Rotate the projected loop and DC axis together so the single cycle starts
+   at a consistent DC point, then scale the projected loop by `1e3` before
+   BEPS fitting.
+6. Save the SHO and BEPS fit-parameter datasets.
 
 There are two valid input styles:
 
@@ -82,8 +85,10 @@ This one sidpy tool performs the full sequence internally:
 3. Save the SHO fit-parameter map as a `sidpy.Dataset` with a DC axis.
 4. Project the fitted SHO amplitude and phase into a piezoresponse loop with
    BGlib `projectLoop`.
-5. Fit the BEPS loop slice from that projected piezoresponse.
-6. Save the BEPS fit-parameter map as a `sidpy.Dataset`.
+5. Rotate the projected loop and the DC axis together so the two branches are
+   contiguous, then scale the projected loop by `1e3`.
+6. Fit the BEPS loop slice from that normalized projected piezoresponse.
+7. Save the BEPS fit-parameter map as a `sidpy.Dataset`.
 
 ## Expected Result
 
@@ -100,5 +105,10 @@ The tool returns:
 - Keep the workflow in this exact order: optionally read with SciFiReaders, then one sidpy workflow call.
 - The SHO metadata stores the loop slice indices used by the derive step.
 - The SHO metadata stores the selected cycle index used by the derive step.
-- The loop fit uses the projected piezoresponse, not the raw loop waveform.
+- The loop fit uses the normalized projected piezoresponse, not the raw loop
+  waveform.
+- The loop fitter only accepts one cycle at a time. If the source file
+  contains multiple cycles, select one cycle first.
+- The derived loop metadata stores the roll used to align the cycle. Reuse
+  that roll when plotting or recomputing R² so the points stay aligned.
 - The saved fit datasets preserve the `X` and `Y` axes from the source file.
diff --git a/sidpy/proc/mcp_server_beps.py b/sidpy/proc/mcp_server_beps.py
@@ -52,9 +52,9 @@
             "goal": (
                 "Read a BEPS HDF5/NSID file with the SciFiReaders MCP server, "
                 "then run one sidpy workflow tool that fits SHO over the full DC "
-                "sweep for one cycle, derives the loop input by projecting the "
-                "fitted SHO amplitude and phase with BGlib projectLoop, and "
-                "saves both fit-parameter maps as sidpy.Datasets."
+                "sweep for one selected cycle only, derives the loop input by "
+                "projecting the fitted SHO amplitude and phase with BGlib "
+                "projectLoop, and saves both fit-parameter maps as sidpy.Datasets."
             ),
             "inputs": {
                 "file_path": "/path/to/PTO_5x5.h5",
@@ -737,6 +737,24 @@ def _project_piezoresponse_loop(
     return dict(results)
 
 
+def _align_loop_trace(
+    dc_voltage: Sequence[float],
+    loop_trace: Sequence[Any],
+    roll_steps: Optional[int] = None,
+) -> tuple[np.ndarray, np.ndarray, int]:
+    """Rotate a loop trace so the cycle starts at the minimum DC voltage."""
+    dc = np.asarray(dc_voltage, dtype=float).reshape(-1)
+    trace = np.asarray(loop_trace, dtype=float)
+    if dc.size != trace.shape[-1]:
+        raise ValueError("dc_voltage and loop_trace must have matching last-axis lengths.")
+
+    if roll_steps is None:
+        roll_steps = -int(np.argmin(dc))
+    else:
+        roll_steps = int(roll_steps)
+    return np.roll(dc, roll_steps), np.roll(trace, roll_steps, axis=-1), roll_steps
+
+
 def get_dataset(dataset_id: str) -> Dict[str, Any]:
     """Return a JSON-friendly summary of a registered dataset."""
     dataset = _get_dataset(dataset_id)
@@ -768,8 +786,9 @@ def derive_loop_input_from_sho_result(sho_dataset_id: str) -> Dict[str, Any]:
         raise ValueError("The SHO result dataset must include amplitude and phase fit parameters.")
 
     dc_voltage = np.asarray(source_dataset._axes[3].values)
-    sho_amplitude = np.asarray(sho_dataset[..., 0], dtype=float)
-    sho_phase = np.asarray(sho_dataset[..., 3], dtype=float)
+    sho_array = np.asarray(sho_dataset)
+    sho_amplitude = np.asarray(sho_array[..., 0], dtype=float)
+    sho_phase = np.asarray(sho_array[..., 3], dtype=float)
 
     projected = np.zeros(sho_amplitude.shape[:2] + (dc_voltage.size,), dtype=float)
     rotation = np.zeros(sho_amplitude.shape[:2] + (2,), dtype=float)
@@ -787,12 +806,16 @@ def derive_loop_input_from_sho_result(sho_dataset_id: str) -> Dict[str, Any]:
             centroid[row, col] = np.asarray(projection["Centroid"], dtype=float)
             geometric_area[row, col] = float(projection["Geometric Area"])
 
-    dc_voltage = np.asarray(source_dataset._axes[3].values)
+    # Rotate the loop so it starts at the minimum DC voltage. That makes the
+    # first half the increasing branch and the second half the decreasing branch.
+    dc_voltage, projected, roll_steps = _align_loop_trace(dc_voltage, projected)
+    projected = projected * 1e3
 
     return {
         "source_dataset_id": source_dataset_id,
         "source_slice": {
             "cycle_index": cycle_index,
+            "loop_roll_steps": roll_steps,
             "signal": "projected_piezoresponse",
             "projection_tool": "BGlib.projectLoop",
         },
@@ -919,6 +942,7 @@ def fit_beps_loops(
     dataset_name: str = "beps_loop_data",
 ) -> Dict[str, Any]:
     """Fit BEPS loop data where the last axis contains the loop trace."""
+    _validate_single_cycle_dc_voltage(dc_voltage)
     dataset = _build_dataset(
         data,
         dc_voltage,
@@ -1127,6 +1151,63 @@ def _r2_score(y_true: Sequence[Any], y_pred: Sequence[Any]) -> float:
     return 1.0 - (ss_res / ss_tot)
 
 
+def _branch_r2_scores(dc_voltage: Sequence[Any], y_true: Sequence[Any], y_pred: Sequence[Any]) -> Dict[str, float]:
+    """Compute R^2 scores for the increasing and decreasing halves of a loop."""
+    dc = np.asarray(dc_voltage, dtype=float).reshape(-1)
+    true = np.asarray(y_true, dtype=float)
+    pred = np.asarray(y_pred, dtype=float)
+    if dc.size != true.shape[-1] or dc.size != pred.shape[-1]:
+        raise ValueError("dc_voltage, y_true, and y_pred must have matching last-axis lengths.")
+    if dc.size < 2:
+        return {
+            "increasing_branch_r2": float("nan"),
+            "decreasing_branch_r2": float("nan"),
+        }
+
+    half = dc.size // 2
+    first_dc = dc[:half]
+    second_dc = dc[half:]
+    first_true = true[..., :half]
+    first_pred = pred[..., :half]
+    second_true = true[..., half:]
+    second_pred = pred[..., half:]
+
+    first_score = _r2_score(first_true, first_pred)
+    second_score = _r2_score(second_true, second_pred)
+
+    first_slope = float(np.mean(np.diff(first_dc))) if first_dc.size > 1 else 0.0
+    second_slope = float(np.mean(np.diff(second_dc))) if second_dc.size > 1 else 0.0
+    if first_slope >= second_slope:
+        increasing_score, decreasing_score = first_score, second_score
+    else:
+        increasing_score, decreasing_score = second_score, first_score
+
+    return {
+        "increasing_branch_r2": float(increasing_score),
+        "decreasing_branch_r2": float(decreasing_score),
+    }
+
+
+def _validate_single_cycle_dc_voltage(dc_voltage: Sequence[Any]) -> None:
+    """Ensure the loop input contains one forward branch and one reverse branch."""
+    dc = np.asarray(dc_voltage, dtype=float).reshape(-1)
+    if dc.size < 3:
+        return
+
+    diffs = np.diff(dc)
+    diffs = diffs[np.abs(diffs) > np.finfo(float).eps]
+    if diffs.size == 0:
+        return
+
+    signs = np.sign(diffs)
+    sign_changes = np.count_nonzero(signs[1:] * signs[:-1] < 0)
+    if sign_changes > 2:
+        raise ValueError(
+            "BEPS loop fitting expects one cycle with a single forward branch and a single "
+            "reverse branch. Provide one cycle at a time."
+        )
+
+
 def _fit_parameter_dataset(
     parameters: Sequence[Any],
     *,
@@ -1326,6 +1407,7 @@ def fit_beps_dataset_from_reader_payload(
     ]).reshape(beps_data.shape)
     beps_quality = {
         "overall_r2": _r2_score(beps_data, beps_pred),
+        **_branch_r2_scores(dc_voltage, beps_data, beps_pred),
     }
     beps_dataset = _fit_parameter_dataset(
         beps_params,
@@ -1491,6 +1573,7 @@ def fit_beps_dataset_from_scifireaders_file(
     ]).reshape(beps_data.shape)
     beps_quality = {
         "overall_r2": _r2_score(beps_data, beps_pred),
+        **_branch_r2_scores(dc_voltage, beps_data, beps_pred),
     }
     beps_dataset = _fit_parameter_dataset(
         beps_params,
diff --git a/tests/proc/test_mcp_server_beps.py b/tests/proc/test_mcp_server_beps.py
@@ -232,6 +232,88 @@ def test_fit_parameter_dataset_preserves_dc_axis_for_4d_sho_maps(self):
         self.assertEqual(result["dimensions"][3]["name"], "fit_parameter")
         self.assertEqual(result["dimensions"][3]["values"], [0, 1, 2, 3, 4])
 
+    def test_align_loop_trace_honors_explicit_roll_steps(self):
+        dc_voltage = np.array([1.0, 2.0, -2.0, -1.0, 0.0])
+        loop_trace = np.array([10.0, 11.0, 12.0, 13.0, 14.0])
+
+        aligned_dc, aligned_trace, applied_roll = mcp_mod._align_loop_trace(
+            dc_voltage,
+            loop_trace,
+            roll_steps=-2,
+        )
+
+        np.testing.assert_array_equal(aligned_dc, np.array([-2.0, -1.0, 0.0, 1.0, 2.0]))
+        np.testing.assert_array_equal(aligned_trace, np.array([12.0, 13.0, 14.0, 10.0, 11.0]))
+        self.assertEqual(applied_roll, -2)
+
+    def test_align_loop_trace_infers_minimum_dc_roll(self):
+        dc_voltage = np.array([1.0, 2.0, -2.0, -1.0, 0.0])
+        loop_trace = np.array([10.0, 11.0, 12.0, 13.0, 14.0])
+
+        aligned_dc, aligned_trace, applied_roll = mcp_mod._align_loop_trace(dc_voltage, loop_trace)
+
+        np.testing.assert_array_equal(aligned_dc, np.array([-2.0, -1.0, 0.0, 1.0, 2.0]))
+        np.testing.assert_array_equal(aligned_trace, np.array([12.0, 13.0, 14.0, 10.0, 11.0]))
+        self.assertEqual(applied_roll, -2)
+
+    def test_real_file_path_workflow_returns_sho_and_beps_datasets(self):
+        from pathlib import Path
+        import shutil
+        import os
+        import tempfile
+
+        os.environ.setdefault("NUMBA_DISABLE_JIT", "1")
+        try:
+            import SciFiReaders as sr
+        except Exception as exc:
+            self.skipTest(f"SciFiReaders is required for the workflow file-path regression test: {exc}")
+
+        file_path = Path("/Users/rvv/Downloads/PTO_5x5.h5")
+        if not file_path.exists():
+            self.skipTest(f"Real BEPS fixture not found: {file_path}")
+
+        fd, tmp_name = tempfile.mkstemp(prefix="pto_5x5_sidpy_", suffix=".h5")
+        os.close(fd)
+        workflow_path = Path(tmp_name)
+        shutil.copy2(file_path, workflow_path)
+
+        reader = sr.NSIDReader(str(workflow_path))
+        data = reader.read()["Channel_000"]
+        array = np.asarray(data)
+
+        result = mcp_mod.fit_beps_dataset_from_scifireaders_file(
+            str(workflow_path),
+            channel_name="Channel_000",
+            dataset_name="PTO_5x5.h5",
+            sho_cycle_index=1,
+            use_kmeans=False,
+            n_clusters=4,
+            return_cov=False,
+            loss="linear",
+            sho_dataset_name="pto_5x5_sho_workflow",
+            beps_dataset_name="pto_5x5_beps_workflow",
+        )
+
+        self.assertIn("sho_dataset_id", result)
+        self.assertIn("beps_dataset_id", result)
+        self.assertEqual(result["sho_dataset"]["shape"], [5, 5, array.shape[3], 4])
+        self.assertEqual(result["beps_dataset"]["shape"], [5, 5, 9])
+        self.assertEqual(
+            result["loop_input"]["source_slice"],
+            {
+                "cycle_index": 1,
+                "signal": "projected_piezoresponse",
+                "projection_tool": "BGlib.projectLoop",
+            },
+        )
+        self.assertEqual(result["sho_dataset"]["metadata"]["fit_kind"], "sho")
+        self.assertEqual(result["beps_dataset"]["metadata"]["fit_kind"], "beps_loop")
+        self.assertGreater(result["sho_dataset"]["metadata"]["fit_quality"]["overall_r2"], 0.8)
+        self.assertGreater(result["sho_dataset"]["metadata"]["fit_quality"]["amplitude_overall_r2"], 0.9)
+        self.assertGreater(result["beps_dataset"]["metadata"]["fit_quality"]["overall_r2"], 0.75)
+        self.assertIn("increasing_branch_r2", result["beps_dataset"]["metadata"]["fit_quality"])
+        self.assertIn("decreasing_branch_r2", result["beps_dataset"]["metadata"]["fit_quality"])
+
     def test_real_file_scifireaders_payload_can_drive_the_executable_workflow(self):
         from pathlib import Path
         import shutil
@@ -348,10 +430,13 @@ def test_real_file_scifireaders_payload_can_drive_the_executable_workflow(self):
         loop_input = mcp_mod.derive_loop_input_from_sho_result(sho_dataset["dataset_id"])
         self.assertIn("beps_data", loop_input)
         self.assertIn("dc_voltage", loop_input)
+        expected_roll_steps = -int(np.argmin(beps_axis))
+        expected_dc_voltage = np.roll(beps_axis, expected_roll_steps)
         self.assertEqual(
             loop_input["source_slice"],
             {
                 "cycle_index": sho_cycle_idx,
+                "loop_roll_steps": expected_roll_steps,
                 "signal": "projected_piezoresponse",
                 "projection_tool": "BGlib.projectLoop",
             },
@@ -367,7 +452,10 @@ def test_real_file_scifireaders_payload_can_drive_the_executable_workflow(self):
                 for col in range(sho_params.shape[1])
             ]
         ).reshape(beps_data.shape)
+        expected_projected = np.roll(expected_projected, expected_roll_steps, axis=-1)
+        expected_projected = expected_projected * 1e3
         np.testing.assert_allclose(np.asarray(loop_input["beps_data"]), expected_projected)
+        np.testing.assert_allclose(np.asarray(loop_input["dc_voltage"]), expected_dc_voltage)
 
         beps_result = mcp_mod.fit_beps_loops(
             loop_input["beps_data"],
@@ -388,6 +476,25 @@ def test_real_file_scifireaders_payload_can_drive_the_executable_workflow(self):
         ]).reshape(np.asarray(loop_input["beps_data"]).shape)
         beps_r2 = r2_score(np.asarray(loop_input["beps_data"]).reshape(-1), beps_pred.reshape(-1))
         self.assertGreater(beps_r2, 0.75)
+        branch_quality = mcp_mod._branch_r2_scores(
+            np.asarray(loop_input["dc_voltage"]),
+            np.asarray(loop_input["beps_data"]),
+            beps_pred,
+        )
+        self.assertIn("increasing_branch_r2", branch_quality)
+        self.assertIn("decreasing_branch_r2", branch_quality)
+
+    def test_fit_beps_loops_rejects_multi_cycle_dc_voltage(self):
+        dc_voltage = np.array([
+            -1.0, -0.5, 0.0, 0.5, 1.0,
+            0.5, 0.0, -0.5, -1.0,
+            -0.5, 0.0, 0.5, 1.0,
+            0.5, 0.0, -0.5, -1.0,
+        ])
+        data = np.zeros((1, 1, dc_voltage.size))
+
+        with self.assertRaises(ValueError):
+            mcp_mod.fit_beps_loops(data, dc_voltage)
 
     def test_real_file_fits_round_trip_to_sidpy_datasets(self):
         from pathlib import Path
@@ -942,6 +1049,7 @@ def test_stdio_mcp_client_can_run_full_beps_dataset_workflow(self):
                 payload["loop_input"]["source_slice"],
                 {
                     "cycle_index": 1,
+                    "loop_roll_steps": expected_roll_steps,
                     "signal": "projected_piezoresponse",
                     "projection_tool": "BGlib.projectLoop",
                 },
@@ -961,6 +1069,9 @@ def test_stdio_mcp_client_can_run_full_beps_dataset_workflow(self):
                     for col in range(sho_params.shape[1])
                 ]
             ).reshape(beps_data.shape)
+            expected_roll_steps = -int(np.argmin(beps_axis))
+            expected_projected = np.roll(expected_projected, expected_roll_steps, axis=-1)
+            expected_projected = expected_projected * 1e3
             np.testing.assert_allclose(np.asarray(payload["loop_input"]["beps_data"]), expected_projected)
 
 
