From f723cf6bf90d50d980cb0cb8e0d09d3face99b5f Mon Sep 17 00:00:00 2001 From: Phil Date: Wed, 5 Aug 2026 20:02:27 -0400 Subject: [PATCH 1/5] fix: bill knock-ins once, and stop them doubling as knockout candidates Two remaining consequences of a reaction sitting in both cost dictionaries at once. Neither can arise for a knockout-only problem, which is why the existing tests do not reach them. SDSolutions._compute_costs_and_bounds walked each cost dictionary in turn and added whatever it found, so a reaction present in both was billed twice. With ki_cost given for four reactions and ko_cost defaulting to all of them, a three-intervention design was reported at cost 6. It now selects the dictionary by the sign of the entry -- positive for an addition, negative for a removal -- which is how SDProblem prices interventions and how filter_sd_maxcost selects a cost. Design sets are untouched; only reported costs change, and only where the dictionaries overlapped. compute_strain_designs defaulted ko_cost to every reaction in the model whenever the caller gave no knockout costs, including reactions the caller had named as knock-in candidates. Knock-ins override knockouts in the MILP, so those entries described a state that could never be solved. Worse, with every reaction named addable, a SUPPRESS problem still answered with knockouts the caller had not asked for. The default now skips the named knock-ins. Reactions the caller did not name stay knockable, so mixed problems -- "R2 may be added, anything else may be removed" -- behave exactly as before. Four regression tests cover the pair, over every installed solver: elementary flux vector enumeration through a PROTECT module with all reactions addable, a negative cost forcing a reaction into every design, single billing of an intervention, and knock-in candidates not being offered as knockouts. Twelve of the sixteen fail without these changes. e_coli_core MCS (growth >= 0.001, max_cost 3) returns the identical 353 designs. Test suite: 372 passed, 2 skipped. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01RFtof9nZXFvCNoXz19po8C --- straindesign/compute_strain_designs.py | 9 ++- straindesign/strainDesignSolutions.py | 25 ++++--- tests/test_05_straindesign.py | 93 ++++++++++++++++++++++++++ 3 files changed, 115 insertions(+), 12 deletions(-) diff --git a/straindesign/compute_strain_designs.py b/straindesign/compute_strain_designs.py index 217e6bf..33ef1f2 100644 --- a/straindesign/compute_strain_designs.py +++ b/straindesign/compute_strain_designs.py @@ -473,8 +473,15 @@ def compute_strain_designs(model: Model, **kwargs: dict) -> SDSolutions: else: kwargs['gene_kos'] = False has_gene_names = False + # With no knockout costs given, every reaction is a knockout candidate -- except those + # the caller named as knock-in candidates. Knock-ins override knockouts in the MILP + # (SDProblem masks ko_cost wherever ki_cost is set), so listing such a reaction in both + # dicts described a state that could never be solved, and left the downstream cost + # lookups to disambiguate it. Reactions the caller did not name stay knockable, which + # keeps mixed problems -- "R2 may be added, anything else may be removed" -- intact. if KOCOST not in kwargs and not kwargs['gene_kos']: - uncmp_ko_cost = {k: 1.0 for k in model.reactions.list_attr('id')} + named_ki = set(kwargs.get(KICOST) or {}) + uncmp_ko_cost = {k: 1.0 for k in model.reactions.list_attr('id') if k not in named_ki} elif KOCOST not in kwargs or not kwargs[KOCOST]: uncmp_ko_cost = {} if KICOST not in kwargs or not kwargs[KICOST]: diff --git a/straindesign/strainDesignSolutions.py b/straindesign/strainDesignSolutions.py index 2889f27..a9fd016 100644 --- a/straindesign/strainDesignSolutions.py +++ b/straindesign/strainDesignSolutions.py @@ -214,18 +214,21 @@ def _compute_costs_and_bounds(cost_sd, reaction_sd, model, sd_setup): Returns: (sd_cost, itv_bounds, has_complex_regul_itv) """ - # compute intervention costs + # compute intervention costs. + # A reaction may appear in both the knockout and the knock-in dict, since ko_cost + # defaults to every reaction whenever only ki_cost is given. Charging each dict + # separately would then bill one intervention twice, so the sign of the entry + # picks the dict: positive is an addition, negative a removal. This matches how + # SDProblem prices interventions and how filter_sd_maxcost selects the cost. sd_cost = [0 for _ in range(len(cost_sd))] - if KOCOST in sd_setup: - for k, v in sd_setup[KOCOST].items(): - for i, s in enumerate(cost_sd): - if k in s and s[k] != 0: - sd_cost[i] += float(v) - if KICOST in sd_setup: - for k, v in sd_setup[KICOST].items(): - for i, s in enumerate(cost_sd): - if k in s and s[k] != 0: - sd_cost[i] += float(v) + ko_c, ki_c = sd_setup.get(KOCOST, {}), sd_setup.get(KICOST, {}) + for i, s in enumerate(cost_sd): + for k, v in s.items(): + if v == 0: + continue + src = (ki_c if k in ki_c else ko_c) if v > 0 else (ko_c if k in ko_c else ki_c) + if k in src: + sd_cost[i] += float(src[k]) if GKOCOST in sd_setup: for k, v in sd_setup[GKOCOST].items(): for i, s in enumerate(cost_sd): diff --git a/tests/test_05_straindesign.py b/tests/test_05_straindesign.py index f24ea1f..b0eb474 100644 --- a/tests/test_05_straindesign.py +++ b/tests/test_05_straindesign.py @@ -427,3 +427,96 @@ def test_lazy_expansion(model_small_example): pass finally: csd.LAZY_EXPANSION_THRESHOLD = orig_threshold + + +def _two_route_network(): + """R1: -> A, R2/R3: A -> B, R4: B ->. + + Two elementary flux vectors, {R1,R2,R4} and {R1,R3,R4}; exactly one uses R2. + Small enough that the expected answer can be written down rather than computed. + """ + from cobra import Model, Reaction, Metabolite + m = Model('two_route') + A, B = Metabolite('A'), Metabolite('B') + for rid, stoich in [('R1', {A: 1.0}), ('R2', {A: -1.0, B: 1.0}), + ('R3', {A: -1.0, B: 1.0}), ('R4', {B: -1.0})]: + r = Reaction(rid) + r.lower_bound, r.upper_bound = 0.0, 10.0 + m.add_reactions([r]) + r.add_metabolites(stoich) + m.objective = 'R4' + return m + + +def _designs(solution): + return sorted(sorted(k for k, v in d.items() if v != 0) for d in solution.get_reaction_sd()) + + +@pytest.mark.timeout(30) +def test_efv_enumeration_via_knockins(curr_solver): + """A PROTECT module with every reaction addable enumerates minimal supports. + + A minimal set of reactions whose presence keeps the region inhabited is an + elementary flux vector, so both routes must come back and neither may be padded. + """ + model = _two_route_network() + sol = sd.compute_strain_designs(model, + sd_modules=[sd.SDModule(model, PROTECT, constraints=['R4 >= 1'])], + max_cost=3, + ki_cost={r: 1.0 for r in ['R1', 'R2', 'R3', 'R4']}, + solution_approach='populate', + solver=curr_solver, + compress=False) + assert _designs(sol) == [['R1', 'R2', 'R4'], ['R1', 'R3', 'R4']] + + +@pytest.mark.timeout(30) +def test_negative_cost_forces_reaction_into_every_design(curr_solver): + """A large negative cost forces a reaction in while the budget still bounds the rest. + + Omitting R2 leaves only positive costs, which cannot reach a negative budget, so + the route through R3 is excluded and only the R2 route survives. + """ + model = _two_route_network() + sol = sd.compute_strain_designs(model, + sd_modules=[sd.SDModule(model, PROTECT, constraints=['R4 >= 1'])], + max_cost=-97, + ki_cost={'R1': 1.0, 'R2': -100.0, 'R3': 1.0, 'R4': 1.0}, + solution_approach='populate', + solver=curr_solver, + compress=False) + assert _designs(sol) == [['R1', 'R2', 'R4']] + assert list(sol.sd_cost) == [-98.0] + + +@pytest.mark.timeout(30) +def test_intervention_costs_are_not_counted_twice(curr_solver): + """A reaction named as a knock-in is billed once, not once per cost dictionary.""" + model = _two_route_network() + sol = sd.compute_strain_designs(model, + sd_modules=[sd.SDModule(model, PROTECT, constraints=['R4 >= 1'])], + max_cost=3, + ki_cost={r: 1.0 for r in ['R1', 'R2', 'R3', 'R4']}, + solution_approach='populate', + solver=curr_solver, + compress=False) + for design, cost in zip(sol.get_reaction_sd(), sol.sd_cost): + assert cost == len([k for k, v in design.items() if v != 0]) + + +@pytest.mark.timeout(30) +def test_knockin_candidates_are_not_also_knockout_candidates(curr_solver): + """Naming a reaction addable must not leave it removable as well. + + Reactions the caller did not name stay knockable, so mixed problems are unaffected; + but when every reaction is a knock-in candidate no design may propose a knockout. + """ + model = _two_route_network() + sol = sd.compute_strain_designs(model, + sd_modules=[sd.SDModule(model, SUPPRESS, constraints=['R4 >= 1'])], + max_cost=3, + ki_cost={r: 1.0 for r in ['R1', 'R2', 'R3', 'R4']}, + solution_approach='populate', + solver=curr_solver, + compress=False) + assert not [k for d in sol.get_reaction_sd() for k, v in d.items() if v < 0] From 7794f36e4796e06eead4300643a3b347b26c61a9 Mon Sep 17 00:00:00 2001 From: Phil Date: Thu, 6 Aug 2026 14:52:20 -0400 Subject: [PATCH 2/5] fix: keep compression sound for intervention candidates, and price designs by cost rather than set inclusion Compression could invent designs and lose them. A parallel lump mixing a candidate with a non-targetable reaction let the lump's binary force flux that no intervention can remove, so a knockout was reported that cuts nothing; and the knock-in cost guard compared original reaction ids against a vector already re-keyed by lump names, so it never matched and a mixed lump became a knock-in instead of a knockout. Reaction targetability now takes part in the parallel key, and the guards read the vectors as they entered the compression step. With compress=False the compressed cost vectors aliased the uncompressed ones, which then had essential and size-1-MCS reactions popped out of them. filter_sd_maxcost priced every knockout at zero from there on, so max_cost was not enforced. Designs were selected by set inclusion, which agrees with cost only while every intervention is expensive. A free or rewarding intervention makes the larger design the cheaper one, so a design is now reported when no comparable design costs strictly less: unchanged for positive costs, both variants kept at zero cost, and the design taking a reward preferred over the one without it. Where a reward provably cannot invalidate a design -- a knock-in with only PROTECT modules, a knockout with only SUPPRESS ones -- it is required by a constraint row instead of being discovered, which keeps the enumeration from walking the designs it dominates. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01RFtof9nZXFvCNoXz19po8C --- straindesign/compression.py | 42 +++++++++---- straindesign/compute_strain_designs.py | 57 +++++++++++++++-- straindesign/networktools.py | 44 +++++++++++-- straindesign/strainDesignMILP.py | 85 ++++++++++++++++++++++++-- tests/test_05_straindesign.py | 54 ++++++++++++++++ 5 files changed, 255 insertions(+), 27 deletions(-) diff --git a/straindesign/compression.py b/straindesign/compression.py index a1ba305..70b8464 100644 --- a/straindesign/compression.py +++ b/straindesign/compression.py @@ -2152,7 +2152,20 @@ def simplify_model_gprs(model, budget=50000): logging.info(' GPR rule simplification: %d rules, %d rewritten.' % (n, nchg)) -def compress_model(model, no_par_compress_reacs=set(), propagate_gpr=False, no_coupled_compress_reacs=set()): +def _rename_lumped(reac_set, reac_map_exp): + """Carry a set of reaction ids through one compression step, in place.""" + if reac_set is None: + return + for new_reac, old_reac_val in reac_map_exp.items(): + old_reacs = [r for r in reac_set if r in old_reac_val] + if old_reacs: + for r in old_reacs: + reac_set.discard(r) + reac_set.add(new_reac) + + +def compress_model(model, no_par_compress_reacs=set(), propagate_gpr=False, no_coupled_compress_reacs=set(), + targetable_rxns=None): """Compress a metabolic model using multiple techniques. Performs blocked reaction removal, conservation relation removal, and @@ -2178,6 +2191,7 @@ def compress_model(model, no_par_compress_reacs=set(), propagate_gpr=False, no_c """ from straindesign.networktools import suppress_lp_context, _is_lp_suppressed no_coupled_compress_reacs = set(no_coupled_compress_reacs) + targetable = None if targetable_rxns is None else set(targetable_rxns) with suppress_lp_context(model): cmp_mapReac = [] LOG.info(' Removing blocked reactions.') @@ -2191,11 +2205,13 @@ def compress_model(model, no_par_compress_reacs=set(), propagate_gpr=False, no_c # 1. Parallel (cheap — hash-based, no RREF) LOG.info(f' Compression {run}: Lumping parallel reactions.') - reac_map_exp = compress_model_parallel(model, no_par_compress_reacs, propagate_gpr=propagate_gpr) + reac_map_exp = compress_model_parallel(model, no_par_compress_reacs, propagate_gpr=propagate_gpr, + targetable_rxns=targetable) parallel_changed = numr > len(reac_map_exp) if parallel_changed: LOG.info(f' Reduced to {len(reac_map_exp)} reactions.') cmp_mapReac.append({"reac_map_exp": reac_map_exp, "parallel": True}) + _rename_lumped(targetable, reac_map_exp) # 2. Conservation relation removal (reduces S rows for RREF) remove_conservation_relations(model) @@ -2212,12 +2228,8 @@ def compress_model(model, no_par_compress_reacs=set(), propagate_gpr=False, no_c numr_pre = len(model.reactions) LOG.info(f' Compression {run}: Lumping coupled reactions.') reac_map_exp = compress_model_coupled(model, propagate_gpr=propagate_gpr, protected_reactions=no_coupled_compress_reacs) - for new_reac, old_reac_val in reac_map_exp.items(): - old_reacs = [r for r in no_par_compress_reacs if r in old_reac_val] - if old_reacs: - for r in old_reacs: - no_par_compress_reacs.remove(r) - no_par_compress_reacs.add(new_reac) + _rename_lumped(no_par_compress_reacs, reac_map_exp) + _rename_lumped(targetable, reac_map_exp) coupled_changed = numr_pre > len(reac_map_exp) if coupled_changed: LOG.info(f' Reduced to {len(reac_map_exp)} reactions.') @@ -2284,7 +2296,7 @@ def compress_model_coupled(model, propagate_gpr=False, protected_reactions=set() return reaction_map -def compress_model_parallel(model, protected_rxns=set(), propagate_gpr=False): +def compress_model_parallel(model, protected_rxns=set(), propagate_gpr=False, targetable_rxns=None): """Compress by lumping parallel reactions. Args: @@ -2292,6 +2304,11 @@ def compress_model_parallel(model, protected_rxns=set(), propagate_gpr=False): protected_rxns: Reactions exempt from parallel compression propagate_gpr: If True, OR-combine GPR rules of lumped reactions (with sympy simplification). Default False. + targetable_rxns: Reactions that carry an intervention binary. When given, + targetable and non-targetable reactions are never lumped together: the + lump's binary would force the non-targetable member's flux to zero, an + intervention the original model cannot perform. Pass None to lump purely + by stoichiometry. Returns: dict: Mapping {compressed_id: {orig_id: factor, ...}} @@ -2316,13 +2333,16 @@ def compress_model_parallel(model, protected_rxns=set(), propagate_gpr=False): # scale factor share a key (e.g. -1 A -> 2 B and -3 A -> 6 B both give ((A, 1), (B, -2))). The # reversibility (fwd/rev) and inhomogeneous-bound (inh) flags are part of the key so only # reactions with matching bound structure are lumped. + tgt = [0] * len(model.reactions) if targetable_rxns is None else \ + [1 if r.id in targetable_rxns else 0 for r in model.reactions] + def _parallel_key(i): cols, vals = stoichmat_T.rows[i], stoichmat_T.data[i] if not vals: - return ((), fwd[i], rev[i], inh[i]) + return ((), fwd[i], rev[i], inh[i], tgt[i]) f0 = float_to_fraction(vals[0]) stoich = tuple((int(c), float_to_fraction(v) / f0) for c, v in zip(cols, vals)) - return (stoich, fwd[i], rev[i], inh[i]) + return (stoich, fwd[i], rev[i], inh[i], tgt[i]) # Find parallel reactions by exact key comparison (hash pre-filter, then full compare) protected = [r.id in protected_rxns for r in model.reactions] diff --git a/straindesign/compute_strain_designs.py b/straindesign/compute_strain_designs.py index 33ef1f2..dcef7b0 100644 --- a/straindesign/compute_strain_designs.py +++ b/straindesign/compute_strain_designs.py @@ -33,6 +33,7 @@ from straindesign.networktools import remove_ext_mets, bound_blocked_or_irrevers_fva, \ extend_model_gpr, extend_model_regulatory, evaluate_gpr_ast, \ compress_model, compress_modules, compress_ki_ko_cost, expand_sd, filter_sd_maxcost, \ + filter_sd_dominated, \ estimate_expansion_size, with_suppressed_lp, _silent_io from straindesign.compression import simplify_model_gprs @@ -505,6 +506,13 @@ def compute_strain_designs(model: Model, **kwargs: dict) -> SDSolutions: with _silent_io(): orig_model = model model = model.copy() + _free = [k for k, v in list(uncmp_ko_cost.items()) + list(uncmp_ki_cost.items()) if v == 0.0] + if _free: + logging.warning('%d intervention(s) entered at zero cost (%s%s). Adding one to a design ' + 'changes no cost, so designs that differ only in these are all reported and ' + 'the result set can grow accordingly. Give them a small positive cost to ' + 'have the smaller design reported alone.' % + (len(_free), ', '.join(sorted(_free)[:5]), ', ...' if len(_free) > 5 else '')) orig_ko_cost = deepcopy(uncmp_ko_cost) orig_ki_cost = deepcopy(uncmp_ki_cost) orig_reg_cost = deepcopy(uncmp_reg_cost) @@ -592,18 +600,32 @@ def compute_strain_designs(model: Model, **kwargs: dict) -> SDSolutions: logging.info(' Reversibility pre-tightening fixed %d reaction directions (%.1fs).' % (_n_tight, time.time() - t0)) logging.info('Compressing Network (' + str(len(cmp_model.reactions)) + ' reactions).') t0 = time.time() + # Reactions carrying an intervention binary must not be lumped in parallel with ones + # that carry none; in gene mode a reaction is targetable through its GPR. + targetable_rxns = set(uncmp_ko_cost) | set(uncmp_ki_cost) + if kwargs['gene_kos']: + targetable_rxns |= {r.id for r in cmp_model.reactions if r.gene_reaction_rule} + # A lump mixing knockout and knock-in candidates keeps only one of the two kinds, which + # is sound only while the dropped one costs something. A free or rewarding intervention + # belongs in designs of its own, so keep those candidates out of both kinds of lump. + _free_cands = {k for k, v in list(uncmp_ko_cost.items()) + list(uncmp_ki_cost.items()) if v <= 0.0} + no_par_compress_reacs |= _free_cands + no_coupled_compress_reacs |= _free_cands cmp_mapReac_1 = compress_model(cmp_model, no_par_compress_reacs, propagate_gpr=True, - no_coupled_compress_reacs=no_coupled_compress_reacs) + no_coupled_compress_reacs=no_coupled_compress_reacs, + targetable_rxns=targetable_rxns) sd_modules = compress_modules(sd_modules, cmp_mapReac_1) # Compress reaction + regulatory costs only (gene costs not yet added) cmp_ko_cost, cmp_ki_cost, cmp_mapReac_1 = compress_ki_ko_cost(uncmp_ko_cost, uncmp_ki_cost, cmp_mapReac_1) logging.info(' Compressed to ' + str(len(cmp_model.reactions)) + ' reactions (%.1fs).' % (time.time() - t0)) else: cmp_mapReac_1 = [] - cmp_ko_cost = uncmp_ko_cost - cmp_ki_cost = uncmp_ki_cost + # copy: the compressed vectors get essential and size-1-MCS reactions popped out of + # them, while the uncompressed ones must keep every cost for filter_sd_maxcost + cmp_ko_cost = dict(uncmp_ko_cost) + cmp_ki_cost = dict(uncmp_ki_cost) # --- FVAs on (possibly compressed) model --- logging.info(' FVA(s) to identify essential reactions.') essential_reacs = set() @@ -675,9 +697,12 @@ def compute_strain_designs(model: Model, **kwargs: dict) -> SDSolutions: logging.info('Compressing after GPR extension (' + str(len(cmp_model.reactions)) + ' reactions).') t0 = time.time() no_par_compress_reacs = _collect_no_par_compress_reacs(sd_modules) + _free_cands = {k for k, v in list(cmp_ko_cost.items()) + list(cmp_ki_cost.items()) if v <= 0.0} cmp_mapReac_2 = compress_model( cmp_model, - no_par_compress_reacs, + no_par_compress_reacs | _free_cands, + targetable_rxns=set(cmp_ko_cost) | set(cmp_ki_cost), + no_coupled_compress_reacs=_free_cands, ) sd_modules = compress_modules(sd_modules, cmp_mapReac_2) cmp_ko_cost, cmp_ki_cost, cmp_mapReac_2 = compress_ki_ko_cost(cmp_ko_cost, cmp_ki_cost, cmp_mapReac_2) @@ -759,6 +784,11 @@ def compute_strain_designs(model: Model, **kwargs: dict) -> SDSolutions: size1_mcs = suppress_essential - essential_reacs # Filter to only knockable reactions (in ko_cost, not ki_cost or regulatory) size1_mcs_knockable = {r for r in size1_mcs if r in cmp_ko_cost} + # A rewarding intervention makes a size-1 MCS that can absorb it the dominated design, + # and this shortcut hands solutions to the caller without consulting the exclusion + # constraints that would settle that. Leave those cases to the MILP. + if any(c < 0.0 for c in list(cmp_ko_cost.values()) + list(cmp_ki_cost.values())): + size1_mcs_knockable = set() if size1_mcs_knockable: cmp_size1_mcs = [{r: -1} for r in size1_mcs_knockable] logging.info(' Found ' + str(len(cmp_size1_mcs)) + ' size-1 MCS via SUPPRESS FVA.') @@ -769,8 +799,12 @@ def compute_strain_designs(model: Model, **kwargs: dict) -> SDSolutions: # They are already found; any larger MCS containing them is non-minimal. # But we only remove pure KO candidates — reactions with regulatory or KI # interventions may still participate in non-KO solutions. - for r in size1_mcs_knockable: - cmp_ko_cost.pop(r, None) + # Larger designs containing them are only redundant while every intervention costs + # something. Once one is free or rewarding, such a design can cost no more than the + # size-1 one, so the candidates have to stay in the MILP for it to be reachable. + if not any(c <= 0.0 for c in list(cmp_ko_cost.values()) + list(cmp_ki_cost.values())): + for r in size1_mcs_knockable: + cmp_ko_cost.pop(r, None) # remove ko-costs (and thus knockability) of essential reactions [cmp_ko_cost.pop(er) for er in essential_reacs if er in cmp_ko_cost] @@ -995,6 +1029,17 @@ def _decompress_solutions(cmp_sd_solution, cmp_mapReac, cmp_size1_mcs, max_cost, if cmp_sd_solution.status not in [OPTIMAL, TIME_LIMIT_W_SOL] and sd: cmp_sd_solution.status = OPTIMAL + # The exclusion constraints already keep the MILP's own designs free of dominated ones when + # every intervention costs something, so this scan is only needed where that breaks down: + # free or rewarding interventions, and size-1 MCS, which are injected past the constraints. + if sd and (cmp_size1_mcs or + any(c <= 0.0 for c in list(uncmp_ko_cost.values()) + list(uncmp_ki_cost.values()))): + kept = filter_sd_dominated(sd, uncmp_ko_cost, uncmp_ki_cost) + if len(kept) < len(sd): + keep_ids = {id(s) for s in kept} + group_map = [g for s, g in zip(sd, group_map) if id(s) in keep_ids] + sd = kept + sd_solutions = SDSolutions(orig_model, sd, cmp_sd_solution.status, setup) sd_solutions.compressed_sd = compressed_sd sd_solutions.compression_map = cmp_mapReac diff --git a/straindesign/networktools.py b/straindesign/networktools.py index 59342c5..2b97c88 100644 --- a/straindesign/networktools.py +++ b/straindesign/networktools.py @@ -1206,14 +1206,18 @@ def compress_ki_ko_cost(kocost, kicost, cmp_mapReac): reac_map_exp = cmp["reac_map_exp"] parallel = cmp["parallel"] cmp.update({KOCOST: kocost, KICOST: kicost}) + # Both blocks below decide by which original reactions a lump contains, so they must + # read the vectors as they entered this step. The KO block rebinds kocost to a vector + # keyed by lump names, against which the original ids in reac_map_exp never match. + kocost_in = kocost if kocost: ko_cost_new = {} for r in reac_map_exp: - if np.any([s in kocost for s in reac_map_exp[r]]): + if np.any([s in kocost_in for s in reac_map_exp[r]]): if not parallel and not np.any([s in kicost for s in reac_map_exp[r]]): - ko_cost_new[r] = np.min([kocost[s] for s in reac_map_exp[r] if s in kocost]) + ko_cost_new[r] = np.min([kocost_in[s] for s in reac_map_exp[r] if s in kocost_in]) elif parallel: - ko_cost_new[r] = np.sum([kocost[s] for s in reac_map_exp[r] if s in kocost]) + ko_cost_new[r] = np.sum([kocost_in[s] for s in reac_map_exp[r] if s in kocost_in]) kocost = ko_cost_new if kicost: ki_cost_new = {} @@ -1221,7 +1225,7 @@ def compress_ki_ko_cost(kocost, kicost, cmp_mapReac): if np.any([s in kicost for s in reac_map_exp[r]]): if not parallel: ki_cost_new[r] = np.sum([kicost[s] for s in reac_map_exp[r] if s in kicost]) - elif parallel and not np.any([s in kocost for s in reac_map_exp[r]]): + elif parallel and not np.any([s in kocost_in for s in reac_map_exp[r]]): ki_cost_new[r] = np.min([kicost[s] for s in reac_map_exp[r] if s in kicost]) kicost = ki_cost_new return kocost, kicost, cmp_mapReac @@ -1350,6 +1354,38 @@ def expand_sd(sd, cmp_mapReac): return sd +def filter_sd_dominated(sd, kocost, kicost): + """Drop strain designs that a comparable, strictly cheaper design dominates. + + A design is reported when no other design that is a subset or a superset of it costs + strictly less. With only positive intervention costs every superset is more expensive, + so this reduces to keeping the inclusion-minimal designs -- which the exclusion + constraints already guarantee, making the filter a no-op. It earns its keep once an + intervention is free or rewarding: then taking it can pay for itself, and the smaller + design is the dominated one. Designs of equal cost are all kept. + + Returns: + (list): The non-dominated strain designs, order preserved. + """ + 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 + + sets = [frozenset(k for k, v in m.items() if v != 0 and k != '**cost**') for m in sd] + costs = [np.sum([itv_cost(k, v) for k, v in m.items() if k != '**cost**']) for m in sd] + keep = [] + for i, si in enumerate(sets): + if not any(costs[j] < costs[i] - 1e-8 and (sj <= si or sj >= si) + for j, sj in enumerate(sets) if j != i): + keep.append(i) + if len(keep) < len(sd): + logging.info(' Discarded %d strain design(s) dominated by a cheaper comparable one.' % (len(sd) - len(keep))) + return [sd[i] for i in keep] + + def filter_sd_maxcost(sd, max_cost, kocost, kicost): """Filter out strain designs that exceed the maximum allowed intervention costs diff --git a/straindesign/strainDesignMILP.py b/straindesign/strainDesignMILP.py index fdc6064..b96ceba 100644 --- a/straindesign/strainDesignMILP.py +++ b/straindesign/strainDesignMILP.py @@ -86,8 +86,38 @@ class SDMILP(SDProblem, MILP_LP): def __init__(self, model: Model, sd_modules: List[SDModule], **kwargs): # Construct problem SDProblem.__init__(self, model, sd_modules, **kwargs) + # A knock-in only enlarges the flux space and a knock-out only shrinks it. So when every + # module is of the kind that survives that direction, a rewarding intervention can never + # invalidate a design, and every design is beaten by the one that also takes it. Fixing + # it here is exact and spares the enumeration from walking the designs it dominates. + _types = {m[MODULE_TYPE] for m in sd_modules} + _forced = [] + for i in range(self.num_z): + if self.z_non_targetable[i] or np.isnan(self.cost[i]) or self.cost[i] >= 0.0: + continue + if (self.z_inverted[i] and _types == {PROTECT}) or \ + (not self.z_inverted[i] and _types == {SUPPRESS}): + _forced.append(i) + logging.info(' Rewarding intervention %s cannot invalidate a design here and is ' + 'taken in every one of them.' % model.reactions[i].id) + if _forced: + # as a row rather than by fixing the bound: a binary pinned by lb == ub sends SCIP's + # populate into a loop where it re-reports the same design indefinitely + rows = sparse.lil_matrix((len(_forced), self.A_ineq.shape[1])) + for k, i in enumerate(_forced): + rows[k, i] = -1.0 + self.A_ineq = sparse.vstack((self.A_ineq, rows.tocsr()), format='csr') + self.b_ineq = list(self.b_ineq) + [-1.0] * len(_forced) # Remove non-knockable z-variables before solver sees them self._trim_z_variables() + # Interventions that do not cost anything to take. Adding one to a design can only keep + # its cost equal or lower, so such a design is not dominated by the smaller one and the + # exclusion constraints must leave it reachable. Empty for the usual all-positive setup, + # where dominance and set inclusion agree and the cuts stay exactly as they were. + self._free_z = [i for i in range(self.num_z) if not self.z_non_targetable[i] and self.cost[i] <= 0.0] + # A rewarding intervention strictly lowers the cost of any design that can absorb it, + # so a design is only worth reporting once none of them can be added while staying valid. + self._rewarding_z = [i for i in self._free_z if self.cost[i] < 0.0] # Build MILP object from constructed problem MILP_LP.__init__(self, c=self.c, @@ -104,6 +134,21 @@ def __init__(self, model: Model, sd_modules: List[SDModule], **kwargs): seed=self.seed, milp_threads=self.milp_threads) + def is_dominated(self, z): + """True if a rewarding intervention can join this design without invalidating it. + + Such a design is strictly cheaper and contains this one, so this one is not worth + reporting. Only reachable when some intervention carries a negative cost. + """ + for j in self._rewarding_z: + if z[0, j]: + continue + z_more = z.tolil() + z_more[0, j] = 1.0 + if all(self.verify_sd(z_more.tocsr())): + return True + return False + def _trim_z_variables(self): """Remove non-knockable (ub=0) z-variables from MILP matrices. @@ -160,16 +205,24 @@ def _expand_z_to_orig(self, z_trimmed): return expanded.tocsr() def add_exclusion_constraints(self, z): - """Exclude binary solution in z and all supersets from MILP""" + """Exclude a design and every superset of it that cannot be cheaper. + + With only positive intervention costs a superset always costs more, so this excludes + all supersets and the rule is plain set inclusion. When some interventions are free or + rewarding (cost <= 0), a superset taking one of those is not dominated by the design + found here, so those literals enter the row negated and such supersets stay reachable. + """ for i in range(z.shape[0]): + free = [j for j in self._free_z if not z[i, j]] # introduce constraint to make MILP infeasible. Some solvers cannot handle empty rows - if z[i].nnz == 0: + if z[i].nnz == 0 and not free: A_ineq = sparse.csr_matrix([1.0] * z[i].shape[1]) A_ineq.resize((1, self.A_ineq.shape[1])) b_ineq = -1 self.add_ineq_constraints(A_ineq, [b_ineq]) - # otherwise, introduce integer cut constraint - elif z[i].nnz == 1: + # a single intervention can be switched off outright, but only while no free + # intervention could join it to form a design that is not dominated by this one + elif z[i].nnz == 1 and not free: interv_idx = int(z[i].indices[0]) self.z_non_targetable[interv_idx] = True self.set_ub([[interv_idx, 0.0]]) @@ -178,6 +231,11 @@ def add_exclusion_constraints(self, z): A_ineq = z[i].copy() A_ineq.resize((1, self.A_ineq.shape[1])) b_ineq = np.sum(z[i]) - 1 + if free: + A_ineq = A_ineq.tolil() + for j in free: + A_ineq[0, j] = -1.0 + A_ineq = A_ineq.tocsr() self.add_ineq_constraints(A_ineq, [b_ineq]) def add_exclusion_constraints_ineq(self, z): @@ -556,10 +614,14 @@ def enumerate(self, **kwargs): if self.show_no_ki is None: self.show_no_ki = True # first check if strain doesn't already fulfill the strain design setup - if self.is_mcs_computation and self.verify_sd(sparse.csr_matrix((1, self.num_z)))[0]: + wt_is_design = self.is_mcs_computation and self.verify_sd(sparse.csr_matrix((1, self.num_z)))[0] + if wt_is_design: logging.warning('The strain already meets the requirements defined in the strain design setup. ' \ 'No interventions are needed.') - return self.build_sd_solution([{}], OPTIMAL, POPULATE) + # Free interventions can still yield designs that cost no more than doing nothing, so + # keep enumerating and let the exclusion constraints carry the empty design forward. + if not self._free_z: + return self.build_sd_solution([{}], OPTIMAL, POPULATE) # otherwise continue if self.solver == 'scip': logging.warning("SCIP does not natively support solution pool generation. "+ \ @@ -574,6 +636,11 @@ def enumerate(self, **kwargs): endtime = time.time() + self.time_limit status = OPTIMAL sols = sparse.csr_matrix((0, self.num_z)) + if wt_is_design: + empty = sparse.csr_matrix((1, self.num_z)) + if not (self._rewarding_z and self.is_dominated(empty)): + sols = sparse.vstack((sols, empty)) + self.add_exclusion_constraints(empty) logging.info('Enumerating strain designs ...') while sols.shape[0] < self.max_solutions and \ status == OPTIMAL and \ @@ -594,6 +661,12 @@ def enumerate(self, **kwargs): for i in range(z.shape[0]): output = [self.sd2dict(z[i])] if all(self.verify_sd(z[i])): + if self._rewarding_z and self.is_dominated(z[i]): + # a cheaper design contains this one; cut only the exact pattern so + # that the design dominating it stays reachable + logging.info('Dominated by a cheaper superset, skipping: ' + str(output)) + self.add_exclusion_constraints_ineq(z[i]) + continue logging.info('Strain designs with cost ' + str(round((z[i] * self.cost)[0], 6)) + ': ' + str(output)) self.add_exclusion_constraints(z[i]) sols = sparse.vstack((sols, z[i])) diff --git a/tests/test_05_straindesign.py b/tests/test_05_straindesign.py index b0eb474..2e534a3 100644 --- a/tests/test_05_straindesign.py +++ b/tests/test_05_straindesign.py @@ -520,3 +520,57 @@ def test_knockin_candidates_are_not_also_knockout_candidates(curr_solver): solver=curr_solver, compress=False) assert not [k for d in sol.get_reaction_sd() for k, v in d.items() if v < 0] + + +@pytest.mark.timeout(30) +def test_rewarding_intervention_beats_the_design_without_it(curr_solver): + """A design is reported only when no comparable design costs less. + + With a negative cost, taking R2 costs less than not taking it, so the design that + takes it dominates the one that does not -- the reverse of the usual case, where a + larger design always costs more and the smaller one wins. + """ + model = _two_route_network() + sol = sd.compute_strain_designs(model, + sd_modules=[sd.SDModule(model, PROTECT, constraints=['R4 >= 1'])], + max_cost=3, + ki_cost={'R1': 1.0, 'R2': -5.0, 'R3': 1.0, 'R4': 1.0}, + solution_approach='populate', + solver=curr_solver, + compress=False) + designs = _designs(sol) + assert designs, 'expected at least one design' + assert all('R2' in d for d in designs), designs + + +@pytest.mark.timeout(30) +def test_zero_cost_intervention_ties_instead_of_dominating(curr_solver): + """At zero cost neither variant is cheaper, so both are reported.""" + model = _two_route_network() + sol = sd.compute_strain_designs(model, + sd_modules=[sd.SDModule(model, PROTECT, constraints=['R4 >= 1'])], + max_cost=3, + ki_cost={'R1': 1.0, 'R2': 0.0, 'R3': 1.0, 'R4': 1.0}, + solution_approach='populate', + solver=curr_solver, + compress=False) + designs = _designs(sol) + assert ['R1', 'R2', 'R4'] in designs, designs + assert ['R1', 'R3', 'R4'] in designs, designs + + +@pytest.mark.timeout(60) +def test_positive_costs_keep_the_smaller_design_only(curr_solver): + """The classical rule is unchanged: no reported design contains another.""" + model = _two_route_network() + sol = sd.compute_strain_designs(model, + sd_modules=[sd.SDModule(model, SUPPRESS, constraints=['R4 >= 1'])], + max_cost=3, + ko_cost={'R1': 1.0, 'R2': 1.0, 'R3': 1.0, 'R4': 1.0}, + solution_approach='populate', + solver=curr_solver, + compress=False) + designs = [set(d) for d in _designs(sol)] + for i, a in enumerate(designs): + for j, b in enumerate(designs): + assert i == j or not a < b, (a, b) From 853b6b392c9ca41c533f8c94938cd73295d7d982 Mon Sep 17 00:00:00 2001 From: Phil Date: Thu, 6 Aug 2026 15:04:45 -0400 Subject: [PATCH 3/5] fix: require essential knock-ins by a row instead of a pinned bound Both ways of saying "this intervention is part of every design" are equivalent to gurobi, cplex and glpk, but with the bound form SCIP's populate re-reports the same design instead of reporting infeasible once the designs are exhausted, so enumeration never terminates. Measured head to head with only the form varying. link_z's step 7 skips inverted z-variables, so nothing downstream relied on the pinned bound. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01RFtof9nZXFvCNoXz19po8C --- straindesign/strainDesignMILP.py | 5 +++-- straindesign/strainDesignProblem.py | 15 ++++++++++++++- 2 files changed, 17 insertions(+), 3 deletions(-) diff --git a/straindesign/strainDesignMILP.py b/straindesign/strainDesignMILP.py index b96ceba..41726ff 100644 --- a/straindesign/strainDesignMILP.py +++ b/straindesign/strainDesignMILP.py @@ -101,8 +101,9 @@ def __init__(self, model: Model, sd_modules: List[SDModule], **kwargs): logging.info(' Rewarding intervention %s cannot invalidate a design here and is ' 'taken in every one of them.' % model.reactions[i].id) if _forced: - # as a row rather than by fixing the bound: a binary pinned by lb == ub sends SCIP's - # populate into a loop where it re-reports the same design indefinitely + # as a row rather than by fixing the bound: with the bound form SCIP's populate + # re-reports the same design instead of reporting infeasible once the designs are + # exhausted, which never terminates (measured; the other solvers accept either form) rows = sparse.lil_matrix((len(_forced), self.A_ineq.shape[1])) for k, i in enumerate(_forced): rows[k, i] = -1.0 diff --git a/straindesign/strainDesignProblem.py b/straindesign/strainDesignProblem.py index 44a0738..d2d2936 100644 --- a/straindesign/strainDesignProblem.py +++ b/straindesign/strainDesignProblem.py @@ -175,7 +175,7 @@ def __init__(self, model: Model, sd_modules: List[SDModule], *args, **kwargs): else: 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.lb = [0.0 for _ in model.reactions] self.ub = [1.0 - float(i) for i in self.z_non_targetable] self.idx_z = [i for i in range(0, numr)] self.c = [0.0] * numr @@ -222,6 +222,19 @@ def __init__(self, model: Model, sd_modules: List[SDModule], *args, **kwargs): # 4. Link LP module to z-variables self.link_z() + # An essential knock-in has to be part of every design. Required by a row rather than by + # pinning the binary's lower bound: with the bound form SCIP's populate re-reports the + # same design instead of reporting infeasible once the designs are exhausted (measured; + # the other three solvers are indifferent to which form is used). link_z's step 7 leaves + # inverted z's alone, so nothing downstream depends on the pinned bound. + _ess = [i for i, r in enumerate(model.reactions) if r.id in self.essential_kis] + if _ess: + _rows = sparse.lil_matrix((len(_ess), self.A_ineq.shape[1])) + for _k, _i in enumerate(_ess): + _rows[_k, _i] = -1.0 + self.A_ineq = sparse.vstack((self.A_ineq, _rows.tocsr()), format='csr') + self.b_ineq = list(self.b_ineq) + [-1.0] * len(_ess) + # if there are only mcs modules, minimize the knockout costs, # otherwise use objective function(s) from modules if all([mod[MODULE_TYPE] in [PROTECT, SUPPRESS, DOUBLEOPT] for mod in sd_modules]): From 7d5a9c8ba091322739a0205868996662a64f0aaa Mon Sep 17 00:00:00 2001 From: Phil Date: Thu, 6 Aug 2026 15:36:40 -0400 Subject: [PATCH 4/5] refactor: share the intervention-cost lookup and select designs by index Both cost filters derived an intervention's price the same way; that lookup is now one function. The dominance scan returns indices so the caller can keep group_map aligned by position rather than by object identity, and it now also runs on the lazy-expansion path, where size-1 MCS are injected just as they are on the eager one. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01RFtof9nZXFvCNoXz19po8C --- straindesign/compute_strain_designs.py | 32 ++++++++++------ straindesign/networktools.py | 53 ++++++++++++-------------- straindesign/strainDesignMILP.py | 2 +- 3 files changed, 46 insertions(+), 41 deletions(-) diff --git a/straindesign/compute_strain_designs.py b/straindesign/compute_strain_designs.py index dcef7b0..c52648c 100644 --- a/straindesign/compute_strain_designs.py +++ b/straindesign/compute_strain_designs.py @@ -33,7 +33,7 @@ from straindesign.networktools import remove_ext_mets, bound_blocked_or_irrevers_fva, \ extend_model_gpr, extend_model_regulatory, evaluate_gpr_ast, \ compress_model, compress_modules, compress_ki_ko_cost, expand_sd, filter_sd_maxcost, \ - filter_sd_dominated, \ + nondominated_sd, \ estimate_expansion_size, with_suppressed_lp, _silent_io from straindesign.compression import simplify_model_gprs @@ -958,6 +958,22 @@ def postprocess_reg_sd(reg_cost, sd): LAZY_EXPANSION_THRESHOLD = 100_000 +def _drop_dominated(sd, group_map, cmp_size1_mcs, uncmp_ko_cost, uncmp_ki_cost): + """Remove designs a comparable, cheaper design dominates, keeping group_map aligned. + + The exclusion constraints already keep the MILP's own designs free of dominated ones while + every intervention costs something, so this scan is only needed where that breaks down: + free or rewarding interventions, and size-1 MCS, which are injected past the constraints. + """ + if not sd or not (cmp_size1_mcs or + any(c <= 0.0 for c in list(uncmp_ko_cost.values()) + list(uncmp_ki_cost.values()))): + return sd, group_map + keep = nondominated_sd(sd, uncmp_ko_cost, uncmp_ki_cost) + if len(keep) == len(sd): + return sd, group_map + return [sd[i] for i in keep], [group_map[i] for i in keep] + + def _decompress_solutions(cmp_sd_solution, cmp_mapReac, cmp_size1_mcs, max_cost, uncmp_ko_cost, uncmp_ki_cost, uncmp_reg_cost, orig_model, setup, gene_kos, orig_gko_cost, orig_gki_cost): """Decompress MILP solutions, using lazy expansion if estimated count exceeds threshold.""" @@ -977,6 +993,9 @@ def _decompress_solutions(cmp_sd_solution, cmp_mapReac, cmp_size1_mcs, max_cost, logging.info(' Estimated %d expanded solutions - using lazy expansion.' % estimated) sd, group_map, compressed_sd = _build_lazy_representatives(cmp_sds, cmp_size1_mcs, cmp_mapReac, max_cost, uncmp_ko_cost, uncmp_ki_cost, uncmp_reg_cost) + # only the representatives can be compared here; the rest of each group is never + # materialised, so a dominated design inside an unexpanded group is not caught + sd, group_map = _drop_dominated(sd, group_map, cmp_size1_mcs, uncmp_ko_cost, uncmp_ki_cost) status = cmp_sd_solution.status if status not in [OPTIMAL, TIME_LIMIT_W_SOL] and sd: @@ -1029,16 +1048,7 @@ def _decompress_solutions(cmp_sd_solution, cmp_mapReac, cmp_size1_mcs, max_cost, if cmp_sd_solution.status not in [OPTIMAL, TIME_LIMIT_W_SOL] and sd: cmp_sd_solution.status = OPTIMAL - # The exclusion constraints already keep the MILP's own designs free of dominated ones when - # every intervention costs something, so this scan is only needed where that breaks down: - # free or rewarding interventions, and size-1 MCS, which are injected past the constraints. - if sd and (cmp_size1_mcs or - any(c <= 0.0 for c in list(uncmp_ko_cost.values()) + list(uncmp_ki_cost.values()))): - kept = filter_sd_dominated(sd, uncmp_ko_cost, uncmp_ki_cost) - if len(kept) < len(sd): - keep_ids = {id(s) for s in kept} - group_map = [g for s, g in zip(sd, group_map) if id(s) in keep_ids] - sd = kept + sd, group_map = _drop_dominated(sd, group_map, cmp_size1_mcs, uncmp_ko_cost, uncmp_ki_cost) sd_solutions = SDSolutions(orig_model, sd, cmp_sd_solution.status, setup) sd_solutions.compressed_sd = compressed_sd diff --git a/straindesign/networktools.py b/straindesign/networktools.py index 2b97c88..db2c4be 100644 --- a/straindesign/networktools.py +++ b/straindesign/networktools.py @@ -1354,36 +1354,41 @@ def expand_sd(sd, cmp_mapReac): return sd -def filter_sd_dominated(sd, kocost, kicost): - """Drop strain designs that a comparable, strictly cheaper design dominates. +def itv_cost(k, v, kocost, kicost): + """Cost of one intervention, picking the dict by the sign of the entry. + + A reaction may appear in both dicts, since kocost defaults to every reaction. The value's + sign says which kind of intervention was made: positive is an addition, negative a removal. + """ + 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 + + +def nondominated_sd(sd, kocost, kicost): + """Indices of the strain designs no comparable, strictly cheaper design dominates. A design is reported when no other design that is a subset or a superset of it costs strictly less. With only positive intervention costs every superset is more expensive, so this reduces to keeping the inclusion-minimal designs -- which the exclusion - constraints already guarantee, making the filter a no-op. It earns its keep once an + constraints already guarantee, making the scan a no-op. It earns its keep once an intervention is free or rewarding: then taking it can pay for itself, and the smaller design is the dominated one. Designs of equal cost are all kept. Returns: - (list): The non-dominated strain designs, order preserved. + (list of int): Indices into sd to keep, in order. """ - 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 - sets = [frozenset(k for k, v in m.items() if v != 0 and k != '**cost**') for m in sd] - costs = [np.sum([itv_cost(k, v) for k, v in m.items() if k != '**cost**']) for m in sd] - keep = [] - for i, si in enumerate(sets): - if not any(costs[j] < costs[i] - 1e-8 and (sj <= si or sj >= si) - for j, sj in enumerate(sets) if j != i): - keep.append(i) + costs = [np.sum([itv_cost(k, v, kocost, kicost) for k, v in m.items() if k != '**cost**']) for m in sd] + keep = [ + i for i, si in enumerate(sets) + if not any(costs[j] < costs[i] - 1e-8 and (sj <= si or sj >= si) for j, sj in enumerate(sets) if j != i) + ] if len(keep) < len(sd): logging.info(' Discarded %d strain design(s) dominated by a cheaper comparable one.' % (len(sd) - len(keep))) - return [sd[i] for i in keep] + return keep def filter_sd_maxcost(sd, max_cost, kocost, kicost): @@ -1397,18 +1402,8 @@ 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([itv_cost(k, v) for k, v in m.items()]) for m in sd] + costs = [np.sum([itv_cost(k, v, kocost, kicost) 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 diff --git a/straindesign/strainDesignMILP.py b/straindesign/strainDesignMILP.py index 41726ff..adada85 100644 --- a/straindesign/strainDesignMILP.py +++ b/straindesign/strainDesignMILP.py @@ -93,7 +93,7 @@ def __init__(self, model: Model, sd_modules: List[SDModule], **kwargs): _types = {m[MODULE_TYPE] for m in sd_modules} _forced = [] for i in range(self.num_z): - if self.z_non_targetable[i] or np.isnan(self.cost[i]) or self.cost[i] >= 0.0: + if self.z_non_targetable[i] or self.cost[i] >= 0.0: continue if (self.z_inverted[i] and _types == {PROTECT}) or \ (not self.z_inverted[i] and _types == {SUPPRESS}): From eadb0943731002f05339a59c00e60bc30c335d7e Mon Sep 17 00:00:00 2001 From: Phil Date: Thu, 6 Aug 2026 16:23:32 -0400 Subject: [PATCH 5/5] fix: correct the gene-cost compression expectations, and narrow what stays out of coupled lumps compress_ki_ko_cost's knock-in guard now works, and that propagates: the knockout branch consults the knock-in vector, so a corrected knock-in vector in one compression step changes the knockout vector in the next. In the GPR test model the second step lumps g7's block in parallel with g5/g9's. A parallel lump carries flux while any member does, so it is a knockout candidate and g7 stops being an addition candidate -- the counts the test pinned described the behaviour of the dead guard. Coupled lumps only discard candidates when a knock-in member turns the lump into an addition candidate, so a lone knock-in no longer has to stay out of them. Knockout candidates still do: the lump's cost is the cheapest member, while expanding it offers each member as its own design at its own price, and a free or rewarding member makes those disagree. Zero-cost warning now covers gene costs too. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01RFtof9nZXFvCNoXz19po8C --- straindesign/compute_strain_designs.py | 47 ++++++++++++++++++++------ tests/test_04_preprocessing.py | 8 +++-- 2 files changed, 43 insertions(+), 12 deletions(-) diff --git a/straindesign/compute_strain_designs.py b/straindesign/compute_strain_designs.py index c52648c..25c158b 100644 --- a/straindesign/compute_strain_designs.py +++ b/straindesign/compute_strain_designs.py @@ -506,7 +506,9 @@ def compute_strain_designs(model: Model, **kwargs: dict) -> SDSolutions: with _silent_io(): orig_model = model model = model.copy() - _free = [k for k, v in list(uncmp_ko_cost.items()) + list(uncmp_ki_cost.items()) if v == 0.0] + _free = [k for k, v in list(uncmp_ko_cost.items()) + list(uncmp_ki_cost.items()) + + list(locals().get('uncmp_gko_cost', {}).items()) + list(locals().get('uncmp_gki_cost', {}).items()) + if v == 0.0] if _free: logging.warning('%d intervention(s) entered at zero cost (%s%s). Adding one to a design ' 'changes no cost, so designs that differ only in these are all reported and ' @@ -605,12 +607,8 @@ def compute_strain_designs(model: Model, **kwargs: dict) -> SDSolutions: targetable_rxns = set(uncmp_ko_cost) | set(uncmp_ki_cost) if kwargs['gene_kos']: targetable_rxns |= {r.id for r in cmp_model.reactions if r.gene_reaction_rule} - # A lump mixing knockout and knock-in candidates keeps only one of the two kinds, which - # is sound only while the dropped one costs something. A free or rewarding intervention - # belongs in designs of its own, so keep those candidates out of both kinds of lump. - _free_cands = {k for k, v in list(uncmp_ko_cost.items()) + list(uncmp_ki_cost.items()) if v <= 0.0} - no_par_compress_reacs |= _free_cands - no_coupled_compress_reacs |= _free_cands + no_par_compress_reacs |= _free_par_reacs(uncmp_ko_cost, uncmp_ki_cost) + no_coupled_compress_reacs |= _free_coupled_reacs(uncmp_ko_cost, uncmp_ki_cost) cmp_mapReac_1 = compress_model(cmp_model, no_par_compress_reacs, propagate_gpr=True, @@ -697,12 +695,11 @@ def compute_strain_designs(model: Model, **kwargs: dict) -> SDSolutions: logging.info('Compressing after GPR extension (' + str(len(cmp_model.reactions)) + ' reactions).') t0 = time.time() no_par_compress_reacs = _collect_no_par_compress_reacs(sd_modules) - _free_cands = {k for k, v in list(cmp_ko_cost.items()) + list(cmp_ki_cost.items()) if v <= 0.0} cmp_mapReac_2 = compress_model( cmp_model, - no_par_compress_reacs | _free_cands, + no_par_compress_reacs | _free_par_reacs(cmp_ko_cost, cmp_ki_cost), targetable_rxns=set(cmp_ko_cost) | set(cmp_ki_cost), - no_coupled_compress_reacs=_free_cands, + no_coupled_compress_reacs=_free_coupled_reacs(cmp_ko_cost, cmp_ki_cost), ) sd_modules = compress_modules(sd_modules, cmp_mapReac_2) cmp_ko_cost, cmp_ki_cost, cmp_mapReac_2 = compress_ki_ko_cost(cmp_ko_cost, cmp_ki_cost, cmp_mapReac_2) @@ -958,6 +955,36 @@ def postprocess_reg_sd(reg_cost, sd): LAZY_EXPANSION_THRESHOLD = 100_000 +def _free_par_reacs(kocost, kicost): + """Candidates that must stay out of parallel lumps because their cost is not positive. + + A parallel lump carries flux while any member does. Switching it off means knocking out + every knockout candidate in it, so its cost is their sum -- which stops being the cheapest + way to leave it on once one of them is free or rewarding, since knocking only that one is + cheaper still and leaves the lump on. The lump also drops its knock-in members, whose own + designs are worth reporting at that point. Both kinds are therefore kept out. + """ + return {k for k, v in kocost.items() if v <= 0.0} | {k for k, v in kicost.items() if v <= 0.0} + + +def _free_coupled_reacs(kocost, kicost): + """The same for coupled lumps, where slightly less has to be excluded. + + A coupled lump dies with any one of its members, so its knockout cost is the cheapest of + them -- but expanding it offers every member as an alternative design, each at its own + price. That is harmless while the members cost the same and the cheapest is representative; + a free or rewarding member breaks it, and the siblings come back at prices the lump never + searched at. Knockout candidates are therefore excluded outright. + + A knock-in member instead makes the lump an addition candidate costing the sum of its + knock-ins, with no such choice at expansion time. One on its own is safe; only a second one + that could join it in the same lump makes exclusion necessary. + """ + ko = {k for k, v in kocost.items() if v <= 0.0} + ki = {k for k, v in kicost.items() if v <= 0.0} if len(kicost) > 1 else set() + return ko | ki + + def _drop_dominated(sd, group_map, cmp_size1_mcs, uncmp_ko_cost, uncmp_ki_cost): """Remove designs a comparable, cheaper design dominates, keeping group_map aligned. diff --git a/tests/test_04_preprocessing.py b/tests/test_04_preprocessing.py index 005e7a5..5df183e 100644 --- a/tests/test_04_preprocessing.py +++ b/tests/test_04_preprocessing.py @@ -38,8 +38,12 @@ def test_gpr_extension_compression1(model_gpr): cmp_map = sd.compress_model(model_gpr) assert (len(model_gpr.reactions) == 16) gkocost, gkicost, cmp_map = sd.compress_ki_ko_cost(gkocost, gkicost, cmp_map) - assert (len(gkocost) == 4) - assert (len(gkicost) == 3) + # The second compression step lumps g7's block in parallel with g5/g9's. A parallel lump + # carries flux while any member does, so it is switched off only by knocking out all its + # knockout candidates, and the knock-in among them cannot turn it on -- it is a knockout + # candidate costing the sum of those, and g7 stops being an addition candidate. + assert (len(gkocost) == 5) + assert (len(gkicost) == 2) @pytest.mark.timeout(15)