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 217e6bf..25c158b 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, \ + nondominated_sd, \ estimate_expansion_size, with_suppressed_lp, _silent_io from straindesign.compression import simplify_model_gprs @@ -473,8 +474,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]: @@ -498,6 +506,15 @@ 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()) + + 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 ' + '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) @@ -585,18 +602,28 @@ 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} + 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, - 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() @@ -670,7 +697,9 @@ def compute_strain_designs(model: Model, **kwargs: dict) -> SDSolutions: no_par_compress_reacs = _collect_no_par_compress_reacs(sd_modules) cmp_mapReac_2 = compress_model( cmp_model, - no_par_compress_reacs, + 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_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) @@ -752,6 +781,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.') @@ -762,8 +796,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] @@ -917,6 +955,52 @@ 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. + + 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.""" @@ -936,6 +1020,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: @@ -988,6 +1075,8 @@ 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 + 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 sd_solutions.compression_map = cmp_mapReac diff --git a/straindesign/networktools.py b/straindesign/networktools.py index 59342c5..db2c4be 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,43 @@ def expand_sd(sd, cmp_mapReac): return sd +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 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 of int): Indices into sd to keep, in order. + """ + 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, 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 keep + + def filter_sd_maxcost(sd, max_cost, kocost, kicost): """Filter out strain designs that exceed the maximum allowed intervention costs @@ -1361,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 fdc6064..adada85 100644 --- a/straindesign/strainDesignMILP.py +++ b/straindesign/strainDesignMILP.py @@ -86,8 +86,39 @@ 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 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: 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 + 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 +135,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 +206,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 +232,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 +615,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 +637,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 +662,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/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]): 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_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) diff --git a/tests/test_05_straindesign.py b/tests/test_05_straindesign.py index f24ea1f..2e534a3 100644 --- a/tests/test_05_straindesign.py +++ b/tests/test_05_straindesign.py @@ -427,3 +427,150 @@ 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] + + +@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)