diff --git a/src/rtichoke/processing/combinations.py b/src/rtichoke/processing/combinations.py
index b7909291..24969e64 100644
--- a/src/rtichoke/processing/combinations.py
+++ b/src/rtichoke/processing/combinations.py
@@ -17,13 +17,17 @@ def create_strata_combinations(stratified_by: str, by: float, breaks) -> pl.Data
     fmt = f"{{:.{decimals}f}}"
 
     if stratified_by == "probability_threshold":
-        upper_bound = breaks[1:]  # breaks
-        lower_bound = breaks[:-1]  # np.roll(upper_bound, 1)
-        # lower_bound[0] = 0.0
+        cutoff_values = np.asarray(breaks)
+        bin_edges = (
+            cutoff_values
+            if cutoff_values[-1] == 1.0
+            else np.append(cutoff_values, 1.0)
+        )
+        upper_bound = bin_edges[1:]
+        lower_bound = bin_edges[:-1]
         mid_point = upper_bound - by / 2
         include_lower_bound = lower_bound > -0.1
-        include_upper_bound = upper_bound == 1.0  # upper_bound != 0.0
-        # chosen_cutoff = upper_bound
+        include_upper_bound = upper_bound == 1.0
         strata = format_strata_column(
             lower_bound=lower_bound,
             upper_bound=upper_bound,
@@ -33,13 +37,13 @@ def create_strata_combinations(stratified_by: str, by: float, breaks) -> pl.Data
         )
 
     elif stratified_by == "ppcr":
+        cutoff_values = np.asarray(breaks)
         strata_mid = breaks[1:]
         lower_bound = strata_mid - by / 2
         upper_bound = strata_mid + by / 2
         mid_point = breaks[1:]
         include_lower_bound = np.ones_like(strata_mid, dtype=bool)
         include_upper_bound = np.zeros_like(strata_mid, dtype=bool)
-        # chosen_cutoff = strata_mid
         strata = np.array([fmt.format(x) for x in strata_mid], dtype=object)
     else:
         raise ValueError(f"Unsupported stratified_by: {stratified_by}")
@@ -52,12 +56,11 @@ def create_strata_combinations(stratified_by: str, by: float, breaks) -> pl.Data
             "mid_point": mid_point,
             "include_lower_bound": include_lower_bound,
             "include_upper_bound": include_upper_bound,
-            # "chosen_cutoff": chosen_cutoff,
             "stratified_by": [stratified_by] * len(strata),
         }
     )
 
-    cutoffs_df = pl.DataFrame({"chosen_cutoff": breaks})
+    cutoffs_df = pl.DataFrame({"chosen_cutoff": cutoff_values})
 
     return bins_df.join(cutoffs_df, how="cross")
 
@@ -92,9 +95,11 @@ def create_breaks_values(probs_vec, stratified_by, by):
     if stratified_by != "probability_threshold":
         breaks = np.quantile(probs_vec, np.linspace(1, 0, int(1 / by) + 1))
     else:
-        breaks = np.round(
-            np.arange(0, 1 + by, by), decimals=len(str(by).split(".")[-1])
-        )
+        decimals = len(str(by).split(".")[-1])
+        n_steps = int(np.floor((1.0 / by) + 1e-12))
+        breaks = np.round(np.arange(n_steps + 1) * by, decimals=decimals)
+        if probs_vec is not None and breaks[-1] != 1.0:
+            breaks = np.append(breaks, 1.0)
     return breaks
 
 
@@ -124,7 +129,6 @@ def _create_aj_data_combinations_binary(
         ]
     )
 
-    # Define values for Cartesian product
     reals_labels = ["real_negatives", "real_positives"]
 
     combinations_frames: list[pl.DataFrame] = [
@@ -152,8 +156,6 @@ def create_aj_data_combinations(
     dfs = [create_strata_combinations(sb, by, breaks) for sb in stratified_by]
     strata_combinations = pl.concat(dfs, how="vertical")
 
-    # strata_enum = pl.Enum(strata_combinations["strata"])
-
     strata_cats = (
         strata_combinations.select(pl.col("strata").unique(maintain_order=True))
         .to_series()
@@ -178,7 +180,6 @@ def create_aj_data_combinations(
         }
     )
 
-    # Define values for Cartesian product
     reals_labels = [
         "real_negatives",
         "real_positives",
diff --git a/tests/test_cutoff_grid_endpoint.py b/tests/test_cutoff_grid_endpoint.py
new file mode 100644
index 00000000..0436c4dd
--- /dev/null
+++ b/tests/test_cutoff_grid_endpoint.py
@@ -0,0 +1,59 @@
+import numpy as np
+import pytest
+
+from rtichoke import prepare_performance_data, prepare_performance_data_times
+from rtichoke.processing.combinations import create_breaks_values
+
+
+@pytest.mark.parametrize(
+    ("by", "expected"),
+    [
+        (0.3, [0.0, 0.3, 0.6, 0.9]),
+        (0.4, [0.0, 0.4, 0.8]),
+        (0.25, [0.0, 0.25, 0.5, 0.75, 1.0]),
+    ],
+)
+def test_probability_threshold_breaks_match_r_seq_semantics(by, expected):
+    breaks = create_breaks_values(None, "probability_threshold", by)
+
+    assert breaks.tolist() == expected
+    assert np.max(breaks) <= 1.0
+
+
+@pytest.mark.parametrize(
+    ("by", "expected_max"),
+    [(0.3, 0.9), (0.4, 0.8), (0.25, 1.0)],
+)
+def test_binary_performance_data_matches_r_cutoff_endpoint(by, expected_max):
+    probs = {"model": np.array([0.05, 0.25, 0.45, 0.65, 0.85, 0.95])}
+    reals = np.array([0, 0, 1, 0, 1, 1])
+
+    result = prepare_performance_data(probs, reals, by=by)
+    cutoffs = result.get_column("chosen_cutoff")
+
+    assert cutoffs.min() == 0.0
+    assert cutoffs.max() == expected_max
+    assert cutoffs.max() <= 1.0
+
+
+@pytest.mark.parametrize(
+    ("by", "expected_max"),
+    [(0.3, 0.9), (0.4, 0.8), (0.25, 1.0)],
+)
+def test_time_performance_data_matches_r_cutoff_endpoint(by, expected_max):
+    probs = {"model": np.array([0.05, 0.25, 0.45, 0.65, 0.85, 0.95])}
+    reals = np.array([0, 0, 1, 0, 1, 1])
+    times = np.array([2.0, 4.0, 3.0, 8.0, 5.0, 6.0])
+
+    result = prepare_performance_data_times(
+        probs,
+        reals,
+        times,
+        fixed_time_horizons=[6.0],
+        by=by,
+    )
+    cutoffs = result.get_column("chosen_cutoff")
+
+    assert cutoffs.min() == 0.0
+    assert cutoffs.max() == expected_max
+    assert cutoffs.max() <= 1.0