From 0d39102204d1aa7df7c87376eb9f1e3e478fc1ed Mon Sep 17 00:00:00 2001 From: Phil Date: Sat, 25 Jul 2026 16:57:37 -0400 Subject: [PATCH] chore: apply yapf and publish artifacts before tagging Two changes to get out of the release's way. The repository had drifted from the yapf style the release workflow enforces, so `yapf -i -r .` would have reformatted 25 files as a side effect of the version bump, burying the release in unrelated churn. Running it now separates the formatting from the release. The diff is cosmetic: comparing token streams before and after, the only non-whitespace change is the removal of statement-separating semicolons, where yapf splits `a = 1; b = 2` onto two lines. e_coli_core gene-MCS gates at 455, set-identical to the reference. The publish workflow created the git tag and the GitHub release before it built or uploaded anything, so any failure in the build, the metadata check or the upload left behind a tag and a release advertising a version that does not exist on PyPI. Since PyPI refuses to accept a version a second time even after it is deleted, that is awkward to unwind: the fix is a new version number, not a retry. Now the distributions are built and pass `twine check` before the repository is touched at all, and the release is cut only after PyPI has accepted the upload. A failure before that point leaves nothing to clean up. Dry run of the reordered sequence at 1.19: both artifacts build, twine check passes on each, and the sdist installs and imports with the expected metadata. Co-Authored-By: Claude Opus 4.8 (1M context) --- .github/workflows/python-publish.yml | 36 +-- straindesign/compression.py | 191 +++++++++------- straindesign/compute_strain_designs.py | 178 ++++++++------- straindesign/cplex_interface.py | 23 +- straindesign/glpk_interface.py | 3 +- straindesign/gurobi_interface.py | 24 +- straindesign/lptools.py | 115 ++++------ straindesign/networktools.py | 59 +++-- straindesign/parse_constr.py | 1 - straindesign/scip_interface.py | 19 +- straindesign/solver_interface.py | 3 +- straindesign/speedy_fva.py | 219 +++++++++++-------- straindesign/strainDesignMILP.py | 2 +- straindesign/strainDesignModule.py | 11 +- straindesign/strainDesignProblem.py | 54 ++--- straindesign/strainDesignSolutions.py | 33 ++- tests/conftest.py | 13 +- tests/test_03_plots.py | 9 +- tests/test_04_preprocessing.py | 103 ++++----- tests/test_05_straindesign.py | 68 ++++-- tests/test_07_compression.py | 28 +-- tests/test_08_gene_ko.py | 8 + tests/test_09_performance.py | 128 +++++------ tests/test_10_gene_design_validity.py | 27 ++- tests/test_11_gene_regulatory_compression.py | 57 +++-- tests/test_12_save_load_embed_model.py | 147 ++++++++----- 26 files changed, 848 insertions(+), 711 deletions(-) diff --git a/.github/workflows/python-publish.yml b/.github/workflows/python-publish.yml index 2cf9fdc..6de21e4 100644 --- a/.github/workflows/python-publish.yml +++ b/.github/workflows/python-publish.yml @@ -38,6 +38,20 @@ jobs: pip install yapf - name: format repository with yapf in google-python format run: yapf -i -r --style ./.github/style.yapf . + # Build and validate the distributions before anything is published or pushed: + # a failure here leaves the repository untouched, with no tag, release or + # version-bump commit to unwind. + - name: Install build and twine + run: | + pip install --upgrade build twine + - name: Clean PyPi build directories + run: rm -rf dist + - name: Build for PyPi + run: | + python -m build + - name: Check the built distributions + run: | + twine check dist/* - name: retrieve name and E-Mail configuration run: | git config --global user.name ${{ secrets.NAME_GITHUB }} @@ -52,6 +66,13 @@ jobs: run: git fetch - name: push code to main run: git push + - name: Publish to PyPI + uses: pypa/gh-action-pypi-publish@release/v1 + with: + password: ${{ secrets.PYPI_API_TOKEN }} + # Tag and release only once PyPI has accepted the upload. A version cannot be + # re-uploaded to PyPI even after deletion, so the release must not advertise a + # version that never made it there. - name: Create Release env: GITHUB_TOKEN: ${{ github.token }} @@ -59,21 +80,6 @@ jobs: gh release create "v${{ github.event.inputs.version }}" \ --title "v${{ github.event.inputs.version }}" \ --notes "${{ github.event.inputs.description }}" - - name: Install build and twine - run: | - pip install --upgrade build twine - - name: Clean PyPi build directories - run: rm -rf dist - - name: Build for PyPi - run: | - python -m build - - name: Check the built distributions - run: | - twine check dist/* - - name: Publish to PyPI - uses: pypa/gh-action-pypi-publish@release/v1 - with: - password: ${{ secrets.PYPI_API_TOKEN }} - name: Install conda dependencies run: conda install anaconda-client conda-build - name: Clean Anaconda build directories diff --git a/straindesign/compression.py b/straindesign/compression.py index 5c08692..a1ba305 100644 --- a/straindesign/compression.py +++ b/straindesign/compression.py @@ -101,7 +101,6 @@ def _lcm_list(numbers: List[int]) -> int: # Rational Matrix with Sparse Storage - _INT64_MAX = (1 << 63) - 1 @@ -372,8 +371,9 @@ def scale_column(self, col: int, scalar_num: int, scalar_den: int) -> None: num_lil, den_lil = self._num_sparse, self._den_sparse num_csc = num_lil.tocsc() if num_lil.format != 'csc' else num_lil den_csc = den_lil.tocsc() if den_lil.format != 'csc' else den_lil - entries = [(num_csc.indices[i], int(num_csc.data[i]), int(den_csc.data[i])) - for i in range(num_csc.indptr[col], num_csc.indptr[col + 1])] + entries = [ + (num_csc.indices[i], int(num_csc.data[i]), int(den_csc.data[i])) for i in range(num_csc.indptr[col], num_csc.indptr[col + 1]) + ] for row, cur_num, cur_den in entries: if cur_num == 0: continue @@ -407,10 +407,9 @@ def to_sparse_csr(self) -> Tuple[csr_matrix, int]: Returns numerator matrix scaled by LCM of denominators, plus the LCM. """ if self.is_bigint(): - raise OverflowError( - "coefficients exceed int64 and cannot be stored in a scipy sparse matrix; " - "use to_coo_exact() / to_sparse_pattern(), or the public sparse_nullspace() helper " - "which returns an ExactCOO in that case.") + raise OverflowError("coefficients exceed int64 and cannot be stored in a scipy sparse matrix; " + "use to_coo_exact() / to_sparse_pattern(), or the public sparse_nullspace() helper " + "which returns an ExactCOO in that case.") # Empty matrix (e.g. a 0-dimensional nullspace) if self._den_sparse is None: return csr_matrix((self._rows, self._cols), dtype=np.int64), 1 @@ -421,8 +420,7 @@ def to_sparse_csr(self) -> Tuple[csr_matrix, int]: common_denom = _lcm_list([int(d) for d in dens if d != 0]) coo_num = self._num_sparse.tocoo() coo_den = self._den_sparse.tocoo() - scaled_data = [int(num) * (common_denom // int(den)) if num != 0 else 0 - for num, den in zip(coo_num.data, coo_den.data)] + scaled_data = [int(num) * (common_denom // int(den)) if num != 0 else 0 for num, den in zip(coo_num.data, coo_den.data)] scaled = csr_matrix((scaled_data, (coo_num.row, coo_num.col)), shape=(self._rows, self._cols), dtype=np.int64) return scaled, common_denom @@ -443,13 +441,17 @@ def to_coo_exact(self) -> 'ExactCOO': if self._dict_frac is not None: for r, cd in self._dict_frac.items(): for c, f in cd.items(): - rows.append(int(r)); cols.append(int(c)); fracs.append(f) + rows.append(int(r)) + cols.append(int(c)) + fracs.append(f) elif self._num_sparse is not None: coo_num = self._num_sparse.tocoo() coo_den = self._den_sparse.tocoo() for r, c, n, d in zip(coo_num.row, coo_num.col, coo_num.data, coo_den.data): if n != 0: - rows.append(int(r)); cols.append(int(c)); fracs.append(Fraction(int(n), int(d))) + rows.append(int(r)) + cols.append(int(c)) + fracs.append(Fraction(int(n), int(d))) denom = _lcm_list([f.denominator for f in fracs]) if fracs else 1 data = [int(f.numerator) * (denom // f.denominator) for f in fracs] return ExactCOO(rows, cols, data, (self._rows, self._cols), denom) @@ -470,7 +472,9 @@ def to_sparse_pattern(self) -> Tuple[csr_matrix, Dict[int, Dict[int, Fraction]]] rd = {} for c, f in cd.items(): if f != 0: - ai.append(int(r)); aj.append(int(c)); rd[int(c)] = f + ai.append(int(r)) + aj.append(int(c)) + rd[int(c)] = f if rd: row_data[int(r)] = rd pattern = csr_matrix(([1] * len(ai), (ai, aj)), shape=(self._rows, self._cols), dtype=np.int8) @@ -695,7 +699,7 @@ def _eliminate(prd, pivot_val, pivot_col, targets, index): _eliminate(prd, pval, pcol, targets, False) holders = pivcol_holders.get(pcol) if holders: - holders.difference_update(above) # pcol is now cleared from those rows + holders.difference_update(above) # pcol is now cleared from those rows # Final GCD reduction of the pivot rows (insurance; rows are already reduced per step). for key in pivot_keys: @@ -706,14 +710,12 @@ def _eliminate(prd, pivot_val, pivot_col, targets, index): row_data[c] //= row_gcd # Translate results back to original column space, keyed by pivot index (rref_data[i] = pivot i). - original_data = {i: {col_order[sc]: v for sc, v in data[key].items()} - for i, key in enumerate(pivot_keys)} + original_data = {i: {col_order[sc]: v for sc, v in data[key].items()} for i, key in enumerate(pivot_keys)} pivot_cols_original = [col_order[p] for p in pivot_cols_sorted] return original_data, rank, pivot_cols_original - def _nullspace_sparse(matrix: RationalMatrix) -> RationalMatrix: """Compute nullspace using integer RREF with row scaling. @@ -940,7 +942,10 @@ def __init__(self, metas: int, reacs: int): class _WorkRecord: """Mutable state during compression algorithm.""" - def __init__(self, stoich: RationalMatrix, meta_names: List[str], reac_names: List[str], + def __init__(self, + stoich: RationalMatrix, + meta_names: List[str], + reac_names: List[str], bounds: Optional[List[Tuple[float, float]]] = None): rows, cols = stoich.get_row_count(), stoich.get_column_count() self.pre = RationalMatrix.identity(rows) @@ -1346,10 +1351,9 @@ def _small_bound(r): return min(fin) if fin else None def _lam(r): - return 1.0 if ratios[r] is None else float(abs(ratios[r])) # master's own ratio is 1 + return 1.0 if ratios[r] is None else float(abs(ratios[r])) # master's own ratio is 1 - bounded = [(b, r) for r in group for b in [_small_bound(r)] - if b is not None and 0.1 <= _lam(r) <= 10] + bounded = [(b, r) for r in group for b in [_small_bound(r)] if b is not None and 0.1 <= _lam(r) <= 10] keep = min(bounded)[1] if bounded else max(group, key=lambda r: nnz[r]) self._restore_group_scale(work, group, ratios, keep) @@ -1363,8 +1367,7 @@ def _lam(r): return contradicting_removed - def _restore_group_scale(self, work: _WorkRecord, group: List[int], - ratios: List[Optional[Fraction]], keep: int) -> None: + def _restore_group_scale(self, work: _WorkRecord, group: List[int], ratios: List[Optional[Fraction]], keep: int) -> None: """Express a merged group in the units of one of its members. A lump's ratios are fixed but its overall scale is free, and merging into ``group[0]`` can @@ -1376,7 +1379,7 @@ def _restore_group_scale(self, work: _WorkRecord, group: List[int], master = group[0] if keep == master: return - lam = abs(ratios[keep]) # |.| so the reaction keeps its orientation + lam = abs(ratios[keep]) # |.| so the reaction keeps its orientation if lam == 0 or lam == 1: return # v_master = ratios[keep] * v_keep, so re-expressing the lump in v_keep multiplies the @@ -1509,11 +1512,12 @@ def remove_conservation_relations(model) -> None: model.remove_metabolites(dependent_mets) -def compress_cobra_model(model, - methods: Optional[List[Union[str, CompressionMethod]]] = None, - in_place: bool = True, - suppressed_reactions: Set[str] = set(), - protected_reactions: Set[str] = set()) -> CompressionResult: +def compress_cobra_model( + model, + methods: Optional[List[Union[str, CompressionMethod]]] = None, + in_place: bool = True, + suppressed_reactions: Set[str] = set(), + protected_reactions: Set[str] = set()) -> CompressionResult: """ Compress a COBRA model using nullspace-based coupling detection. @@ -1561,7 +1565,8 @@ def compress_cobra_model(model, # Run compression compressor = StoichMatrixCompressor(*compression_methods) bounds = [(float(r.lower_bound), float(r.upper_bound)) for r in model.reactions] - compression_record = compressor.compress(stoich_matrix, metabolite_names, reaction_names, suppressed_reactions, bounds, protected_reactions) + compression_record = compressor.compress(stoich_matrix, metabolite_names, reaction_names, suppressed_reactions, bounds, + protected_reactions) # Apply to model (uses direct manipulation, bypasses solver) reaction_map = _apply_compression_to_model(model, compression_record, reaction_names) @@ -1691,10 +1696,7 @@ def _apply_compression_to_model(model, compression_record, original_reaction_nam # Update stored objective through compression factors obj_dict = getattr(model, '_suppressed_obj', None) if obj_dict is not None: - merged_obj = sum( - obj_dict.pop(original_reaction_names[idx], 0.0) * float(coeff) - for idx, coeff in contributing - ) + merged_obj = sum(obj_dict.pop(original_reaction_names[idx], 0.0) * float(coeff) for idx, coeff in contributing) if merged_obj != 0: obj_dict[main_rxn.id] = merged_obj @@ -1784,12 +1786,12 @@ def stoichmat_coeff_to_fraction(model) -> None: for rxn in model.reactions: for met, coeff in rxn._metabolites.items(): if isinstance(coeff, Fraction): - continue # already exact + continue # already exact elif isinstance(coeff, (float, int)): - rxn._metabolites[met] = float_to_fraction(coeff) # -> Fraction - elif hasattr(coeff, 'p'): # sympy.Rational -> Fraction + rxn._metabolites[met] = float_to_fraction(coeff) # -> Fraction + elif hasattr(coeff, 'p'): # sympy.Rational -> Fraction rxn._metabolites[met] = Fraction(int(coeff.p), int(coeff.q)) - elif hasattr(coeff, 'numerator'): # other Rational -> Fraction + elif hasattr(coeff, 'numerator'): # other Rational -> Fraction rxn._metabolites[met] = Fraction(coeff.numerator, coeff.denominator) else: raise TypeError(f"Unsupported coefficient type: {type(coeff)}") @@ -1848,8 +1850,10 @@ def _expr_to_gpr_string(expr): return expr op, args = expr other = 'or' if op == 'and' else 'and' - parts = [f'({s})' if isinstance(a, tuple) and a[0] == other else s - for a, s in sorted(((a, _expr_to_gpr_string(a)) for a in args), key=lambda p: p[1])] + parts = [ + f'({s})' if isinstance(a, tuple) and a[0] == other else s + for a, s in sorted(((a, _expr_to_gpr_string(a)) for a in args), key=lambda p: p[1]) + ] return f' {op} '.join(parts) @@ -1868,7 +1872,6 @@ def _combine_gprs(gpr_bodies, op): # High-Level Compression API - # Monotone (positive-unate) GPR-rule simplification # # Pipeline: parse -> minimal SOP (DNF + absorption) -> algebraic factoring. @@ -1896,35 +1899,63 @@ def _gpr_parse(s): Accepts both ``and``/``or`` and ``*``/``+`` operators, and is robust to any gene id, including digit-leading or dotted names. """ - toks = list(_gpr_tokenize(s)); pos = 0 - def peek(): return toks[pos] if pos < len(toks) else None + toks = list(_gpr_tokenize(s)) + pos = 0 + + def peek(): + return toks[pos] if pos < len(toks) else None + def eat(): - nonlocal pos; t = toks[pos]; pos += 1; return t + nonlocal pos + t = toks[pos] + pos += 1 + return t + def p_or(): n = [p_and()] - while peek() in ('or', '+'): eat(); n.append(p_and()) + while peek() in ('or', '+'): + eat() + n.append(p_and()) return ('OR', n) if len(n) > 1 else n[0] + def p_and(): n = [p_atom()] - while peek() in ('and', '*'): eat(); n.append(p_atom()) + while peek() in ('and', '*'): + eat() + n.append(p_atom()) return ('AND', n) if len(n) > 1 else n[0] + def p_atom(): - if peek() == '(': eat(); e = p_or(); eat(); return e + if peek() == '(': + eat() + e = p_or() + eat() + return e return ('VAR', eat()) + return p_or() # variable <-> bit mapping (reset per rule via simplify_gpr_string) -_GPR_VMAP = {}; _GPR_VINV = [] +_GPR_VMAP = {} +_GPR_VINV = [] + + def _gpr_bit(v): i = _GPR_VMAP.get(v) if i is None: - i = len(_GPR_VINV); _GPR_VMAP[v] = i; _GPR_VINV.append(v) + i = len(_GPR_VINV) + _GPR_VMAP[v] = i + _GPR_VINV.append(v) return 1 << i + + def _gpr_lits_of(mask): out = [] while mask: - l = mask & -mask; out.append(('VAR', _GPR_VINV[l.bit_length() - 1])); mask ^= l + l = mask & -mask + out.append(('VAR', _GPR_VINV[l.bit_length() - 1])) + mask ^= l return out @@ -1945,11 +1976,12 @@ def _gpr_absorb(cubes): def _gpr_to_dnf(node): t = node[0] - if t == 'VAR': return [_gpr_bit(node[1])] + if t == 'VAR': return [_gpr_bit(node[1])] if t == 'CONST': return [] if not node[1] else [0] if t == 'OR': cov = [] - for ch in node[1]: cov += _gpr_to_dnf(ch) + for ch in node[1]: + cov += _gpr_to_dnf(ch) return _gpr_absorb(cov) if t == 'AND': cov = [0] @@ -1961,8 +1993,10 @@ def _gpr_to_dnf(node): def _gpr_common(cubes): - it = iter(cubes); c = next(it) - for x in it: c &= x + it = iter(cubes) + c = next(it) + for x in it: + c &= x return c @@ -1971,7 +2005,9 @@ def _gpr_lit_counts(F): for c in F: m = c while m: - l = m & -m; cnt[l] = cnt.get(l, 0) + 1; m ^= l + l = m & -m + cnt[l] = cnt.get(l, 0) + 1 + m ^= l return cnt @@ -1991,18 +2027,21 @@ def _gpr_candidate_divisors(F): if len(F) < 2: return [] cnt = _gpr_lit_counts(F) reps = sorted((x for x, n in cnt.items() if n >= 2), key=lambda x: -cnt[x]) - seen = set(); out = [] + seen = set() + out = [] for l in reps: K = tuple(sorted(_gpr_one_kernel(F, l))) if len(K) >= 2 and K not in seen: - seen.add(K); out.append(list(K)) + seen.add(K) + out.append(list(K)) return out def _gpr_divide(F, D): """Exact algebraic division: (Q, R) with D*Q disjoint-union R == F (correctness guaranteed regardless of divisor quality -- a quotient cube is accepted only if D*Q stays inside F).""" - Fs = set(F); quo = None + Fs = set(F) + quo = None for d in D: vd = {c & ~d for c in F if (c & d) == d} quo = vd if quo is None else (quo & vd) @@ -2015,8 +2054,8 @@ def _gpr_divide(F, D): def _gpr_factor(F): F = _gpr_absorb(F) - if not F: return ('CONST', False) - if F == [0]: return ('CONST', True) + if not F: return ('CONST', False) + if F == [0]: return ('CONST', True) if len(F) == 1: lits = _gpr_lits_of(F[0]) return lits[0] if len(lits) == 1 else ('AND', lits) @@ -2044,9 +2083,9 @@ def _gpr_factor(F): def _gpr_est_cubes(node): """Upper bound on DNF cube count (product across ANDs, sum across ORs); cheap, no expansion.""" t = node[0] - if t == 'VAR': return 1 + if t == 'VAR': return 1 if t == 'CONST': return 1 - if t == 'OR': return sum(_gpr_est_cubes(c) for c in node[1]) + if t == 'OR': return sum(_gpr_est_cubes(c) for c in node[1]) if t == 'AND': p = 1 for c in node[1]: @@ -2056,6 +2095,8 @@ def _gpr_est_cubes(node): _GPR_WARN = [] + + def _gpr_factor_auto(node, budget=50000): """Global factoring within budget; AND-split above it. Never splits an OR unless one single OR-block alone exceeds budget (logged as a last resort -- raise the budget to avoid).""" @@ -2083,7 +2124,9 @@ def simplify_gpr_string(rule, budget=50000): """Return a leaf-minimized, boolean-equivalent monotone GPR string ('' passes through).""" if not rule or not rule.strip(): return rule - _GPR_VMAP.clear(); _GPR_VINV.clear(); _GPR_WARN.clear() + _GPR_VMAP.clear() + _GPR_VINV.clear() + _GPR_WARN.clear() return _gpr_to_string(_gpr_factor_auto(_gpr_parse(rule), budget)) @@ -2102,14 +2145,14 @@ def simplify_model_gprs(model, budget=50000): try: new = simplify_gpr_string(s, budget) if new and new != s: - r.gene_reaction_rule = new; nchg += 1 + r.gene_reaction_rule = new + nchg += 1 except Exception as e: logging.warning('gpr_simplify: kept original GPR for %s (%s)' % (r.id, type(e).__name__)) 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 compress_model(model, no_par_compress_reacs=set(), propagate_gpr=False, no_coupled_compress_reacs=set()): """Compress a metabolic model using multiple techniques. Performs blocked reaction removal, conservation relation removal, and @@ -2148,8 +2191,7 @@ def compress_model(model, no_par_compress_reacs=set(), propagate_gpr=False, # 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) parallel_changed = numr > len(reac_map_exp) if parallel_changed: LOG.info(f' Reduced to {len(reac_map_exp)} reactions.') @@ -2169,9 +2211,7 @@ def compress_model(model, no_par_compress_reacs=set(), propagate_gpr=False, # 4. Coupled (expensive — nullspace/RREF) 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) + 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: @@ -2228,8 +2268,7 @@ def compress_model_coupled(model, propagate_gpr=False, protected_reactions=set() for r in model.reactions: r.gene_reaction_rule = '' - result = compress_cobra_model(model, methods=CompressionMethod.standard(), in_place=True, - protected_reactions=protected_reactions) + result = compress_cobra_model(model, methods=CompressionMethod.standard(), in_place=True, protected_reactions=protected_reactions) reaction_map = result.reaction_map # Propagate GPR rules: AND-combine contributing reactions' GPR ASTs @@ -2270,8 +2309,7 @@ def compress_model_parallel(model, protected_rxns=set(), propagate_gpr=False): ub = [float(r.upper_bound) for r in model.reactions] fwd = [1 if (isinf(u) and f > 0 or isinf(l) and f < 0) else 0 for f, l, u in zip(factor, lb, ub)] rev = [1 if (isinf(l) and f > 0 or isinf(u) and f < 0) else 0 for f, l, u in zip(factor, lb, ub)] - inh = [i + 1 if not ((isinf(ub[i]) or ub[i] == 0) and (isinf(lb[i]) or lb[i] == 0)) else 0 - for i in range(len(model.reactions))] + inh = [i + 1 if not ((isinf(ub[i]) or ub[i] == 0) and (isinf(lb[i]) or lb[i] == 0)) else 0 for i in range(len(model.reactions))] # Canonical scale-invariant key per reaction: normalize the stoichiometry row by its first # nonzero coefficient in exact rational arithmetic, so reactions parallel up to any rational @@ -2350,10 +2388,7 @@ def _parallel_key(i): else: scales = [abs(factor[j]) for j in group] total = sum(Fraction(s).limit_denominator(1000) for s in scales) - rational_map[model.reactions[i].id] = { - old_reac_ids[j]: Fraction(abs(factor[j])).limit_denominator(1000) / total - for j in group - } + rational_map[model.reactions[i].id] = {old_reac_ids[j]: Fraction(abs(factor[j])).limit_denominator(1000) / total for j in group} return rational_map diff --git a/straindesign/compute_strain_designs.py b/straindesign/compute_strain_designs.py index 6db62c6..4f5c167 100644 --- a/straindesign/compute_strain_designs.py +++ b/straindesign/compute_strain_designs.py @@ -67,8 +67,9 @@ def _essentials_from_limits(flux_limits): ``flux_limits`` is an FVA result over the module's constrained polytope. Both bounds must share a sign and stay clear of zero, so the reaction carries flux in every point of the module. """ - return {reac_id for reac_id, limits in flux_limits.iterrows() - if np.min(abs(limits)) > _ESSENTIALITY_TOL and np.prod(np.sign(limits)) > 0} + return { + reac_id for reac_id, limits in flux_limits.iterrows() if np.min(abs(limits)) > _ESSENTIALITY_TOL and np.prod(np.sign(limits)) > 0 + } def reduce_model_gprs(model, essential_reacs, gkis, gkos): @@ -265,7 +266,6 @@ def is_gene_essential_to_reaction_ast(reaction, gene_id): @with_suppressed_lp - def compute_strain_designs(model: Model, **kwargs: dict) -> SDSolutions: """Computes strain designs for a user-defined strain design problem @@ -577,20 +577,21 @@ def compute_strain_designs(model: Model, **kwargs: dict) -> SDSolutions: for r in cmp_model.reactions: can_fwd, can_rev = _rev[r.id] if not can_fwd and float(r._upper_bound) > 0.0: - r._upper_bound = min(0.0, float(r._upper_bound)); _n_tight += 1 + r._upper_bound = min(0.0, float(r._upper_bound)) + _n_tight += 1 if not can_rev and float(r._lower_bound) < 0.0: - r._lower_bound = max(0.0, float(r._lower_bound)); _n_tight += 1 - logging.info(' Reversibility pre-tightening fixed %d reaction directions (%.1fs).' - % (_n_tight, time.time() - t0)) + r._lower_bound = max(0.0, float(r._lower_bound)) + _n_tight += 1 + 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() - cmp_mapReac_1 = compress_model(cmp_model, no_par_compress_reacs, - propagate_gpr=True, - no_coupled_compress_reacs=no_coupled_compress_reacs) + cmp_mapReac_1 = compress_model(cmp_model, + no_par_compress_reacs, + propagate_gpr=True, + no_coupled_compress_reacs=no_coupled_compress_reacs) 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) + 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 = [] @@ -627,8 +628,7 @@ def compute_strain_designs(model: Model, **kwargs: dict) -> SDSolutions: len([True for r in cmp_model.reactions if r.gene_reaction_rule]) < num_gpr): num_genes = len(cmp_model.genes) num_gpr = len([True for r in cmp_model.reactions if r.gene_reaction_rule]) - logging.info(' Simplified to ' + str(num_genes) + ' genes and ' + - str(num_gpr) + ' gpr rules (%.1fs).' % (time.time() - t_gpr)) + logging.info(' Simplified to ' + str(num_genes) + ' genes and ' + str(num_gpr) + ' gpr rules (%.1fs).' % (time.time() - t_gpr)) logging.info(' Extending metabolic network with gpr associations.') reac_map = extend_model_gpr(cmp_model, has_gene_names) for i, m in enumerate(sd_modules): @@ -668,11 +668,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) - cmp_mapReac_2 = compress_model(cmp_model, no_par_compress_reacs, -) + cmp_mapReac_2 = compress_model( + cmp_model, + no_par_compress_reacs, + ) 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) + cmp_ko_cost, cmp_ki_cost, cmp_mapReac_2 = compress_ki_ko_cost(cmp_ko_cost, cmp_ki_cost, cmp_mapReac_2) cmp_mapReac = cmp_mapReac_1 + cmp_mapReac_2 logging.info(' Compressed to ' + str(len(cmp_model.reactions)) + ' reactions (%.1fs).' % (time.time() - t0)) else: @@ -684,14 +685,15 @@ def compute_strain_designs(model: Model, **kwargs: dict) -> SDSolutions: t0 = time.time() # Save pre-FVA bounds for dump_preprocessed (bound config experiments) pre_fva_bounds = {r.id: (r.lower_bound, r.upper_bound) for r in cmp_model.reactions} - # Subset-scope the bound-stripping FVA off reactions already at (0,+inf). - # FVA on such a reaction can only ever tighten a genuinely blocked one - # to (0,0); leaving it at (0,+inf) changes no feasible flux (the stoichiometry + # Subset-scope the bound-stripping FVA off reactions already at (0,+inf). + # FVA on such a reaction can only ever tighten a genuinely blocked one + # to (0,0); leaving it at (0,+inf) changes no feasible flux (the stoichiometry # already forces it to 0) and a blocked knockable can never belong to a minimal cut set. - _fva_scope = [r.id for r in cmp_model.reactions - if not (float(r.lower_bound) == 0.0 - and np.isinf(float(r.upper_bound)) - and float(r.upper_bound) > 0)] + _fva_scope = [ + r.id + for r in cmp_model.reactions + if not (float(r.lower_bound) == 0.0 and np.isinf(float(r.upper_bound)) and float(r.upper_bound) > 0) + ] essential_reacs = set() suppress_essential = set() cmp_size1_mcs = [] @@ -700,19 +702,17 @@ def compute_strain_designs(model: Model, **kwargs: dict) -> SDSolutions: # With exactly one classical module, one FVA over the constrained module polytope can serve both # model-bound tightening and module essentiality. This is only sound for a single module: applying # one module's tighter ranges to the shared model could otherwise alter another module's polytope. - fold_module_fva = ( - len(sd_modules) == 1 - and sd_modules[0][MODULE_TYPE] in [SUPPRESS, PROTECT] - and sd_modules[0][INNER_OBJECTIVE] is None - ) + fold_module_fva = (len(sd_modules) == 1 and sd_modules[0][MODULE_TYPE] in [SUPPRESS, PROTECT] and + sd_modules[0][INNER_OBJECTIVE] is None) if fold_module_fva: module = sd_modules[0] fold_scope = sorted(set(_fva_scope) | set(knockable_ids)) - flux_limits = bound_blocked_or_irrevers_fva( - cmp_model, solver=kwargs[SOLVER], constraints=module[CONSTRAINTS], - compress=False, reaction_list=fold_scope) - module_limits = flux_limits.loc[ - [reac_id for reac_id in knockable_ids if reac_id in flux_limits.index]] + flux_limits = bound_blocked_or_irrevers_fva(cmp_model, + solver=kwargs[SOLVER], + constraints=module[CONSTRAINTS], + compress=False, + reaction_list=fold_scope) + module_limits = flux_limits.loc[[reac_id for reac_id in knockable_ids if reac_id in flux_limits.index]] module['fva_bounds'] = module_limits essentials_in_module = _essentials_from_limits(module_limits) if module[MODULE_TYPE] == SUPPRESS: @@ -721,8 +721,7 @@ def compute_strain_designs(model: Model, **kwargs: dict) -> SDSolutions: essential_reacs.update(essentials_in_module) logging.info(' Folded model/module FVA done (%.1fs).' % (time.time() - t0)) else: - bound_blocked_or_irrevers_fva( - cmp_model, solver=kwargs[SOLVER], compress=False, reaction_list=_fva_scope) + bound_blocked_or_irrevers_fva(cmp_model, solver=kwargs[SOLVER], compress=False, reaction_list=_fva_scope) logging.info(' FVA done (%.1fs).' % (time.time() - t0)) # FVA to identify essential reactions and size-1 MCS before building MILP @@ -733,8 +732,11 @@ def compute_strain_designs(model: Model, **kwargs: dict) -> SDSolutions: # module and flows to SDMILP via sd_modules. Scoping to knockable reactions keeps # the LP count down (and only knockable reactions carry z-links to tighten anyway). for module in sd_modules: - flux_limits = fva(cmp_model, solver=kwargs[SOLVER], constraints=module[CONSTRAINTS], - compress=False, reaction_list=knockable_ids) + flux_limits = fva(cmp_model, + solver=kwargs[SOLVER], + constraints=module[CONSTRAINTS], + compress=False, + reaction_list=knockable_ids) module['fva_bounds'] = flux_limits essentials_in_module = _essentials_from_limits(flux_limits) if module[MODULE_TYPE] == SUPPRESS: @@ -814,34 +816,35 @@ def compute_strain_designs(model: Model, **kwargs: dict) -> SDSolutions: import os, pickle as _pickle dump_path = dump_preprocessed with open(dump_path, 'wb') as f: - _pickle.dump({ - 'cmp_model': cmp_model, - 'sd_modules': sd_modules, - 'kwargs_milp': kwargs_milp, - 'kwargs_computation': kwargs_computation, - 'solution_approach': solution_approach, - 'cmp_mapReac': cmp_mapReac, - # Expansion/filtering data - 'uncmp_ko_cost': uncmp_ko_cost, - 'uncmp_ki_cost': uncmp_ki_cost, - 'uncmp_reg_cost': uncmp_reg_cost, - 'orig_model': orig_model, - 'orig_sd_modules': orig_sd_modules, - 'orig_ko_cost': orig_ko_cost, - 'orig_ki_cost': orig_ki_cost, - 'orig_reg_cost': orig_reg_cost, - 'gene_kos': kwargs['gene_kos'], - 'orig_gko_cost': locals().get('orig_gko_cost'), - 'orig_gki_cost': locals().get('orig_gki_cost'), - 'max_cost': kwargs[MAX_COST], - 'cmp_size1_mcs': cmp_size1_mcs, - 'pre_fva_bounds': pre_fva_bounds, - }, f) + _pickle.dump( + { + 'cmp_model': cmp_model, + 'sd_modules': sd_modules, + 'kwargs_milp': kwargs_milp, + 'kwargs_computation': kwargs_computation, + 'solution_approach': solution_approach, + 'cmp_mapReac': cmp_mapReac, + # Expansion/filtering data + 'uncmp_ko_cost': uncmp_ko_cost, + 'uncmp_ki_cost': uncmp_ki_cost, + 'uncmp_reg_cost': uncmp_reg_cost, + 'orig_model': orig_model, + 'orig_sd_modules': orig_sd_modules, + 'orig_ko_cost': orig_ko_cost, + 'orig_ki_cost': orig_ki_cost, + 'orig_reg_cost': orig_reg_cost, + 'gene_kos': kwargs['gene_kos'], + 'orig_gko_cost': locals().get('orig_gko_cost'), + 'orig_gki_cost': locals().get('orig_gki_cost'), + 'max_cost': kwargs[MAX_COST], + 'cmp_size1_mcs': cmp_size1_mcs, + 'pre_fva_bounds': pre_fva_bounds, + }, + f) logging.info('Preprocessed data saved to ' + dump_path) logging.info(' Resume with:') logging.info(' from straindesign import compute_strain_designs_from_preprocessed') - logging.info(" sol = compute_strain_designs_from_preprocessed('%s', seed=42)" % - dump_path.replace('\\', '\\\\')) + logging.info(" sol = compute_strain_designs_from_preprocessed('%s', seed=42)" % dump_path.replace('\\', '\\\\')) # Return early with size-1 MCS only (or empty) setup = deepcopy(cmp_sd_solution.sd_setup) if 'cmp_sd_solution' in dir() else {MODEL_ID: orig_model.id} @@ -885,15 +888,14 @@ def compute_strain_designs(model: Model, **kwargs: dict) -> SDSolutions: if kwargs['gene_kos']: setup.update({GKOCOST: orig_gko_cost, GKICOST: orig_gki_cost}) - sd_solutions = _decompress_solutions( - cmp_sd_solution, cmp_mapReac, cmp_size1_mcs, - kwargs[MAX_COST], uncmp_ko_cost, uncmp_ki_cost, uncmp_reg_cost, - orig_model, setup, kwargs['gene_kos'], - locals().get('orig_gko_cost'), locals().get('orig_gki_cost')) + sd_solutions = _decompress_solutions(cmp_sd_solution, cmp_mapReac, cmp_size1_mcs, kwargs[MAX_COST], uncmp_ko_cost, uncmp_ki_cost, + uncmp_reg_cost, orig_model, setup, kwargs['gene_kos'], + locals().get('orig_gko_cost'), + locals().get('orig_gki_cost')) sd_solutions._cmp_model = cmp_model - logging.info(str(sd_solutions.get_num_materialized()) + ' solutions found' - + (' (lazy, estimated %d total).' % sd_solutions.get_num_sols() - if sd_solutions.is_lazy else '.')) + logging.info( + str(sd_solutions.get_num_materialized()) + ' solutions found' + + (' (lazy, estimated %d total).' % sd_solutions.get_num_sols() if sd_solutions.is_lazy else '.')) return sd_solutions @@ -915,9 +917,8 @@ def postprocess_reg_sd(reg_cost, sd): LAZY_EXPANSION_THRESHOLD = 100_000 -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): +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.""" logging.info(' Decompressing.') @@ -933,9 +934,8 @@ def _decompress_solutions(cmp_sd_solution, cmp_mapReac, cmp_size1_mcs, if estimated > LAZY_EXPANSION_THRESHOLD: 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) + 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) status = cmp_sd_solution.status if status not in [OPTIMAL, TIME_LIMIT_W_SOL] and sd: @@ -995,8 +995,7 @@ def _decompress_solutions(cmp_sd_solution, cmp_mapReac, cmp_size1_mcs, return sd_solutions -def _build_lazy_representatives(cmp_sds, cmp_size1_mcs, cmp_mapReac, max_cost, - uncmp_ko_cost, uncmp_ki_cost, uncmp_reg_cost): +def _build_lazy_representatives(cmp_sds, cmp_size1_mcs, cmp_mapReac, max_cost, uncmp_ko_cost, uncmp_ki_cost, uncmp_reg_cost): """Build one representative expanded solution per compressed group. Returns (sd, group_map, compressed_sd). @@ -1033,9 +1032,7 @@ def _build_lazy_representatives(cmp_sds, cmp_size1_mcs, cmp_mapReac, max_cost, return sd, group_map, compressed_sd -def compute_strain_designs_from_preprocessed(dump, seed=None, solver=None, - solution_approach=None, max_solutions=None, - time_limit=None): +def compute_strain_designs_from_preprocessed(dump, seed=None, solver=None, solution_approach=None, max_solutions=None, time_limit=None): """Load preprocessed model and run MILP solve with optional overrides. Args: @@ -1094,8 +1091,7 @@ def compute_strain_designs_from_preprocessed(dump, seed=None, solver=None, from straindesign.networktools import suppress_lp_context with suppress_lp_context(cmp_model): logging.info('Loading preprocessed data from ' + (dump if isinstance(dump, str) else 'dict input.')) - logging.info(' Seed: %s, Solver: %s, Approach: %s' % ( - kwargs_milp.get(SEED), kwargs_milp.get(SOLVER), sol_approach)) + logging.info(' Seed: %s, Solver: %s, Approach: %s' % (kwargs_milp.get(SEED), kwargs_milp.get(SOLVER), sol_approach)) t0 = time.time() sd_milp = SDMILP(cmp_model, sd_modules, **kwargs_milp) @@ -1116,13 +1112,11 @@ def compute_strain_designs_from_preprocessed(dump, seed=None, solver=None, if gene_kos: setup.update({GKOCOST: orig_gko_cost, GKICOST: orig_gki_cost}) - sd_solutions = _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) + sd_solutions = _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) sd_solutions._cmp_model = cmp_model - logging.info(str(sd_solutions.get_num_materialized()) + ' solutions found' - + (' (lazy, estimated %d total).' % sd_solutions.get_num_sols() - if sd_solutions.is_lazy else '.')) + logging.info( + str(sd_solutions.get_num_materialized()) + ' solutions found' + + (' (lazy, estimated %d total).' % sd_solutions.get_num_sols() if sd_solutions.is_lazy else '.')) return sd_solutions diff --git a/straindesign/cplex_interface.py b/straindesign/cplex_interface.py index 11453d4..e064b1e 100644 --- a/straindesign/cplex_interface.py +++ b/straindesign/cplex_interface.py @@ -95,7 +95,18 @@ class Cplex_MILP_LP(Cplex): A CPLEX MILP/LP interface class. """ - def __init__(self, c=None, A_ineq=None, b_ineq=None, A_eq=None, b_eq=None, lb=None, ub=None, vtype=None, indic_constr=None, seed=None, milp_threads=None): + def __init__(self, + c=None, + A_ineq=None, + b_ineq=None, + A_eq=None, + b_eq=None, + lb=None, + ub=None, + vtype=None, + indic_constr=None, + seed=None, + milp_threads=None): super().__init__() self.objective.set_sense(self.objective.sense.minimize) try: @@ -331,14 +342,12 @@ def set_lp_method(self, method): method: LP_METHOD_AUTO, LP_METHOD_PRIMAL, LP_METHOD_DUAL, or LP_METHOD_BARRIER """ # CPLEX lpmethod: 0=auto, 1=primal, 2=dual, 4=barrier - _map = {LP_METHOD_AUTO: 0, LP_METHOD_PRIMAL: 1, - LP_METHOD_DUAL: 2, LP_METHOD_BARRIER: 4} + _map = {LP_METHOD_AUTO: 0, LP_METHOD_PRIMAL: 1, LP_METHOD_DUAL: 2, LP_METHOD_BARRIER: 4} self.parameters.lpmethod.set(_map.get(method, 0)) def get_lp_method(self): """Return the current LP method as a solver-neutral string.""" - _rmap = {0: LP_METHOD_AUTO, 1: LP_METHOD_PRIMAL, - 2: LP_METHOD_DUAL, 4: LP_METHOD_BARRIER} + _rmap = {0: LP_METHOD_AUTO, 1: LP_METHOD_PRIMAL, 2: LP_METHOD_DUAL, 4: LP_METHOD_BARRIER} return _rmap.get(self.parameters.lpmethod.get(), LP_METHOD_AUTO) def get_basis(self): @@ -362,9 +371,7 @@ def set_basis(self, basis): """ if basis is None: return - self.start.set_start( - col_status=basis['vbasis'], row_status=basis['cbasis'], - col_primal=[], row_primal=[], col_dual=[], row_dual=[]) + self.start.set_start(col_status=basis['vbasis'], row_status=basis['cbasis'], col_primal=[], row_primal=[], col_dual=[], row_dual=[]) def set_time_limit(self, t): """Set the computation time limit (in seconds)""" diff --git a/straindesign/glpk_interface.py b/straindesign/glpk_interface.py index 90ebaab..80a3e1c 100644 --- a/straindesign/glpk_interface.py +++ b/straindesign/glpk_interface.py @@ -387,8 +387,7 @@ def set_lp_method(self, method): # GLPK meth: 1=primal, 2=dual, 3=dual+pricing if method == LP_METHOD_BARRIER: logging.warning('GLPK does not support barrier method, falling back to dual simplex.') - _map = {LP_METHOD_AUTO: 1, LP_METHOD_PRIMAL: 1, - LP_METHOD_DUAL: 2, LP_METHOD_BARRIER: 2} + _map = {LP_METHOD_AUTO: 1, LP_METHOD_PRIMAL: 1, LP_METHOD_DUAL: 2, LP_METHOD_BARRIER: 2} self.lp_params.meth = _map.get(method, 1) def get_lp_method(self): diff --git a/straindesign/gurobi_interface.py b/straindesign/gurobi_interface.py index a9aa730..3087d09 100644 --- a/straindesign/gurobi_interface.py +++ b/straindesign/gurobi_interface.py @@ -31,6 +31,7 @@ # Shared quiet environment — avoids creating a new license session per model _quiet_env = None + def _get_quiet_env(): """Return a shared Gurobi environment with OutputFlag=0.""" global _quiet_env @@ -107,7 +108,18 @@ class Gurobi_MILP_LP(gp.Model): A Gurobi MILP/LP interface class. """ - def __init__(self, c=None, A_ineq=None, b_ineq=None, A_eq=None, b_eq=None, lb=None, ub=None, vtype=None, indic_constr=None, seed=None, milp_threads=None): + def __init__(self, + c=None, + A_ineq=None, + b_ineq=None, + A_eq=None, + b_eq=None, + lb=None, + ub=None, + vtype=None, + indic_constr=None, + seed=None, + milp_threads=None): super().__init__(env=_get_quiet_env()) try: numvars = A_ineq.shape[1] @@ -142,8 +154,8 @@ def __init__(self, c=None, A_ineq=None, b_ineq=None, A_eq=None, b_eq=None, lb=No cols = A_csr.indices[start:end] vals = A_csr.data[start:end] lhs = gp.quicksum(float(val) * x[int(col)] for col, val in zip(cols, vals)) - self.addGenConstrIndicator(x[indic_constr.binv[i]], bool(indic_constr.indicval[i]), - lhs, '=' if indic_constr.sense[i] == 'E' else '<', indic_constr.b[i]) + self.addGenConstrIndicator(x[indic_constr.binv[i]], bool(indic_constr.indicval[i]), lhs, + '=' if indic_constr.sense[i] == 'E' else '<', indic_constr.b[i]) # set parameters self.params.OutputFlag = 0 @@ -379,14 +391,12 @@ def set_lp_method(self, method): Args: method: LP_METHOD_AUTO, LP_METHOD_PRIMAL, LP_METHOD_DUAL, or LP_METHOD_BARRIER """ - _map = {LP_METHOD_AUTO: -1, LP_METHOD_PRIMAL: 0, - LP_METHOD_DUAL: 1, LP_METHOD_BARRIER: 2} + _map = {LP_METHOD_AUTO: -1, LP_METHOD_PRIMAL: 0, LP_METHOD_DUAL: 1, LP_METHOD_BARRIER: 2} self.params.Method = _map.get(method, -1) def get_lp_method(self): """Return the current LP method as a solver-neutral string.""" - _rmap = {-1: LP_METHOD_AUTO, 0: LP_METHOD_PRIMAL, - 1: LP_METHOD_DUAL, 2: LP_METHOD_BARRIER} + _rmap = {-1: LP_METHOD_AUTO, 0: LP_METHOD_PRIMAL, 1: LP_METHOD_DUAL, 2: LP_METHOD_BARRIER} return _rmap.get(self.params.Method, LP_METHOD_AUTO) def get_basis(self): diff --git a/straindesign/lptools.py b/straindesign/lptools.py index 1ba04fe..9da5494 100644 --- a/straindesign/lptools.py +++ b/straindesign/lptools.py @@ -319,7 +319,10 @@ def fva_legacy(model, **kwargs) -> DataFrame: if status not in [OPTIMAL, UNBOUNDED]: logging.info('FVA problem not feasible.') return DataFrame( - {"minimum": [nan] * numr, "maximum": [nan] * numr}, + { + "minimum": [nan] * numr, + "maximum": [nan] * numr + }, index=reaction_ids, ) @@ -328,18 +331,14 @@ def fva_legacy(model, **kwargs) -> DataFrame: x = [nan] * 2 * numr if processes > 1 and numr > 300: - with SDPool(processes, initializer=fva_worker_init, - initargs=(A_ineq, b_ineq, A_eq, b_eq, lb, ub, solver)) as pool: + with SDPool(processes, initializer=fva_worker_init, initargs=(A_ineq, b_ineq, A_eq, b_eq, lb, ub, solver)) as pool: chunk_size = len(reaction_ids) // processes - for i, value in pool.imap_unordered(fva_worker_compute, range(2 * numr), - chunksize=chunk_size): + for i, value in pool.imap_unordered(fva_worker_compute, range(2 * numr), chunksize=chunk_size): x[i] = value elif processes > 1 and numr > 500 and solver == GLPK: - with SDPool(processes, initializer=fva_worker_init_glpk, - initargs=(A_ineq, b_ineq, A_eq, b_eq, lb, ub)) as pool: + with SDPool(processes, initializer=fva_worker_init_glpk, initargs=(A_ineq, b_ineq, A_eq, b_eq, lb, ub)) as pool: chunk_size = len(reaction_ids) // processes - for i, value in pool.imap_unordered(fva_worker_compute_glpk, range(2 * numr), - chunksize=chunk_size): + for i, value in pool.imap_unordered(fva_worker_compute_glpk, range(2 * numr), chunksize=chunk_size): x[i] = value else: fva_worker_init(A_ineq, b_ineq, A_eq, b_eq, lb, ub, solver) @@ -352,8 +351,7 @@ def fva_legacy(model, **kwargs) -> DataFrame: logging.warning(f'FVA: {len(nan_remaining)}/{2*numr} LP solves returned NaN, re-solving.') _BATCH = 50 while nan_remaining: - lp_retry = MILP_LP(A_ineq=A_ineq, b_ineq=b_ineq, A_eq=A_eq, b_eq=b_eq, - lb=lb, ub=ub, solver=solver) + lp_retry = MILP_LP(A_ineq=A_ineq, b_ineq=b_ineq, A_eq=A_eq, b_eq=b_eq, lb=lb, ub=ub, solver=solver) prev_retry = 0 for i in nan_remaining[:_BATCH]: C = idx2c(i, prev_retry) @@ -374,8 +372,7 @@ def fva_legacy(model, **kwargs) -> DataFrame: sig = sign(mod(i, 2) - 0.5) c_vec = [0.0] * numr c_vec[col] = sig - lp_last = MILP_LP(c=c_vec, A_ineq=A_ineq, b_ineq=b_ineq, A_eq=A_eq, b_eq=b_eq, - lb=lb, ub=ub, solver=solver) + lp_last = MILP_LP(c=c_vec, A_ineq=A_ineq, b_ineq=b_ineq, A_eq=A_eq, b_eq=b_eq, lb=lb, ub=ub, solver=solver) x[i] = lp_last.slim_solve() nan_remaining = [i for i in nan_remaining if isnan(x[i])] if nan_remaining: @@ -383,8 +380,10 @@ def fva_legacy(model, **kwargs) -> DataFrame: x = [v if abs(v) >= 1e-11 else 0.0 for v in x] return DataFrame( - {"minimum": [x[i] for i in range(1, 2 * numr, 2)], - "maximum": [-x[i] for i in range(0, 2 * numr, 2)]}, + { + "minimum": [x[i] for i in range(1, 2 * numr, 2)], + "maximum": [-x[i] for i in range(0, 2 * numr, 2)] + }, index=reaction_ids, ) @@ -623,8 +622,7 @@ def slim_fba_via_cmp(model, cmp_model, cmp_map, **kwargs) -> float: if isinstance(obj, str): obj = linexpr2dict(obj, orig_reaction_ids) else: - obj = {r.id: r.objective_coefficient for r in model.reactions - if r.objective_coefficient != 0} + obj = {r.id: r.objective_coefficient for r in model.reactions if r.objective_coefficient != 0} # Trace each objective reaction through compression obj_cmp = {} @@ -673,8 +671,7 @@ def slim_fba_via_cmp(model, cmp_model, cmp_map, **kwargs) -> float: lb = [float(v.lower_bound) for v in cmp_model.reactions] ub = [float(v.upper_bound) for v in cmp_model.reactions] - fba_prob = MILP_LP(c=c, A_ineq=A_ineq, b_ineq=b_ineq, - A_eq=A_eq, b_eq=b_eq, lb=lb, ub=ub, solver=solver) + fba_prob = MILP_LP(c=c, A_ineq=A_ineq, b_ineq=b_ineq, A_eq=A_eq, b_eq=b_eq, lb=lb, ub=ub, solver=solver) opt_cx = fba_prob.slim_solve() if isnan(opt_cx): @@ -906,6 +903,7 @@ def yopt(model, **kwargs) -> Solution: else: status = INFEASIBLE + def expand_fluxes(fluxes_cmp, cmp_map, orig_reaction_ids): """Expand a compressed flux vector to the full (uncompressed) model. @@ -944,6 +942,7 @@ def expand_fluxes(fluxes_cmp, cmp_map, orig_reaction_ids): fluxes[rid] = 0.0 return fluxes + def _make_fix_constraint(axes, ax_idx, ax_type, value): """Create an equality constraint fixing axis ax_idx to value.""" if ax_type == 'rate': @@ -960,8 +959,7 @@ def _optimize_axis(model, ax_idx, axes, ax_type, constraints, solver, sense): if ax_type == 'rate': sol = fba(model, obj=axes[ax_idx][0], constraints=constraints, solver=solver, obj_sense=sense) else: # yield - sol = yopt(model, obj_num=axes[ax_idx][0], obj_den=axes[ax_idx][1], - constraints=constraints, solver=solver, obj_sense=sense) + sol = yopt(model, obj_num=axes[ax_idx][0], obj_den=axes[ax_idx][1], constraints=constraints, solver=solver, obj_sense=sense) return sol @@ -1000,8 +998,7 @@ def _fba_project(obj_dict): # Step 1: Find 4 extremes extremes = [] - for coeff, sense in [(ax0_coeff, 'maximize'), (ax0_coeff, 'minimize'), - (ax1_coeff, 'maximize'), (ax1_coeff, 'minimize')]: + for coeff, sense in [(ax0_coeff, 'maximize'), (ax0_coeff, 'minimize'), (ax1_coeff, 'maximize'), (ax1_coeff, 'minimize')]: sol = fba(model, obj=coeff, constraints=constraints, solver=solver, obj_sense=sense) if sol.status == OPTIMAL: x0 = sum(c * sol.fluxes.get(r, 0) for r, c in ax0_coeff.items()) @@ -1027,16 +1024,13 @@ def _fba_project(obj_dict): unique_pts.sort(key=lambda p: arctan2(p[1] - cy, p[0] - cx)) # Step 4: Recursive edge refinement - diameter = max( - ((a[0]-b[0])**2 + (a[1]-b[1])**2)**0.5 - for a in unique_pts for b in unique_pts - ) + diameter = max(((a[0] - b[0])**2 + (a[1] - b[1])**2)**0.5 for a in unique_pts for b in unique_pts) min_edge = 1e-10 * diameter if diameter > 0 else 1e-10 def _refine(vi, vj, depth): if depth > 50: return [vi] - edge_len = ((vi[0]-vj[0])**2 + (vi[1]-vj[1])**2)**0.5 + edge_len = ((vi[0] - vj[0])**2 + (vi[1] - vj[1])**2)**0.5 if edge_len < min_edge: return [vi] # Outward normal (perpendicular to edge, pointing outward from centroid) @@ -1063,8 +1057,8 @@ def _refine(vi, vj, depth): dist /= norm_len if dist > tol and edge_len > min_edge: # Check p_new is not a duplicate of vi or vj - if ((abs(p_new[0]-vi[0]) < tol and abs(p_new[1]-vi[1]) < tol) or - (abs(p_new[0]-vj[0]) < tol and abs(p_new[1]-vj[1]) < tol)): + if ((abs(p_new[0] - vi[0]) < tol and abs(p_new[1] - vi[1]) < tol) or + (abs(p_new[0] - vj[0]) < tol and abs(p_new[1] - vj[1]) < tol)): return [vi] left = _refine(vi, p_new, depth + 1) right = _refine(p_new, vj, depth + 1) @@ -1091,6 +1085,7 @@ def _trace_boundary_adaptive(model, axes, ax_types, constraints, solver, max_dep Returns (upper_boundary, lower_boundary) as sorted lists of (x, y) tuples. Uses recursive midpoint refinement where linear interpolation error exceeds tolerance. """ + def _fix_and_opt(x_val, sense): constr = constraints.copy() constr.append(_make_fix_constraint(axes, 0, ax_types[0], x_val)) @@ -1100,11 +1095,8 @@ def _fix_and_opt(x_val, sense): return nan # Step 1: axis-0 range endpoints (already known from val_limits, but we need y values) - x_min, x_max = ceil_dec( - _optimize_axis(model, 0, axes, ax_types[0], constraints, solver, 'minimize').objective_value, 8 - ), floor_dec( - _optimize_axis(model, 0, axes, ax_types[0], constraints, solver, 'maximize').objective_value, 8 - ) + x_min, x_max = ceil_dec(_optimize_axis(model, 0, axes, ax_types[0], constraints, solver, 'minimize').objective_value, + 8), floor_dec(_optimize_axis(model, 0, axes, ax_types[0], constraints, solver, 'maximize').objective_value, 8) # Step 2: y values at endpoints y_min_at_xmin = _fix_and_opt(x_min, 'minimize') @@ -1168,10 +1160,7 @@ def _trace_polytope_3d_rate(model, axes, constraints, solver): tol = 1e-8 def _project(sol): - return tuple( - ceil_dec(sum(c * sol.fluxes.get(r, 0) for r, c in coeffs[j].items()), 9) - for j in range(3) - ) + return tuple(ceil_dec(sum(c * sol.fluxes.get(r, 0) for r, c in coeffs[j].items()), 9) for j in range(3)) def _is_dup(pt, pts): return any(all(abs(pt[k] - q[k]) < tol for k in range(3)) for q in pts) @@ -1248,10 +1237,8 @@ def _hull_face_polygons(hull): # Build orthonormal basis in the face plane ref = array([1, 0, 0]) if abs(normal[0]) < 0.9 else array([0, 1, 0]) u = ref - normal * normal.dot(ref) - u = u / (u.dot(u) ** 0.5) - v = array([normal[1]*u[2] - normal[2]*u[1], - normal[2]*u[0] - normal[0]*u[2], - normal[0]*u[1] - normal[1]*u[0]]) + u = u / (u.dot(u)**0.5) + v = array([normal[1] * u[2] - normal[2] * u[1], normal[2] * u[0] - normal[0] * u[2], normal[0] * u[1] - normal[1] * u[0]]) # Sort by angle in face plane angles = [arctan2((pt - centroid).dot(v), (pt - centroid).dot(u)) for pt in pts] order = sorted(range(len(verts)), key=lambda i: angles[i]) @@ -1420,14 +1407,13 @@ def plot_flux_space(model, axes, **kwargs) -> Tuple[list, list, list]: variable contains information about which datapoints need to be connected in triangles to render a closed surface. The last variable contains the matplotlib object. """ - + cmp_model = kwargs.pop('cmp_model', None) cmp_map = kwargs.pop('cmp_map', None) _orig_axes = None # store original axis names for labelling if cmp_model is not None and cmp_map is not None: - from straindesign.networktools import (resolve_gene_constraints, - compress_constraints, _build_cmp_reverse_map) + from straindesign.networktools import (resolve_gene_constraints, compress_constraints, _build_cmp_reverse_map) # Resolve gene constraints on the original model, then compress if CONSTRAINTS in kwargs and kwargs[CONSTRAINTS]: kwargs[CONSTRAINTS] = resolve_gene_constraints(model, kwargs[CONSTRAINTS]) @@ -1455,7 +1441,7 @@ def plot_flux_space(model, axes, **kwargs) -> Tuple[list, list, list]: axes[i] = [reverse.get(a, a) if isinstance(a, str) else a for a in ax] # Switch to compressed model model = cmp_model - + reaction_ids = model.reactions.list_attr("id") if CONSTRAINTS in kwargs: @@ -1592,8 +1578,7 @@ def plot_flux_space(model, axes, **kwargs) -> Tuple[list, list, list]: vertices = _trace_polygon_rate_rate(model, axes, kwargs[CONSTRAINTS], solver) else: adapt_depth = max(5, int(log2(max(points, 2)))) - upper, lower = _trace_boundary_adaptive( - model, axes, ax_type, kwargs[CONSTRAINTS], solver, max_depth=adapt_depth) + upper, lower = _trace_boundary_adaptive(model, axes, ax_type, kwargs[CONSTRAINTS], solver, max_depth=adapt_depth) # Build polygon from upper (left-to-right) + reversed lower (right-to-left) if upper and lower: vertices = upper + list(reversed(lower)) @@ -1634,9 +1619,7 @@ def plot_flux_space(model, axes, **kwargs) -> Tuple[list, list, list]: x0, y0 = unique_verts[0] dx, dy = unique_verts[1][0] - x0, unique_verts[1][1] - y0 span = max(abs(dx), abs(dy), 1e-12) - is_collinear = all( - abs((v[0] - x0) * dy - (v[1] - y0) * dx) / span < 1e-6 - for v in unique_verts[2:]) + is_collinear = all(abs((v[0] - x0) * dy - (v[1] - y0) * dx) / span < 1e-6 for v in unique_verts[2:]) if is_collinear and len(unique_verts) <= 1: # Collapsed to a point @@ -1786,8 +1769,8 @@ def plot_flux_space(model, axes, **kwargs) -> Tuple[list, list, list]: elif n_yields == 1: # 1 yield + 2 rate: slice along yield, trace rate-rate polygon per slice - datapoints, triang, slice_outlines = _trace_3d_slice_polygon( - model, axes, ax_type, val_limits, kwargs[CONSTRAINTS], solver, points) + datapoints, triang, slice_outlines = _trace_3d_slice_polygon(model, axes, ax_type, val_limits, kwargs[CONSTRAINTS], solver, + points) else: # 2+ yields: grid-based scanning (fallback) @@ -1807,8 +1790,7 @@ def plot_flux_space(model, axes, **kwargs) -> Tuple[list, list, list]: if abs(y_hi - y_lo) < 1e-10: y_space = [y_lo] else: - n_pts = max(3, int(points * abs(y_hi - y_lo) / max(1e-10, - max(abs(val_limits[1][1] - val_limits[1][0]), 1e-10)))) + n_pts = max(3, int(points * abs(y_hi - y_lo) / max(1e-10, max(abs(val_limits[1][1] - val_limits[1][0]), 1e-10)))) n_pts = min(n_pts, points) y_space = linspace(y_lo, y_hi, n_pts).tolist() upper_slice = [] @@ -1843,8 +1825,7 @@ def plot_flux_space(model, axes, **kwargs) -> Tuple[list, list, list]: triang = [] for s in range(len(slices) - 1): _triangulate_strips(datapoints_top[s], datapoints_top[s + 1], datapoints, triang) - _triangulate_strips(datapoints_bottom[s], datapoints_bottom[s + 1], datapoints, triang, - flip_winding=True) + _triangulate_strips(datapoints_bottom[s], datapoints_bottom[s + 1], datapoints, triang, flip_winding=True) front_top = [t[0] for t in datapoints_top] front_bot = [b[0] for b in datapoints_bottom] _triangulate_strips(front_top, front_bot, datapoints, triang, flip_winding=True) @@ -1852,8 +1833,7 @@ def plot_flux_space(model, axes, **kwargs) -> Tuple[list, list, list]: back_bot = [b[-1] for b in datapoints_bottom] _triangulate_strips(back_top, back_bot, datapoints, triang) _triangulate_strips(datapoints_top[0], datapoints_bottom[0], datapoints, triang) - _triangulate_strips(datapoints_top[-1], datapoints_bottom[-1], datapoints, triang, - flip_winding=True) + _triangulate_strips(datapoints_top[-1], datapoints_bottom[-1], datapoints, triang, flip_winding=True) if not datapoints: raise Exception('No feasible points found. Problem may be infeasible.') @@ -1882,10 +1862,8 @@ def _normal_color(face_pts): return 0.5 e1 = pts[1] - pts[0] e2 = pts[2] - pts[0] - normal = array([e1[1]*e2[2] - e1[2]*e2[1], - e1[2]*e2[0] - e1[0]*e2[2], - e1[0]*e2[1] - e1[1]*e2[0]]) - length = (normal.dot(normal)) ** 0.5 + normal = array([e1[1] * e2[2] - e1[2] * e2[1], e1[2] * e2[0] - e1[0] * e2[2], e1[0] * e2[1] - e1[1] * e2[0]]) + length = (normal.dot(normal))**0.5 if length > 0: normal = normal / length # Map normal direction to scalar: use spherical angles @@ -1899,8 +1877,7 @@ def _normal_color(face_pts): rng = mx - mn if mx > mn else 1 face_colors = plt.get_cmap(cmap)([(c - mn) / rng for c in color_vals]) face_colors[:, 3] = 1.0 - collection = Poly3DCollection(poly_verts, facecolors=face_colors, - edgecolors='black', linewidths=1.0) + collection = Poly3DCollection(poly_verts, facecolors=face_colors, edgecolors='black', linewidths=1.0) ax3.add_collection3d(collection) plot1 = collection elif triang: @@ -1911,16 +1888,14 @@ def _normal_color(face_pts): rng = mx - mn if mx > mn else 1 face_colors = plt.get_cmap(cmap)([(c - mn) / rng for c in color_vals]) face_colors[:, 3] = 1.0 - collection = Poly3DCollection(tri_verts, facecolors=face_colors, - edgecolors='none', linewidths=0) + collection = Poly3DCollection(tri_verts, facecolors=face_colors, edgecolors='none', linewidths=0) ax3.add_collection3d(collection) # Draw slice polygon outlines (exact contour at each yield level) if slice_outlines: for idx_list in slice_outlines: pts = array([datapoints[i] for i in idx_list]) pts = array(list(pts) + [pts[0]]) # close the loop - ax3.plot(pts[:, 0], pts[:, 1], pts[:, 2], - color='gray', linewidth=0.4) + ax3.plot(pts[:, 0], pts[:, 1], pts[:, 2], color='gray', linewidth=0.4) plot1 = collection else: plot1 = ax3.scatter(x, y, z, s=20) diff --git a/straindesign/networktools.py b/straindesign/networktools.py index 060dc88..0c5027e 100644 --- a/straindesign/networktools.py +++ b/straindesign/networktools.py @@ -53,6 +53,7 @@ def _sb_noop(self, lb, ub): # -- Cobra-level suppression replacements ------------------------------------ + def _suppressed_set_id(self, value): """Bypass solver variable rename: write _id and update DictList index.""" old_id = self._id @@ -158,7 +159,7 @@ def _suppressed_copy(model): extension's ``add_metabolites`` calls run without building (or pushing constraints into) a live solver. """ - iface = model.solver.interface # captured before stubbing + iface = model.solver.interface # captured before stubbing saved = model._solver orig_copy = next(o for cls, attr, o in _ORIG_COBRA if cls is Model and attr == 'copy') try: @@ -169,6 +170,7 @@ def _suppressed_copy(model): new._solver = _CarrierSolver(iface) return new + _ORIG_CONTAINER_GETITEM = None # saved Container.__getitem__ @@ -272,10 +274,10 @@ def _suppressed_add_metabolites(self, metabolite_list): # -- Saved originals (None = not suppressed) ---------------------------------- -_ORIG_SLC = None # (cls, method) for Constraint.set_linear_coefficients -_ORIG_SB = None # (cls, method) for Variable.set_bounds +_ORIG_SLC = None # (cls, method) for Constraint.set_linear_coefficients +_ORIG_SB = None # (cls, method) for Variable.set_bounds _ORIG_OSLC = None # (cls, method) for Objective.set_linear_coefficients -_ORIG_COBRA = [] # list of (cls, attr_name, original) for cobra-level patches +_ORIG_COBRA = [] # list of (cls, attr_name, original) for cobra-level patches def _suppress_lp_updates(model): @@ -438,8 +440,7 @@ def suppress_lp_context(model): del model._suppressed_obj if current_ids != _pre_ids: if model.groups: - kept = {c.id for c in - list(model.reactions) + list(model.metabolites) + list(model.genes)} + kept = {c.id for c in list(model.reactions) + list(model.metabolites) + list(model.genes)} for grp in model.groups: stale = [m for m in grp.members if m.id not in kept] if stale: @@ -461,11 +462,13 @@ def with_suppressed_lp(func): The first positional argument must be a cobra Model. """ + @wraps(func) def wrapper(*args, **kwargs): model = args[0] with suppress_lp_context(model): return func(*args, **kwargs) + return wrapper @@ -473,6 +476,7 @@ def wrapper(*args, **kwargs): # I/O Suppression # ============================================================================= + @contextmanager def _silent_io(): """Suppress stdout, stderr and logging.""" @@ -608,6 +612,7 @@ def gene_kos_to_constraints(model, gene_kos): return [[{r: 1}, '=', 0] for r in sorted(knocked_out_reactions)] + def _build_cmp_reverse_map(cmp_map): """Build reverse lookup: original reaction ID -> final compressed ID. @@ -668,11 +673,10 @@ def compress_constraints(constraints, cmp_mapReac): coeff_dict = c[0] lumped = [k for k in coeff_dict if k in old_reac_val] if lumped: - coeff_dict[new_reac] = sum( - coeff_dict.pop(k) * old_reac_val[k] for k in lumped - ) + coeff_dict[new_reac] = sum(coeff_dict.pop(k) * old_reac_val[k] for k in lumped) return constraints + def resolve_gene_constraints(model, constraints): """Scan constraints for gene IDs/names and replace with reaction constraints. @@ -1358,8 +1362,7 @@ def filter_sd_maxcost(sd, max_cost, kocost, kicost): # 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. 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] + 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] # sort strain designs by intervention costs [s.update({'**cost**': c}) for s, c in zip(sd, costs)] @@ -1425,14 +1428,12 @@ def bound_blocked_or_irrevers_fva(model, **kwargs): from operator import attrgetter as _attrgetter from fractions import Fraction as _Fraction import math as _math -from cobra.io.dict import (_metabolite_to_dict, _metabolite_from_dict, - _gene_to_dict, gene_from_dict, _fix_type, _update_optional, - _OPTIONAL_REACTION_ATTRIBUTES, _ORDERED_OPTIONAL_REACTION_KEYS, - _OPTIONAL_MODEL_ATTRIBUTES, _ORDERED_OPTIONAL_MODEL_KEYS) +from cobra.io.dict import (_metabolite_to_dict, _metabolite_from_dict, _gene_to_dict, gene_from_dict, _fix_type, _update_optional, + _OPTIONAL_REACTION_ATTRIBUTES, _ORDERED_OPTIONAL_REACTION_KEYS, _OPTIONAL_MODEL_ATTRIBUTES, + _ORDERED_OPTIONAL_MODEL_KEYS) from cobra.util.solver import set_objective as _set_objective -_INF_TOKENS = {'inf': _math.inf, 'infinity': _math.inf, '+inf': _math.inf, - '-inf': -_math.inf, '-infinity': -_math.inf} +_INF_TOKENS = {'inf': _math.inf, 'infinity': _math.inf, '+inf': _math.inf, '-inf': -_math.inf, '-infinity': -_math.inf} def _num_to_json(v): @@ -1442,13 +1443,13 @@ def _num_to_json(v): if isinstance(v, bool): return v if isinstance(v, _Fraction): - return str(v) # exact: '1/3', '10', '-1/2' + return str(v) # exact: '1/3', '10', '-1/2' num, den = getattr(v, 'numerator', None), getattr(v, 'denominator', None) if num is not None and den is not None and not isinstance(v, int): - return '%d/%d' % (int(num), int(den)) # sympy Rational etc. - v = _fix_type(v) # numpy -> native python + return '%d/%d' % (int(num), int(den)) # sympy Rational etc. + v = _fix_type(v) # numpy -> native python if isinstance(v, float) and (_math.isinf(v) or _math.isnan(v)): - return str(v) # 'inf', '-inf', 'nan' + return str(v) # 'inf', '-inf', 'nan' return v @@ -1461,7 +1462,7 @@ def _num_from_json(v): return _INF_TOKENS[low] if low == 'nan': return _math.nan - return _Fraction(v) # 'num/den' or integer string + return _Fraction(v) # 'num/den' or integer string return v @@ -1473,8 +1474,7 @@ def _reaction_to_json(reaction, obj_coeff=0): new['id'] = reaction.id new['name'] = reaction.name new['metabolites'] = _OrderedDict( - (str(met), _num_to_json(reaction.metabolites[met])) - for met in sorted(reaction.metabolites, key=_attrgetter('id'))) + (str(met), _num_to_json(reaction.metabolites[met])) for met in sorted(reaction.metabolites, key=_attrgetter('id'))) new['lower_bound'] = _num_to_json(reaction.lower_bound) new['upper_bound'] = _num_to_json(reaction.upper_bound) new['gene_reaction_rule'] = reaction.gene_reaction_rule @@ -1495,9 +1495,8 @@ def _reaction_from_json(reaction, model): if k in {'objective_coefficient', 'reversibility', 'reaction'}: continue elif k == 'metabolites': - new_reaction.add_metabolites(_OrderedDict( - (model.metabolites.get_by_id(str(met)), _num_from_json(coeff)) - for met, coeff in v.items())) + new_reaction.add_metabolites( + _OrderedDict((model.metabolites.get_by_id(str(met)), _num_from_json(coeff)) for met, coeff in v.items())) elif k in {'lower_bound', 'upper_bound'}: setattr(new_reaction, k, _num_from_json(v)) else: @@ -1552,10 +1551,8 @@ def model_from_dict(obj): model.add_metabolites([_metabolite_from_dict(m) for m in obj['metabolites']]) model.genes.extend([gene_from_dict(g) for g in obj['genes']]) model.add_reactions([_reaction_from_json(r, model) for r in obj['reactions']]) - objective_reactions = [r for r in obj['reactions'] - if _num_from_json(r.get('objective_coefficient', 0)) != 0] - coefficients = {model.reactions.get_by_id(r['id']): _num_from_json(r['objective_coefficient']) - for r in objective_reactions} + objective_reactions = [r for r in obj['reactions'] if _num_from_json(r.get('objective_coefficient', 0)) != 0] + coefficients = {model.reactions.get_by_id(r['id']): _num_from_json(r['objective_coefficient']) for r in objective_reactions} _set_objective(model, coefficients) for k, v in obj.items(): if k in {'id', 'name', 'notes', 'compartments', 'annotation'}: diff --git a/straindesign/parse_constr.py b/straindesign/parse_constr.py index ce628e2..81843f7 100644 --- a/straindesign/parse_constr.py +++ b/straindesign/parse_constr.py @@ -327,4 +327,3 @@ def linexprdict2str(D): return expr else: return "" - diff --git a/straindesign/scip_interface.py b/straindesign/scip_interface.py index 8798918..7a7fc8c 100644 --- a/straindesign/scip_interface.py +++ b/straindesign/scip_interface.py @@ -94,7 +94,18 @@ class SCIP_MILP(pso.Model): A SCIP MILP interface class. """ - def __init__(self, c=None, A_ineq=None, b_ineq=None, A_eq=None, b_eq=None, lb=None, ub=None, vtype=None, indic_constr=None, seed=None, milp_threads=None): + def __init__(self, + c=None, + A_ineq=None, + b_ineq=None, + A_eq=None, + b_eq=None, + lb=None, + ub=None, + vtype=None, + indic_constr=None, + seed=None, + milp_threads=None): super().__init__() # uncomment to forward SCIP output to python terminal # self.redirectOutput() @@ -371,16 +382,14 @@ def set_lp_method(self, method): method: LP_METHOD_AUTO, LP_METHOD_PRIMAL, LP_METHOD_DUAL, or LP_METHOD_BARRIER """ # SCIP lp/initalgorithm + lp/resolvealgorithm: 's'=auto, 'p'=primal, 'd'=dual, 'b'=barrier - _map = {LP_METHOD_AUTO: 's', LP_METHOD_PRIMAL: 'p', - LP_METHOD_DUAL: 'd', LP_METHOD_BARRIER: 'b'} + _map = {LP_METHOD_AUTO: 's', LP_METHOD_PRIMAL: 'p', LP_METHOD_DUAL: 'd', LP_METHOD_BARRIER: 'b'} algo = _map.get(method, 's') self.setParam('lp/initalgorithm', algo) self.setParam('lp/resolvealgorithm', algo) def get_lp_method(self): """Return the current LP method as a solver-neutral string.""" - _rmap = {'s': LP_METHOD_AUTO, 'p': LP_METHOD_PRIMAL, - 'd': LP_METHOD_DUAL, 'b': LP_METHOD_BARRIER} + _rmap = {'s': LP_METHOD_AUTO, 'p': LP_METHOD_PRIMAL, 'd': LP_METHOD_DUAL, 'b': LP_METHOD_BARRIER} return _rmap.get(self.getParam('lp/initalgorithm'), LP_METHOD_AUTO) def get_basis(self): diff --git a/straindesign/solver_interface.py b/straindesign/solver_interface.py index b89a074..0462cec 100644 --- a/straindesign/solver_interface.py +++ b/straindesign/solver_interface.py @@ -102,7 +102,8 @@ class MILP_LP(object): def __init__(self, **kwargs): allowed_keys = { - 'c', 'A_ineq', 'b_ineq', 'A_eq', 'b_eq', 'lb', 'ub', 'vtype', 'indic_constr', 'M', SOLVER, 'skip_checks', 'tlim', SEED, MILP_THREADS + 'c', 'A_ineq', 'b_ineq', 'A_eq', 'b_eq', 'lb', 'ub', 'vtype', 'indic_constr', 'M', SOLVER, 'skip_checks', 'tlim', SEED, + MILP_THREADS } # set all keys passed in kwargs for key, value in kwargs.items(): diff --git a/straindesign/speedy_fva.py b/straindesign/speedy_fva.py index 50be8f0..b4190f2 100644 --- a/straindesign/speedy_fva.py +++ b/straindesign/speedy_fva.py @@ -47,8 +47,12 @@ from cobra import Configuration from straindesign.lptools import ( - select_solver, idx2c, fva_worker_init, fva_worker_compute, - fva_worker_init_glpk, fva_worker_compute_glpk, + select_solver, + idx2c, + fva_worker_init, + fva_worker_compute, + fva_worker_init_glpk, + fva_worker_compute_glpk, ) from straindesign.solver_interface import MILP_LP from straindesign.pool import SDPool @@ -56,21 +60,24 @@ from straindesign.names import CONSTRAINTS, SOLVER, OPTIMAL, UNBOUNDED, INFEASIBLE, GLPK, LP_METHOD_DUAL from straindesign.networktools import suppress_lp_context from straindesign.compression import ( - compress_cobra_model, CompressionMethod, remove_conservation_relations, - stoichmat_coeff_to_fraction, stoichmat_coeff2float, remove_blocked_reactions, + compress_cobra_model, + CompressionMethod, + remove_conservation_relations, + stoichmat_coeff_to_fraction, + stoichmat_coeff2float, + remove_blocked_reactions, ) - # Tunable threshold at which parallel FVA kicks in. Is compared against the number of # LPs that are left to be solved by Phase 2. Empirically, parallel is only worthwhile # at high LP counts. _PARALLEL_PHASE2_MIN = 4000 - # --------------------------------------------------------------------------- # Compression helpers # --------------------------------------------------------------------------- + def _compress_for_fva(model): """Copy and compress model for FVA (single-pass coupled compression + conservation removal). @@ -95,9 +102,7 @@ def _compress_for_fva(model): # Single-pass coupled compression (NULLSPACE only, no RECURSIVE iteration) for r in cmp_model.reactions: r.gene_reaction_rule = '' - result = compress_cobra_model(cmp_model, - methods=[CompressionMethod.NULLSPACE], - in_place=True) + result = compress_cobra_model(cmp_model, methods=[CompressionMethod.NULLSPACE], in_place=True) rmap = result.reaction_map if len(rmap) < n_before: cmp_maps.append(rmap) @@ -120,9 +125,7 @@ def _map_constraints(parsed_constraints, cmp_maps, cmp_reaction_ids): coeff_dict = constraint[0] lumped = [k for k in coeff_dict if k in orig_map] if lumped: - coeff_dict[cmp_id] = sum( - coeff_dict.pop(k) * float(orig_map[k]) for k in lumped - ) + coeff_dict[cmp_id] = sum(coeff_dict.pop(k) * float(orig_map[k]) for k in lumped) # Remove references to reactions not in compressed model (e.g., blocked) for constraint in parsed_constraints: for k in list(constraint[0].keys()): @@ -166,8 +169,10 @@ def _expand_fva(fva_cmp, cmp_maps, orig_reaction_ids): result_max[rxn_id] = 0.0 return DataFrame( - {"minimum": [result_min[r] for r in orig_reaction_ids], - "maximum": [result_max[r] for r in orig_reaction_ids]}, + { + "minimum": [result_min[r] for r in orig_reaction_ids], + "maximum": [result_max[r] for r in orig_reaction_ids] + }, index=orig_reaction_ids, ) @@ -176,6 +181,7 @@ def _expand_fva(fva_cmp, cmp_maps, orig_reaction_ids): # Global scan LP helper # --------------------------------------------------------------------------- + def _build_abssum_lp(S_eq, b_eq, A_ineq, b_ineq, lb, ub, solver, BIG=1000.0): """Build LP for min sum(|x|) via variable splitting. @@ -193,9 +199,9 @@ def _build_abssum_lp(S_eq, b_eq, A_ineq, b_ineq, lb, ub, solver, BIG=1000.0): tol = 1e-9 # Classify reactions - fwd = lb >= -tol # forward-only or fixed - bwd = ub <= tol # backward-only or fixed - rev = (~fwd) & (~bwd) # truly reversible + fwd = lb >= -tol # forward-only or fixed + bwd = ub <= tol # backward-only or fixed + rev = (~fwd) & (~bwd) # truly reversible n_rev = int(rev.sum()) rev_idx = np.where(rev)[0] @@ -210,9 +216,9 @@ def _build_abssum_lp(S_eq, b_eq, A_ineq, b_ineq, lb, ub, solver, BIG=1000.0): elif bwd[j]: c[j] = -1.0 # |x_j| = -x_j else: - c[j] = 1.0 # |x_j| = x_j (fwd or fixed) + c[j] = 1.0 # |x_j| = x_j (fwd or fixed) for k, j in enumerate(rev_idx): - c[n + k] = 1.0 # p_k + c[n + k] = 1.0 # p_k c[n + n_rev + k] = 1.0 # n_k # Equality constraints: original S*x = 0 (+ extras) + splitting equalities @@ -233,8 +239,7 @@ def _build_abssum_lp(S_eq, b_eq, A_ineq, b_ineq, lb, ub, solver, BIG=1000.0): rows_all = np.concatenate([rows_split, rows_split, rows_split]) cols_all = np.concatenate([cols_x, cols_p, cols_n]) data_all = np.concatenate([data_x, data_p, data_n]) - A_split = sparse.csr_matrix((data_all, (rows_all, cols_all)), - shape=(n_rev, n_ext)) + A_split = sparse.csr_matrix((data_all, (rows_all, cols_all)), shape=(n_rev, n_ext)) # Extend original equalities to n_ext columns S_ext = sparse.hstack([S_eq, sparse.csr_matrix((S_eq.shape[0], 2 * n_rev))]) @@ -270,9 +275,14 @@ def _build_abssum_lp(S_eq, b_eq, A_ineq, b_ineq, lb, ub, solver, BIG=1000.0): lb_ext[n + n_rev + k] = 0.0 ub_ext[n + n_rev + k] = min(-lb[j], BIG) - lp = MILP_LP(c=c.tolist(), A_ineq=A_ineq_ext, b_ineq=list(b_ineq), - A_eq=A_eq_full, b_eq=b_eq_full, - lb=lb_ext.tolist(), ub=ub_ext.tolist(), solver=solver) + lp = MILP_LP(c=c.tolist(), + A_ineq=A_ineq_ext, + b_ineq=list(b_ineq), + A_eq=A_eq_full, + b_eq=b_eq_full, + lb=lb_ext.tolist(), + ub=ub_ext.tolist(), + solver=solver) return lp, n @@ -280,6 +290,7 @@ def _build_abssum_lp(S_eq, b_eq, A_ineq, b_ineq, lb, ub, solver, BIG=1000.0): # Main entry point # --------------------------------------------------------------------------- + def speedy_fva(model, **kwargs): """Accelerated FVA using global scan LPs and KKT-based optimality propagation. @@ -345,10 +356,8 @@ def speedy_fva(model, **kwargs): kwargs[CONSTRAINTS] = resolve_gene_constraints(orig_model, kwargs[CONSTRAINTS]) kwargs[CONSTRAINTS] = parse_constraints(kwargs[CONSTRAINTS], orig_reaction_ids) if cmp_maps: - kwargs[CONSTRAINTS] = _map_constraints( - kwargs[CONSTRAINTS], cmp_maps, reaction_ids) - A_ineq_extra, b_ineq_extra, A_eq_extra, b_eq_extra = lineqlist2mat( - kwargs[CONSTRAINTS], reaction_ids) + kwargs[CONSTRAINTS] = _map_constraints(kwargs[CONSTRAINTS], cmp_maps, reaction_ids) + A_ineq_extra, b_ineq_extra, A_eq_extra, b_eq_extra = lineqlist2mat(kwargs[CONSTRAINTS], reaction_ids) if SOLVER not in kwargs: kwargs[SOLVER] = None @@ -376,14 +385,16 @@ def speedy_fva(model, **kwargs): lb = np.array([v.lower_bound for v in model.reactions], dtype=np.float64) ub = np.array([v.upper_bound for v in model.reactions], dtype=np.float64) - lp = MILP_LP(A_ineq=A_ineq, b_ineq=b_ineq, A_eq=A_eq, b_eq=b_eq, - lb=lb.tolist(), ub=ub.tolist(), solver=solver) + lp = MILP_LP(A_ineq=A_ineq, b_ineq=b_ineq, A_eq=A_eq, b_eq=b_eq, lb=lb.tolist(), ub=ub.tolist(), solver=solver) _, _, status = lp.solve() if status not in [OPTIMAL, UNBOUNDED]: logging.info('FVA problem not feasible.') return DataFrame( - {"minimum": [nan] * n_orig, "maximum": [nan] * n_orig}, + { + "minimum": [nan] * n_orig, + "maximum": [nan] * n_orig + }, index=reaction_ids, ) @@ -438,8 +449,7 @@ def _bound_scan(x_local): # ------------------------------------------------------------------ _tp1 = _time.perf_counter() # 1a: v=0 feasibility check — free resolutions, no LP needed - v0_feasible = (not np.any(lb > tol_bound) and not np.any(ub < -tol_bound) - and not has_constraints) + v0_feasible = (not np.any(lb > tol_bound) and not np.any(ub < -tol_bound) and not has_constraints) if v0_feasible: zero_lb = np.abs(lb) < tol_bound newly_min = zero_lb & (~res_min) @@ -456,8 +466,7 @@ def _bound_scan(x_local): if precheck: # 1b: min(sum(|x|)) scan — pushes reactions toward zero # Effective for resolving lb=0 / ub=0 bounds in one shot. - scan_lp, n_scan = _build_abssum_lp( - A_eq, b_eq, A_ineq, b_ineq, lb, ub, solver) + scan_lp, n_scan = _build_abssum_lp(A_eq, b_eq, A_ineq, b_ineq, lb, ub, solver) scan_lp.set_lp_method(LP_METHOD_DUAL) n_ext = len(scan_lp.c) @@ -478,9 +487,8 @@ def _bound_scan(x_local): if verbose: n_done_iter = int(res_max.sum() + res_min.sum()) - logging.debug( - f" Phase 1 min|x|: +{resolved_absmin} " - f"({n_done_iter}/{2*n_orig} resolved)") + logging.debug(f" Phase 1 min|x|: +{resolved_absmin} " + f"({n_done_iter}/{2*n_orig} resolved)") # 1c: Iterative push-to-bounds — directed per-reaction objectives # Push unresolved-max reactions toward ub, unresolved-min toward lb. @@ -510,7 +518,7 @@ def _bound_scan(x_local): x_scan = np.array(x_list_scan[:n_scan], dtype=np.float64) before = int(res_max.sum() + res_min.sum()) _bound_scan(x_scan) - + np.maximum(incumbent_max, x_scan, out=incumbent_max) np.minimum(incumbent_min, x_scan, out=incumbent_min) resolved_push_ub = int(res_max.sum() + res_min.sum()) - before @@ -535,7 +543,7 @@ def _bound_scan(x_local): x_scan = np.array(x_list_scan[:n_scan], dtype=np.float64) before = int(res_max.sum() + res_min.sum()) _bound_scan(x_scan) - + np.maximum(incumbent_max, x_scan, out=incumbent_max) np.minimum(incumbent_min, x_scan, out=incumbent_min) resolved_push_lb = int(res_max.sum() + res_min.sum()) - before @@ -543,10 +551,9 @@ def _bound_scan(x_local): if verbose: n_done_iter = int(res_max.sum() + res_min.sum()) - logging.debug( - f" Phase 1 push {push_iter}: " - f"ub +{resolved_push_ub}, lb +{resolved_push_lb} " - f"({n_done_iter}/{2*n_orig} resolved)") + logging.debug(f" Phase 1 push {push_iter}: " + f"ub +{resolved_push_ub}, lb +{resolved_push_lb} " + f"({n_done_iter}/{2*n_orig} resolved)") push_iter += 1 if resolved_this_round < 5: @@ -570,29 +577,22 @@ def _bound_scan(x_local): unresolved = [] for j in range(n_orig): if not res_max[j]: - unresolved.append(2 * j) # even = max + unresolved.append(2 * j) # even = max if not res_min[j]: unresolved.append(2 * j + 1) # odd = min x_par = [nan] * (2 * n_orig) t0 = _time.perf_counter() if solver == GLPK: - with SDPool(threads, initializer=fva_worker_init_glpk, - initargs=(A_ineq, b_ineq, A_eq, b_eq, - lb.tolist(), ub.tolist())) as pool: + with SDPool(threads, initializer=fva_worker_init_glpk, initargs=(A_ineq, b_ineq, A_eq, b_eq, lb.tolist(), ub.tolist())) as pool: chunk_size = max(1, len(unresolved) // threads) - for i, value in pool.imap_unordered( - fva_worker_compute_glpk, unresolved, - chunksize=chunk_size): + for i, value in pool.imap_unordered(fva_worker_compute_glpk, unresolved, chunksize=chunk_size): x_par[i] = value else: with SDPool(threads, initializer=fva_worker_init, - initargs=(A_ineq, b_ineq, A_eq, b_eq, - lb.tolist(), ub.tolist(), solver)) as pool: + initargs=(A_ineq, b_ineq, A_eq, b_eq, lb.tolist(), ub.tolist(), solver)) as pool: chunk_size = max(1, len(unresolved) // threads) - for i, value in pool.imap_unordered( - fva_worker_compute, unresolved, - chunksize=chunk_size): + for i, value in pool.imap_unordered(fva_worker_compute, unresolved, chunksize=chunk_size): x_par[i] = value t_solve += _time.perf_counter() - t0 lps_solved += len(unresolved) @@ -602,9 +602,7 @@ def _bound_scan(x_local): if nan_idx: _BATCH = 50 while nan_idx: - lp_retry = MILP_LP(A_ineq=A_ineq, b_ineq=b_ineq, - A_eq=A_eq, b_eq=b_eq, - lb=lb.tolist(), ub=ub.tolist(), solver=solver) + lp_retry = MILP_LP(A_ineq=A_ineq, b_ineq=b_ineq, A_eq=A_eq, b_eq=b_eq, lb=lb.tolist(), ub=ub.tolist(), solver=solver) prev_retry = 0 for i in nan_idx[:_BATCH]: C = idx2c(i, prev_retry) @@ -638,8 +636,7 @@ def _bound_scan(x_local): def _rebuild_lp(): nonlocal lp, prev_col - lp = MILP_LP(A_ineq=A_ineq, b_ineq=b_ineq, A_eq=A_eq, b_eq=b_eq, - lb=lb.tolist(), ub=ub.tolist(), solver=solver) + lp = MILP_LP(A_ineq=A_ineq, b_ineq=b_ineq, A_eq=A_eq, b_eq=b_eq, lb=lb.tolist(), ub=ub.tolist(), solver=solver) prev_col = -1 _rebuild_lp() @@ -738,7 +735,10 @@ def _rebuild_lp(): incumbent_min[np.abs(incumbent_min) < 1e-11] = 0.0 fva_result = DataFrame( - {"minimum": incumbent_min, "maximum": incumbent_max}, + { + "minimum": incumbent_min, + "maximum": incumbent_max + }, index=reaction_ids, ) @@ -746,12 +746,10 @@ def _rebuild_lp(): cmp_msg = "" if cmp_maps: cmp_msg = f" (compressed {n_original}→{n_orig} rxns)" - logging.debug( - f" speedy_fva done{cmp_msg}: {lps_solved} LPs, " - f"{total_bound_resolved} bound-resolved, " - f"{2*n_orig} total objectives") - logging.debug( - f" timing: solve={t_solve:.2f}s, threads={threads}") + logging.debug(f" speedy_fva done{cmp_msg}: {lps_solved} LPs, " + f"{total_bound_resolved} bound-resolved, " + f"{2*n_orig} total objectives") + logging.debug(f" timing: solve={t_solve:.2f}s, threads={threads}") t_phase['phase2'] = _time.perf_counter() - _tp2 t_phase['total'] = _time.perf_counter() - _tp0 @@ -777,10 +775,10 @@ def _rebuild_lp(): # Fast exact reversibility (sign-only FVA) for pre-compression tightening # --------------------------------------------------------------------------- -_REV_SCAN_TOL = 1e-3 # co-option certifies only on flux comfortably above solver noise +_REV_SCAN_TOL = 1e-3 # co-option certifies only on flux comfortably above solver noise _REV_REBUILD_EVERY = 200 -_ZERO_SNAP = 1e-11 # |flux| below this is solver noise, not a direction; 0 disables snapping -_DEGEN_TOL = 1e-6 # warm-start guard: fresh optimum must not fall below a known-achievable incumbent +_ZERO_SNAP = 1e-11 # |flux| below this is solver noise, not a direction; 0 disables snapping +_DEGEN_TOL = 1e-6 # warm-start guard: fresh optimum must not fall below a known-achievable incumbent def _rev_structural_sweep(model): @@ -797,22 +795,28 @@ def _rev_structural_sweep(model): for e in met_rx.values(): prod = [(i, 'f') for i, c in e if c > 0 and af[i]] + [(i, 'r') for i, c in e if c < 0 and ar[i]] cons = [(i, 'f') for i, c in e if c < 0 and af[i]] + [(i, 'r') for i, c in e if c > 0 and ar[i]] - prx = {i for i, _ in prod}; crx = {i for i, _ in cons} + prx = {i for i, _ in prod} + crx = {i for i, _ in cons} def kill(i, d): nonlocal changed if d == 'f' and af[i]: - af[i] = False; changed = True + af[i] = False + changed = True if d == 'r' and ar[i]: - ar[i] = False; changed = True + ar[i] = False + changed = True + if not prod: - for i, d in cons: kill(i, d) + for i, d in cons: + kill(i, d) elif len(prx) == 1: s = next(iter(prx)) for i, d in cons: if i == s: kill(i, d) if not cons: - for i, d in prod: kill(i, d) + for i, d in prod: + kill(i, d) elif len(crx) == 1: s = next(iter(crx)) for i, d in prod: @@ -860,33 +864,42 @@ def fast_reversibility(model, solver=None, compress=True): lb = np.array([float(r.lower_bound) for r in m.reactions]) ub = np.array([float(r.upper_bound) for r in m.reactions]) S = sparse.csr_matrix(create_stoichiometric_matrix(m)) - A_ineq = sparse.csr_matrix((0, n)); b_ineq = [] - A_eq = S; b_eq = [0.0] * S.shape[0] + A_ineq = sparse.csr_matrix((0, n)) + b_ineq = [] + A_eq = S + b_eq = [0.0] * S.shape[0] def build(): - return MILP_LP(A_ineq=A_ineq, b_ineq=b_ineq, A_eq=A_eq, b_eq=b_eq, - lb=lb.tolist(), ub=ub.tolist(), solver=solver) + return MILP_LP(A_ineq=A_ineq, b_ineq=b_ineq, A_eq=A_eq, b_eq=b_eq, lb=lb.tolist(), ub=ub.tolist(), solver=solver) + lp = build() - incumbent_max = np.full(n, -np.inf); incumbent_min = np.full(n, np.inf) + incumbent_max = np.full(n, -np.inf) + incumbent_min = np.full(n, np.inf) # Bound signs are exact model data, not solver output, so they decide on the exact sign: a # direction counts as blocked only if the bound itself forbids it. A tolerance here would # discard a direction whose achievable flux is merely small (a max of 1e-8 is still forward). - res_max = ub <= 0.0 # fwd already blocked by bounds (sweep/original) + res_max = ub <= 0.0 # fwd already blocked by bounds (sweep/original) res_min = lb >= 0.0 incumbent_max[res_max] = ub[res_max] incumbent_min[res_min] = lb[res_min] fixed = ub == lb - res_max[fixed] = True; res_min[fixed] = True - incumbent_max[fixed] = ub[fixed]; incumbent_min[fixed] = lb[fixed] + res_max[fixed] = True + res_min[fixed] = True + incumbent_max[fixed] = ub[fixed] + incumbent_min[fixed] = lb[fixed] def scan(x): - nm = (~res_max) & (x > _REV_SCAN_TOL); res_max[nm] = True + nm = (~res_max) & (x > _REV_SCAN_TOL) + res_max[nm] = True np.maximum(incumbent_max, x, out=incumbent_max) - nn = (~res_min) & (x < -_REV_SCAN_TOL); res_min[nn] = True + nn = (~res_min) & (x < -_REV_SCAN_TOL) + res_min[nn] = True np.minimum(incumbent_min, x, out=incumbent_min) - n_lp = 0; prev_col = -1; seq = 0 + n_lp = 0 + prev_col = -1 + seq = 0 def solve_dir(j, direction): nonlocal prev_col, seq, n_lp @@ -897,14 +910,17 @@ def solve_dir(j, direction): else: lp.set_objective_idx(C) prev_col = j - r = lp.solve(); n_lp += 1; seq += 1 + r = lp.solve() + n_lp += 1 + seq += 1 return r # Feasibility preflight: one zero-objective solve proves the polytope is non-empty (so a later # infeasible status can only come from the objective change, not the model) and its flux vector # seeds every incumbent, so no subsequent warm-start optimum can silently contradict a flux # already witnessed as achievable. - x_feas, _, status_feas = lp.solve(); n_lp += 1 + x_feas, _, status_feas = lp.solve() + n_lp += 1 if status_feas == INFEASIBLE: raise ValueError('fast_reversibility: the model has no steady-state flux distribution.') if status_feas == OPTIMAL and x_feas: @@ -914,15 +930,20 @@ def unknown(j, direction): """Record 'direction not determined': the incumbent goes to +/-inf so the direction is reported as achievable. Tightening on an unproven bound could delete a feasible pathway; reporting a spurious direction only forgoes tightening.""" - if direction == 1: res_max[j] = True; incumbent_max[j] = np.inf - else: res_min[j] = True; incumbent_min[j] = -np.inf + if direction == 1: + res_max[j] = True + incumbent_max[j] = np.inf + else: + res_min[j] = True + incumbent_min[j] = -np.inf for j in range(n): for direction in (1, -1): if (direction == 1 and res_max[j]) or (direction == -1 and res_min[j]): continue if seq > 0 and seq % _REV_REBUILD_EVERY == 0: - lp = build(); prev_col = -1 + lp = build() + prev_col = -1 x_list, obj_val, status = solve_dir(j, direction) if status != OPTIMAL: # UNBOUNDED is a proven infinite direction, every other nonoptimal status (time @@ -934,14 +955,19 @@ def unknown(j, direction): degen = (direction == 1 and np.isfinite(inc) and val < inc - _DEGEN_TOL * (1 + abs(inc))) or \ (direction == -1 and np.isfinite(inc) and val > inc + _DEGEN_TOL * (1 + abs(inc))) if degen: - lp = build(); prev_col = -1 + lp = build() + prev_col = -1 x_list, obj_val, status = solve_dir(j, direction) if status != OPTIMAL: unknown(j, direction) continue val = -obj_val if direction == 1 else obj_val - if direction == 1: res_max[j] = True; incumbent_max[j] = max(incumbent_max[j], val) - else: res_min[j] = True; incumbent_min[j] = min(incumbent_min[j], val) + if direction == 1: + res_max[j] = True + incumbent_max[j] = max(incumbent_max[j], val) + else: + res_min[j] = True + incumbent_min[j] = min(incumbent_min[j], val) scan(np.array(x_list[:n], dtype=np.float64)) # (4) expand compressed min/max back to original reactions @@ -952,5 +978,4 @@ def unknown(j, direction): df = DataFrame({"minimum": incumbent_min, "maximum": incumbent_max}, index=cmp_rid) if cmp_maps: df = _expand_fva(df, cmp_maps, orig_rid) - return {r: (float(df.at[r, 'maximum']) > 0.0, float(df.at[r, 'minimum']) < 0.0) - for r in orig_rid} + return {r: (float(df.at[r, 'maximum']) > 0.0, float(df.at[r, 'minimum']) < 0.0) for r in orig_rid} diff --git a/straindesign/strainDesignMILP.py b/straindesign/strainDesignMILP.py index 325e77e..f7c603a 100644 --- a/straindesign/strainDesignMILP.py +++ b/straindesign/strainDesignMILP.py @@ -116,7 +116,7 @@ def _trim_z_variables(self): self._z_orig_indices = None # no trimming needed return - self._z_orig_indices = keep_z # trimmed_idx -> orig_idx + self._z_orig_indices = keep_z # trimmed_idx -> orig_idx self._orig_num_z = self.num_z n_cont = len(self.c) - self.num_z diff --git a/straindesign/strainDesignModule.py b/straindesign/strainDesignModule.py index df24d3f..e7cd16a 100644 --- a/straindesign/strainDesignModule.py +++ b/straindesign/strainDesignModule.py @@ -231,8 +231,8 @@ def __init__(self, model, module_type, *args, **kwargs): self[MODEL_ID] = model.id self[MODULE_TYPE] = module_type allowed_keys = { - CONSTRAINTS, INNER_OBJECTIVE, INNER_OPT_SENSE, OUTER_OBJECTIVE, OUTER_OPT_SENSE, - INNER_OPT_TOL, OUTER_OPT_TOL, PROD_ID, 'skip_checks', MIN_GCP, 'reac_ids' + CONSTRAINTS, INNER_OBJECTIVE, INNER_OPT_SENSE, OUTER_OBJECTIVE, OUTER_OPT_SENSE, INNER_OPT_TOL, OUTER_OPT_TOL, PROD_ID, + 'skip_checks', MIN_GCP, 'reac_ids' } # set all keys passed in kwargs as properties of the SD_Module object for key, value in kwargs.items(): @@ -262,8 +262,8 @@ def __init__(self, model, module_type, *args, **kwargs): elif self[INNER_OPT_SENSE] not in [MINIMIZE, MAXIMIZE] or self[OUTER_OPT_SENSE] not in [MINIMIZE, MAXIMIZE]: raise Exception('Inner and outer optimization sense must be "' + MINIMIZE + '" or "' + MAXIMIZE + '" (default).') if ((self[INNER_OBJECTIVE] == None) or (self[OUTER_OBJECTIVE] == None)): - raise Exception('When module type is "' + OPTKNOCK + '", "' + ROBUSTKNOCK + - '" or "' + DOUBLEOPT + '", an inner and outer objective function must be provided.') + raise Exception('When module type is "' + OPTKNOCK + '", "' + ROBUSTKNOCK + '" or "' + DOUBLEOPT + + '", an inner and outer objective function must be provided.') elif (self[MODULE_TYPE] == OPTCOUPLE): if self[INNER_OPT_SENSE] is None: self[INNER_OPT_SENSE] = MAXIMIZE @@ -277,8 +277,7 @@ def __init__(self, model, module_type, *args, **kwargs): raise Exception('When module type is "' + OPTCOUPLE + '", the production reaction id must be provided.') if self[MODULE_TYPE] in [PROTECT, SUPPRESS] and self[OUTER_OBJECTIVE] is not None: if self[INNER_OBJECTIVE] is None: - raise Exception('When outer_objective is set for "' + self[MODULE_TYPE] + - '", inner_objective must also be provided.') + raise Exception('When outer_objective is set for "' + self[MODULE_TYPE] + '", inner_objective must also be provided.') if self[OUTER_OPT_SENSE] is None: self[OUTER_OPT_SENSE] = MAXIMIZE if self[OUTER_OPT_SENSE] not in [MINIMIZE, MAXIMIZE]: diff --git a/straindesign/strainDesignProblem.py b/straindesign/strainDesignProblem.py index d994c5e..6bfe031 100644 --- a/straindesign/strainDesignProblem.py +++ b/straindesign/strainDesignProblem.py @@ -278,9 +278,9 @@ def _module_bound_override(self, sd_module): logging.warning(' Module FVA range for %s is inconsistent (min %g > max %g), no bound ' 'override applied.' % (rid, lim.minimum, lim.maximum)) continue - lo = 0.0 if lim.minimum >= tol else None # convincingly non-negative in-module - hi = 0.0 if lim.maximum <= -tol else None # convincingly non-positive in-module - if lo is not None or hi is not None: # both sides fixed == blocked in-module + lo = 0.0 if lim.minimum >= tol else None # convincingly non-negative in-module + hi = 0.0 if lim.maximum <= -tol else None # convincingly non-positive in-module + if lo is not None or hi is not None: # both sides fixed == blocked in-module override[rid] = (lo, hi) return override @@ -335,23 +335,19 @@ def addModule(self, sd_module): # Inequality: c_v · x_v + tol · c_dual · d <= 0 (actual >= tol * optimal) n_v, n_ref, n_d = len(c_v), len(c_inner), len(c_inner_dual) c_dual_scaled = [inner_opt_tol * c for c in c_inner_dual] - ineq_link = sparse.hstack((sparse.csr_matrix(c_v), sparse.csr_matrix((1, n_ref)), - sparse.csr_matrix(c_dual_scaled))) - eq_link = sparse.hstack((sparse.csr_matrix((1, n_v)), sparse.csr_matrix(c_inner), - sparse.csr_matrix(c_inner_dual))) - A_ineq_p = sparse.vstack((sparse.block_diag((A_ineq_v, A_ineq_inner, A_ineq_dual)), - ineq_link)).tocsr() + ineq_link = sparse.hstack((sparse.csr_matrix(c_v), sparse.csr_matrix((1, n_ref)), sparse.csr_matrix(c_dual_scaled))) + eq_link = sparse.hstack((sparse.csr_matrix((1, n_v)), sparse.csr_matrix(c_inner), sparse.csr_matrix(c_inner_dual))) + A_ineq_p = sparse.vstack((sparse.block_diag((A_ineq_v, A_ineq_inner, A_ineq_dual)), ineq_link)).tocsr() b_ineq_p = b_ineq_v + b_ineq_inner + b_ineq_dual + [0.0] - A_eq_p = sparse.vstack((sparse.block_diag((A_eq_v, A_eq_inner, A_eq_dual)), - eq_link)).tocsr() + A_eq_p = sparse.vstack((sparse.block_diag((A_eq_v, A_eq_inner, A_eq_dual)), eq_link)).tocsr() b_eq_p = b_eq_v + b_eq_inner + b_eq_dual + [0.0] lb_p = lb_v + lb_inner + lb_dual ub_p = ub_v + ub_inner + ub_dual z_map_vars_p = sparse.hstack((z_map_vars_v, z_map_vars_inner, z_map_vars_dual)) - z_map_constr_ineq_p = sparse.hstack((z_map_constr_ineq_v, z_map_constr_ineq_inner, - z_map_constr_ineq_dual, sparse.csc_matrix((self.num_z, 1)))) - z_map_constr_eq_p = sparse.hstack((z_map_constr_eq_v, z_map_constr_eq_inner, - z_map_constr_eq_dual, sparse.csc_matrix((self.num_z, 1)))) + z_map_constr_ineq_p = sparse.hstack( + (z_map_constr_ineq_v, z_map_constr_ineq_inner, z_map_constr_ineq_dual, sparse.csc_matrix((self.num_z, 1)))) + z_map_constr_eq_p = sparse.hstack( + (z_map_constr_eq_v, z_map_constr_eq_inner, z_map_constr_eq_dual, sparse.csc_matrix((self.num_z, 1)))) else: # Exact inner optimality (original code) A_ineq_p = sparse.block_diag((A_ineq_v, A_ineq_dual)).tocsr() @@ -679,12 +675,10 @@ def addModule(self, sd_module): A_eq_base = sparse.block_diag((A_eq_p, A_eq_p, A_eq_dl)) b_eq_ext = b_eq_p + b_eq_p + b_eq_dl # Equality anchor: c_out · x_ref + c_dl · d_out = 0 - eq_link = sparse.hstack((sparse.csr_matrix((1, n_p)), - sparse.csr_matrix(c_out_in_p), sparse.csr_matrix(c_dl))) + eq_link = sparse.hstack((sparse.csr_matrix((1, n_p)), sparse.csr_matrix(c_out_in_p), sparse.csr_matrix(c_dl))) # Relaxed inequality: c_out · x_actual + tol * c_dl · d_out <= 0 c_dl_scaled = [outer_opt_tol * v for v in c_dl] - iq_link = sparse.hstack((sparse.csr_matrix(c_out_in_p), - sparse.csr_matrix((1, n_p)), sparse.csr_matrix(c_dl_scaled))) + iq_link = sparse.hstack((sparse.csr_matrix(c_out_in_p), sparse.csr_matrix((1, n_p)), sparse.csr_matrix(c_dl_scaled))) A_eq_p = sparse.vstack((A_eq_base, eq_link)).tocsr() b_eq_p = b_eq_ext + [0.0] A_ineq_p = sparse.vstack((A_ineq_ext, iq_link)).tocsr() @@ -692,16 +686,16 @@ def addModule(self, sd_module): lb_p = lb_p + lb_p + lb_dl ub_p = ub_p + ub_p + ub_dl z_map_vars_p = sparse.hstack((z_map_vars_p, z_map_vars_p, z_map_vars_dl)) - z_map_constr_ineq_p = sparse.hstack((z_map_constr_ineq_p, z_map_constr_ineq_p, - z_map_constr_ineq_dl, sparse.csc_matrix((self.num_z, 1)))) - z_map_constr_eq_p = sparse.hstack((z_map_constr_eq_p, z_map_constr_eq_p, - z_map_constr_eq_dl, sparse.csc_matrix((self.num_z, 1)))) + z_map_constr_ineq_p = sparse.hstack( + (z_map_constr_ineq_p, z_map_constr_ineq_p, z_map_constr_ineq_dl, sparse.csc_matrix((self.num_z, 1)))) + z_map_constr_eq_p = sparse.hstack( + (z_map_constr_eq_p, z_map_constr_eq_p, z_map_constr_eq_dl, sparse.csc_matrix((self.num_z, 1)))) else: # Exact outer optimality (original code) A_ineq_p = sparse.block_diag((A_ineq_p, A_ineq_dl)).tocsr() b_ineq_p = b_ineq_p + b_ineq_dl - A_eq_p = sparse.vstack((sparse.block_diag((A_eq_p, A_eq_dl)), - sparse.hstack((sparse.csr_matrix(c_out_in_p), sparse.csr_matrix(c_dl))))).tocsr() + A_eq_p = sparse.vstack((sparse.block_diag( + (A_eq_p, A_eq_dl)), sparse.hstack((sparse.csr_matrix(c_out_in_p), sparse.csr_matrix(c_dl))))).tocsr() b_eq_p = b_eq_p + b_eq_dl + [0.0] lb_p = lb_p + lb_dl ub_p = ub_p + ub_dl @@ -817,7 +811,7 @@ def link_z(self): # constraint (gurobi/cplex) or the constant self.M (glpk/user-M). knockable_constr_ineq = np.unique(self.z_map_constr_ineq.nonzero()[1]) - _idxz = set(self.idx_z) # O(1) membership in the scan below + _idxz = set(self.idx_z) # O(1) membership in the scan below cont_vars = [i not in _idxz for i in range(0, numvars)] M_lb = [self.lb[i] for i in np.nonzero(cont_vars)[0]] M_ub = [self.ub[i] for i in np.nonzero(cont_vars)[0]] @@ -849,7 +843,7 @@ def link_z(self): max_Ax[i] = np.inf n_multi += 1 logging.info(' Bounding MILP: %d constraints (%d zero, %d single-var, %d multi-var->indicator/M).' % - (num_Ms, n_zero, n_single, n_multi)) + (num_Ms, n_zero, n_single, n_multi)) # round Ms up to 5 digits Ms = [np.ceil(M * 1e5) / 1e5 if not isinf(M) else self.M for M in max_Ax] @@ -1038,9 +1032,9 @@ def link_z(self): continue col = Aic.getcol(z).tocoo() if any((r not in budget_rows) and v != 0 for r, v in zip(col.row.tolist(), col.data.tolist())): - continue # finite-M z-link -> controls a row + continue # finite-M z-link -> controls a row if Aec is not None and Aec.getcol(z).nnz: - continue # equality z-link + continue # equality z-link self.ub[z] = 0.0 n_free += 1 if n_free: @@ -1126,7 +1120,7 @@ def build_primal_from_cbm(model, V_ineq=None, v_ineq=None, V_eq=None, v_eq=None, lb[i] = max(lb[i], float(lo)) if hi is not None: ub[i] = min(ub[i], float(hi)) - if lb[i] > ub[i]: # numeric guard: never emit an inconsistent block + if lb[i] > ub[i]: # numeric guard: never emit an inconsistent block lb[i], ub[i] = float(lo), float(hi) z_map_vars = sparse.identity(numr, 'd', format="csc") z_map_constr_eq = sparse.csc_matrix((numr, A_eq.shape[0])) diff --git a/straindesign/strainDesignSolutions.py b/straindesign/strainDesignSolutions.py index 2d86e51..2889f27 100644 --- a/straindesign/strainDesignSolutions.py +++ b/straindesign/strainDesignSolutions.py @@ -185,13 +185,13 @@ def _translate_genes_to_reactions(sd_list, model): candidate_reacs.update(r.id for r in model.genes.get_by_id(g).reactions) for r in candidate_reacs: gpr_r = rxn_gpr[r] - if gpr_r.eval(ko_off): # reaction still possible under the interventions - if not gpr_r.eval(all_off): # ... only because of a knock-in + if gpr_r.eval(ko_off): # reaction still possible under the interventions + if not gpr_r.eval(all_off): # ... only because of a knock-in reac_ki.add(r) - else: # reaction dead under the interventions - if gpr_r.eval(noki_off): # ... the knock-out is what killed it + else: # reaction dead under the interventions + if gpr_r.eval(noki_off): # ... the knock-out is what killed it reac_ko.add(r) - else: # ... dead regardless (e.g. an un-made knock-in) + else: # ... dead regardless (e.g. an un-made knock-in) reac_no_ki.add(r) reaction_sd[i].update({k: -1.0 for k in reac_ko}) reaction_sd[i].update({k: 1.0 for k in reac_ki}) @@ -461,8 +461,7 @@ def expand_group(self, grp_idx): # Expand expanded = expand_sd([compressed_sd[grp_idx].copy()], cmp_mapReac) - expanded = filter_sd_maxcost(expanded, meta['max_cost'], - meta['uncmp_ko_cost'], meta['uncmp_ki_cost']) + expanded = filter_sd_maxcost(expanded, meta['max_cost'], meta['uncmp_ko_cost'], meta['uncmp_ki_cost']) # Postprocess regulatory interventions (inline to avoid circular import) reg_cost = meta.get('uncmp_reg_cost', {}) for s in expanded: @@ -475,12 +474,11 @@ def expand_group(self, grp_idx): # GPR translation + costs/bounds if self._model is None: - raise RuntimeError( - 'This SDSolutions was loaded without a model, so compressed ' - 'groups cannot be expanded. Reload with ' - 'SDSolutions.load(file, model=True) to rebuild the embedded ' - 'model, or call attach_model(model) with the original model ' - 'that was passed to compute_strain_designs.') + raise RuntimeError('This SDSolutions was loaded without a model, so compressed ' + 'groups cannot be expanded. Reload with ' + 'SDSolutions.load(file, model=True) to rebuild the embedded ' + 'model, or call attach_model(model) with the original model ' + 'that was passed to compute_strain_designs.') model = self._model if self.is_gene_sd: reaction_sd_exp, gene_sd_exp = self._translate_genes_to_reactions(expanded, model) @@ -490,8 +488,7 @@ def expand_group(self, grp_idx): gene_sd_exp = None cost_sd = expanded - sd_cost_exp, itv_bounds_exp, has_complex = self._compute_costs_and_bounds( - cost_sd, reaction_sd_exp, model, self.sd_setup) + sd_cost_exp, itv_bounds_exp, has_complex = self._compute_costs_and_bounds(cost_sd, reaction_sd_exp, model, self.sd_setup) if has_complex: self.has_complex_regul_itv = True @@ -676,11 +673,11 @@ def load(cls, filename, model=None, cmp_model=None): obj = pickle.load(f) def _resolve(arg, embedded): - if arg is True: # rebuild the embedded snapshot + if arg is True: # rebuild the embedded snapshot return model_from_dict(embedded) if embedded is not None else None - if arg is None or arg is False: # attach nothing + if arg is None or arg is False: # attach nothing return None - return arg # an explicit cobra model + return arg # an explicit cobra model obj._model = _resolve(model, getattr(obj, '_embedded_model_dict', None)) obj._cmp_model = _resolve(cmp_model, getattr(obj, '_embedded_cmp_model_dict', None)) diff --git a/tests/conftest.py b/tests/conftest.py index f4556fc..23d3ff5 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -3,15 +3,15 @@ from cobra import Configuration from straindesign.names import * - # --------------------------------------------------------------------------- # Custom CLI flags for test_performance.py tiered benchmarks # --------------------------------------------------------------------------- + def pytest_addoption(parser): for name, help_text in [ ("--medium", "Run iMLcore genome-scale benchmarks (~4 min total)."), - ("--large", "Run iML1515 large-model benchmarks (several min/solver)."), + ("--large", "Run iML1515 large-model benchmarks (several min/solver)."), ]: try: parser.addoption(name, action="store_true", default=False, help=help_text) @@ -22,14 +22,12 @@ def pytest_addoption(parser): def pytest_configure(config): for marker, desc in [ ("medium", "genome-scale benchmark; enable with --medium"), - ("large", "large-model benchmark; enable with --large"), + ("large", "large-model benchmark; enable with --large"), ]: config.addinivalue_line("markers", f"{marker}: {desc}") # Suppress known third-party warnings - config.addinivalue_line("filterwarnings", - "ignore:FigureCanvasTemplate is non-interactive:UserWarning") - config.addinivalue_line("filterwarnings", - "ignore:builtin type Swig.*has no __module__ attribute:DeprecationWarning") + config.addinivalue_line("filterwarnings", "ignore:FigureCanvasTemplate is non-interactive:UserWarning") + config.addinivalue_line("filterwarnings", "ignore:builtin type Swig.*has no __module__ attribute:DeprecationWarning") def pytest_collection_modifyitems(config, items): @@ -40,6 +38,7 @@ def pytest_collection_modifyitems(config, items): if marker in item.keywords: item.add_marker(skip) + cobra_conf = Configuration() bound_thres = max((abs(cobra_conf.lower_bound), abs(cobra_conf.upper_bound))) diff --git a/tests/test_03_plots.py b/tests/test_03_plots.py index 9c1e1b5..0805974 100644 --- a/tests/test_03_plots.py +++ b/tests/test_03_plots.py @@ -28,8 +28,7 @@ def test_plot_3d_space(curr_solver, model_weak_coupling): def test_plot_2d_point(curr_solver, model_weak_coupling): """0D: fix both axes to single values, should show a dot.""" constr = ['r4 = 0', 'r7 = 0', 'r9 = 0', 'r_P = 2', 'r_BM = 4'] - dp, tri, _ = sd.plot_flux_space(model_weak_coupling, ('r_P', 'r_BM'), - constraints=constr, plt_backend='template') + dp, tri, _ = sd.plot_flux_space(model_weak_coupling, ('r_P', 'r_BM'), constraints=constr, plt_backend='template') assert len(dp) == 1 @@ -37,8 +36,7 @@ def test_plot_2d_point(curr_solver, model_weak_coupling): def test_plot_2d_line(curr_solver, model_weak_coupling): """1D: fix one axis, line along the other.""" constr = ['r4 = 0', 'r7 = 0', 'r9 = 0', 'r_P = 2'] - dp, tri, _ = sd.plot_flux_space(model_weak_coupling, ('r_P', 'r_BM'), - constraints=constr, plt_backend='template') + dp, tri, _ = sd.plot_flux_space(model_weak_coupling, ('r_P', 'r_BM'), constraints=constr, plt_backend='template') assert len(dp) >= 2 @@ -46,6 +44,5 @@ def test_plot_2d_line(curr_solver, model_weak_coupling): def test_plot_3d_point(curr_solver, model_weak_coupling): """0D in 3D: all axes fixed, should show a dot.""" constr = ['r4 = 0', 'r7 = 0', 'r9 = 0', 'r_BM = 2', 'r_P = 2', 'r_Q = 2'] - dp, tri, _ = sd.plot_flux_space(model_weak_coupling, ('r_P', 'r_BM', 'r_Q'), - constraints=constr, plt_backend='template') + dp, tri, _ = sd.plot_flux_space(model_weak_coupling, ('r_P', 'r_BM', 'r_Q'), constraints=constr, plt_backend='template') assert len(dp) == 1 diff --git a/tests/test_04_preprocessing.py b/tests/test_04_preprocessing.py index a0f60d9..005e7a5 100644 --- a/tests/test_04_preprocessing.py +++ b/tests/test_04_preprocessing.py @@ -20,9 +20,9 @@ from cobra.core.gene import GPR from sympy import simplify_logic - # ── GPR extension + compression ────────────────────────────────────── + @pytest.mark.timeout(15) def test_gpr_extension_compression1(model_gpr): gkocost = { @@ -61,33 +61,33 @@ def test_gpr_extension_compression2(model_gpr): # ── GPR propagation helper unit tests ──────────────────────────────── + class TestGprAstToExpr: # None is special-cased (identity check); the rest share one assertion shape. def test_none_returns_none(self): assert _gpr_ast_to_expr(None) is None - @pytest.mark.parametrize("node,expected", [ - (ast.Name(id='g1'), 'g1'), - (ast.BoolOp(op=ast.And(), values=[ast.Name(id='g1'), ast.Name(id='g2')]), - ('and', ('g1', 'g2'))), - (ast.BoolOp(op=ast.Or(), values=[ast.Name(id='g1'), ast.Name(id='g2')]), - ('or', ('g1', 'g2'))), - # (g1 and g2) or g3 - (ast.BoolOp(op=ast.Or(), values=[ - ast.BoolOp(op=ast.And(), values=[ast.Name(id='g1'), ast.Name(id='g2')]), - ast.Name(id='g3')]), - ('or', (('and', ('g1', 'g2')), 'g3'))), - # g1 and (g2 and g3) -> flattened - (ast.BoolOp(op=ast.And(), values=[ - ast.Name(id='g1'), - ast.BoolOp(op=ast.And(), values=[ast.Name(id='g2'), ast.Name(id='g3')])]), - ('and', ('g1', 'g2', 'g3'))), - ], ids=['single_gene', 'and', 'or', 'nested', 'nested_same_op_flattened']) + @pytest.mark.parametrize( + "node,expected", + [ + (ast.Name(id='g1'), 'g1'), + (ast.BoolOp(op=ast.And(), values=[ast.Name(id='g1'), ast.Name(id='g2')]), ('and', ('g1', 'g2'))), + (ast.BoolOp(op=ast.Or(), values=[ast.Name(id='g1'), ast.Name(id='g2')]), ('or', ('g1', 'g2'))), + # (g1 and g2) or g3 + (ast.BoolOp(op=ast.Or(), values=[ast.BoolOp(op=ast.And(), values=[ast.Name(id='g1'), ast.Name(id='g2')]), + ast.Name(id='g3')]), ('or', (('and', ('g1', 'g2')), 'g3'))), + # g1 and (g2 and g3) -> flattened + (ast.BoolOp(op=ast.And(), values=[ast.Name(id='g1'), + ast.BoolOp(op=ast.And(), values=[ast.Name(id='g2'), ast.Name(id='g3')])]), + ('and', ('g1', 'g2', 'g3'))), + ], + ids=['single_gene', 'and', 'or', 'nested', 'nested_same_op_flattened']) def test_gpr_ast_to_expr(self, node, expected): assert _gpr_ast_to_expr(node) == expected class TestExprToGprString: + @pytest.mark.parametrize("expr,expected", [ (None, ''), ('g1', 'g1'), @@ -95,7 +95,8 @@ class TestExprToGprString: (('or', ['g1', 'g2']), 'g1 or g2'), (('or', [('and', ['g1', 'g2']), 'g3']), '(g1 and g2) or g3'), (('and', [('or', ['g1', 'g2']), 'g3']), '(g1 or g2) and g3'), - ], ids=['none', 'single_gene', 'and', 'or', 'nested_and_in_or', 'nested_or_in_and']) + ], + ids=['none', 'single_gene', 'and', 'or', 'nested_and_in_or', 'nested_or_in_and']) def test_expr_to_gpr_string(self, expr, expected): assert _expr_to_gpr_string(expr) == expected @@ -112,38 +113,42 @@ def test_roundtrips_through_cobra(self): class TestCombineGprAnd: - @pytest.mark.parametrize("nodes,expected", [ - ([None, None], ''), - ([ast.BoolOp(op=ast.And(), values=[ast.Name(id='g1'), ast.Name(id='g2')])], - 'g1 and g2'), - # Empty GPR (None) should be skipped in AND combination. - ([ast.BoolOp(op=ast.And(), values=[ast.Name(id='g1'), ast.Name(id='g2')]), None, None], - 'g1 and g2'), - ([ast.BoolOp(op=ast.And(), values=[ast.Name(id='g1'), ast.Name(id='g2')]), - ast.Name(id='g3')], - 'g1 and g2 and g3'), - # AND of overlapping expressions should simplify: (g1 and g2) AND (g1 and g3). - ([ast.BoolOp(op=ast.And(), values=[ast.Name(id='g1'), ast.Name(id='g2')]), - ast.BoolOp(op=ast.And(), values=[ast.Name(id='g1'), ast.Name(id='g3')])], - 'g1 and g2 and g3'), - ([], ''), - ], ids=['all_empty', 'single_non_empty', 'skip_empty', 'two_non_empty', - 'simplification', 'empty_list']) + + @pytest.mark.parametrize( + "nodes,expected", + [ + ([None, None], ''), + ([ast.BoolOp(op=ast.And(), values=[ast.Name(id='g1'), ast.Name(id='g2')])], 'g1 and g2'), + # Empty GPR (None) should be skipped in AND combination. + ([ast.BoolOp(op=ast.And(), values=[ast.Name(id='g1'), ast.Name(id='g2')]), None, None], 'g1 and g2'), + ([ast.BoolOp(op=ast.And(), values=[ast.Name(id='g1'), ast.Name(id='g2')]), + ast.Name(id='g3')], 'g1 and g2 and g3'), + # AND of overlapping expressions should simplify: (g1 and g2) AND (g1 and g3). + ([ + ast.BoolOp(op=ast.And(), values=[ast.Name(id='g1'), ast.Name(id='g2')]), + ast.BoolOp(op=ast.And(), values=[ast.Name(id='g1'), ast.Name(id='g3')]) + ], 'g1 and g2 and g3'), + ([], ''), + ], + ids=['all_empty', 'single_non_empty', 'skip_empty', 'two_non_empty', 'simplification', 'empty_list']) def test_combine_gprs_and(self, nodes, expected): assert _combine_gprs(nodes, 'and') == expected class TestCombineGprOr: - @pytest.mark.parametrize("nodes,expected", [ - # If any reaction has empty GPR (always active), result is empty. - ([ast.Name(id='g1'), None], ''), - ([None, None], ''), - ([ast.Name(id='g1'), ast.Name(id='g2')], 'g1 or g2'), - # OR with duplicate terms should deduplicate: g1 OR g1 -> g1. - ([ast.Name(id='g1'), ast.Name(id='g1')], 'g1'), - ([], ''), - ], ids=['any_empty_returns_empty', 'all_empty', 'two_non_empty', - 'deduplication', 'empty_list']) + + @pytest.mark.parametrize( + "nodes,expected", + [ + # If any reaction has empty GPR (always active), result is empty. + ([ast.Name(id='g1'), None], ''), + ([None, None], ''), + ([ast.Name(id='g1'), ast.Name(id='g2')], 'g1 or g2'), + # OR with duplicate terms should deduplicate: g1 OR g1 -> g1. + ([ast.Name(id='g1'), ast.Name(id='g1')], 'g1'), + ([], ''), + ], + ids=['any_empty_returns_empty', 'all_empty', 'two_non_empty', 'deduplication', 'empty_list']) def test_combine_gprs_or(self, nodes, expected): assert _combine_gprs(nodes, 'or') == expected @@ -158,12 +163,14 @@ def test_no_absorption(self): # ── GPR propagation integration tests (model_gpr.xml) ──────────────── + @pytest.fixture def gpr_model(): return read_sbml_model('tests/model_gpr.xml') class TestModelGprCompression: + def test_coupled_compression_propagates_gpr(self, gpr_model): """Coupled compression should AND-combine GPR rules, skipping empty ones.""" remove_blocked_reactions(gpr_model) @@ -234,8 +241,7 @@ def test_coupled_group_r4_r5_r6_rdex(self, gpr_model): from cobra.core.gene import GPR result_sympy = GPR.from_string(target_rxn.gene_reaction_rule).as_symbolic() - expected = GPR.from_string( - '(g1 and g4 and g7 and g8) or (g1 and g4 and g5 and g8 and g9)').as_symbolic() + expected = GPR.from_string('(g1 and g4 and g7 and g8) or (g1 and g4 and g5 and g8 and g9)').as_symbolic() assert simplify_logic(result_sympy ^ expected) == False, \ f"GPR mismatch. Got: {result_sympy}, expected: {expected}" @@ -267,4 +273,3 @@ def test_coupled_group_r3_rpex(self, gpr_model): assert simplify_logic(result_sympy ^ expected) == False, \ f"GPR mismatch. Got: {result_sympy}, expected: {expected}" - diff --git a/tests/test_05_straindesign.py b/tests/test_05_straindesign.py index e3381a5..f24ea1f 100644 --- a/tests/test_05_straindesign.py +++ b/tests/test_05_straindesign.py @@ -259,32 +259,47 @@ def test_doubleopt(curr_solver, model_doubleopt): need to achieve 60% of their optima. """ # Knockable: organism-internal reactions + shared sinks (enable cross-feeding KOs) - ko_reacs = [r.id for r in model_doubleopt.reactions if r.id.startswith('A_R') or r.id.startswith('B_R') - or r.id in ['A_10', 'B_10', 'shared_R_D', 'shared_R_C']] + ko_reacs = [ + r.id + for r in model_doubleopt.reactions + if r.id.startswith('A_R') or r.id.startswith('B_R') or r.id in ['A_10', 'B_10', 'shared_R_D', 'shared_R_C'] + ] kocost = {r: 1 for r in ko_reacs} # --- Exact DOUBLEOPT: tight coupling via shared sink KOs --- - modules_exact = [sd.SDModule(model_doubleopt, DOUBLEOPT, - inner_objective='A_BM', - outer_objective='B_BM', - constraints=['A_BM >= 0.1', 'B_BM >= 0.1'])] + modules_exact = [ + sd.SDModule(model_doubleopt, DOUBLEOPT, inner_objective='A_BM', outer_objective='B_BM', constraints=['A_BM >= 0.1', 'B_BM >= 0.1']) + ] sol_exact = sd.compute_strain_designs(model_doubleopt.copy(), - sd_modules=modules_exact, max_cost=6, max_solutions=inf, - solution_approach=POPULATE, ko_cost=kocost, solver=curr_solver, compress=False) + sd_modules=modules_exact, + max_cost=6, + max_solutions=inf, + solution_approach=POPULATE, + ko_cost=kocost, + solver=curr_solver, + compress=False) assert len(sol_exact.reaction_sd) > 0, \ "Exact DOUBLEOPT should find solutions when shared sinks are knockable" min_cost_exact = min(sum(abs(v) for v in s.values()) for s in sol_exact.reaction_sd) # --- Relaxed DOUBLEOPT: 60% optimality tolerance finds cheaper solutions --- - modules_relaxed = [sd.SDModule(model_doubleopt, DOUBLEOPT, - inner_objective='A_BM', - outer_objective='B_BM', - inner_opt_tol=0.6, - outer_opt_tol=0.6, - constraints=['A_BM >= 0.1', 'B_BM >= 0.1'])] + modules_relaxed = [ + sd.SDModule(model_doubleopt, + DOUBLEOPT, + inner_objective='A_BM', + outer_objective='B_BM', + inner_opt_tol=0.6, + outer_opt_tol=0.6, + constraints=['A_BM >= 0.1', 'B_BM >= 0.1']) + ] sol_relaxed = sd.compute_strain_designs(model_doubleopt.copy(), - sd_modules=modules_relaxed, max_cost=3, max_solutions=inf, - solution_approach=POPULATE, ko_cost=kocost, solver=curr_solver, compress=False) + sd_modules=modules_relaxed, + max_cost=3, + max_solutions=inf, + solution_approach=POPULATE, + ko_cost=kocost, + solver=curr_solver, + compress=False) assert len(sol_relaxed.reaction_sd) > 0, \ "Relaxed DOUBLEOPT should find solutions with opt_tol=0.6" min_cost_relaxed = min(sum(abs(v) for v in s.values()) for s in sol_relaxed.reaction_sd) @@ -348,9 +363,14 @@ def test_dump_preprocessed(model_small_example, tmp_path): # Step 1: Dump preprocessed data (should return early without solving) sol_dump = sd.compute_strain_designs(model_small_example, - sd_modules=modules, max_cost=inf, max_solutions=inf, - solution_approach='any', ki_cost=kicost, solver=solver, - compress=True, dump_preprocessed=dump_path) + sd_modules=modules, + max_cost=inf, + max_solutions=inf, + solution_approach='any', + ki_cost=kicost, + solver=solver, + compress=True, + dump_preprocessed=dump_path) import os assert os.path.exists(dump_path), "Dump file should exist" @@ -387,9 +407,13 @@ def test_lazy_expansion(model_small_example): csd.LAZY_EXPANSION_THRESHOLD = 1 # Force lazy mode try: sol = sd.compute_strain_designs(model_small_example, - sd_modules=modules, max_cost=inf, max_solutions=inf, - solution_approach='any', ki_cost=kicost, solver=solver, - compress=True) + sd_modules=modules, + max_cost=inf, + max_solutions=inf, + solution_approach='any', + ki_cost=kicost, + solver=solver, + compress=True) if sol.is_lazy: assert sol.get_num_materialized() < sol.get_num_sols(), \ "Lazy mode: materialized < estimated total" diff --git a/tests/test_07_compression.py b/tests/test_07_compression.py index a60961b..d5d8aae 100644 --- a/tests/test_07_compression.py +++ b/tests/test_07_compression.py @@ -150,11 +150,9 @@ def test_fba_optimum_recovered_through_map(): model = load_model("e_coli_core") cmp_maps = nt.compress_model(model) cmp_id, factor = _trace_lump(cmp_maps, biomass) - assert cmp_id in [r.id for r in model.reactions], ( - f"compression map names {cmp_id}, which is not in the compressed model") + assert cmp_id in [r.id for r in model.reactions], (f"compression map names {cmp_id}, which is not in the compressed model") val = sd.fba(model, obj={cmp_id: 1}, obj_sense='maximize').objective_value - assert abs(factor * val - ref) < 1e-6, ( - f"recovered optimum {factor * val} != uncompressed {ref}") + assert abs(factor * val - ref) < 1e-6, (f"recovered optimum {factor * val} != uncompressed {ref}") def test_cobra_optimize_after_compression(): @@ -188,8 +186,7 @@ def test_cobra_optimize_after_compression(): sol = model_cmp.optimize() assert sol.status == 'optimal', f"Expected optimal solution, got {sol.status}" val_expanded = sol.objective_value * biomass_coeff - assert abs(val_orig - val_expanded) < 1e-6, ( - f"Expanded objective mismatch: original={val_orig}, expanded={val_expanded}") + assert abs(val_orig - val_expanded) < 1e-6, (f"Expanded objective mismatch: original={val_orig}, expanded={val_expanded}") # ============================================================================= @@ -252,13 +249,12 @@ def test_mcs_e_coli_core(): solver = SCIP if SCIP in strong_solvers else next(iter(strong_solvers)) model = load_model('e_coli_core') modules = [sd.SDModule(model, SUPPRESS, constraints='BIOMASS_Ecoli_core_w_GAM >= 0.001')] - sols = sd.compute_strain_designs(model, - sd_modules=modules, - solution_approach=POPULATE, - max_cost=3, - gene_kos=True, - solver=solver, -) - assert len(sols.reaction_sd) == 455, ( - f"Expected 455 MCS for e_coli_core, got {len(sols.reaction_sd)}") - + sols = sd.compute_strain_designs( + model, + sd_modules=modules, + solution_approach=POPULATE, + max_cost=3, + gene_kos=True, + solver=solver, + ) + assert len(sols.reaction_sd) == 455, (f"Expected 455 MCS for e_coli_core, got {len(sols.reaction_sd)}") diff --git a/tests/test_08_gene_ko.py b/tests/test_08_gene_ko.py index b28d1f6..b1925b1 100644 --- a/tests/test_08_gene_ko.py +++ b/tests/test_08_gene_ko.py @@ -17,7 +17,9 @@ def wt_growth(ecoli_core): # ── gene_kos_to_constraints unit tests ─────────────────────────────── + class TestGeneKosToConstraints: + def test_single_gene_by_id(self, ecoli_core): """b0727 (sucB) is in AND rule for AKGDH → single KO kills AKGDH.""" c = gene_kos_to_constraints(ecoli_core, ['b0727']) @@ -62,7 +64,9 @@ def test_empty_input(self, ecoli_core): # ── resolve_gene_constraints unit tests ────────────────────────────── + class TestResolveGeneConstraints: + def test_string_format(self, ecoli_core): """Gene KO as string 'b0727 = 0' should produce reaction constraints.""" c = resolve_gene_constraints(ecoli_core, 'b0727 = 0') @@ -116,7 +120,9 @@ def test_sd_grammar_mixed(self, ecoli_core): # ── FBA integration tests ──────────────────────────────────────────── + class TestFbaWithGeneKO: + def test_nonessential_single_ko(self, ecoli_core, wt_growth): """b0727 (sucB) KO reduces growth slightly but is not lethal.""" sol = sd.fba(ecoli_core, constraints='b0727 = 0') @@ -160,7 +166,9 @@ def test_sd_grammar_minus_one_in_fba(self, ecoli_core, wt_growth): # ── FVA integration tests ──────────────────────────────────────────── + class TestFvaWithGeneKO: + def test_knocked_reaction_fixed_at_zero(self, ecoli_core): """FVA after b0727 KO: AKGDH must have min=max=0.""" fva_r = sd.fva(ecoli_core, constraints='b0727 = 0') diff --git a/tests/test_09_performance.py b/tests/test_09_performance.py index 715052d..9e6252c 100644 --- a/tests/test_09_performance.py +++ b/tests/test_09_performance.py @@ -88,7 +88,8 @@ def _git_sha() -> str: try: return subprocess.check_output( ["git", "rev-parse", "--short", "HEAD"], - cwd=TESTS_DIR, stderr=subprocess.DEVNULL, + cwd=TESTS_DIR, + stderr=subprocess.DEVNULL, ).decode().strip() except Exception: return "unknown" @@ -107,8 +108,7 @@ def _solver_ver(solver: str) -> str: return "n/a" -def record(name: str, solver: str, model_id: str, - elapsed: float, n_sol: int, status: str) -> None: +def record(name: str, solver: str, model_id: str, elapsed: float, n_sol: int, status: str) -> None: entry = { "test": name, "solver": solver, @@ -128,11 +128,11 @@ def record(name: str, solver: str, model_id: str, STRONG_SOLVERS = [s for s in [CPLEX, GUROBI] if s in sd.avail_solvers] - # --------------------------------------------------------------------------- # Model fixtures (session-scoped → loaded once per run) # --------------------------------------------------------------------------- + @pytest.fixture(scope="session") def model_core(): return load_model("e_coli_core") @@ -152,6 +152,7 @@ def model_imlcore(): # Quick suite — e_coli_core (MCS correctness + speed, ~8 s / solver) # =========================================================================== + @pytest.mark.parametrize("solver", STRONG_SOLVERS) @pytest.mark.timeout(180) def test_ecoli_core_mcs_455(solver, model_core): @@ -164,18 +165,15 @@ def test_ecoli_core_mcs_455(solver, model_core): t0 = time.perf_counter() sol = sd.compute_strain_designs( m, - sd_modules=[sd.SDModule(m, SUPPRESS, - constraints="BIOMASS_Ecoli_core_w_GAM >= 0.001")], + sd_modules=[sd.SDModule(m, SUPPRESS, constraints="BIOMASS_Ecoli_core_w_GAM >= 0.001")], solution_approach=POPULATE, max_cost=3, gene_kos=True, solver=solver, ) elapsed = time.perf_counter() - t0 - record("mcs_455", solver, "e_coli_core", elapsed, - len(sol.reaction_sd), sol.status) - assert len(sol.reaction_sd) == 455, ( - f"[{solver}] Expected 455 MCS, got {len(sol.reaction_sd)}") + record("mcs_455", solver, "e_coli_core", elapsed, len(sol.reaction_sd), sol.status) + assert len(sol.reaction_sd) == 455, (f"[{solver}] Expected 455 MCS, got {len(sol.reaction_sd)}") # =========================================================================== @@ -189,8 +187,17 @@ def test_ecoli_core_mcs_455(solver, model_core): # Shared KO / KI costs for model_weak_coupling (from test_05_straindesign.py) _WEAK_KO = { - 'r1': 1, 'r2': 1, 'r4': 1.1, 'r5': 0.75, 'r7': 0.8, 'r8': 1, - 'r9': 1, 'r_S': 1.0, 'r_P': 1, 'r_BM': 1, 'r_Q': 1.5, + 'r1': 1, + 'r2': 1, + 'r4': 1.1, + 'r5': 0.75, + 'r7': 0.8, + 'r8': 1, + 'r9': 1, + 'r_S': 1.0, + 'r_P': 1, + 'r_BM': 1, + 'r_Q': 1.5, } _WEAK_KI = {'r3': 0.6, 'r6': 1.0} _WEAK_REG = {'r6 >= 4.5': 1.2} @@ -208,9 +215,7 @@ def test_weak_mcs_wgcp(solver, model_weak): """ m = model_weak.copy() modules = [ - sd.SDModule(m, SUPPRESS, - inner_objective="r_BM", - constraints=["r_P - 0.4 r_S <= 0", "r_S >= 0.1"]), + sd.SDModule(m, SUPPRESS, inner_objective="r_BM", constraints=["r_P - 0.4 r_S <= 0", "r_S >= 0.1"]), sd.SDModule(m, PROTECT, constraints=["r_BM >= 0.2"]), ] t0 = time.perf_counter() @@ -222,14 +227,12 @@ def test_weak_mcs_wgcp(solver, model_weak): max_solutions=inf, ko_cost=_WEAK_KO, ki_cost=_WEAK_KI, - reg_cost=dict(_WEAK_REG), # fresh copy: extend_model_regulatory mutates its arg + reg_cost=dict(_WEAK_REG), # fresh copy: extend_model_regulatory mutates its arg solver=solver, ) elapsed = time.perf_counter() - t0 - record("mcs_wgcp", solver, "weak_coupling", elapsed, - len(sol.reaction_sd), sol.status) - assert len(sol.reaction_sd) == 3, ( - f"[{solver}] Expected 3 wGCP MCS solutions, got {len(sol.reaction_sd)}") + record("mcs_wgcp", solver, "weak_coupling", elapsed, len(sol.reaction_sd), sol.status) + assert len(sol.reaction_sd) == 3, (f"[{solver}] Expected 3 wGCP MCS solutions, got {len(sol.reaction_sd)}") @pytest.mark.parametrize("solver", STRONG_SOLVERS) @@ -242,10 +245,7 @@ def test_weak_optknock(solver, model_weak): """ m = model_weak.copy() modules = [ - sd.SDModule(m, OPTKNOCK, - outer_objective="r_P", - inner_objective="r_BM", - constraints="r_BM >= 1"), + sd.SDModule(m, OPTKNOCK, outer_objective="r_P", inner_objective="r_BM", constraints="r_BM >= 1"), ] t0 = time.perf_counter() sol = sd.compute_strain_designs( @@ -256,14 +256,12 @@ def test_weak_optknock(solver, model_weak): max_solutions=3, ko_cost=_WEAK_KO, ki_cost=_WEAK_KI, - reg_cost=dict(_WEAK_REG), # fresh copy: extend_model_regulatory mutates its arg + reg_cost=dict(_WEAK_REG), # fresh copy: extend_model_regulatory mutates its arg solver=solver, ) elapsed = time.perf_counter() - t0 - record("optknock", solver, "weak_coupling", elapsed, - len(sol.reaction_sd), sol.status) - assert len(sol.reaction_sd) == 3, ( - f"[{solver}] Expected 3 OptKnock solutions, got {len(sol.reaction_sd)}") + record("optknock", solver, "weak_coupling", elapsed, len(sol.reaction_sd), sol.status) + assert len(sol.reaction_sd) == 3, (f"[{solver}] Expected 3 OptKnock solutions, got {len(sol.reaction_sd)}") @pytest.mark.parametrize("solver", STRONG_SOLVERS) @@ -276,10 +274,9 @@ def test_weak_robustknock(solver, model_weak): """ m = model_weak.copy() modules = [ - sd.SDModule(m, ROBUSTKNOCK, - outer_objective="r_P", - inner_objective="r_BM", - constraints=[[{"r_BM": 1.0}, ">=", 1.0]]), + sd.SDModule(m, ROBUSTKNOCK, outer_objective="r_P", inner_objective="r_BM", constraints=[[{ + "r_BM": 1.0 + }, ">=", 1.0]]), ] t0 = time.perf_counter() sol = sd.compute_strain_designs( @@ -290,14 +287,12 @@ def test_weak_robustknock(solver, model_weak): max_solutions=3, ko_cost=_WEAK_KO, ki_cost=_WEAK_KI, - reg_cost=dict(_WEAK_REG), # fresh copy: extend_model_regulatory mutates its arg + reg_cost=dict(_WEAK_REG), # fresh copy: extend_model_regulatory mutates its arg solver=solver, ) elapsed = time.perf_counter() - t0 - record("robustknock", solver, "weak_coupling", elapsed, - len(sol.reaction_sd), sol.status) - assert len(sol.reaction_sd) >= 2, ( - f"[{solver}] Expected ≥2 RobustKnock solutions, got {len(sol.reaction_sd)}") + record("robustknock", solver, "weak_coupling", elapsed, len(sol.reaction_sd), sol.status) + assert len(sol.reaction_sd) >= 2, (f"[{solver}] Expected ≥2 RobustKnock solutions, got {len(sol.reaction_sd)}") @pytest.mark.parametrize("solver", STRONG_SOLVERS) @@ -310,10 +305,7 @@ def test_weak_optcouple(solver, model_weak): """ m = model_weak.copy() modules = [ - sd.SDModule(m, OPTCOUPLE, - prod_id="r_P", - inner_objective="r_BM", - min_gcp=1.0), + sd.SDModule(m, OPTCOUPLE, prod_id="r_P", inner_objective="r_BM", min_gcp=1.0), ] t0 = time.perf_counter() sol = sd.compute_strain_designs( @@ -324,20 +316,19 @@ def test_weak_optcouple(solver, model_weak): max_solutions=3, ko_cost=_WEAK_KO, ki_cost=_WEAK_KI, - reg_cost=dict(_WEAK_REG), # fresh copy: extend_model_regulatory mutates its arg + reg_cost=dict(_WEAK_REG), # fresh copy: extend_model_regulatory mutates its arg solver=solver, ) elapsed = time.perf_counter() - t0 - record("optcouple", solver, "weak_coupling", elapsed, - len(sol.reaction_sd), sol.status) - assert len(sol.reaction_sd) == 2, ( - f"[{solver}] Expected 2 OptCouple solutions, got {len(sol.reaction_sd)}") + record("optcouple", solver, "weak_coupling", elapsed, len(sol.reaction_sd), sol.status) + assert len(sol.reaction_sd) == 2, (f"[{solver}] Expected 2 OptCouple solutions, got {len(sol.reaction_sd)}") # =========================================================================== # Standard suite (--medium): iMLcore gene-level MCS (~47 s / solver each) # =========================================================================== + @pytest.mark.medium @pytest.mark.parametrize("solver", STRONG_SOLVERS) @pytest.mark.timeout(90) @@ -350,11 +341,8 @@ def test_imlcore_mcs_ethanol(solver, model_imlcore): """ m = model_imlcore.copy() modules = [ - sd.SDModule(m, SUPPRESS, - constraints=["EX_etoh_e <= 1.0", - "BIOMASS_Ec_iML1515_core_75p37M >= 0.14"]), - sd.SDModule(m, PROTECT, - constraints=["BIOMASS_Ec_iML1515_core_75p37M >= 0.15"]), + sd.SDModule(m, SUPPRESS, constraints=["EX_etoh_e <= 1.0", "BIOMASS_Ec_iML1515_core_75p37M >= 0.14"]), + sd.SDModule(m, PROTECT, constraints=["BIOMASS_Ec_iML1515_core_75p37M >= 0.15"]), ] t0 = time.perf_counter() sol = sd.compute_strain_designs( @@ -366,10 +354,8 @@ def test_imlcore_mcs_ethanol(solver, model_imlcore): solver=solver, ) elapsed = time.perf_counter() - t0 - record("imlcore_ethanol", solver, "iMLcore", elapsed, - len(sol.reaction_sd), sol.status) - assert len(sol.reaction_sd) > 0, ( - f"[{solver}] Expected ≥1 MCS for iMLcore ethanol scenario, got 0") + record("imlcore_ethanol", solver, "iMLcore", elapsed, len(sol.reaction_sd), sol.status) + assert len(sol.reaction_sd) > 0, (f"[{solver}] Expected ≥1 MCS for iMLcore ethanol scenario, got 0") @pytest.mark.medium @@ -383,8 +369,7 @@ def test_imlcore_mcs_growth(solver, model_imlcore): """ m = model_imlcore.copy() modules = [ - sd.SDModule(m, SUPPRESS, - constraints="BIOMASS_Ec_iML1515_core_75p37M >= 0.001"), + sd.SDModule(m, SUPPRESS, constraints="BIOMASS_Ec_iML1515_core_75p37M >= 0.001"), ] t0 = time.perf_counter() sol = sd.compute_strain_designs( @@ -396,16 +381,15 @@ def test_imlcore_mcs_growth(solver, model_imlcore): solver=solver, ) elapsed = time.perf_counter() - t0 - record("imlcore_growth", solver, "iMLcore", elapsed, - len(sol.reaction_sd), sol.status) - assert len(sol.reaction_sd) > 0, ( - f"[{solver}] Expected ≥1 MCS for iMLcore growth scenario, got 0") + record("imlcore_growth", solver, "iMLcore", elapsed, len(sol.reaction_sd), sol.status) + assert len(sol.reaction_sd) > 0, (f"[{solver}] Expected ≥1 MCS for iMLcore growth scenario, got 0") # =========================================================================== # Large suite (--large): iML1515 — known answer 393 # =========================================================================== + @pytest.mark.large @pytest.mark.parametrize("solver", [GUROBI]) # CPLEX is very slow on this one; skip to save time def test_iml1515_mcs_393(solver): @@ -421,24 +405,22 @@ def test_iml1515_mcs_393(solver): t0 = time.perf_counter() sol = sd.compute_strain_designs( m, - sd_modules=[sd.SDModule(m, SUPPRESS, - constraints="BIOMASS_Ec_iML1515_core_75p37M >= 0.001")], + sd_modules=[sd.SDModule(m, SUPPRESS, constraints="BIOMASS_Ec_iML1515_core_75p37M >= 0.001")], solution_approach=POPULATE, max_cost=3, gene_kos=True, solver=solver, ) elapsed = time.perf_counter() - t0 - record("iml1515_393", solver, "iML1515", elapsed, - len(sol.reaction_sd), sol.status) - assert len(sol.reaction_sd) == 393, ( - f"[{solver}] Expected 393 MCS for iML1515, got {len(sol.reaction_sd)}") + record("iml1515_393", solver, "iML1515", elapsed, len(sol.reaction_sd), sol.status) + assert len(sol.reaction_sd) == 393, (f"[{solver}] Expected 393 MCS for iML1515, got {len(sol.reaction_sd)}") # =========================================================================== # Session teardown: write JSON + print comparison table # =========================================================================== + @pytest.fixture(scope="session", autouse=True) def _write_results(): yield @@ -449,7 +431,9 @@ def _write_results(): "git_sha": _git_sha(), "platform": platform.platform(), "python": platform.python_version(), - "solver_versions": {s: _solver_ver(s) for s in STRONG_SOLVERS}, + "solver_versions": { + s: _solver_ver(s) for s in STRONG_SOLVERS + }, "results": _RESULTS, } out_file = RESULTS_DIR / f"{_SESSION_TS}.json" @@ -460,14 +444,12 @@ def _write_results(): tests_ran = sorted({(r["test"], r["model"]) for r in _RESULTS}) col = 16 - header = (f"{'Test':<22} {'Model':<16}" - + "".join(f" {s:>{col}}" for s in solvers_ran)) + header = (f"{'Test':<22} {'Model':<16}" + "".join(f" {s:>{col}}" for s in solvers_ran)) print("=== Solver Performance Summary ===") print(header) print("-" * len(header)) for (test, model) in tests_ran: - row = {r["solver"]: r for r in _RESULTS - if r["test"] == test and r["model"] == model} + row = {r["solver"]: r for r in _RESULTS if r["test"] == test and r["model"] == model} line = f"{test:<22} {model:<16}" for s in solvers_ran: if s in row: diff --git a/tests/test_10_gene_design_validity.py b/tests/test_10_gene_design_validity.py index 8f144ca..25fd85e 100644 --- a/tests/test_10_gene_design_validity.py +++ b/tests/test_10_gene_design_validity.py @@ -9,8 +9,7 @@ from math import inf from cobra.io import read_sbml_model import straindesign as sd -from straindesign.names import (MODULES, MAX_COST, MAX_SOLUTIONS, SOLUTION_APPROACH, - KOCOST, GKOCOST, SOLVER) +from straindesign.names import (MODULES, MAX_COST, MAX_SOLUTIONS, SOLUTION_APPROACH, KOCOST, GKOCOST, SOLVER) from straindesign import SUPPRESS, PROTECT TOL = 1e-6 @@ -37,11 +36,25 @@ def _apply_gene_design(model, design): def _gpr_mcs_setup(model, gko_cost, solver, approach): - modules = [sd.SDModule(model, SUPPRESS, constraints=["1.0 rd_ex >= 1.0 "]), - sd.SDModule(model, PROTECT, constraints=[[{'r_bm': 1.0}, '>=', 1.0]])] - return {MODULES: modules, MAX_COST: 3, MAX_SOLUTIONS: inf, SOLUTION_APPROACH: approach, - KOCOST: {'rs_up': 1.0, 'rd_ex': 1.0, 'rp_ex': 1.1}, - GKOCOST: gko_cost, SOLVER: solver} + modules = [ + sd.SDModule(model, SUPPRESS, constraints=["1.0 rd_ex >= 1.0 "]), + sd.SDModule(model, PROTECT, constraints=[[{ + 'r_bm': 1.0 + }, '>=', 1.0]]) + ] + return { + MODULES: modules, + MAX_COST: 3, + MAX_SOLUTIONS: inf, + SOLUTION_APPROACH: approach, + KOCOST: { + 'rs_up': 1.0, + 'rd_ex': 1.0, + 'rp_ex': 1.1 + }, + GKOCOST: gko_cost, + SOLVER: solver + } @pytest.fixture diff --git a/tests/test_11_gene_regulatory_compression.py b/tests/test_11_gene_regulatory_compression.py index 54fc60e..8049c45 100644 --- a/tests/test_11_gene_regulatory_compression.py +++ b/tests/test_11_gene_regulatory_compression.py @@ -34,12 +34,23 @@ def _coupled_chain(): m = Model("chain") mets = {x: Metabolite(x) for x in "ABCD"} m.add_metabolites(list(mets.values())) - vin = Reaction("vin"); vin.add_metabolites({mets["A"]: 1}); vin.bounds = (0, 10) - r1 = Reaction("r1"); r1.add_metabolites({mets["A"]: -1, mets["B"]: 1}); r1.bounds = (0, 1000) - r2 = Reaction("r2"); r2.add_metabolites({mets["B"]: -1, mets["C"]: 1}); r2.bounds = (0, 1000) - r3 = Reaction("r3"); r3.add_metabolites({mets["C"]: -1, mets["D"]: 1}); r3.bounds = (0, 1000) - vout = Reaction("vout"); vout.add_metabolites({mets["D"]: -1}); vout.bounds = (0, 1000) - m.add_reactions([vin, r1, r2, r3, vout]); m.objective = "vout" + vin = Reaction("vin") + vin.add_metabolites({mets["A"]: 1}) + vin.bounds = (0, 10) + r1 = Reaction("r1") + r1.add_metabolites({mets["A"]: -1, mets["B"]: 1}) + r1.bounds = (0, 1000) + r2 = Reaction("r2") + r2.add_metabolites({mets["B"]: -1, mets["C"]: 1}) + r2.bounds = (0, 1000) + r3 = Reaction("r3") + r3.add_metabolites({mets["C"]: -1, mets["D"]: 1}) + r3.bounds = (0, 1000) + vout = Reaction("vout") + vout.add_metabolites({mets["D"]: -1}) + vout.bounds = (0, 1000) + m.add_reactions([vin, r1, r2, r3, vout]) + m.objective = "vout" return m @@ -48,7 +59,8 @@ def test_coupled_exemption_keeps_reaction_and_merges_rest(): and flux space is preserved.""" base = _fba(_coupled_chain(), {"vout": 1}) - m_full = _coupled_chain(); compress_model(m_full) + m_full = _coupled_chain() + compress_model(m_full) assert len(m_full.reactions) == 1, "unprotected coupled chain should fully collapse" m_prot = _coupled_chain() @@ -62,20 +74,31 @@ def test_coupled_exemption_keeps_reaction_and_merges_rest(): def _gene_model(): """g1 controls r1 (A->2B) and r2 (B->C), which are coupled 1:2.""" - m = Model("g"); A, B, C = Metabolite("A"), Metabolite("B"), Metabolite("C") + m = Model("g") + A, B, C = Metabolite("A"), Metabolite("B"), Metabolite("C") m.add_metabolites([A, B, C]) - vS = Reaction("vS"); vS.add_metabolites({A: 1}); vS.bounds = (0, 10) - r1 = Reaction("r1"); r1.add_metabolites({A: -1, B: 2}); r1.bounds = (0, 1000); r1.gene_reaction_rule = "g1" - r2 = Reaction("r2"); r2.add_metabolites({B: -1, C: 1}); r2.bounds = (0, 1000); r2.gene_reaction_rule = "g1" - vC = Reaction("vC"); vC.add_metabolites({C: -1}); vC.bounds = (0, 1000) - m.add_reactions([vS, r1, r2, vC]); m.objective = "vC" + vS = Reaction("vS") + vS.add_metabolites({A: 1}) + vS.bounds = (0, 10) + r1 = Reaction("r1") + r1.add_metabolites({A: -1, B: 2}) + r1.bounds = (0, 1000) + r1.gene_reaction_rule = "g1" + r2 = Reaction("r2") + r2.add_metabolites({B: -1, C: 1}) + r2.bounds = (0, 1000) + r2.gene_reaction_rule = "g1" + vC = Reaction("vC") + vC.add_metabolites({C: -1}) + vC.bounds = (0, 1000) + m.add_reactions([vS, r1, r2, vC]) + m.objective = "vC" return m def _max_C_under_g1_le_1(protect): m = _gene_model() - cmap = compress_model(m, no_coupled_compress_reacs=({"r1", "r2"} if protect else set()), - propagate_gpr=True) + cmap = compress_model(m, no_coupled_compress_reacs=({"r1", "r2"} if protect else set()), propagate_gpr=True) extend_model_gpr(m, use_names=False) m.reactions.g1.upper_bound = 1.0 return _fba(m, _remap_obj({"vC": 1}, cmap)) @@ -84,7 +107,9 @@ def _max_C_under_g1_le_1(protect): def test_gene_regulatory_multiplicity_preserved_under_compression(): """With g1 controlling coupled r1,r2, the regulatory bound g1<=1 must give the same max product compressed (with protection) as uncompressed; without protection it is wrong (3x).""" - mu = _gene_model(); extend_model_gpr(mu, use_names=False); mu.reactions.g1.upper_bound = 1.0 + mu = _gene_model() + extend_model_gpr(mu, use_names=False) + mu.reactions.g1.upper_bound = 1.0 ref = _fba(mu, {"vC": 1}) assert abs(ref - 2.0 / 3.0) < 1e-4, "reference should be 2/3" diff --git a/tests/test_12_save_load_embed_model.py b/tests/test_12_save_load_embed_model.py index 176ace6..e3ef907 100644 --- a/tests/test_12_save_load_embed_model.py +++ b/tests/test_12_save_load_embed_model.py @@ -12,8 +12,7 @@ # imported into the package), so reach the module via sys.modules. csd = sys.modules["straindesign.compute_strain_designs"] from straindesign import SUPPRESS, PROTECT -from straindesign.names import (MODULES, MAX_COST, MAX_SOLUTIONS, - SOLUTION_APPROACH, KOCOST, GKOCOST, SOLVER, SEED) +from straindesign.names import (MODULES, MAX_COST, MAX_SOLUTIONS, SOLUTION_APPROACH, KOCOST, GKOCOST, SOLVER, SEED) GPR = os.path.join(os.path.dirname(__file__), "model_gpr.xml") TOL = 1e-6 @@ -36,12 +35,15 @@ def _compute(model, threshold=None): if threshold is not None: old, csd.LAZY_EXPANSION_THRESHOLD = csd.LAZY_EXPANSION_THRESHOLD, threshold try: - return sd.compute_strain_designs( - model, - sd_modules=[sd.SDModule(model, SUPPRESS, - constraints="Biomass_Ecoli_core >= 0.1")], - gene_kos=True, max_cost=3, max_solutions=8, - solution_approach="any", solver=_solver(), seed=1, compress=True) + return sd.compute_strain_designs(model, + sd_modules=[sd.SDModule(model, SUPPRESS, constraints="Biomass_Ecoli_core >= 0.1")], + gene_kos=True, + max_cost=3, + max_solutions=8, + solution_approach="any", + solver=_solver(), + seed=1, + compress=True) finally: if threshold is not None: csd.LAZY_EXPANSION_THRESHOLD = old @@ -53,7 +55,7 @@ def test_embed_and_restore_roundtrip(model, tmp_path): ref_rsd = sols.get_reaction_sd() ref_gsd = sols.get_gene_sd() f = str(tmp_path / "sd.pkl") - sols.save(f) # embed_model=True by default + sols.save(f) # embed_model=True by default # default load: model NOT rebuilt, but solutions are there loaded = sd.SDSolutions.load(f) @@ -67,14 +69,13 @@ def test_embed_and_restore_roundtrip(model, tmp_path): assert len(restored._model.reactions) == len(model.reactions) assert len(restored._model.genes) == len(model.genes) r = "AKGDH" - assert (restored._model.reactions.get_by_id(r).gene_reaction_rule - == model.reactions.get_by_id(r).gene_reaction_rule) + assert (restored._model.reactions.get_by_id(r).gene_reaction_rule == model.reactions.get_by_id(r).gene_reaction_rule) def test_explicit_model_takes_precedence(model, tmp_path): sols = _compute(model) f = str(tmp_path / "sd.pkl") - sols.save(f, embed_model=False) # leaner file, no snapshot + sols.save(f, embed_model=False) # leaner file, no snapshot loaded = sd.SDSolutions.load(f) assert loaded._embedded_model_dict is None # nothing to restore from, stays model-less @@ -88,19 +89,19 @@ def test_live_model_and_solver_never_pickled(model, tmp_path): f = str(tmp_path / "sd.pkl") sols.save(f) with open(f, "rb") as fh: - raw = pickle.load(fh) # must not raise (no live solver in pickle) + raw = pickle.load(fh) # must not raise (no live solver in pickle) assert raw._model is None assert raw._embedded_model_dict is not None # ── lazy round-trip: save must not force-expand; restore then expand ───── def test_lazy_save_no_expand_then_restore_and_expand(model, tmp_path): - sols = _compute(model, threshold=1) # force lazy expansion + sols = _compute(model, threshold=1) # force lazy expansion assert sols.is_lazy materialized = sols.get_num_materialized() f = str(tmp_path / "lazy.pkl") - sols.save(f) # must NOT expand_all (no hang) - assert sols.is_lazy # still lazy after save + sols.save(f) # must NOT expand_all (no hang) + assert sols.is_lazy # still lazy after save # default load: lazy, materialized reps available, expand errors clearly loaded = sd.SDSolutions.load(f) @@ -123,6 +124,7 @@ def test_lazy_save_no_expand_then_restore_and_expand(model, tmp_path): # on compressed representatives AND after expansion, then the whole # computation is REPRODUCED from the artifact alone (embedded model + sd_setup). + def _max_flux(model, rid): with model: model.objective = rid @@ -157,18 +159,32 @@ def _validate_designs(reaction_designs): def _keyset(reaction_designs): """Canonical set of designs = set of frozensets of knocked-out reaction ids.""" - return {frozenset(k for k, v in d.items() if v in (-1, -1.0, False)) - for d in reaction_designs} + return {frozenset(k for k, v in d.items() if v in (-1, -1.0, False)) for d in reaction_designs} def _combined_setup(model, solver, approach): - modules = [sd.SDModule(model, SUPPRESS, constraints=["1.0 rd_ex >= 1.0 "]), - sd.SDModule(model, PROTECT, constraints=[[{"r_bm": 1.0}, ">=", 1.0]])] - return {MODULES: modules, MAX_COST: 3, MAX_SOLUTIONS: inf, - SOLUTION_APPROACH: approach, - KOCOST: {"rs_up": 1.0, "rd_ex": 1.0, "rp_ex": 1.1}, # reaction KOs - GKOCOST: {g.id: 1.0 for g in model.genes}, # gene KOs - SOLVER: solver, SEED: 7} + modules = [ + sd.SDModule(model, SUPPRESS, constraints=["1.0 rd_ex >= 1.0 "]), + sd.SDModule(model, PROTECT, constraints=[[{ + "r_bm": 1.0 + }, ">=", 1.0]]) + ] + return { + MODULES: modules, + MAX_COST: 3, + MAX_SOLUTIONS: inf, + SOLUTION_APPROACH: approach, + KOCOST: { + "rs_up": 1.0, + "rd_ex": 1.0, + "rp_ex": 1.1 + }, # reaction KOs + GKOCOST: { + g.id: 1.0 for g in model.genes + }, # gene KOs + SOLVER: solver, + SEED: 7 + } @pytest.mark.parametrize("approach", ["any", "best", "populate"]) @@ -192,12 +208,12 @@ def test_full_reproducibility_roundtrip(tmp_path, approach, force_lazy): ref_gene_sd = sol.get_gene_sd() ref_keys = _keyset(sol.get_reaction_sd()) assert ref_gene_sd - _validate_designs(sol.get_reaction_sd()) # sanity on the fresh result + _validate_designs(sol.get_reaction_sd()) # sanity on the fresh result # save a self-contained artifact (embed_model=True by default) f = str(tmp_path / "repro.pkl") sol.save(f) - assert sol.is_lazy == force_lazy # save never expanded a lazy result + assert sol.is_lazy == force_lazy # save never expanded a lazy result # (A) reload WITHOUT a model: designs preserved; validate the COMPRESSED # representatives directly via SUPPRESS/PROTECT FBA @@ -210,11 +226,11 @@ def test_full_reproducibility_roundtrip(tmp_path, approach, force_lazy): # expanded design set via FBA restored = sd.SDSolutions.load(f, model=True, cmp_model=True) assert restored._model is not None - restored.expand_all() # no-op if already non-lazy + restored.expand_all() # no-op if already non-lazy assert not restored.is_lazy _validate_designs(restored.get_reaction_sd()) - expanded_keys = _keyset(restored.get_reaction_sd()) # full design set - assert expanded_keys >= ref_keys # reps ⊆ full set + expanded_keys = _keyset(restored.get_reaction_sd()) # full design set + assert expanded_keys >= ref_keys # reps ⊆ full set # (C) REPRODUCE the computation from the artifact ALONE: the embedded # (uncompressed) model + the stored sd_setup must regenerate the exact @@ -229,11 +245,27 @@ def test_full_reproducibility_roundtrip(tmp_path, approach, force_lazy): def _combined_gpr(model, solver, approach="any"): - modules = [sd.SDModule(model, SUPPRESS, constraints=["1.0 rd_ex >= 1.0 "]), - sd.SDModule(model, PROTECT, constraints=[[{"r_bm": 1.0}, ">=", 1.0]])] - return {MODULES: modules, MAX_COST: 3, SOLUTION_APPROACH: approach, - KOCOST: {"rs_up": 1.0, "rd_ex": 1.0, "rp_ex": 1.1}, - GKOCOST: {g.id: 1.0 for g in model.genes}, SOLVER: solver, SEED: 7} + modules = [ + sd.SDModule(model, SUPPRESS, constraints=["1.0 rd_ex >= 1.0 "]), + sd.SDModule(model, PROTECT, constraints=[[{ + "r_bm": 1.0 + }, ">=", 1.0]]) + ] + return { + MODULES: modules, + MAX_COST: 3, + SOLUTION_APPROACH: approach, + KOCOST: { + "rs_up": 1.0, + "rd_ex": 1.0, + "rp_ex": 1.1 + }, + GKOCOST: { + g.id: 1.0 for g in model.genes + }, + SOLVER: solver, + SEED: 7 + } def test_networktools_model_dict_rational_exact(): @@ -243,18 +275,21 @@ def test_networktools_model_dict_rational_exact(): from straindesign.networktools import model_to_dict, model_from_dict import json, math m = Model("rat") - A = Metabolite("A_c", compartment="c"); B = Metabolite("B_c", compartment="c") - r = Reaction("r1"); m.add_reactions([r]) + A = Metabolite("A_c", compartment="c") + B = Metabolite("B_c", compartment="c") + r = Reaction("r1") + m.add_reactions([r]) r.add_metabolites({A: Fraction(1, 3), B: Fraction(-7, 3)}) - r.lower_bound = Fraction(1, 3); r.upper_bound = float("inf") + r.lower_bound = Fraction(1, 3) + r.upper_bound = float("inf") d = model_to_dict(m) - rt = model_from_dict(json.loads(json.dumps(d))) # force through JSON primitives + rt = model_from_dict(json.loads(json.dumps(d))) # force through JSON primitives lb = rt.reactions.r1.lower_bound coeff = rt.reactions.r1.metabolites[rt.metabolites.A_c] - assert isinstance(lb, Fraction) and lb == Fraction(1, 3) # exact, still rational + assert isinstance(lb, Fraction) and lb == Fraction(1, 3) # exact, still rational assert isinstance(coeff, Fraction) and coeff == Fraction(1, 3) - assert float(Fraction(1, 3)) != Fraction(1, 3) # sanity: float would differ - assert math.isinf(rt.reactions.r1.upper_bound) # inf round-trips + assert float(Fraction(1, 3)) != Fraction(1, 3) # sanity: float would differ + assert math.isinf(rt.reactions.r1.upper_bound) # inf round-trips def test_networktools_model_dict_float_matches_cobra(model): @@ -263,7 +298,7 @@ def test_networktools_model_dict_float_matches_cobra(model): from straindesign.networktools import model_to_dict, model_from_dict import json d = model_to_dict(model) - json.dumps(d) # must be JSON-clean + json.dumps(d) # must be JSON-clean rt = model_from_dict(d) assert len(rt.reactions) == len(model.reactions) assert len(rt.genes) == len(model.genes) @@ -278,12 +313,12 @@ def test_compressed_model_embedded_by_default_restore_optin(tmp_path): assert sol._cmp_model is not None and len(sol._cmp_model.reactions) < len(m.reactions) cm0 = sol._cmp_model f = str(tmp_path / "sd.pkl") - sol.save(f) # embeds BOTH by default + sol.save(f) # embeds BOTH by default # default load restores neither model o0 = sd.SDSolutions.load(f) assert o0.get_model() is None and o0.get_compressed_model() is None - assert o0._embedded_cmp_model_dict is not None # ...but it IS embedded + assert o0._embedded_cmp_model_dict is not None # ...but it IS embedded # model=True, cmp_model=True rebuild both; compressed model is exact + usable o = sd.SDSolutions.load(f, model=True, cmp_model=True) @@ -296,18 +331,19 @@ def test_compressed_model_embedded_by_default_restore_optin(tmp_path): assert r0.lower_bound == r1.lower_bound and r0.upper_bound == r1.upper_bound for a in ("lower_bound", "upper_bound"): if isinstance(getattr(r0, a), Fraction): - assert isinstance(getattr(r1, a), Fraction) # rational not float-ified + assert isinstance(getattr(r1, a), Fraction) # rational not float-ified for mt, c0 in r0.metabolites.items(): c1 = r1.metabolites[cm.metabolites.get_by_id(mt.id)] assert c0 == c1 and (not isinstance(c0, Fraction) or isinstance(c1, Fraction)) assert any(isinstance(b, Fraction) for r in cm.reactions for b in (r.lower_bound, r.upper_bound)) cm.objective = cm.reactions[0].id - assert cm.optimize().status == "optimal" # a fully working model + assert cm.optimize().status == "optimal" # a fully working model # explicit overrides take precedence assert sd.SDSolutions.load(f, cmp_model=cm0).get_compressed_model() is cm0 # embed_model=False -> no compressed snapshot - f2 = str(tmp_path / "lean.pkl"); sol.save(f2, embed_model=False) + f2 = str(tmp_path / "lean.pkl") + sol.save(f2, embed_model=False) assert sd.SDSolutions.load(f2, model=True, cmp_model=True).get_compressed_model() is None @@ -317,14 +353,15 @@ def test_compressed_solutions_analyzable_in_restored_cmp_model(tmp_path): solver = _solver() m = read_sbml_model(GPR) sol = sd.compute_strain_designs(m, sd_setup=_combined_gpr(m, solver)) - f = str(tmp_path / "sd.pkl"); sol.save(f) + f = str(tmp_path / "sd.pkl") + sol.save(f) o = sd.SDSolutions.load(f, model=True, cmp_model=True) cm = o.get_compressed_model() cm_rxn_ids = {r.id for r in cm.reactions} non_empty = [cs for cs in o.compressed_sd if cs] assert non_empty, "expected at least one non-empty compressed solution" for cs in non_empty: - assert set(cs).issubset(cm_rxn_ids) # analysable in the compressed model + assert set(cs).issubset(cm_rxn_ids) # analysable in the compressed model # apply one compressed KO set and confirm FBA still runs in the small model cs = non_empty[0] with cm: @@ -339,8 +376,12 @@ def test_networktools_preserves_objective_direction(): """networktools serializers keep the objective SENSE, which cobra drops.""" from cobra import Model, Metabolite, Reaction from straindesign.networktools import model_to_dict, model_from_dict - m = Model("obj"); A = Metabolite("A_c", compartment="c") - r = Reaction("r1", lower_bound=0, upper_bound=10); m.add_reactions([r]) - r.add_metabolites({A: 1.0}); m.objective = "r1"; m.objective.direction = "min" + m = Model("obj") + A = Metabolite("A_c", compartment="c") + r = Reaction("r1", lower_bound=0, upper_bound=10) + m.add_reactions([r]) + r.add_metabolites({A: 1.0}) + m.objective = "r1" + m.objective.direction = "min" rt = model_from_dict(model_to_dict(m)) assert rt.objective.direction == "min"