From 104760a5b4809b0f6791f42f8127ddfc204dfc18 Mon Sep 17 00:00:00 2001 From: Phil Date: Wed, 5 Aug 2026 17:53:57 -0400 Subject: [PATCH] fix: price knock-ins by their own cost, and stop rejecting negative budgets Three defects on the same path. Each is invisible on its own, because a knockout-only problem never touches any of them. 1. filter_sd_maxcost priced an intervention by dict membership, knockout dict first, while SDProblem prices it knock-in over knockout. A reaction present in both dicts was therefore charged its knockout cost after the MILP had already priced it as a knock-in, and the design was discarded. The overlap is the normal case rather than an edge case: compute_strain_designs defaults ko_cost to every reaction whenever the caller passes only ki_cost. This one is not specific to negative costs. With ko_cost={'R2': 50}, ki_cost={all: 1} and max_cost=3 on a four-reaction network, the MILP emits {R1,R2,R4} at cost 3 and the filter re-prices it at 52, so the design is lost silently. Gene-based problems escape it because ko_cost starts empty there. Dispatch on the sign of the design value instead, which is the value semantics the function already documents. Where ki_cost is empty the new expression is identical to the old term for term, so knockout behaviour cannot change. The same expression zipped the filtered designs against the unfiltered cost list when stamping '**cost**', misaligning labels and the resulting sort whenever the filter removed anything. Fixed alongside, since correct stamping is required once anything is filtered. 2. cont_MILP kept the three leading rows that constrain z alone. Stripped of their z columns they read 0 <= max_cost, which is vacuous for a non-negative budget but makes verify_sd reject every design once max_cost is negative. Neutralise their right-hand sides in the continuous copy rather than dropping the rows, whose positions z_map_constr_ineq columns are aligned to. 3. Row 0 asserted total cost >= 0. That is implied whenever costs are non-negative, and z is binary so the objective is bounded without it. Its only effect was to make negative costs infeasible by construction. Left open. Together these restore negative intervention costs as a way of forcing a reaction into every solution: give it a large negative cost and set max_cost so that omitting it is unaffordable, while the budget still bounds the rest. On a network with elementary flux vectors {R1,R2,R4} and {R1,R3,R4}, cost(R2)=-100 with max_cost=-97 now returns {R1,R2,R4} alone. e_coli_core MCS (growth >= 0.001, max_cost 3) returns the identical 353 designs, compared as sets. Test suite: 368 passed, 2 skipped. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01RFtof9nZXFvCNoXz19po8C --- straindesign/networktools.py | 17 ++++++++++++++--- straindesign/strainDesignProblem.py | 26 +++++++++++++++++++++----- 2 files changed, 35 insertions(+), 8 deletions(-) diff --git a/straindesign/networktools.py b/straindesign/networktools.py index 0c5027e..59342c5 100644 --- a/straindesign/networktools.py +++ b/straindesign/networktools.py @@ -1361,11 +1361,22 @@ def filter_sd_maxcost(sd, max_cost, kocost, kicost): # introduced KIs and KOs carry values of +1.0 and -1.0 respectively # non-made KIs are marked by 0.0 and non-made KOs don't appear. # We count costs of interventions made, which are marked by v != 0. + # A reaction may occur in both cost dicts (kocost defaults to all reactions). + # The value's sign says which kind of intervention was made, so it selects + # the cost dict, matching how SDProblem prices interventions. + def itv_cost(k, v): + if v == 0: + return 0 + if v > 0: # knock-in + return kicost[k] if k in kicost else kocost.get(k, 0) + return kocost[k] if k in kocost else kicost.get(k, 0) # knock-out + if max_cost: - costs = [np.sum([(kocost[k] if k in kocost else kicost.get(k, 0)) if v != 0 else 0 for k, v in m.items()]) for m in sd] - sd = [sd[i] for i in range(len(sd)) if costs[i] <= max_cost + 1e-8] + costs = [np.sum([itv_cost(k, v) for k, v in m.items()]) for m in sd] + keep = [i for i in range(len(sd)) if costs[i] <= max_cost + 1e-8] + sd = [sd[i] for i in keep] # sort strain designs by intervention costs - [s.update({'**cost**': c}) for s, c in zip(sd, costs)] + [sd[j].update({'**cost**': costs[i]}) for j, i in enumerate(keep)] sd.sort(key=lambda x: x.pop('**cost**')) return sd diff --git a/straindesign/strainDesignProblem.py b/straindesign/strainDesignProblem.py index 626a581..44a0738 100644 --- a/straindesign/strainDesignProblem.py +++ b/straindesign/strainDesignProblem.py @@ -154,19 +154,26 @@ def __init__(self, model: Model, sd_modules: List[SDModule], *args, **kwargs): for i in [i for i, x in enumerate(self.cost) if np.isnan(x)]: self.cost[i] = 0.0 # Top 3 fixed rows of A_ineq (with the b_ineq values assembled just below): - # row 0 (idx_row_maxcost): -cost . z <= 0 -> total intervention cost >= 0 + # row 0 (idx_row_maxcost): -cost . z <= inf -> left open, see below # row 1 (idx_row_mincost): cost . z <= max_cost -> total intervention cost <= max_cost (budget cap) # row 2 (idx_row_obj): objective placeholder row (set later via fixObjective) # NOTE: the maxcost/mincost names are historical and read counter-intuitively -- row 0 is the - # >= 0 lower bracket, row 1 is the <= max_cost upper cap. + # lower bracket on total cost, row 1 is the <= max_cost upper cap. + # Row 0 used to read `total cost >= 0`. That is implied whenever every cost is + # non-negative, and z is binary, so the objective is bounded without it. Its only + # effect was on negative intervention costs, which are a documented way of forcing + # a reaction into every solution: give it a large negative cost and set max_cost + # accordingly, so any solution omitting it is unaffordable. The bracket made that + # infeasible by construction. The row is kept so that row indices -- and the + # z_map_constr_ineq columns aligned to them -- stay put. self.idx_row_maxcost = 0 self.idx_row_mincost = 1 self.idx_row_obj = 2 self.A_ineq = sparse.csr_matrix([[-i for i in self.cost], self.cost, [0 for _ in range(self.num_z)]]) if self.max_cost is None: - self.b_ineq = [0.0, float(np.sum(np.abs(self.cost))), np.inf] + self.b_ineq = [np.inf, float(np.sum(np.abs(self.cost))), np.inf] else: - self.b_ineq = [0.0, float(self.max_cost), np.inf] + self.b_ineq = [np.inf, float(self.max_cost), np.inf] self.z_map_constr_ineq = sparse.csc_matrix((numr, 3)) self.lb = [1.0 if r.id in self.essential_kis else 0.0 for r in model.reactions] self.ub = [1.0 - float(i) for i in self.z_non_targetable] @@ -199,7 +206,16 @@ def __init__(self, model: Model, sd_modules: List[SDModule], *args, **kwargs): # Save continous part of MILP for easy strain design validation cont_vars = [i for i in range(0, self.A_ineq.shape[1]) if not i in self.idx_z] - self.cont_MILP = ContMILP(self.A_ineq[:, cont_vars], self.b_ineq.copy(), self.A_eq[:, cont_vars], self.b_eq.copy(), + # The three leading rows constrain z alone (both cost brackets and the objective + # placeholder). Once the z columns are dropped they read 0 <= rhs, which says + # nothing about the continuous system but turns every verification infeasible as + # soon as a cost is negative enough to make max_cost negative. Neutralise their + # right-hand sides here rather than dropping the rows, since z_map_constr_ineq + # columns are indexed by row position. + cont_b_ineq = self.b_ineq.copy() + for r in (self.idx_row_maxcost, self.idx_row_mincost, self.idx_row_obj): + cont_b_ineq[r] = np.inf + self.cont_MILP = ContMILP(self.A_ineq[:, cont_vars], cont_b_ineq, self.A_eq[:, cont_vars], self.b_eq.copy(), [self.lb[i] for i in cont_vars], [self.ub[i] for i in cont_vars], [self.c[i] for i in cont_vars], self.z_map_constr_ineq.tocoo(), self.z_map_constr_eq.tocoo(), self.z_map_vars[:, cont_vars].tocoo())