diff --git a/docs/kombine/06_compare_thomas_grunkemeier.md b/docs/kombine/06_compare_thomas_grunkemeier.md
index f84ddac..301c8ef 100644
--- a/docs/kombine/06_compare_thomas_grunkemeier.md
+++ b/docs/kombine/06_compare_thomas_grunkemeier.md
@@ -249,7 +249,7 @@ def compare_tg_vs_kombine(  # pylint: disable=too-many-arguments,too-many-positi
   datacard_path,
   conf_level=0.95,
   atol_surv=1e-6,
-  atol_ci=3e-5,
+  atol_ci=2e-4,  # brentq on MINLP 2NLL uses xtol/rtol 1e-4
 ):
   """Compare TG (R) and KoMbine (Python), return merged table and summary."""
   kombine_df, datacard = kombine_binomial_only(datacard_path, conf_level=conf_level)
diff --git a/docs/kombine/compile_km_plots.py b/docs/kombine/compile_km_plots.py
index 4f9c2a0..077c187 100755
--- a/docs/kombine/compile_km_plots.py
+++ b/docs/kombine/compile_km_plots.py
@@ -23,6 +23,7 @@
 import re
 import subprocess
 import sys
+import time
 import warnings
 
 import matplotlib.lines
@@ -838,14 +839,22 @@ def main():
     print("TESTING MODE: Using smaller datasets")
     print("=" * 60)
 
+  t_all = time.perf_counter()
+
   if generate_all or args.km_example:
+    t0 = time.perf_counter()
     plot_km_example(testing=args.testing)
+    print(f"  [km_example] elapsed {time.perf_counter() - t0:.1f}s", flush=True)
 
   if generate_all or args.greenwood:
+    t0 = time.perf_counter()
     plot_compare_to_greenwood(testing=args.testing)
+    print(f"  [greenwood] elapsed {time.perf_counter() - t0:.1f}s", flush=True)
 
   if generate_all or args.p_value:
+    t0 = time.perf_counter()
     plot_compare_p_value(testing=args.testing)
+    print(f"  [p_value] elapsed {time.perf_counter() - t0:.1f}s", flush=True)
 
   if generate_all or args.lung:
     # Handle lung plot testing options
@@ -867,10 +876,15 @@ def main():
             f"Testing panels: {', '.join(p for p, t in zip(panel_names, lung_testing) if t)}")
       print("=" * 60)
 
+    t0 = time.perf_counter()
     plot_lung_dataset(testing=lung_testing, survival_type=args.lung_survival_type)
+    print(f"  [lung] elapsed {time.perf_counter() - t0:.1f}s", flush=True)
 
   print("=" * 60)
-  print("Plot generation complete!")
+  print(
+    f"Plot generation complete! total_elapsed={time.perf_counter() - t_all:.1f}s",
+    flush=True,
+  )
 
 
 if __name__ == "__main__":
diff --git a/kombine/discrete_optimization.py b/kombine/discrete_optimization.py
index e54a716..01e28c2 100644
--- a/kombine/discrete_optimization.py
+++ b/kombine/discrete_optimization.py
@@ -7,6 +7,7 @@
 import typing
 
 import numpy as np
+import scipy.optimize
 
 from .utilities import InspectableCache
 
@@ -351,3 +352,177 @@ def minimize_discrete_single_minimum( #pylint: disable=too-many-locals, too-many
     print(f"{i_min + left:3d} {candidates[i_min]:6.3f} {values[i_min]:9.5g}")
 
   return candidates[i_min], values[i_min]
+
+
+def _innermost_bracket(
+  cache: dict[float, float],
+  level: float,
+  x_outer: float,
+  x_inner: float,
+) -> tuple[float, float]:
+  """Return ``(last_inside, first_outside)`` for the MLE-connected crossing.
+
+  ``first_outside`` is the cached outside point closest to ``x_inner``.
+  ``last_inside`` is ``x_inner`` (or an inside point still closer to the MLE
+  than that outside point).  Do not walk through other inside points: a
+  MINLP profile can dip back below the cut in a disconnected valley.
+  Endpoints must be in ``cache`` with ``f(x_outer) >= level >= f(x_inner)``.
+  """
+  lo = min(x_outer, x_inner)
+  hi = max(x_outer, x_inner)
+  outside = [x for x, fx in cache.items() if lo <= x <= hi and fx >= level]
+  inside = [x for x, fx in cache.items() if lo <= x <= hi and fx < level]
+  if not outside or not inside:
+    raise ValueError(
+      f"No cached sign change for level {level} on [{x_outer}, {x_inner}]"
+    )
+  if x_outer < x_inner:
+    first_outside = max(outside)
+    inner_side = [x for x in inside if x > first_outside]
+    last_inside = max(inner_side) if inner_side else x_inner
+  else:
+    first_outside = min(outside)
+    inner_side = [x for x in inside if x < first_outside]
+    last_inside = min(inner_side) if inner_side else x_inner
+  if last_inside == first_outside:
+    raise ValueError(
+      f"Cached bracket collapsed for level {level}: {last_inside}"
+    )
+  return last_inside, first_outside
+
+
+def _adjacent_outside_bracket(  # pylint: disable=too-many-arguments, too-many-positional-arguments
+  f_cached: collections.abc.Callable[[float], float],
+  level: float,
+  last_inside: float,
+  first_outside: float,
+  xtol: float,
+  rtol: float,
+) -> tuple[float, float]:
+  """Double from ``last_inside`` toward ``first_outside`` until ``f >= level``.
+
+  Returns a left-to-right pair for ``brentq``.
+  """
+  if _is_close(last_inside, first_outside, xtol, rtol):
+    return min(last_inside, first_outside), max(last_inside, first_outside)
+  direction = 1.0 if first_outside > last_inside else -1.0
+  remaining = abs(first_outside - last_inside)
+  step = min(
+    max(
+      xtol,
+      rtol * max(abs(last_inside), abs(first_outside), 1e-12),
+      0.05 * remaining,
+    ),
+    remaining,
+  )
+  x = last_inside
+  while True:
+    x_next = x + direction * step
+    if (x_next - first_outside) * direction >= 0.0:
+      x_next = first_outside
+    if f_cached(x_next) >= level:
+      return (x, x_next) if x < x_next else (x_next, x)
+    if x_next == first_outside:
+      return (x, x_next) if x < x_next else (x_next, x)
+    x = x_next
+    step = min(step * 2.0, abs(first_outside - x))
+    if step <= 0.0:
+      return (x, first_outside) if x < first_outside else (first_outside, x)
+
+
+def cached_level_crossings(  # pylint: disable=too-many-locals, too-many-arguments, too-many-positional-arguments
+  func: collections.abc.Callable[[float], float],
+  x_outer: float,
+  x_inner: float,
+  levels: collections.abc.Sequence[float],
+  *,
+  xtol: float = 1e-4,
+  rtol: float = 1e-4,
+) -> list[float]:
+  """Find where a 1D profile crosses several levels, sharing evaluations.
+
+  ``func`` should be nonnegative and smaller toward ``x_inner`` (e.g.
+  ``2NLL(p) - 2NLL_min`` on one side of the MLE).  ``levels`` are the
+  target heights (e.g. χ² cuts).  After every new evaluation, each
+  remaining cut walks cached points from ``x_inner`` outward and takes
+  the first exit (closest outside point, last inside point), then doubles
+  from that inside point until the cut is exceeded so a thin near-MLE
+  crossing is not skipped.
+
+  If ``func(x_outer) <= level``, the whole interval is inside that cut and
+  the crossing is ``x_outer``.  Remaining cuts are solved widest-bracket
+  first (typically the outermost level).
+  """
+  if not levels:
+    return []
+  if x_outer == x_inner:
+    return [x_inner] * len(levels)
+
+  cache: dict[float, float] = {}
+
+  def f_cached(x: float) -> float:
+    if x not in cache:
+      cache[x] = float(func(x))
+    return cache[x]
+
+  f_cached(x_outer)
+  f_cached(x_inner)
+
+  results: list[float | None] = [None] * len(levels)
+  pending = list(range(len(levels)))
+
+  def remaining_width(index: int) -> float:
+    level = levels[index]
+    if f_cached(x_outer) <= level:
+      return -1.0
+    if f_cached(x_inner) >= level:
+      return -1.0
+    left, right = _innermost_bracket(cache, level, x_outer, x_inner)
+    return abs(right - left)
+
+  while pending:
+    pending.sort(key=lambda i: (remaining_width(i), levels[i]), reverse=True)
+    index = pending.pop(0)
+    level = levels[index]
+    f_outer = f_cached(x_outer)
+    f_inner = f_cached(x_inner)
+    if f_outer <= level:
+      results[index] = x_outer
+      continue
+    if f_inner >= level:
+      results[index] = x_inner
+      continue
+    last_inside, first_outside = _innermost_bracket(cache, level, x_outer, x_inner)
+    left, right = _adjacent_outside_bracket(
+      f_cached, level, last_inside, first_outside, xtol, rtol,
+    )
+    if left == right or _is_close(left, right, xtol, rtol):
+      closer = left if abs(f_cached(left) - level) <= abs(f_cached(right) - level) else right
+      results[index] = closer
+      continue
+    g_left = f_cached(left) - level
+    g_right = f_cached(right) - level
+    if g_left == 0.0:
+      results[index] = left
+      continue
+    if g_right == 0.0:
+      results[index] = right
+      continue
+    results[index] = typing.cast(
+      float,
+      scipy.optimize.brentq(
+        lambda x, _level=level: f_cached(x) - _level,
+        left,
+        right,
+        xtol=xtol,
+        rtol=np.float64(rtol),
+        full_output=False,
+      ),
+    )
+
+  crossings: list[float] = []
+  for x in results:
+    if x is None:
+      raise RuntimeError("cached_level_crossings left a level unsolved")
+    crossings.append(x)
+  return crossings
diff --git a/kombine/kaplan_meier_MINLP.py b/kombine/kaplan_meier_MINLP.py
index f8807d0..ef66c87 100644
--- a/kombine/kaplan_meier_MINLP.py
+++ b/kombine/kaplan_meier_MINLP.py
@@ -38,6 +38,36 @@ def n_choose_d_term_table(n_patients) -> dict[tuple[int, int], float]:
   return table
 
 
+def km_survival_from_risk_counts(
+  r_vals: npt.ArrayLike,
+  s_vals: npt.ArrayLike,
+  *,
+  log_zero_epsilon: float = LOG_ZERO_EPSILON_DEFAULT,
+  cumulative: bool = False,
+) -> float | list[float]:
+  """
+  Kaplan-Meier survival from at-risk and survived counts.
+
+  Matches the Gurobi encoding: skip r == 0; otherwise use
+  log(s + eps) - log(r + eps). If ``cumulative``, return the survival
+  probability after each interval; otherwise the product over all intervals.
+  """
+  r_vals = np.asarray(r_vals, dtype=float)
+  s_vals = np.asarray(s_vals, dtype=float)
+  if r_vals.shape != s_vals.shape:
+    raise ValueError("r_vals and s_vals must have the same shape")
+  log_cum = 0.0
+  probs: list[float] = []
+  for r, s in zip(r_vals.tolist(), s_vals.tolist()):
+    if r > 0.5:
+      log_cum += math.log(s + log_zero_epsilon) - math.log(r + log_zero_epsilon)
+    if cumulative:
+      probs.append(float(math.exp(log_cum)))
+  if cumulative:
+    return probs
+  return float(math.exp(log_cum))
+
+
 class KaplanMeierPatientNLL(KaplanMeierPatientBase):
   """
   A patient with a time and a parameter.
@@ -488,7 +518,11 @@ def __init__(  # pylint: disable=too-many-arguments
     endpoint_epsilon: float = 1e-6,
     log_zero_epsilon: float = LOG_ZERO_EPSILON_DEFAULT, # New parameter for log arguments
     collapse_consecutive_deaths: bool = True,
+    binomial_only: bool = False,
+    patient_wise_only: bool = False,
   ):
+    if binomial_only and patient_wise_only:
+      raise ValueError("binomial_only and patient_wise_only cannot both be True")
     self.__all_patients = all_patients
     self.__parameter_min = parameter_min
     self.__parameter_max = parameter_max
@@ -496,14 +530,113 @@ def __init__(  # pylint: disable=too-many-arguments
     self.__endpoint_epsilon = endpoint_epsilon
     self.__log_zero_epsilon = log_zero_epsilon # Store the epsilon
     self.__collapse_consecutive_deaths = collapse_consecutive_deaths
+    self.__binomial_only = binomial_only
+    self.__patient_wise_only = patient_wise_only
     self.__expected_probability_constraint = None
-    self.__binomial_penalty_constraint = None
-    self.__patient_constraints_for_binomial_only = None
+    self.__risk_set_r_vars = None
+    self.__risk_set_s_vars = None
+    self.__profile_p_died = None
+    self.__profile_p_survived = None
+    self.__profile_log_p_died = None
+    self.__profile_log_p_survived = None
+    # MIP starts from the previous nearby expected_probability solve.
+    self.__mip_start_a: dict[int, float] | None = None
+    self.__mip_start_mode: tuple[bool, bool] | None = None
+    self.__mip_start_profile: dict[str, list[float]] | None = None
     if not np.isfinite(self.__parameter_min and self.__parameter_min != -np.inf):
       raise ValueError("parameter_min must be finite or -inf")
     if not np.isfinite(self.__parameter_max and self.__parameter_max != np.inf):
       raise ValueError("parameter_max must be finite or inf")
 
+  def seed_assignment_starts(self, starts: dict[int, float] | list[int]) -> None:
+    """
+    Seed binary assignment Starts from a previous solve (e.g. another time point).
+
+    ``starts`` may be a dict of patient index -> 0/1, or a list of selected indices.
+    """
+    if isinstance(starts, dict):
+      self.__mip_start_a = {int(k): float(v) for k, v in starts.items()}
+    else:
+      selected = set(int(j) for j in starts)
+      self.__mip_start_a = {
+        j: 1.0 if j in selected else 0.0 for j in range(self.n_patients)
+      }
+    # Mode unknown when seeded externally; allow first solve to use them.
+    self.__mip_start_mode = None
+
+  def _apply_assignment_mip_starts(
+    self,
+    a,
+    *,
+    binomial_only: bool,
+    patient_wise_only: bool,
+  ) -> None:
+    """Apply cached assignment Starts when the constraint mode matches."""
+    mode = (binomial_only, patient_wise_only)
+    if self.__mip_start_mode is not None and self.__mip_start_mode != mode:
+      self.__mip_start_a = None
+      self.__mip_start_mode = None
+      self.__mip_start_profile = None
+    if self.__mip_start_a is None:
+      for j in range(self.n_patients):
+        a[j].Start = GRB.UNDEFINED
+      return
+    for j, value in self.__mip_start_a.items():
+      if j < self.n_patients:
+        a[j].Start = value
+
+  def _store_assignment_mip_starts(
+    self,
+    a,
+    *,
+    binomial_only: bool,
+    patient_wise_only: bool,
+  ) -> None:
+    """Cache assignment incumbents for the next nearby expected_probability."""
+    self.__mip_start_a = {j: float(a[j].X) for j in range(self.n_patients)}
+    self.__mip_start_mode = (binomial_only, patient_wise_only)
+
+  def export_assignment_mip_starts(self) -> dict[int, float] | None:
+    """Return the cached assignment Starts, if any."""
+    if self.__mip_start_a is None:
+      return None
+    return dict(self.__mip_start_a)
+
+  def _apply_profile_mip_starts(self) -> None:
+    """Start continuous profile probabilities from the previous nearby solve."""
+    starts = self.__mip_start_profile
+    if starts is None:
+      return
+    var_map = {
+      "p_died": self.__profile_p_died,
+      "p_survived": self.__profile_p_survived,
+      "log_p_died": self.__profile_log_p_died,
+      "log_p_survived": self.__profile_log_p_survived,
+    }
+    for name, values in starts.items():
+      vars_by_time = var_map.get(name)
+      if vars_by_time is None:
+        continue
+      for i, value in enumerate(values):
+        if i < len(vars_by_time):
+          vars_by_time[i].Start = value
+
+  def _store_profile_mip_starts(self) -> None:
+    """Cache continuous profile incumbents for the next nearby expected_probability."""
+    p_died = self.__profile_p_died
+    p_survived = self.__profile_p_survived
+    log_p_died = self.__profile_log_p_died
+    log_p_survived = self.__profile_log_p_survived
+    if p_died is None or p_survived is None or log_p_died is None or log_p_survived is None:
+      self.__mip_start_profile = None
+      return
+    self.__mip_start_profile = {
+      "p_died": [float(p_died[i].X) for i in range(self.n_times_to_consider)],
+      "p_survived": [float(p_survived[i].X) for i in range(self.n_times_to_consider)],
+      "log_p_died": [float(log_p_died[i].X) for i in range(self.n_times_to_consider)],
+      "log_p_survived": [float(log_p_survived[i].X) for i in range(self.n_times_to_consider)],
+    }
+
   @property
   def all_patients(self) -> list[KaplanMeierPatientNLL]:
     """
@@ -957,6 +1090,53 @@ def n_choose_d_term_table(self) -> dict[tuple[int, int], float]:
     """
     return n_choose_d_term_table(n_patients=self.n_patients)
 
+  def _fixed_selected_counts(self) -> tuple[list[int], list[int], list[int], list[int]]:
+    """
+    At-risk, died, survived, and collapsed sub-death counts for the nominal
+    in-range assignment (used when ``binomial_only`` fixes ``a``).
+    """
+    selected = self.parameter_in_range
+    r_vals: list[int] = []
+    d_vals: list[int] = []
+    s_vals: list[int] = []
+    sub_d_vals: list[int] = []
+    for dt in self.times_to_consider:
+      r_value = int(np.count_nonzero(self.patient_still_at_risk(dt) & selected))
+      d_value = int(np.count_nonzero(self.patient_died(dt) & selected))
+      r_vals.append(r_value)
+      d_vals.append(d_value)
+      s_vals.append(r_value - d_value)
+      for collapsed_time in self._collapsed_time_groups[dt]:
+        sub_d_vals.append(int(np.count_nonzero(
+          self.patient_died(collapsed_time, collapse_consecutive_deaths=False)
+          & selected
+        )))
+    return r_vals, d_vals, s_vals, sub_d_vals
+
+  def _rd_pair_is_feasible(self, i: int, r_value: int, d_value: int) -> bool:
+    """Whether (r, d) at death-time index i can occur given risk-set flow.
+
+    Survivors at time i cannot exceed the next risk-set maximum plus all
+    censoring in ``[times[i], times[i+1])``, and r cannot exceed the previous
+    risk-set maximum.  These bounds are valid for collapsed times because
+    ``times_to_consider`` are group ends and intervening censoring is what
+    splits groups.
+    """
+    if (
+      r_value > self.n_at_risk_max[i]
+      or d_value > self.n_died_max[i]
+      or d_value > r_value
+    ):
+      return False
+    s_value = r_value - d_value
+    if i + 1 < self.n_times_to_consider:
+      cmax = int(self.n_censored_between_times_max[i])
+      if s_value > int(self.n_at_risk_max[i + 1]) + cmax:
+        return False
+    if i > 0 and r_value > int(self.n_at_risk_max[i - 1]):
+      return False
+    return True
+
   def add_counter_variables_and_constraints(
     self,
     model: gp.Model,
@@ -1193,7 +1373,11 @@ def add_binomial_penalty(  # pylint: disable=too-many-locals, too-many-statement
     To complicate things, we only know the overall expected survival probability,
     not the probability of survival in each group.
     So we need to profile those.
+
+    Survived counts enter through ``r - d`` on the choose-(r, d) indicators;
+    the ``s`` tupledict is unused here but kept for call-site compatibility.
     """
+    _ = s  # survived encoded as r - d on choose-(r, d) indicators
 
     #p_i = probability of dying at death time i
     p_died = model.addVars(
@@ -1229,7 +1413,6 @@ def add_binomial_penalty(  # pylint: disable=too-many-locals, too-many-statement
       lb=log_p_bounds[0],
       ub=log_p_bounds[1],
     )
-    sub_d_counter = -1
     for i in range(self.n_times_to_consider):
       model.addGenConstrExp(log_p_died[i], p_died[i], name=f"log_p_died_constr_{i}")
       model.addGenConstrExp(log_p_survived[i], p_survived[i], name=f"log_p_survived_constr_{i}")
@@ -1261,184 +1444,121 @@ def add_binomial_penalty(  # pylint: disable=too-many-locals, too-many-statement
       log_expected_probability == log_p_survived.sum(),
       name="overall_expected_probability_constraint",
     )
+    self.__profile_p_died = p_died
+    self.__profile_p_survived = p_survived
+    self.__profile_log_p_died = log_p_died
+    self.__profile_log_p_survived = log_p_survived
 
-    #Binomial terms
-    #binomial probability = (n_at_risk choose n_died)
-    #                       * dying probability ^ n_died
-    #                       * surviving probability ^ n_survived
-    #  ==> log likelihood = log(n_at_risk choose n_died)
-    #                       + n_died * log(dying probability)
-    #                       + (n_at_risk - n_died) * log(surviving probability)
-    #                     = log(n_at_risk choose n_died)
-    #                       + n_died * log_p_died
-    #                       + (n_at_risk - n_died) * log_p_survived
-
-    use_binomial_penalty_indicator = model.addVar(
-      vtype=GRB.BINARY,
-      name="use_binomial_penalty_indicator",
-    )
-
-    #n_at_risk choose n_died term
     n_choose_d_table = self.n_choose_d_term_table
-    all_n_choose_d_indicator_vars = model.addVars(
-      self.n_times_to_consider, len(n_choose_d_table),
-      vtype=GRB.BINARY,
-      name="n_choose_d_indicator",
-    )
-    n_choose_d_indicator_vars = {
-      (i, n, d): all_n_choose_d_indicator_vars[i, idx]
-      for i in range(self.n_times_to_consider)
-      for idx, (n, d) in enumerate(n_choose_d_table.keys())
-    }
-    n_died_indicator_vars = model.addVars(
-      self.n_times_to_consider, int(max(self.n_died_max)+1),
-      vtype=GRB.BINARY,
-      name="n_died_indicator",
-    )
-    n_survived_indicator_vars = model.addVars(
-      self.n_times_to_consider, self.n_patients + 1,
-      #could probably have somewhat fewer of these: the maximum is n_patients,
-      #but the minimum is not 0.
-      vtype=GRB.BINARY,
-      name="n_survived_indicator",
-    )
     binomial_terms = []
-    for i, time in enumerate(self.times_to_consider):
-      for (r_value, d_value), penalty in n_choose_d_table.items():
-        indicator = n_choose_d_indicator_vars[i, r_value, d_value]
-        model.addGenConstrIndicator(
-          indicator,
-          True,
-          r[i],
-          GRB.EQUAL,
-          r_value,
-          name=f"n_choose_d_indicator_r_{i}_{r_value}_{d_value}",
-        )
-        model.addGenConstrIndicator(
-          indicator,
-          True,
-          d[i],
-          GRB.EQUAL,
-          d_value,
-          name=f"n_choose_d_indicator_d_{i}_{r_value}_{d_value}",
-        )
-        if (
-          r_value > self.n_at_risk_max[i]
-          or d_value > self.n_died_max[i]
-          or d_value > r_value
-        ):
-          model.addConstr(
-            indicator == 0,
-            name=f"n_choose_d_indicator_impossible_{i}_{r_value}_{d_value}",
-          )
-        binomial_terms.append(-penalty * indicator)
-      # Ensure that exactly one n_choose_d_indicator is selected for each death time
-      indicators = [
-        all_n_choose_d_indicator_vars[i, idx]
-        for idx in range(len(n_choose_d_table))
-      ]
-      model.addConstr(
-        gp.quicksum(indicators) == 1,
-        name=f"one_n_choose_d_indicator_per_death_time_{i}",
-      )
 
-      for d_value in range(max(self.n_died_max) + 1):
-        model.addGenConstrIndicator(
-          n_died_indicator_vars[i, d_value],
-          True,
-          d[i],
-          GRB.EQUAL,
-          d_value,
-          name=f"n_died_indicator_{i}_{d_value}",
-        )
-        if d_value > self.n_died_max[i]:
-          model.addConstr(
-            n_died_indicator_vars[i, d_value] == 0,
-            name=f"n_died_indicator_impossible_{i}_{d_value}",
-          )
-        binomial_terms.append(
-          -d_value * log_p_died[i] * n_died_indicator_vars[i, d_value]
-        )
-      # Ensure that exactly one n_died_indicator is selected for each group
-      model.addConstr(
-        gp.quicksum(
-          n_died_indicator_vars[i, d_value]
-          for d_value in range(max(self.n_died_max)+1)
-        ) == 1,
-        name=f"one_n_died_indicator_per_death_time_{i}",
-      )
-
-      for s_value in range(self.n_patients + 1):
-        model.addGenConstrIndicator(
-          n_survived_indicator_vars[i, s_value],
-          True,
-          s[i],
-          GRB.EQUAL,
-          s_value,
-          name=f"n_survived_indicator_{i}_{s_value}",
-        )
-        if s_value > self.n_at_risk_max[i]:
-          model.addConstr(
-            n_survived_indicator_vars[i, s_value] == 0,
-            name=f"n_survived_indicator_impossible_{i}_{s_value}",
+    if self.__binomial_only:
+      _ = (r, d, sub_d)
+      r_vals, d_vals, _s_vals, sub_d_vals = self._fixed_selected_counts()
+      sub_d_offset = 0
+      for i, time in enumerate(self.times_to_consider):
+        r_value = r_vals[i]
+        d_value = d_vals[i]
+        s_value = r_value - d_value
+        penalty = n_choose_d_table[(r_value, d_value)]
+        binomial_terms.append(-penalty)
+        binomial_terms.append(-d_value * log_p_died[i])
+        binomial_terms.append(-s_value * log_p_survived[i])
+        if d_value > 0:
+          binomial_terms.append(
+            -(math.lgamma(d_value + 1) - d_value * np.log(d_value))
           )
-        binomial_terms.append(
-          -s_value * log_p_survived[i] * n_survived_indicator_vars[i, s_value]
-        )
-
-      # Ensure that exactly one n_survived_indicator is selected for each group
-      model.addConstr(
-        gp.quicksum(
-          n_survived_indicator_vars[i, s_value]
-          for s_value in range(self.n_patients + 1)
-        ) == 1,
-        name=f"one_n_survived_indicator_per_death_time_{i}",
-      )
-
-      # Additional term needed when collapsing consecutive deaths
-      # See \ref{sec:collapsing-consecutive-deaths} in the paper
-      # Note that if collapse_consecutive_deaths is False (or if there's
-      # only one death time in the group), we add and subtract the same thing.
-      for sub_d_counter, collapsed_time in enumerate(
-        self._collapsed_time_groups[time],
-        start=sub_d_counter+1
-      ):
-        sub_d_var = sub_d[sub_d_counter]
-        max_sub_d = np.count_nonzero(
-          self.patient_died(collapsed_time, collapse_consecutive_deaths=False)
-        )
-        sub_d_indicators = []
-        for sub_d_value in range(max_sub_d + 1):
-          sub_d_indicator = model.addVar(
+        n_sub = len(self._collapsed_time_groups[time])
+        for sub_d_value in sub_d_vals[sub_d_offset:sub_d_offset + n_sub]:
+          if sub_d_value > 0:
+            binomial_terms.append(
+              math.lgamma(sub_d_value + 1) - sub_d_value * np.log(sub_d_value)
+            )
+        sub_d_offset += n_sub
+    else:
+      # Binomial terms: one SOS1-style choose-(r,d) encoding per death time.
+      # Only flow-feasible (r,d) indicators are created; d*log p_d and s*log p_s
+      # are charged on the same indicator.
+      sub_d_counter = -1
+      for i, time in enumerate(self.times_to_consider):
+        feasible_indicators = []
+        for (r_value, d_value), penalty in n_choose_d_table.items():
+          if not self._rd_pair_is_feasible(i, r_value, d_value):
+            continue
+          indicator = model.addVar(
             vtype=GRB.BINARY,
-            name=f"sub_d_indicator_{i}_{sub_d_counter}_{sub_d_value}",
+            name=f"n_choose_d_indicator_{i}_{r_value}_{d_value}",
           )
-          sub_d_indicators.append(sub_d_indicator)
+          feasible_indicators.append(indicator)
           model.addGenConstrIndicator(
-            sub_d_indicator,
+            indicator,
             True,
-            sub_d_var,
+            r[i],
             GRB.EQUAL,
-            sub_d_value,
-            name=f"sub_d_indicator_constr_{i}_{sub_d_counter}_{sub_d_value}",
+            r_value,
+            name=f"n_choose_d_indicator_r_{i}_{r_value}_{d_value}",
           )
-          if sub_d_value > 0:
+          model.addGenConstrIndicator(
+            indicator,
+            True,
+            d[i],
+            GRB.EQUAL,
+            d_value,
+            name=f"n_choose_d_indicator_d_{i}_{r_value}_{d_value}",
+          )
+          s_value = r_value - d_value
+          binomial_terms.append(-penalty * indicator)
+          binomial_terms.append(-d_value * log_p_died[i] * indicator)
+          binomial_terms.append(-s_value * log_p_survived[i] * indicator)
+          if d_value > 0:
             binomial_terms.append(
-              sub_d_indicator * (
-                math.lgamma(sub_d_value + 1) - sub_d_value * np.log(sub_d_value)
+              -indicator * (
+                math.lgamma(d_value + 1) - d_value * np.log(d_value)
               )
             )
+
+        if not feasible_indicators:
+          raise RuntimeError(
+            f"No feasible (r, d) pairs for death-time index {i} "
+            f"(n_at_risk_max={self.n_at_risk_max[i]}, n_died_max={self.n_died_max[i]})"
+          )
         model.addConstr(
-          gp.quicksum(sub_d_indicators) == 1,
-          name=f"one_sub_d_indicator_per_sub_death_time_{i}_{sub_d_counter}",
+          gp.quicksum(feasible_indicators) == 1,
+          name=f"one_n_choose_d_indicator_per_death_time_{i}",
         )
 
-      for d_value in range(max(self.n_died_max) + 1):
-        if d_value > 0:
-          binomial_terms.append(
-            -n_died_indicator_vars[i, d_value] * (
-              math.lgamma(d_value + 1) - d_value * np.log(d_value)
+        for sub_d_counter, collapsed_time in enumerate(
+          self._collapsed_time_groups[time],
+          start=sub_d_counter+1
+        ):
+          sub_d_var = sub_d[sub_d_counter]
+          max_sub_d = np.count_nonzero(
+            self.patient_died(collapsed_time, collapse_consecutive_deaths=False)
+          )
+          sub_d_indicators = []
+          for sub_d_value in range(max_sub_d + 1):
+            sub_d_indicator = model.addVar(
+              vtype=GRB.BINARY,
+              name=f"sub_d_indicator_{i}_{sub_d_counter}_{sub_d_value}",
             )
+            sub_d_indicators.append(sub_d_indicator)
+            model.addGenConstrIndicator(
+              sub_d_indicator,
+              True,
+              sub_d_var,
+              GRB.EQUAL,
+              sub_d_value,
+              name=f"sub_d_indicator_constr_{i}_{sub_d_counter}_{sub_d_value}",
+            )
+            if sub_d_value > 0:
+              binomial_terms.append(
+                sub_d_indicator * (
+                  math.lgamma(sub_d_value + 1) - sub_d_value * np.log(sub_d_value)
+                )
+              )
+          model.addConstr(
+            gp.quicksum(sub_d_indicators) == 1,
+            name=f"one_sub_d_indicator_per_sub_death_time_{i}_{sub_d_counter}",
           )
 
     binom_penalty_expr = gp.quicksum(binomial_terms)
@@ -1446,6 +1566,19 @@ def add_binomial_penalty(  # pylint: disable=too-many-locals, too-many-statement
       vtype=GRB.CONTINUOUS,
       name="binom_penalty",
     )
+    if self.__binomial_only:
+      model.addConstr(
+        binom_penalty == binom_penalty_expr,
+        name="binomial_penalty_definition",
+      )
+      return binom_penalty, expected_probability_var, None
+
+    # Full NLL: bilinear terms (indicator * log p). Keep the indicator + big-M
+    # sandwich so Gurobi does not treat this as a quadratic equality (status 13).
+    use_binomial_penalty_indicator = model.addVar(
+      vtype=GRB.BINARY,
+      name="use_binomial_penalty_indicator",
+    )
     model.addGenConstrIndicator(
       use_binomial_penalty_indicator,
       False,
@@ -1454,10 +1587,8 @@ def add_binomial_penalty(  # pylint: disable=too-many-locals, too-many-statement
       0.0,
       name="binomial_penalty_inactive",
     )
-
-    #big M constraint to ensure binomial penalty is only used when the indicator is set
     max_penalty_term = max(
-      abs(penalty) for penalty in self.n_choose_d_term_table.values()
+      abs(penalty) for penalty in n_choose_d_table.values()
     )
     max_d = max(self.n_died_max)
     max_s = self.n_patients
@@ -1476,7 +1607,6 @@ def add_binomial_penalty(  # pylint: disable=too-many-locals, too-many-statement
       binom_penalty >= binom_penalty_expr - big_M * (1 - use_binomial_penalty_indicator),
       name="binomial_penalty_expr_lower_bound"
     )
-
     return binom_penalty, expected_probability_var, use_binomial_penalty_indicator
 
   def add_patient_wise_penalty(
@@ -1542,24 +1672,48 @@ def _make_gurobi_model(self):  #pylint: disable=too-many-locals
       a=a,
     )
 
-    # Add Kaplan-Meier probability variables and constraints (replaces trajectory logic)
-    km_probability_var = self.add_kaplan_meier_probability_variables_and_constraints(
-      model=model,
-      r=r,
-      s=s,
-    )
+    km_probability_var = None
+    expected_probability_var = None
+    use_binomial_penalty_indicator = None
+    binom_penalty = 0.0
 
-    (
-      binom_penalty,
-      expected_probability_var,
-      use_binomial_penalty_indicator,
-    ) = self.add_binomial_penalty(
-      model=model,
-      r=r,
-      d=d,
-      sub_d=sub_d,
-      s=s,
-    )
+    if self.__patient_wise_only:
+      km_probability_var = self.add_kaplan_meier_probability_variables_and_constraints(
+        model=model,
+        r=r,
+        s=s,
+      )
+    else:
+      (
+        binom_penalty,
+        expected_probability_var,
+        use_binomial_penalty_indicator,
+      ) = self.add_binomial_penalty(
+        model=model,
+        r=r,
+        d=d,
+        sub_d=sub_d,
+        s=s,
+      )
+      if use_binomial_penalty_indicator is not None:
+        model.addConstr(
+          use_binomial_penalty_indicator == 1,
+          name="use_binomial_penalty",
+        )
+      if self.__binomial_only:
+        for j in range(self.n_patients):
+          if self.parameter_in_range[j]:
+            assert self.nll_penalty_for_patient_in_range[j] <= 0
+            model.addConstr(
+              a[j] == 1,
+              name=f"patient_{j}_must_be_selected_binomial_only",
+            )
+          else:
+            assert self.nll_penalty_for_patient_in_range[j] >= 0
+            model.addConstr(
+              a[j] == 0,
+              name=f"patient_{j}_must_not_be_selected_binomial_only",
+            )
 
     patient_penalty = self.add_patient_wise_penalty(
       model=model,
@@ -1572,6 +1726,8 @@ def _make_gurobi_model(self):  #pylint: disable=too-many-locals
       GRB.MINIMIZE,
     )
     model.update()
+    self.__risk_set_r_vars = r
+    self.__risk_set_s_vars = s
 
     return (
       model,
@@ -1594,58 +1750,30 @@ def update_model_with_expected_probability( # pylint: disable=too-many-arguments
     *,
     model: gp.Model,
     expected_probability: float | None,
-    patient_wise_only: bool,
-    binomial_only: bool,
-    a: gp.tupledict[int, gp.Var],
-    km_probability_var: gp.Var,
-    use_binomial_penalty_indicator: gp.Var,
-    expected_probability_var: gp.Var,
+    km_probability_var: gp.Var | None,
+    expected_probability_var: gp.Var | None,
   ):
     """
     Update the Gurobi model with the expected probability constraint.
     This is the only thing that changes between runs of the MINLP.
     """
-    #drop the previous constraints if they exist
     if self.__expected_probability_constraint is not None:
       model.remove(self.__expected_probability_constraint)
       self.__expected_probability_constraint = None
-    if self.__binomial_penalty_constraint is not None:
-      model.remove(self.__binomial_penalty_constraint)
-      self.__binomial_penalty_constraint = None
-    if self.__patient_constraints_for_binomial_only is not None:
-      for constr in self.__patient_constraints_for_binomial_only:
-        model.remove(constr)
-      self.__patient_constraints_for_binomial_only = None
-
-    if not patient_wise_only:
-      # ---------------------------
-      # Binomial penalty is active, constrain its expected probability
-      # ---------------------------
-      self.__binomial_penalty_constraint = model.addConstr(
-        use_binomial_penalty_indicator == 1,
-        name="use_binomial_penalty"
-      )
+
+    if not self.__patient_wise_only:
       if expected_probability is not None:
+        assert expected_probability_var is not None
         self.__expected_probability_constraint = model.addConstr(
           expected_probability_var == expected_probability,
           name="expected_probability_constraint",
         )
     else:
-      #no binomial penalty means there's nothing to constrain the observed
-      #probability to the expected probability.  In that case, what does
-      #it mean to get an NLL for the expected probability?
-      #Instead, we constrain the observed probability to be at least as
-      #far from the nominal observed probability as the expected
-      #and find the minimum patient-wise NLL.
-      self.__binomial_penalty_constraint = model.addConstr(
-        use_binomial_penalty_indicator == 0,
-        name="use_binomial_penalty"
-      )
-
       # Constrain the KM probability based on the expected_probability
       # If expected > observed, then KM_prob >= expected_probability
       # If expected < observed, then KM_prob <= expected_probability
       # If expected == observed or is None, then KM_prob is unconstrained
+      assert km_probability_var is not None
       if expected_probability is None:
         pass
       elif expected_probability > self.observed_KM_probability:
@@ -1658,31 +1786,9 @@ def update_model_with_expected_probability( # pylint: disable=too-many-arguments
           km_probability_var <= expected_probability + self.__endpoint_epsilon,
           name="km_prob_le_expected"
         )
-      else: # expected_probability == self.observed_KM_probability
+      else:
         assert expected_probability == self.observed_KM_probability
 
-    if binomial_only:
-      self.__patient_constraints_for_binomial_only = []
-      for j in range(self.n_patients):
-        if self.parameter_in_range[j]:
-          assert self.nll_penalty_for_patient_in_range[j] <= 0
-          #the patient must be selected
-          self.__patient_constraints_for_binomial_only.append(
-            model.addConstr(
-              a[j] == 1,
-              name=f"patient_{j}_must_be_selected_binomial_only",
-            )
-          )
-        else:
-          assert self.nll_penalty_for_patient_in_range[j] >= 0
-          #the patient must not be selected
-          self.__patient_constraints_for_binomial_only.append(
-            model.addConstr(
-              a[j] == 0,
-              name=f"patient_{j}_must_not_be_selected_binomial_only",
-            )
-          )
-
     model.update()
 
   def run_MINLP( # pylint: disable=too-many-locals, too-many-statements, too-many-branches, too-many-arguments
@@ -1708,13 +1814,17 @@ def run_MINLP( # pylint: disable=too-many-locals, too-many-statements, too-many-
         "Running MINLP for expected probability ", expected_probability,
         " at time point ", self.time_point, " at time ", datetime.datetime.now()
       )
+    if binomial_only != self.__binomial_only or patient_wise_only != self.__patient_wise_only:
+      raise ValueError(
+        "run_MINLP mode must match MINLPForKM construction "
+        f"(constructed binomial_only={self.__binomial_only}, "
+        f"patient_wise_only={self.__patient_wise_only})"
+      )
     if expected_probability is not None:
       if not patient_wise_only and (expected_probability <= 0 or expected_probability >= 1):
         raise ValueError(f"expected_probability={expected_probability} must be in (0, 1) or None")
       if expected_probability < 0 or expected_probability > 1:
         raise ValueError(f"expected_probability={expected_probability} must be in [0, 1] or None")
-    if binomial_only and patient_wise_only:
-      raise ValueError("binomial_only and patient_wise_only cannot both be True")
 
     if MIPGap is None:
       MIPGap = self.__default_MIPGap
@@ -1728,27 +1838,31 @@ def run_MINLP( # pylint: disable=too-many-locals, too-many-statements, too-many-
       a,
       km_probability_var,
       expected_probability_var,
-      use_binomial_penalty_indicator,
+      _use_binomial_penalty_indicator,
     ) = self.gurobi_model
     self.update_model_with_expected_probability(
       model=model,
-      a=a,
       km_probability_var=km_probability_var,
       expected_probability=expected_probability,
-      patient_wise_only=patient_wise_only,
-      binomial_only=binomial_only,
       expected_probability_var=expected_probability_var,
-      use_binomial_penalty_indicator=use_binomial_penalty_indicator,
     )
 
-    # Initial Gurobi parameters
+    self._apply_assignment_mip_starts(
+      a,
+      binomial_only=binomial_only,
+      patient_wise_only=patient_wise_only,
+    )
+    self._apply_profile_mip_starts()
+
+    # Initial Gurobi parameters. FuncPieces starts at 1000; fallbacks may raise it.
+    # Cuts=-1 restores the default after a fallback that set Cuts=2.
     initial_gurobi_params = {
       'OutputFlag': 1 if verbose else 0,
       'DisplayInterval': 1,
       'MIPGap': MIPGap,
       'MIPGapAbs': MIPGapAbs,
       'NonConvex': 2,
-      'NumericFocus': 3 if patient_wise_only else 0,
+      'NumericFocus': 0,
       'Seed': 123456,
       'TimeLimit': TimeLimit,
       'Threads': Threads,
@@ -1759,91 +1873,42 @@ def run_MINLP( # pylint: disable=too-many-locals, too-many-statements, too-many-
     if LogFile is not None:
       initial_gurobi_params['LogFile'] = os.fspath(LogFile)
 
-    # Define fallback strategies
+    # Recovery fallbacks, then a weaker MIPGap certificate, then verbose debug.
+    # Do not accept SUBOPTIMAL; TimeLimit retries are handled in _optimize_with_fallbacks.
     fallback_strategies = []
-
-    # Fallback 1: Try MIPFocus 2 if initial was suboptimal and not already 2
-    if MIPFocus != 2: # Only add this fallback if MIPFocus wasn't already 2
+    if MIPFocus != 2:
       fallback_strategies.append(
         ({'MIPFocus': 2}, "MIPFocus set to 2 (optimality focus)")
       )
-
-    # Fallback 2: Increase TimeLimit if it was set and still suboptimal
     if TimeLimit is not None:
       fallback_strategies.append(
         ({'TimeLimit': TimeLimit * 1.5}, "Increased TimeLimit by 50%")
       )
-
-    # Fallback 3: Increase FuncPieces
-    current_func_pieces = initial_gurobi_params.get('FuncPieces', 0)
-    if current_func_pieces < 2000: # Arbitrary upper limit to prevent excessive FuncPieces
-      fallback_strategies.append(
-        ({'FuncPieces': max(2000, int(current_func_pieces * 2))}, "Increased FuncPieces (doubled)")
-      )
-    if current_func_pieces < 5000:
-      fallback_strategies.append(
-        (
-          {'FuncPieces': max(5000, int(current_func_pieces * 2.5)), 'FuncPieceRatio': 0.75},
-          "Increased FuncPieces and adjusted FuncPieceRatio"
-        )
-      )
-
-    # Fallback 4: Try different NumericFocus
-    fallback_strategies.append(
-      ({'NumericFocus': 2}, "Changed NumericFocus to 2 (accuracy)")
-    )
-
-    # Fallback 5: Experiment with Cuts (more aggressive)
     fallback_strategies.append(
-      ({'Cuts': 2}, "Aggressive cut generation")
+      ({'FuncPieces': 2000, 'FuncPieceRatio': 0.5}, "Increased FuncPieces to 2000")
     )
-
-    # Fallback 6: Experiment with Heuristics (less aggressive)
     fallback_strategies.append(
-      ({'Heuristics': 0.5}, "Less aggressive heuristics")
+      ({'FuncPieces': 5000, 'FuncPieceRatio': 0.75}, "Increased FuncPieces to 5000")
     )
-
-    #Fallback 7: Tighten barrier convergence tolerance
     fallback_strategies.append(
-      ({'BarConvTol': 1e-8}, "Tightened barrier convergence tolerance")
+      ({'NumericFocus': 3}, "NumericFocus set to 3 (highest precision)")
     )
-
-    #Fallback 8: Tighten feasibility tolerance
     fallback_strategies.append(
-      ({'FeasibilityTol': 1e-8}, "Tightened feasibility tolerance")
-    )
-
-    # Fallback 9: Tighten optimality tolerance
-    fallback_strategies.append(
-      ({'OptimalityTol': 1e-8}, "Tightened optimality tolerance")
-    )
-
-    # Fallback 10: Use barrier method
-    fallback_strategies.append(
-      ({'Method': 2}, "Switched to barrier method")
-    )
-
-    #Fallback 11: Avoid switching back to simplex
-    fallback_strategies.append(
-      ({'CrossOver': 0}, "Avoided switching back to simplex method")
-    )
-
-    # Fallback 12: Barrier at all nodes
-    fallback_strategies.append(
-      ({'NodeMethod': 2}, "Used barrier method at all nodes")
+      ({'Cuts': 2}, "Aggressive cut generation")
     )
-
-    # Last fallback: turn verbose output on
-    # This is not going to work, but it will help us debug the issue
+    if MIPGap < 1e-3:
+      fallback_strategies.append(
+        ({'MIPGap': 1e-3}, "Relaxed MIPGap to 1e-3")
+      )
+    if MIPGap < 1e-2:
+      fallback_strategies.append(
+        ({'MIPGap': 1e-2}, "Relaxed MIPGap to 1e-2")
+      )
     if not verbose:
       fallback_strategies.append(
         (
-          {
-            'OutputFlag': 1,
-            'DisplayInterval': 1,
-            'InfUnbdInfo': 1,
-          },
-          "Turned verbose output on"
+          {'OutputFlag': 1, 'DisplayInterval': 1},
+          "Enabled Gurobi output for debug after failed recovery",
         )
       )
 
@@ -1877,6 +1942,12 @@ def run_MINLP( # pylint: disable=too-many-locals, too-many-statements, too-many-
       )
 
     assert all(var is not None for var in a)
+    self._store_assignment_mip_starts(
+      a,
+      binomial_only=binomial_only,
+      patient_wise_only=patient_wise_only,
+    )
+    self._store_profile_mip_starts()
     selected = [j for j in range(self.n_patients) if a[j].X > 0.5]
     n_total_val = sum(selected)
     n_alive_val = sum(
@@ -1892,14 +1963,27 @@ def run_MINLP( # pylint: disable=too-many-locals, too-many-statements, too-many-
       ) for j in range(self.n_patients)
       if np.isfinite(nll_penalty_for_patient_in_range[j])
     )
-    binom_penalty_var = model.getVarByName("binom_penalty")
-    assert binom_penalty_var is not None
-    binomial_penalty_val = binom_penalty_var.X
-    p_survived_val = []
-    for i in range(self.n_times_to_consider):
-      var = model.getVarByName(f"p_survived[{i}]")
-      assert var is not None
-      p_survived_val.append(var.X)
+    if patient_wise_only:
+      binomial_penalty_val = 0.0
+      p_survived_val = [np.nan] * self.n_times_to_consider
+      assert km_probability_var is not None
+      km_probability_val = km_probability_var.X
+    else:
+      binom_penalty_var = model.getVarByName("binom_penalty")
+      assert binom_penalty_var is not None
+      binomial_penalty_val = binom_penalty_var.X
+      p_survived_val = []
+      for i in range(self.n_times_to_consider):
+        var = model.getVarByName(f"p_survived[{i}]")
+        assert var is not None
+        p_survived_val.append(var.X)
+      assert self.__risk_set_r_vars is not None
+      assert self.__risk_set_s_vars is not None
+      km_probability_val = km_survival_from_risk_counts(
+        [self.__risk_set_r_vars[i].X for i in range(self.n_times_to_consider)],
+        [self.__risk_set_s_vars[i].X for i in range(self.n_times_to_consider)],
+        log_zero_epsilon=self.__log_zero_epsilon,
+      )
     if verbose:
       print("Selected patients:", selected)
       print("n_total:          ", int(n_total_val))
@@ -1918,5 +2002,5 @@ def run_MINLP( # pylint: disable=too-many-locals, too-many-statements, too-many-
       patient_penalties=nll_penalty_for_patient_in_range,
       selected=selected,
       model=model,
-      km_probability=km_probability_var.X,
+      km_probability=km_probability_val,
     )
diff --git a/kombine/kaplan_meier_hazard_ratio_MINLP.py b/kombine/kaplan_meier_hazard_ratio_MINLP.py
index 190554c..fee80ba 100644
--- a/kombine/kaplan_meier_hazard_ratio_MINLP.py
+++ b/kombine/kaplan_meier_hazard_ratio_MINLP.py
@@ -69,6 +69,65 @@ def __init__(  # pylint: disable=too-many-arguments
       tie_handling=tie_handling,
       log_hazard_ratio_bounds=log_hazard_ratio_bounds,
     )
+    # MIP starts from the previous nearby-HR solve on this calculator.
+    self.__mip_start_assignments: dict[tuple[int, int], float] | None = None
+    self.__mip_start_cox_only: bool | None = None
+    self.__mip_start_log_hr: float | None = None
+    # Only reuse starts when |Δ log H| is below this (avoids trapping on distant jumps).
+    self.__mip_start_max_log_hr_delta = 1.0
+
+  def _clear_hazard_ratio_mip_starts(self, a, beta_var=None) -> None:
+    """Drop cached starts and reset Gurobi Start attributes."""
+    self.__mip_start_assignments = None
+    self.__mip_start_cox_only = None
+    self.__mip_start_log_hr = None
+    for var in a.values():
+      var.Start = GRB.UNDEFINED
+    if beta_var is not None:
+      beta_var.Start = GRB.UNDEFINED
+
+  def _apply_hazard_ratio_mip_starts(
+    self,
+    a,
+    beta_var,
+    log_hazard_ratio: float,
+    cox_only: bool,
+  ) -> None:
+    """
+    Seed the next solve from the previous incumbent when cox_only matches
+    and the hazard ratio is close in log space.
+    """
+    if self.__mip_start_cox_only is not None and self.__mip_start_cox_only != cox_only:
+      self._clear_hazard_ratio_mip_starts(a, beta_var)
+      return
+
+    if (
+      self.__mip_start_assignments is None
+      or self.__mip_start_log_hr is None
+      or abs(log_hazard_ratio - self.__mip_start_log_hr)
+        > self.__mip_start_max_log_hr_delta
+    ):
+      for var in a.values():
+        var.Start = GRB.UNDEFINED
+      beta_var.Start = GRB.UNDEFINED
+      return
+
+    for key, value in self.__mip_start_assignments.items():
+      a[key].Start = value
+    beta_var.Start = log_hazard_ratio
+
+  def _store_hazard_ratio_mip_starts(
+    self,
+    a,
+    cox_only: bool,
+    log_hazard_ratio: float,
+  ) -> None:
+    """Cache assignment incumbents for the next nearby-HR solve."""
+    self.__mip_start_assignments = {
+      key: float(var.X) for key, var in a.items()
+    }
+    self.__mip_start_cox_only = cox_only
+    self.__mip_start_log_hr = log_hazard_ratio
 
   def compute_2nll_at_hazard_ratio(  # pylint: disable=too-many-locals, too-many-arguments
     self,
@@ -144,8 +203,6 @@ def compute_2nll_at_hazard_ratio(  # pylint: disable=too-many-locals, too-many-a
       model,
       null_hypothesis_indicator,
       a,
-      km_probability_at_time_low,
-      km_probability_at_time_high,
       beta,
       use_cox_penalty_indicator,
     ) = self.gurobi_model
@@ -183,6 +240,8 @@ def compute_2nll_at_hazard_ratio(  # pylint: disable=too-many-locals, too-many-a
       use_cox_penalty_indicator=use_cox_penalty_indicator,
     )
 
+    self._apply_hazard_ratio_mip_starts(a, beta_var, log_hazard_ratio, cox_only)
+
     # Setup and optimize with standard parameters and fallback strategies
     # pylint: disable=duplicate-code
     model = self._setup_and_optimize(
@@ -200,6 +259,8 @@ def compute_2nll_at_hazard_ratio(  # pylint: disable=too-many-locals, too-many-a
       raise ValueError(f"Optimization failed with status {model.status}")
     # pylint: enable=duplicate-code
 
+    self._store_hazard_ratio_mip_starts(a, cox_only, log_hazard_ratio)
+
     # Extract results
     twonll = model.ObjVal
     patients_low, patients_high = self._extract_patients_per_curve(a)
@@ -209,9 +270,7 @@ def compute_2nll_at_hazard_ratio(  # pylint: disable=too-many-locals, too-many-a
     # Extract curve statistics
     (n_total_low, n_alive_low, km_prob_low,
      n_total_high, n_alive_high, km_prob_high) = (
-      self._extract_curve_statistics(
-        model, km_probability_at_time_low, km_probability_at_time_high
-      )
+      self._extract_curve_statistics(model)
     )
 
     result = scipy.optimize.OptimizeResult(
diff --git a/kombine/kaplan_meier_likelihood.py b/kombine/kaplan_meier_likelihood.py
index 5edb5eb..5092ce3 100644
--- a/kombine/kaplan_meier_likelihood.py
+++ b/kombine/kaplan_meier_likelihood.py
@@ -20,7 +20,7 @@
 import scipy.optimize
 import scipy.stats
 
-from .discrete_optimization import binary_search_sign_change
+from .discrete_optimization import binary_search_sign_change, cached_level_crossings
 from .kaplan_meier import (
   KaplanMeierBase,
   KaplanMeierInstance,
@@ -279,9 +279,12 @@ def nominalkm(self) -> KaplanMeierInstance:
   def minlp_for_km(
     self,
     time_point: float,
+    *,
+    binomial_only: bool = False,
+    patient_wise_only: bool = False,
   ):
     """
-    Get the MINLP for the given time point.
+    Get the MINLP for the given time point and likelihood mode.
     """
     return MINLPForKM(
       all_patients=self.all_patients,
@@ -291,6 +294,8 @@ def minlp_for_km(
       endpoint_epsilon=self.__endpoint_epsilon,
       log_zero_epsilon=self.__log_zero_epsilon,
       collapse_consecutive_deaths=self.__collapse_consecutive_deaths,
+      binomial_only=binomial_only,
+      patient_wise_only=patient_wise_only,
     )
 
   def get_twoNLL_function( # pylint: disable=too-many-arguments
@@ -304,9 +309,11 @@ def get_twoNLL_function( # pylint: disable=too-many-arguments
     MIPGap=None,
     MIPGapAbs=None,
     rerun_until_convergence=False,
+    assignment_starts=None,
   ) -> tuple[
-    collections.abc.Callable[[float | None], scipy.optimize.OptimizeResult],
-    collections.abc.Callable[[float | None], float],
+    InspectableCache[float | None, scipy.optimize.OptimizeResult],
+    InspectableCache[float | None, float],
+    MINLPForKM,
   ]:
     """
     Get the twoNLL function for the given time point.
@@ -321,7 +328,14 @@ def get_twoNLL_function( # pylint: disable=too-many-arguments
     if MIPGapAbs is None:
       MIPGapAbs = self.__default_MIPGapAbs
 
-    minlp = self.minlp_for_km(time_point=time_point)
+    minlp = self.minlp_for_km(
+      time_point=time_point,
+      binomial_only=binomial_only,
+      patient_wise_only=patient_wise_only,
+    )
+    if assignment_starts is not None:
+      minlp.seed_assignment_starts(assignment_starts)
+
     @InspectableCache
     def run_MINLP(expected_probability: float | None) -> scipy.optimize.OptimizeResult:
       """
@@ -388,7 +402,7 @@ def twoNLL(expected_probability: float | None) -> float:
       if not result.success:
         return np.inf
       return result.x
-    return run_MINLP, twoNLL
+    return run_MINLP, twoNLL, minlp
 
   def calculate_possible_probabilities(self, time_point: float) -> np.ndarray:
     """
@@ -513,7 +527,7 @@ def survival_probabilities_likelihood( # pylint: disable=too-many-locals, too-ma
         )
       survival_probabilities_time_point = []
       survival_probabilities.append(survival_probabilities_time_point)
-      run_MINLP, twoNLL = self.get_twoNLL_function(
+      run_MINLP, twoNLL, _minlp = self.get_twoNLL_function(
         time_point=t,
         binomial_only=binomial_only,
         patient_wise_only=patient_wise_only,
@@ -544,34 +558,26 @@ def survival_probabilities_likelihood( # pylint: disable=too-many-locals, too-ma
           1 - self.__endpoint_epsilon,
         )
 
-      for CL in CLs:
-        if patient_wise_only and t < min(self.patient_death_times):
-          # If the time point is outside the range of patient times, we cannot
-          # calculate a patient-wise survival probability.
-          survival_probabilities_time_point.append((1, 1))
-          continue
-
-        d2NLLcut = scipy.stats.chi2.ppf(CL, 1).item()
-        def objective_function(
-          expected_probability: float,
-          twoNLL=twoNLL, twoNLL_min=twoNLL_min, d2NLLcut=d2NLLcut
-        ) -> float:
-          return twoNLL(expected_probability) - twoNLL_min - d2NLLcut
-        if best_prob == best_prob_clipped:
-          #only do this if it's not too close to the edge
-          #to avoid edge effects
-          np.testing.assert_allclose(
-            objective_function(
-              best_prob
-            ),
-            -d2NLLcut,
-            atol=1e-2,
-          )
+      if patient_wise_only:
+        for CL in CLs:
+          if t < min(self.patient_death_times):
+            survival_probabilities_time_point.append((1, 1))
+            continue
+
+          d2NLLcut = scipy.stats.chi2.ppf(CL, 1).item()
+          def objective_function(
+            expected_probability: float,
+            twoNLL=twoNLL, twoNLL_min=twoNLL_min, d2NLLcut=d2NLLcut
+          ) -> float:
+            return twoNLL(expected_probability) - twoNLL_min - d2NLLcut
+          if best_prob == best_prob_clipped:
+            np.testing.assert_allclose(
+              objective_function(best_prob),
+              -d2NLLcut,
+              atol=1e-2,
+            )
 
-        if patient_wise_only:
           probs = self.possible_probabilities(time_point=t)
-          #Explicitly add best_prob to the probabilities
-          #It should be there already, but sometimes isn't due to numerical issues
           if best_prob not in probs:
             probs = np.append(probs, best_prob)
             probs = np.sort(probs)
@@ -582,7 +588,6 @@ def objective_function(
             err_msg=f"Best probability {best_prob} not found in possible probabilities {probs}",
           )
 
-          # Check edge case: upper bound
           if objective_function(probs[-1]) < 0:
             upper_bound = 1
           else:
@@ -599,7 +604,6 @@ def objective_function(
               raise RuntimeError("No upper sign change found")
             upper_bound = upper
 
-          # Check edge case: lower bound
           if objective_function(probs[0]) < 0:
             lower_bound = 0
           else:
@@ -615,40 +619,66 @@ def objective_function(
             if lower is None:
               raise RuntimeError("No lower sign change found")
             lower_bound = lower
+          survival_probabilities_time_point.append((lower_bound, upper_bound))
+      else:
+        levels = [scipy.stats.chi2.ppf(cl, 1).item() for cl in CLs]
 
-        else:
-          if (
-            best_prob <= self.__endpoint_epsilon
-            or objective_function(self.__endpoint_epsilon) < 0
-          ):
-            lower_bound = 0
-          elif objective_function(best_prob_clipped) >= 0:
-            lower_bound = best_prob_clipped
-          else:
-            lower_bound = scipy.optimize.brentq(
-              objective_function,
-              self.__endpoint_epsilon,
-              best_prob_clipped,
-              xtol=MIPGapAbs,
-              rtol=np.float64(MIPGap),
-            )
-          if (
-            best_prob >= 1 - self.__endpoint_epsilon
-            or objective_function(1 - self.__endpoint_epsilon) < 0
-          ):
-            upper_bound = 1
-          elif objective_function(best_prob_clipped) >= 0:
-            upper_bound = best_prob_clipped
-          else:
-            upper_bound = scipy.optimize.brentq(
-              objective_function,
-              best_prob_clipped,
-              1 - self.__endpoint_epsilon,
-              xtol=MIPGapAbs,
-              rtol=np.float64(MIPGap),
-            )
+        def profile_excess(
+          expected_probability: float,
+          _twoNLL=twoNLL,
+          _twoNLL_min=twoNLL_min,
+        ) -> float:
+          return _twoNLL(expected_probability) - _twoNLL_min
+
+        if best_prob == best_prob_clipped:
+          np.testing.assert_allclose(profile_excess(best_prob), 0.0, atol=1e-2)
+
+        brentq_xtol = 1e-4
+        brentq_rtol = 1e-4
+        eps = self.__endpoint_epsilon
+        one_minus_eps = 1.0 - eps
 
-        survival_probabilities_time_point.append((lower_bound, upper_bound))
+        if best_prob <= eps:
+          lowers = [0.0] * len(CLs)
+        else:
+          lowers = cached_level_crossings(
+            profile_excess,
+            eps,
+            best_prob_clipped,
+            levels,
+            xtol=brentq_xtol,
+            rtol=brentq_rtol,
+          )
+          f_eps = profile_excess(eps)
+          lowers = [
+            0.0 if f_eps <= level else float(x)
+            for x, level in zip(lowers, levels, strict=True)
+          ]
+
+        if best_prob >= one_minus_eps:
+          uppers = [1.0] * len(CLs)
+        else:
+          uppers = cached_level_crossings(
+            profile_excess,
+            one_minus_eps,
+            best_prob_clipped,
+            levels,
+            xtol=brentq_xtol,
+            rtol=brentq_rtol,
+          )
+          f_hi = profile_excess(one_minus_eps)
+          uppers = [
+            1.0 if f_hi <= level else float(x)
+            for x, level in zip(uppers, levels, strict=True)
+          ]
+
+        if print_progress:
+          print(
+            f"  CL crossings unique p evaluations={len(twoNLL.cache)} "
+            f"for time point {t}"
+          )
+        for lower_bound, upper_bound in zip(lowers, uppers, strict=True):
+          survival_probabilities_time_point.append((lower_bound, upper_bound))
     return np.array(best_probabilities), np.array(survival_probabilities)
 
   def plot(self, config: KaplanMeierPlotConfig | None = None, **kwargs) -> dict:
diff --git a/kombine/kaplan_meier_p_value_MINLP.py b/kombine/kaplan_meier_p_value_MINLP.py
index fc7823d..3c8fde3 100644
--- a/kombine/kaplan_meier_p_value_MINLP.py
+++ b/kombine/kaplan_meier_p_value_MINLP.py
@@ -17,7 +17,11 @@
 import scipy.optimize
 import scipy.stats
 
-from .kaplan_meier_MINLP import KaplanMeierPatientNLL, n_choose_d_term_table
+from .kaplan_meier_MINLP import (
+  KaplanMeierPatientNLL,
+  n_choose_d_term_table,
+  km_survival_from_risk_counts,
+)
 from .utilities import (
   LOG_ZERO_EPSILON_DEFAULT,
   GurobiOptimizerMixin,
@@ -54,7 +58,6 @@ def __init__( # pylint: disable=too-many-arguments
     self.__tie_handling = tie_handling
     self.__log_hazard_ratio_bounds = log_hazard_ratio_bounds
     self.__null_hypothesis_constraint = None
-    self.__patient_constraints_for_cox_only = None
     self.__patient_wise_only_constraint = None
     self.__cox_penalty_constraint = None
 
@@ -90,6 +93,13 @@ def parameter_max(self) -> float:
     """
     return self.__parameter_max
 
+  @property
+  def log_zero_epsilon(self) -> float:
+    """
+    Small epsilon used to avoid log(0) in the Gurobi model.
+    """
+    return self.__log_zero_epsilon
+
   @property
   def log_hazard_ratio_bounds(self) -> tuple[float, float]:
     """
@@ -652,9 +662,7 @@ def _extract_patients_per_curve(self, a: gp.tupledict[tuple[int, ...], gp.Var]):
 
     return patients_low, patients_high
 
-  def _extract_curve_statistics( #pylint: disable=too-many-locals
-    self, model: gp.Model, km_probability_at_time_low, km_probability_at_time_high
-  ):
+  def _extract_curve_statistics(self, model: gp.Model):  # pylint: disable=too-many-locals
     """
     Extract statistics for each curve from the optimized model.
 
@@ -662,39 +670,45 @@ def _extract_curve_statistics( #pylint: disable=too-many-locals
         tuple: (n_total_low, n_alive_low, km_prob_low,
                 n_total_high, n_alive_high, km_prob_high)
     """
-    # Extract KM probabilities for each curve
-    km_prob_low = [
-      km_prob.X for _, km_prob in km_probability_at_time_low.items()
-    ]
-    km_prob_high = [
-      km_prob.X for _, km_prob in km_probability_at_time_high.items()
-    ]
-
-    # For n_total and n_alive, we need to look at the r and d variables per curve
-    # These represent at-risk and died counts for each curve
+    n_times = len(self.all_death_times)
+    r_low = []
+    r_high = []
+    s_low = []
+    s_high = []
     n_total_low = 0
     n_alive_low = 0
     n_total_high = 0
     n_alive_high = 0
-
-    for i in range(len(self.all_death_times)):
-      r_low = model.getVarByName(f"r[{i},0]")
-      r_high = model.getVarByName(f"r[{i},1]")
-      d_low = model.getVarByName(f"d[{i},0]")
-      d_high = model.getVarByName(f"d[{i},1]")
-      assert r_low is not None
-      assert r_high is not None
-      assert d_low is not None
-      assert d_high is not None
-
-      if i == 0:  # Use first time point as representative
-        n_total_low = int(np.rint(r_low.X))
-      if i == 0:  # Use first time point as representative
-        n_total_high = int(np.rint(r_high.X))
-
-      n_alive_low = n_total_low - int(np.rint(d_low.X))
-      n_alive_high = n_total_high - int(np.rint(d_high.X))
-
+    for i in range(n_times):
+      r0 = model.getVarByName(f"r[{i},0]")
+      r1 = model.getVarByName(f"r[{i},1]")
+      d0 = model.getVarByName(f"d[{i},0]")
+      d1 = model.getVarByName(f"d[{i},1]")
+      s0 = model.getVarByName(f"n_survived[{i},0]")
+      s1 = model.getVarByName(f"n_survived[{i},1]")
+      assert r0 is not None and r1 is not None
+      assert d0 is not None and d1 is not None
+      assert s0 is not None and s1 is not None
+      r_low.append(r0.X)
+      r_high.append(r1.X)
+      s_low.append(s0.X)
+      s_high.append(s1.X)
+      if i == 0:
+        n_total_low = int(np.rint(r0.X))
+        n_total_high = int(np.rint(r1.X))
+      n_alive_low = n_total_low - int(np.rint(d0.X))
+      n_alive_high = n_total_high - int(np.rint(d1.X))
+
+    km_prob_low = km_survival_from_risk_counts(
+      r_low, s_low,
+      log_zero_epsilon=self.log_zero_epsilon,
+      cumulative=True,
+    )
+    km_prob_high = km_survival_from_risk_counts(
+      r_high, s_high,
+      log_zero_epsilon=self.log_zero_epsilon,
+      cumulative=True,
+    )
     return (
       n_total_low, n_alive_low, km_prob_low,
       n_total_high, n_alive_high, km_prob_high,
@@ -741,15 +755,7 @@ def _make_gurobi_model(self):
     #Binary decision variables: a[j, k] = 1 if patient j is in curve k
     a = model.addVars(self.n_patients, 2, vtype=gp.GRB.BINARY, name="a")
 
-    r, d, n_survived = self.add_counter_variables_and_constraints(model, a)
-
-    # Add Kaplan-Meier probability variables and constraints
-    (
-      km_probability_at_time_low,
-      km_probability_at_time_high
-    ) = self.add_kaplan_meier_probability_variables_and_constraints(
-      model, r, n_survived
-    )
+    r, d, _n_survived = self.add_counter_variables_and_constraints(model, a)
 
     (
       cox_penalty,
@@ -773,8 +779,6 @@ def _make_gurobi_model(self):
       model,
       null_hypothesis_indicator,
       a,
-      km_probability_at_time_low,
-      km_probability_at_time_high,
       beta,
       use_cox_penalty_indicator,
     )
@@ -818,36 +822,20 @@ def update_model_with_cox_only_constraints(
     cox_only: bool,
   ):
     """
-    Update the model with cox_only constraints.
-    If cox_only is True, we add constraints for a[i, j] to be either 0 or 1,
-    based on parameter_in_range.
+    Fix or unfix assignment variables for cox_only mode.
+
+    When cox_only is True, assignments are locked to the nominal in-range
+    groups via bounds so presolve can drop the assignment MIP.
     """
-    # Remove existing constraints if they exist
-    if self.__patient_constraints_for_cox_only is not None:
-      for constr in self.__patient_constraints_for_cox_only:
-        model.remove(constr)
-      self.__patient_constraints_for_cox_only = None
-
-    if cox_only:
-      self.__patient_constraints_for_cox_only = []
-      for i in range(self.n_patients):
-        for j in range(2):  # j=0 for low curve, j=1 for high curve
-          if self.parameter_in_range[i, j]:
-            # The patient must be selected for this curve
-            self.__patient_constraints_for_cox_only.append(
-              model.addConstr(
-                a[i, j] == 1,
-                name=f"patient_{i}_must_be_selected_curve_{j}_cox_only",
-              )
-            )
-          else:
-            # The patient must not be selected for this curve
-            self.__patient_constraints_for_cox_only.append(
-              model.addConstr(
-                a[i, j] == 0,
-                name=f"patient_{i}_must_not_be_selected_curve_{j}_cox_only",
-              )
-            )
+    for i in range(self.n_patients):
+      for j in range(2):
+        if cox_only:
+          value = 1.0 if self.parameter_in_range[i, j] else 0.0
+          a[i, j].LB = value
+          a[i, j].UB = value
+        else:
+          a[i, j].LB = 0.0
+          a[i, j].UB = 1.0
 
     model.update()
 
@@ -939,7 +927,7 @@ def update_model_with_patient_wise_only_constraint( #pylint: disable=too-many-ar
 
     model.update()
 
-  def solve_and_pvalue( # pylint: disable=too-many-locals, too-many-arguments
+  def solve_and_pvalue( # pylint: disable=too-many-locals, too-many-arguments, too-many-branches, too-many-statements
     self,
     *,
     cox_only: bool = False,
@@ -987,8 +975,6 @@ def solve_and_pvalue( # pylint: disable=too-many-locals, too-many-arguments
       model,
       null_hypothesis_indicator,
       a,
-      km_probability_at_time_low,
-      km_probability_at_time_high,
       beta,
       use_cox_penalty_indicator,
     ) = self.gurobi_model
@@ -1033,9 +1019,7 @@ def solve_and_pvalue( # pylint: disable=too-many-locals, too-many-arguments
     # Extract curve statistics for null hypothesis
     (n_total_low_null, n_alive_low_null, km_prob_low_null,
      n_total_high_null, n_alive_high_null, km_prob_high_null) = (
-      self._extract_curve_statistics(
-        model, km_probability_at_time_low, km_probability_at_time_high
-      )
+      self._extract_curve_statistics(model)
     )
 
     # Extract hazard ratio for null hypothesis (should be 1.0)
@@ -1061,9 +1045,25 @@ def solve_and_pvalue( # pylint: disable=too-many-locals, too-many-arguments
       model=model,
     )
 
+    start_a = None
+    start_beta = None
+    if not cox_only:
+      start_a = {
+        (i, j): float(a[i, j].X)
+        for i in range(self.n_patients)
+        for j in range(2)
+      }
+      if beta is not None:
+        start_beta = float(beta.X)
+
     if print_progress or verbose:
       print("Solving for alternative hypothesis...")
     self.update_model_for_null_hypothesis_or_not(model, null_hypothesis_indicator, False)
+    if start_a is not None:
+      for (i, j), value in start_a.items():
+        a[i, j].Start = value
+      if start_beta is not None:
+        beta.Start = start_beta
     # pylint: disable=duplicate-code
     model = self._setup_and_optimize(
       model,
@@ -1088,9 +1088,7 @@ def solve_and_pvalue( # pylint: disable=too-many-locals, too-many-arguments
     # Extract curve statistics for alternative hypothesis
     (n_total_low_alt, n_alive_low_alt, km_prob_low_alt,
      n_total_high_alt, n_alive_high_alt, km_prob_high_alt) = (
-      self._extract_curve_statistics(
-        model, km_probability_at_time_low, km_probability_at_time_high
-      )
+      self._extract_curve_statistics(model)
     )
 
     # Extract hazard ratio for alternative hypothesis (can be any value)
diff --git a/kombine/utilities.py b/kombine/utilities.py
index eb33023..9d65bb0 100644
--- a/kombine/utilities.py
+++ b/kombine/utilities.py
@@ -295,7 +295,7 @@ def _optimize_with_fallbacks(
   ):
     """
     Attempts to optimize the Gurobi model, applying fallback strategies
-    if the initial optimization is suboptimal.
+    if the initial optimization is suboptimal or hits the time limit.
 
     Args:
         model: The Gurobi model to optimize.
@@ -318,11 +318,14 @@ def _optimize_with_fallbacks(
       print("Attempting initial optimization...")
     model.optimize()
 
-    # Check for suboptimal status and apply fallbacks
-    if model.status == GRB.SUBOPTIMAL:
+    needs_fallback = model.status in (GRB.SUBOPTIMAL, GRB.TIME_LIMIT)
+    if needs_fallback:
       for i, (fallback_params, description) in enumerate(fallback_strategies):
         if verbose:
-          print(f"Model returned suboptimal solution. Applying fallback {i+1}: {description}")
+          print(
+            f"Model returned status {model.status}. "
+            f"Applying fallback {i+1}: {description}"
+          )
           print(f"  New parameters: {fallback_params}")
         self._set_gurobi_params(model, fallback_params)
         model.optimize()
diff --git a/test/kombine/test_discrete_optimization.py b/test/kombine/test_discrete_optimization.py
index 46ed9c9..af24444 100644
--- a/test/kombine/test_discrete_optimization.py
+++ b/test/kombine/test_discrete_optimization.py
@@ -6,10 +6,12 @@
 import unittest
 
 import numpy as np
+import scipy.optimize
 
 from kombine.discrete_optimization import (
   minimize_discrete_single_minimum,
   binary_search_sign_change,
+  cached_level_crossings,
 )
 
 class TestMinimizeDiscreteSingleMinimum(unittest.TestCase):
@@ -439,5 +441,103 @@ def func(x):
     self.assertLessEqual(result, 0.2)
 
 
+class TestCachedLevelCrossings(unittest.TestCase):
+  """Shared-evaluation multi-level brentq for a 1D profile."""
+
+  @staticmethod
+  def _quadratic(x):
+    return (x - 0.5) ** 2
+
+  def test_quadratic_both_sides_match_analytic(self):
+    """Crossings of (x-0.5)^2 match the analytic sqrt(level) points."""
+    inner = 0.01
+    outer = 0.04
+    levels = [inner, outer]
+    left = cached_level_crossings(self._quadratic, 0.0, 0.5, levels, xtol=1e-8, rtol=1e-8)
+    right = cached_level_crossings(self._quadratic, 1.0, 0.5, levels, xtol=1e-8, rtol=1e-8)
+    np.testing.assert_allclose(left, [0.5 - inner ** 0.5, 0.5 - outer ** 0.5], atol=1e-7)
+    np.testing.assert_allclose(right, [0.5 + inner ** 0.5, 0.5 + outer ** 0.5], atol=1e-7)
+
+  def test_level_order_does_not_change_results(self):
+    """Returned roots follow the requested level order."""
+    inner = 0.01
+    outer = 0.04
+    forward = cached_level_crossings(
+      self._quadratic, 0.0, 0.5, [inner, outer], xtol=1e-8, rtol=1e-8
+    )
+    backward = cached_level_crossings(
+      self._quadratic, 0.0, 0.5, [outer, inner], xtol=1e-8, rtol=1e-8
+    )
+    np.testing.assert_allclose(forward, [backward[1], backward[0]], atol=1e-10)
+
+  def test_fewer_evals_than_independent_brentq(self):
+    """Shared cache uses fewer func evals than independent brentq per level."""
+    inner = 0.01
+    outer = 0.04
+    levels = [inner, outer]
+    xtol = 1e-8
+    rtol = 1e-8
+
+    independent = []
+    def counting_independent(x):
+      independent.append(x)
+      return self._quadratic(x)
+
+    for level in levels:
+      scipy.optimize.brentq(
+        lambda x, _level=level: counting_independent(x) - _level,
+        0.0, 0.5, xtol=xtol, rtol=np.float64(rtol),
+      )
+      scipy.optimize.brentq(
+        lambda x, _level=level: counting_independent(x) - _level,
+        0.5, 1.0, xtol=xtol, rtol=np.float64(rtol),
+      )
+
+    shared = []
+    def counting_shared(x):
+      shared.append(x)
+      return self._quadratic(x)
+
+    cached_level_crossings(counting_shared, 0.0, 0.5, levels, xtol=xtol, rtol=rtol)
+    cached_level_crossings(counting_shared, 1.0, 0.5, levels, xtol=xtol, rtol=rtol)
+
+    self.assertLess(len(shared), len(independent))
+    # Early probes on the full [0, 0.5] interval must be reused for the outer cut.
+    self.assertTrue(any(x < 0.5 - inner ** 0.5 - 1e-6 for x in shared))
+
+  def test_outer_already_inside_returns_boundary(self):
+    """If the outer endpoint is already inside the cut, return that endpoint."""
+    # f(0) = 0.25, so a cut of 0.4 has no crossing on [0, 0.5].
+    result = cached_level_crossings(self._quadratic, 0.0, 0.5, [0.4], xtol=1e-6, rtol=1e-6)
+    self.assertEqual(result, [0.0])
+
+  def test_innermost_crossing_when_profile_wiggles(self):
+    """A disconnected outer valley must not steal the MLE-connected crossing."""
+    knots = np.array([0.0, 0.2, 0.3, 0.4, 0.5])
+    values = np.array([0.2, 0.05, 0.02, 0.06, 0.0])
+
+    def wiggle(x):
+      return float(np.interp(x, knots, values))
+
+    # From the MLE at 0.5, 0.04 is first exceeded on (0.4, 0.5).
+    result = cached_level_crossings(wiggle, 0.0, 0.5, [0.04], xtol=1e-8, rtol=1e-8)
+    expected = 0.4 + 0.1 * (0.06 - 0.04) / 0.06
+    np.testing.assert_allclose(result, [expected], atol=1e-6)
+
+  def test_default_tolerances_stop_before_scipy_defaults(self):
+    """Default 1e-4 xtol/rtol stop with fewer evals than scipy's tight defaults."""
+    tight = []
+    loose = []
+    def f_tight(x):
+      tight.append(x)
+      return self._quadratic(x)
+    def f_loose(x):
+      loose.append(x)
+      return self._quadratic(x)
+    cached_level_crossings(f_tight, 0.0, 0.5, [0.04], xtol=2e-12, rtol=2e-15)
+    cached_level_crossings(f_loose, 0.0, 0.5, [0.04], xtol=1e-4, rtol=1e-4)
+    self.assertLess(len(loose), len(tight))
+
+
 if __name__ == "__main__":
   unittest.main()
diff --git a/test/kombine/test_km_likelihood.py b/test/kombine/test_km_likelihood.py
index a71f12d..3ef18c0 100644
--- a/test/kombine/test_km_likelihood.py
+++ b/test/kombine/test_km_likelihood.py
@@ -3,6 +3,7 @@
 """
 
 import argparse
+import itertools
 import pathlib
 import json
 import math
@@ -12,6 +13,7 @@
 
 import kombine.datacard
 from kombine.kaplan_meier_MINLP import KaplanMeierPatientNLL, MINLPForKM
+from kombine.kaplan_meier_likelihood import KaplanMeierLikelihood
 from ..utility_testing_functions import format_value_for_json, Tolerance
 
 warnings.simplefilter("error")
@@ -19,6 +21,30 @@
 here = pathlib.Path(__file__).parent
 datacards = here / "datacards" / "simple_examples"
 
+
+def _nll_patient(time, censored=False):
+  return KaplanMeierPatientNLL(
+    time=time,
+    censored=censored,
+    parameter_nll=lambda x: 0,
+    observed_parameter=0.5,
+  )
+
+
+def _minlp_for_rd(patients, time_point, *, collapse):
+  return MINLPForKM(
+    patients,
+    parameter_min=0.1,
+    parameter_max=0.9,
+    time_point=time_point,
+    collapse_consecutive_deaths=collapse,
+  )
+
+
+def _rd_feasible(minlp, i, r_value, d_value):
+  return minlp._rd_pair_is_feasible(i, r_value, d_value)  # pylint: disable=protected-access
+
+
 def runtest(
   censoring=False,
 ): #pylint: disable=too-many-locals, too-many-statements, too-many-branches
@@ -446,6 +472,143 @@ def runtest(
       )
     raise
 
+def test_rd_pair_is_feasible_per_time_bounds():
+  """Per-time (r, d) bounds: r and d caps, d ≤ r, and the empty assignment."""
+  patients = [
+    _nll_patient(1.0),
+    _nll_patient(1.0),
+    _nll_patient(1.5, censored=True),
+    _nll_patient(2.0),
+    _nll_patient(2.0),
+    _nll_patient(3.0),
+  ]
+  minlp = _minlp_for_rd(patients, 3.0, collapse=False)
+  np.testing.assert_allclose(minlp.times_to_consider, [1.0, 2.0, 3.0])
+  np.testing.assert_array_equal(minlp.n_at_risk_max, [6, 3, 1])
+  np.testing.assert_array_equal(minlp.n_died_max, [2, 2, 1])
+
+  assert _rd_feasible(minlp, 0, 0, 0)
+  assert _rd_feasible(minlp, 0, 6, 2)
+  assert not _rd_feasible(minlp, 0, 7, 0)
+  assert not _rd_feasible(minlp, 0, 6, 3)
+  assert not _rd_feasible(minlp, 0, 2, 3)
+  assert not _rd_feasible(minlp, 2, 2, 0)
+
+
+def test_rd_pair_is_feasible_flow_and_censoring_slack():
+  """Survivors cannot exceed the next risk-set max plus intervening censoring."""
+  patients = [
+    _nll_patient(1.0),
+    _nll_patient(1.0),
+    _nll_patient(1.5, censored=True),
+    _nll_patient(2.0),
+    _nll_patient(2.0),
+    _nll_patient(3.0),
+  ]
+  minlp = _minlp_for_rd(patients, 3.0, collapse=False)
+  np.testing.assert_array_equal(minlp.n_censored_between_times_max, [1, 0])
+
+  # i=0: s <= n_at_risk_max[1] + cmax = 3 + 1
+  assert _rd_feasible(minlp, 0, 6, 2)  # s=4 on the boundary
+  assert not _rd_feasible(minlp, 0, 6, 1)  # s=5
+  assert not _rd_feasible(minlp, 0, 6, 0)  # s=6
+  # cmax slack: s can exceed the next risk-set max, but not by more than cmax
+  assert _rd_feasible(minlp, 0, 4, 0)  # s=4 == 3 + 1
+  assert not _rd_feasible(minlp, 0, 5, 0)  # s=5
+
+  # i=1: no censoring in [2, 3), so s <= n_at_risk_max[2] = 1
+  assert _rd_feasible(minlp, 1, 3, 2)
+  assert _rd_feasible(minlp, 1, 1, 0)
+  assert not _rd_feasible(minlp, 1, 3, 1)
+  assert not _rd_feasible(minlp, 1, 3, 0)
+
+
+def test_rd_pair_is_feasible_last_time_and_single_group():
+  """No next-interval s-rule at the last time, or when collapse leaves one group."""
+  patients = [
+    _nll_patient(1.0),
+    _nll_patient(1.0),
+    _nll_patient(1.5, censored=True),
+    _nll_patient(2.0),
+    _nll_patient(2.0),
+    _nll_patient(3.0),
+  ]
+  minlp = _minlp_for_rd(patients, 3.0, collapse=False)
+  # Last time: s-rule does not apply, so (r=1, d=0) is allowed even though
+  # there is no later risk set.
+  assert _rd_feasible(minlp, 2, 1, 0)
+  assert _rd_feasible(minlp, 2, 1, 1)
+  assert not _rd_feasible(minlp, 2, 1, 2)
+
+  collapsed = _minlp_for_rd(
+    [_nll_patient(float(t)) for t in (1, 2, 3, 4)],
+    4.0,
+    collapse=True,
+  )
+  assert collapsed.n_times_to_consider == 1
+  assert collapsed.n_censored_between_times_max.size == 0
+  # Would be flow-infeasible if a later empty risk set existed.
+  assert _rd_feasible(collapsed, 0, 4, 0)
+  assert _rd_feasible(collapsed, 0, 4, 4)
+  assert not _rd_feasible(collapsed, 0, 5, 0)
+
+
+def test_rd_pair_is_feasible_previous_risk_injected():
+  """r cannot exceed the previous risk-set maximum (redundant on nested risk sets)."""
+  minlp = _minlp_for_rd(
+    [_nll_patient(1.0), _nll_patient(2.0)],
+    2.0,
+    collapse=False,
+  )
+  assert minlp.n_times_to_consider == 2
+  minlp.__dict__["n_at_risk_max"] = np.array([2, 5])
+  minlp.__dict__["n_died_max"] = np.array([2, 2])
+  minlp.__dict__["n_censored_between_times_max"] = np.array([0])
+
+  assert not _rd_feasible(minlp, 1, 3, 0)
+  assert _rd_feasible(minlp, 1, 2, 0)
+  # First time has no previous-r rule.
+  assert _rd_feasible(minlp, 0, 2, 0)
+
+
+def test_rd_pair_is_feasible_realizable_assignments_kept():
+  """Every 0/1 assignment's (r, d) stays feasible on the simple N=12 cards."""
+  cases = (
+    ("poisson_ratio_km.txt", 8.0),
+    ("poisson_ratio_km.txt", 8.8),
+    ("poisson_ratio_km_censoring.txt", 8.8),
+  )
+  for dcname, time_point in cases:
+    datacard = kombine.datacard.Datacard.parse_datacard(datacards / dcname)
+    kml = datacard.km_likelihood(
+      parameter_min=0.19,
+      parameter_max=0.79,
+      endpoint_epsilon=1e-4,
+    )
+    for collapse in (True, False):
+      minlp = MINLPForKM(
+        kml.all_patients,
+        parameter_min=0.19,
+        parameter_max=0.79,
+        time_point=time_point,
+        endpoint_epsilon=1e-4,
+        collapse_consecutive_deaths=collapse,
+      )
+      n_at = [minlp.patient_still_at_risk(t) for t in minlp.times_to_consider]
+      died = [minlp.patient_died(t) for t in minlp.times_to_consider]
+      n_patients = minlp.n_patients
+      assert n_patients <= 16
+      for bits in itertools.product((False, True), repeat=n_patients):
+        selected = np.array(bits, dtype=bool)
+        for i in range(minlp.n_times_to_consider):
+          r_value = int(np.count_nonzero(n_at[i] & selected))
+          d_value = int(np.count_nonzero(died[i] & selected))
+          assert _rd_feasible(minlp, i, r_value, d_value), (
+            f"{dcname} t={time_point} collapse={collapse} i={i} "
+            f"(r, d)=({r_value}, {d_value})"
+          )
+
+
 def test_times_to_consider_collapse_logic(): # pylint: disable=too-many-locals
   """
   Comprehensive tests for the times_to_consider collapsing logic.
@@ -644,7 +807,8 @@ def test_collapse_equivalence_with_tied_death_and_censoring():
       ci_collapse[finite_ci],
       ci_no_collapse[finite_ci],
       rtol=0.0,
-      atol=3e-5,
+      # CI edges come from brentq over MINLP 2NLL; allow ~MIPGap (1e-4) noise.
+      atol=1e-4,
     )
 
 def test_collapse_equivalence_nonfixed_observable():
@@ -703,7 +867,8 @@ def test_collapse_equivalence_nonfixed_observable():
       ci_collapse[finite_ci],
       ci_no_collapse[finite_ci],
       rtol=0.0,
-      atol=3e-5,
+      # CI edges come from brentq over MINLP 2NLL; allow ~MIPGap (1e-4) noise.
+      atol=1e-4,
       err_msg=f"CI mismatch for mode {mode_kwargs}",
     )
 
@@ -730,11 +895,37 @@ def main(args=None):
   g.add_argument(
     "--collapse-logic-only",
     action="store_true",
-    help="Test only the collapse logic.",
+    help="Test only the collapse logic (short; ~2 min).",
+  )
+  p.add_argument(
+    "--print-progress",
+    action="store_true",
+    help="Print progress during MINLP / likelihood walks.",
+  )
+  p.add_argument(
+    "--gurobi-verbose",
+    action="store_true",
+    help="Enable Gurobi OutputFlag logging during solves.",
   )
   args = p.parse_args(args)
 
+  if args.print_progress or args.gurobi_verbose:
+    _orig = KaplanMeierLikelihood.survival_probabilities_likelihood
+    def _wrapped(self, *a, **kw):
+      if args.print_progress:
+        kw["print_progress"] = True
+      if args.gurobi_verbose:
+        kw["gurobi_verbose"] = True
+        kw["optimize_verbose"] = True
+      return _orig(self, *a, **kw)
+    KaplanMeierLikelihood.survival_probabilities_likelihood = _wrapped
+
   if args.collapse_logic_only:
+    test_rd_pair_is_feasible_per_time_bounds()
+    test_rd_pair_is_feasible_flow_and_censoring_slack()
+    test_rd_pair_is_feasible_last_time_and_single_group()
+    test_rd_pair_is_feasible_previous_risk_injected()
+    test_rd_pair_is_feasible_realizable_assignments_kept()
     test_times_to_consider_collapse_logic()
     test_collapse_equivalence_with_tied_death_and_censoring()
     test_collapse_equivalence_nonfixed_observable()
@@ -749,6 +940,11 @@ def main(args=None):
     runtest(censoring=True)
 
   # Run additional tests
+  test_rd_pair_is_feasible_per_time_bounds()
+  test_rd_pair_is_feasible_flow_and_censoring_slack()
+  test_rd_pair_is_feasible_last_time_and_single_group()
+  test_rd_pair_is_feasible_previous_risk_injected()
+  test_rd_pair_is_feasible_realizable_assignments_kept()
   test_times_to_consider_collapse_logic()
   test_collapse_equivalence_with_tied_death_and_censoring()
   test_collapse_equivalence_nonfixed_observable()