diff --git a/dingo/LooplessFluxSampler.py b/dingo/LooplessFluxSampler.py new file mode 100644 index 0000000..1d91fff --- /dev/null +++ b/dingo/LooplessFluxSampler.py @@ -0,0 +1,405 @@ +# dingo : a python library for metabolic networks sampling and analysis +# dingo is part of GeomScale project + +# Copyright (c) 2024 + +# Licensed under GNU LGPL.3, see LICENCE file + +# ============================================================================= +# LooplessFluxSampler: Non-convex sampling for thermodynamically feasible fluxes +# ============================================================================= +# +# Flux sampling in metabolic networks can produce thermodynamically infeasible +# solutions: internal cycles that continuously consume/generate energy without +# net effect. These "Type III pathways" violate the second law of thermodynamics. +# +# This module implements: +# 1. Detection of internal cycles via nullspace analysis +# 2. Rejection-based loopless sampling (sample + filter) +# 3. A penalty-based approach that biases sampling away from loopy solutions +# +# References: +# [1] Chalkis et al. - dingo: Python package for metabolic flux sampling +# [2] De Martino - Scales and multimodal flux distributions via thermodynamics +# [3] Saa et al. - LooplessFluxSampler +# ============================================================================= + +import numpy as np +from scipy.linalg import null_space +from dingo.MetabolicNetwork import MetabolicNetwork +from dingo.PolytopeSampler import PolytopeSampler + + +class LooplessFluxSampler: + """ + A sampler that generates thermodynamically feasible (loopless) steady states + from a metabolic network by detecting and rejecting internal flux loops. + + Internal loops are Type III pathways: cycles of internal reactions that carry + flux without being driven by any exchange reaction. They violate the second + law of thermodynamics and are biologically meaningless. + + Strategy (rejection-based): + 1. Identify internal reactions (non-exchange, non-biomass). + 2. Compute the nullspace of the internal stoichiometric sub-matrix. + Vectors in this nullspace represent potential internal cycles. + 3. Sample from the full flux polytope using dingo's standard samplers. + 4. For each sample, project the internal flux vector onto the cycle + nullspace. If the projection norm exceeds a tolerance, the sample + contains an active loop and is rejected. + + Parameters + ---------- + metabolic_network : MetabolicNetwork + A dingo MetabolicNetwork object. + loop_tolerance : float, optional + Maximum allowed norm of the internal-cycle projection (default 1e-5). + """ + + def __init__(self, metabolic_network, loop_tolerance=1e-5): + if not isinstance(metabolic_network, MetabolicNetwork): + raise TypeError("Expected a MetabolicNetwork object.") + + self._network = metabolic_network + self._loop_tolerance = loop_tolerance + + # Identify internal vs exchange reactions + self._internal_indices = [] + self._exchange_indices = [] + self._classify_reactions() + + # Compute the internal cycle nullspace + self._cycle_nullspace = None + self._compute_cycle_nullspace() + + # ------------------------------------------------------------------ + # Reaction classification + # ------------------------------------------------------------------ + def _classify_reactions(self): + """ + Classify reactions as internal or exchange. + + Exchange reactions are identified as those involving only one + metabolite (single non-zero entry in the stoichiometric column), + or reactions listed in the model's exchanges property. + """ + S = self._network.S + reactions = self._network.reactions + exchanges = set(self._network.exchanges) if self._network.exchanges else set() + n_reactions = S.shape[1] + + for j in range(n_reactions): + col = S[:, j] + n_nonzero = np.count_nonzero(col) + + # A reaction is considered exchange if: + # 1) It has a single nonzero stoichiometric coefficient, or + # 2) It is listed in the model's exchange reactions + if n_nonzero <= 1 or reactions[j] in exchanges: + self._exchange_indices.append(j) + else: + self._internal_indices.append(j) + + # ------------------------------------------------------------------ + # Internal cycle detection via nullspace + # ------------------------------------------------------------------ + def _compute_cycle_nullspace(self): + """ + Compute the nullspace of the internal stoichiometric sub-matrix. + + For internal reactions only, we extract S_int (the sub-matrix of S + with columns corresponding to internal reactions). The nullspace of + S_int gives directions in which internal reactions can carry flux + while still satisfying S_int * v_int = 0 – these are internal cycles. + """ + S = self._network.S + int_idx = self._internal_indices + + if len(int_idx) == 0: + self._cycle_nullspace = np.empty((0, 0)) + return + + # Extract the internal sub-matrix + S_int = S[:, int_idx] + + # Compute the nullspace of S_int + # Each column of N_int is a potential internal cycle direction + N_int = null_space(S_int) + + self._cycle_nullspace = N_int + + def _has_active_loop(self, flux_vector): + """ + Check whether a flux vector contains an active internal loop. + + Projects the internal flux sub-vector onto the cycle nullspace. + If the projection has significant magnitude, the sample carries + an internal cycle. + + Parameters + ---------- + flux_vector : ndarray + A complete flux vector (dimension = number of reactions). + + Returns + ------- + bool + True if an active internal loop is detected. + float + The norm of the loop projection (for diagnostics). + """ + if self._cycle_nullspace.size == 0: + return False, 0.0 + + # Extract internal fluxes + v_int = flux_vector[self._internal_indices] + + # Project onto the cycle nullspace: proj = N * N^T * v_int + N = self._cycle_nullspace + projection = N @ (N.T @ v_int) + loop_norm = np.linalg.norm(projection) + + return loop_norm > self._loop_tolerance, loop_norm + + # ------------------------------------------------------------------ + # Loopless sampling (rejection-based) + # ------------------------------------------------------------------ + def sample_loopless( + self, + n_samples=500, + max_attempts_factor=10, + method="billiard_walk", + burn_in=100, + thinning=2, + opt_percentage=None, + ): + """ + Generate loopless steady-state flux samples using rejection sampling. + + Samples are drawn from the flux polytope and those containing + internal loops are discarded. + + Parameters + ---------- + n_samples : int + Desired number of loopless samples. + max_attempts_factor : int + Maximum total samples to draw = n_samples * max_attempts_factor. + method : str + MCMC sampling method (default: 'billiard_walk'). + burn_in : int + Number of burn-in samples (default: 100). + thinning : int + Thinning factor for the MCMC chain (default: 2). + opt_percentage : int or None + If not None, set the model's opt_percentage before sampling. + + Returns + ------- + loopless_samples : ndarray + Steady states without internal loops (shape: n_reactions x n_accepted). + rejection_stats : dict + Statistics about the rejection process. + """ + if opt_percentage is not None: + self._network.set_opt_percentage(opt_percentage) + + sampler = PolytopeSampler(self._network) + + max_total = n_samples * max_attempts_factor + + accepted = [] + total_sampled = 0 + total_rejected = 0 + loop_norms = [] + + # Draw all samples at once (more efficient than batching) + n_draw = min(max_total, max(n_samples * 5, 1000)) + steady_states = sampler.generate_steady_states_no_multiphase( + method=method, n=n_draw, burn_in=burn_in, thinning=thinning + ) + + total_sampled = n_draw + + # Filter samples for loops + for j in range(steady_states.shape[1]): + flux = steady_states[:, j] + has_loop, norm = self._has_active_loop(flux) + loop_norms.append(norm) + + if not has_loop: + accepted.append(flux) + if len(accepted) >= n_samples: + break + else: + total_rejected += 1 + + # Stack accepted samples + if len(accepted) > 0: + loopless_samples = np.column_stack(accepted[:n_samples]) + else: + loopless_samples = np.empty((self._network.num_of_reactions(), 0)) + + rejection_stats = { + "total_sampled": total_sampled, + "total_accepted": len(accepted), + "total_rejected": total_rejected, + "acceptance_rate": len(accepted) / max(total_sampled, 1), + "mean_loop_norm": float(np.mean(loop_norms)) if loop_norms else 0.0, + "max_loop_norm": float(np.max(loop_norms)) if loop_norms else 0.0, + "loop_tolerance": self._loop_tolerance, + } + + return loopless_samples, rejection_stats + + # ------------------------------------------------------------------ + # Penalty-weighted sampling + # ------------------------------------------------------------------ + def sample_penalized( + self, + n_samples=500, + penalty_weight=10.0, + method="billiard_walk", + burn_in=100, + thinning=2, + opt_percentage=None, + ): + """ + Generate flux samples where loopy samples are down-weighted + using importance sampling with a loop-penalty function. + + Instead of hard rejection, each sample receives a weight: + w_i = exp(-penalty_weight * ||projection_i||^2) + + This produces a weighted sample set that favours loopless solutions. + + Parameters + ---------- + n_samples : int + Number of samples to draw (all are kept, but weighted). + penalty_weight : float + Strength of the penalty (higher = more bias against loops). + method : str + MCMC sampling method. + burn_in : int + Burn-in samples. + thinning : int + Thinning factor. + opt_percentage : int or None + If not None, set the model's opt_percentage before sampling. + + Returns + ------- + samples : ndarray + All steady state samples (shape: n_reactions x n_samples). + weights : ndarray + Importance weights for each sample (shape: n_samples,). + penalty_stats : dict + Statistics about the penalty distribution. + """ + if opt_percentage is not None: + self._network.set_opt_percentage(opt_percentage) + + sampler = PolytopeSampler(self._network) + steady_states = sampler.generate_steady_states_no_multiphase( + method=method, n=n_samples, burn_in=burn_in, thinning=thinning + ) + + weights = np.zeros(steady_states.shape[1]) + loop_norms = np.zeros(steady_states.shape[1]) + + for j in range(steady_states.shape[1]): + flux = steady_states[:, j] + _, norm = self._has_active_loop(flux) + loop_norms[j] = norm + + # Compute weights in log-space to avoid underflow + # w_i = exp(-penalty_weight * norm_i^2) + log_weights = -penalty_weight * loop_norms ** 2 + # Shift so the maximum log-weight is 0 (prevents all-zero underflow) + log_weights -= np.max(log_weights) + weights = np.exp(log_weights) + + # Normalise weights + weight_sum = np.sum(weights) + if weight_sum > 0: + weights /= weight_sum + else: + # Uniform fallback if all weights are zero + weights[:] = 1.0 / len(weights) + + penalty_stats = { + "mean_loop_norm": float(np.mean(loop_norms)), + "max_loop_norm": float(np.max(loop_norms)), + "n_effectively_loopless": int(np.sum(loop_norms < self._loop_tolerance)), + "effective_sample_size": float(1.0 / np.sum(weights ** 2)), + "penalty_weight": penalty_weight, + } + + return steady_states, weights, penalty_stats + + # ------------------------------------------------------------------ + # Diagnostics + # ------------------------------------------------------------------ + def analyze_loops(self, steady_states): + """ + Analyze the loop content of a set of steady-state samples. + + Parameters + ---------- + steady_states : ndarray + Steady states (shape: n_reactions x n_samples). + + Returns + ------- + dict + Loop analysis results including per-sample norms, fraction + with loops, and the most common loop directions. + """ + n_samples = steady_states.shape[1] + loop_norms = np.zeros(n_samples) + has_loops = np.zeros(n_samples, dtype=bool) + + for j in range(n_samples): + has_loop, norm = self._has_active_loop(steady_states[:, j]) + loop_norms[j] = norm + has_loops[j] = has_loop + + return { + "loop_norms": loop_norms, + "fraction_with_loops": float(np.mean(has_loops)), + "mean_loop_norm": float(np.mean(loop_norms)), + "median_loop_norm": float(np.median(loop_norms)), + "max_loop_norm": float(np.max(loop_norms)), + "n_internal_reactions": len(self._internal_indices), + "n_exchange_reactions": len(self._exchange_indices), + "cycle_nullspace_dim": self._cycle_nullspace.shape[1] + if self._cycle_nullspace.size > 0 else 0, + } + + # ------------------------------------------------------------------ + # Properties + # ------------------------------------------------------------------ + @property + def internal_indices(self): + return self._internal_indices + + @property + def exchange_indices(self): + return self._exchange_indices + + @property + def cycle_nullspace(self): + return self._cycle_nullspace + + @property + def loop_tolerance(self): + return self._loop_tolerance + + @loop_tolerance.setter + def loop_tolerance(self, value): + self._loop_tolerance = value + + @property + def network(self): + return self._network diff --git a/dingo/__init__.py b/dingo/__init__.py index 067c8b9..fe30ba5 100644 --- a/dingo/__init__.py +++ b/dingo/__init__.py @@ -26,6 +26,7 @@ from dingo.parser import dingo_args from dingo.MetabolicNetwork import MetabolicNetwork from dingo.PolytopeSampler import PolytopeSampler +from dingo.LooplessFluxSampler import LooplessFluxSampler from dingo.pyoptinterface_based_impl import fba, fva, inner_ball, remove_redundant_facets, set_default_solver diff --git a/dingo/nullspace.py b/dingo/nullspace.py index f6ca33e..fd9501a 100644 --- a/dingo/nullspace.py +++ b/dingo/nullspace.py @@ -8,7 +8,11 @@ import numpy as np from scipy import linalg -import sparseqr +try: + import sparseqr + _HAS_SPARSEQR = True +except ImportError: + _HAS_SPARSEQR = False import scipy.sparse.linalg @@ -43,6 +47,15 @@ def nullspace_sparse(Aeq, beq): """ N_shift = np.linalg.lstsq(Aeq, beq, rcond=None)[0] + + # Fallback when sparseqr (PySPQR / SuiteSparse) is not installed: + # use scipy's dense null_space instead + if not _HAS_SPARSEQR: + N = linalg.null_space(Aeq) + N = np.asarray(N, dtype="float") + N = np.ascontiguousarray(N, dtype="float") + return N, N_shift + Aeq = Aeq.T Aeq = scipy.sparse.csc_matrix(Aeq) diff --git a/dingo/pyoptinterface_based_impl.py b/dingo/pyoptinterface_based_impl.py index 6c83bf7..8a2b708 100644 --- a/dingo/pyoptinterface_based_impl.py +++ b/dingo/pyoptinterface_based_impl.py @@ -3,6 +3,13 @@ import numpy as np import sys +# Ensure HiGHS shared library is loaded (needed on some platforms) +try: + from pyoptinterface._src.highs import autoload_library + autoload_library() +except Exception: + pass + default_solver = "highs" def set_default_solver(solver_name): @@ -14,7 +21,7 @@ def get_solver(solver_name): if solver_name in solvers: return solvers[solver_name] else: - raise Exception("An unknown solver {solver_name} is requested.") + raise Exception(f"An unknown solver {solver_name} is requested.") def dot(c, x): return poi.quicksum(c[i] * x[i] for i in range(len(x)) if abs(c[i]) > 1e-12) @@ -79,12 +86,9 @@ def fba(lb, ub, S, c, solver_name=None): optimum_sol[i] = model.get_value(v[i]) return optimum_sol, optimum_value - except poi.TerminationStatusCode.NUMERICAL_ERROR as e: - print(f"A numerical error occurred: {e}") - except poi.TerminationStatusCode.OTHER_ERROR as e: - print(f"An error occurred: {e}") except Exception as e: print(f"An unexpected error occurred: {e}") + raise def fva(lb, ub, S, c, opt_percentage=100, solver_name=None): @@ -186,12 +190,9 @@ def fva(lb, ub, S, c, opt_percentage=100, solver_name=None): max_biomass_objective, ) - except poi.TerminationStatusCode.NUMERICAL_ERROR as e: - print(f"A numerical error occurred: {e}") - except poi.TerminationStatusCode.OTHER_ERROR as e: - print(f"An error occurred: {e}") except Exception as e: print(f"An unexpected error occurred: {e}") + raise def inner_ball(A, b, solver_name=None): @@ -488,9 +489,6 @@ def remove_redundant_facets(lb, ub, S, c, opt_percentage=100, solver_name=None): A_res = np.ascontiguousarray(A_res, dtype="float") return A_res, b_res, Aeq_res, beq_res - except poi.TerminationStatusCode.NUMERICAL_ERROR as e: - print(f"A numerical error occurred: {e}") - except poi.TerminationStatusCode.OTHER_ERROR as e: - print(f"An error occurred: {e}") except Exception as e: print(f"An unexpected error occurred: {e}") + raise diff --git a/tests/flux_sampling_ecoli.py b/tests/flux_sampling_ecoli.py new file mode 100644 index 0000000..4d74115 --- /dev/null +++ b/tests/flux_sampling_ecoli.py @@ -0,0 +1,252 @@ +# dingo : a python library for metabolic networks sampling and analysis +# dingo is part of GeomScale project + +# Copyright (c) 2024 + +# Licensed under GNU LGPL.3, see LICENCE file + +# ============================================================================= +# Easy Test: Flux Sampling Analysis on the E. coli Core Model +# ============================================================================= +# This script performs flux sampling under three biomass constraint scenarios: +# i. Optimal biomass growth (opt_percentage=100) +# ii. At least half of the optimal (opt_percentage=50) +# iii. Setting biomass free (opt_percentage=0) +# +# For each scenario we: +# - Load the e_coli_core model +# - Configure the biomass constraint +# - Sample steady states from the flux polytope +# - Validate sample shapes, non-zero fluxes, and biomass bounds +# ============================================================================= + +import unittest +import os +import sys +import numpy as np +from dingo import MetabolicNetwork, PolytopeSampler +from dingo.pyoptinterface_based_impl import set_default_solver + + +class TestFluxSamplingEcoli(unittest.TestCase): + """Flux sampling on e_coli_core under three biomass growth scenarios.""" + + MODEL_PATH = os.path.join(os.getcwd(), "ext_data", "e_coli_core.json") + + # ---------------------------------------------------------------- + # Helper utilities + # ---------------------------------------------------------------- + def _load_model(self): + """Load the E. coli core model and return it along with FBA optimum.""" + model = MetabolicNetwork.from_json(self.MODEL_PATH) + fba_res = model.fba() + max_biomass_objective = fba_res[1] + return model, max_biomass_objective + + def _sample_steady_states(self, model, method="billiard_walk", n_samples=500): + """Create a PolytopeSampler and generate steady states.""" + sampler = PolytopeSampler(model) + steady_states = sampler.generate_steady_states_no_multiphase( + method=method, n=n_samples, burn_in=100, thinning=2 + ) + return steady_states, sampler + + # ---------------------------------------------------------------- + # Scenario i: Optimal biomass growth (100%) + # ---------------------------------------------------------------- + def test_optimal_biomass_growth(self): + """ + Scenario i: Sample the flux space requiring optimal biomass growth. + All sampled steady states should achieve close to 100% of the FBA optimum. + """ + print("\n" + "=" * 70) + print("SCENARIO i: Optimal biomass growth (opt_percentage = 100)") + print("=" * 70) + + model, max_biomass = self._load_model() + + # Keep the default 100% optimality constraint + model.set_opt_percentage(100) + print(f" FBA max biomass objective: {max_biomass:.6f}") + + steady_states, sampler = self._sample_steady_states(model) + + # -- Validate shape: 95 reactions for e_coli_core -- + self.assertEqual(steady_states.shape[0], 95, + "Expected 95 reactions in e_coli_core model") + self.assertTrue(steady_states.shape[1] > 0, + "Expected at least some samples") + + # -- Validate non-zero fluxes -- + self.assertTrue(np.any(np.abs(steady_states) > 1e-12), + "Steady states should contain non-zero flux values") + + # -- Validate biomass constraint -- + biomass_idx = model.biomass_index + biomass_fluxes = steady_states[biomass_idx, :] + min_biomass_sample = np.min(biomass_fluxes) + + print(f" Biomass index: {biomass_idx}") + print(f" Biomass flux — mean: {np.mean(biomass_fluxes):.6f}, " + f"min: {min_biomass_sample:.6f}, max: {np.max(biomass_fluxes):.6f}") + print(f" Sample shape: {steady_states.shape}") + + # With 100% opt constraint, biomass should be close to max + self.assertTrue(min_biomass_sample >= max_biomass * 0.95, + f"Min biomass ({min_biomass_sample:.4f}) should be " + f">= 95% of max ({max_biomass:.4f})") + + # Save results for downstream analysis + np.save(os.path.join(os.getcwd(), "tests", "results_optimal.npy"), + steady_states) + + print(" [OK] Scenario i PASSED") + + # ---------------------------------------------------------------- + # Scenario ii: At least half-optimal biomass (50%) + # ---------------------------------------------------------------- + def test_half_optimal_biomass(self): + """ + Scenario ii: Sample the flux space requiring at least 50% of the + optimal biomass growth. + """ + print("\n" + "=" * 70) + print("SCENARIO ii: At least half-optimal biomass (opt_percentage = 50)") + print("=" * 70) + + model, max_biomass = self._load_model() + + # Set 50% optimality constraint + model.set_opt_percentage(50) + print(f" FBA max biomass objective: {max_biomass:.6f}") + print(f" Required minimum biomass: {max_biomass * 0.50:.6f}") + + steady_states, sampler = self._sample_steady_states(model) + + # -- Validate shape -- + self.assertEqual(steady_states.shape[0], 95) + self.assertTrue(steady_states.shape[1] > 0) + + # -- Validate non-zero fluxes -- + self.assertTrue(np.any(np.abs(steady_states) > 1e-12)) + + # -- Validate biomass constraint -- + biomass_idx = model.biomass_index + biomass_fluxes = steady_states[biomass_idx, :] + min_biomass_sample = np.min(biomass_fluxes) + + print(f" Biomass index: {biomass_idx}") + print(f" Biomass flux — mean: {np.mean(biomass_fluxes):.6f}, " + f"min: {min_biomass_sample:.6f}, max: {np.max(biomass_fluxes):.6f}") + print(f" Sample shape: {steady_states.shape}") + + # With 50% opt constraint, biomass should be at least half + self.assertTrue(min_biomass_sample >= max_biomass * 0.45, + f"Min biomass ({min_biomass_sample:.4f}) should be " + f">= 45% of max ({max_biomass:.4f}) (with tolerance)") + + # Save results + np.save(os.path.join(os.getcwd(), "tests", "results_half_optimal.npy"), + steady_states) + + print(" [OK] Scenario ii PASSED") + + # ---------------------------------------------------------------- + # Scenario iii: Biomass free (0% constraint) + # ---------------------------------------------------------------- + def test_biomass_free(self): + """ + Scenario iii: Sample the flux space with biomass unconstrained. + The biomass reaction can take any feasible value, including zero. + """ + print("\n" + "=" * 70) + print("SCENARIO iii: Biomass free (no biomass constraint)") + print("=" * 70) + + model, max_biomass = self._load_model() + + # Remove biomass constraint by setting opt_percentage to 0 + # This effectively removes the objective function lower bound + n = model.num_of_reactions() + model.set_opt_percentage(0) + + print(f" FBA max biomass objective: {max_biomass:.6f}") + print(f" Biomass constraint: NONE (free)") + + steady_states, sampler = self._sample_steady_states(model) + + # -- Validate shape -- + self.assertEqual(steady_states.shape[0], 95) + self.assertTrue(steady_states.shape[1] > 0) + + # -- Validate non-zero fluxes -- + self.assertTrue(np.any(np.abs(steady_states) > 1e-12)) + + # -- Examine biomass distribution -- + biomass_idx = model.biomass_index + biomass_fluxes = steady_states[biomass_idx, :] + + print(f" Biomass index: {biomass_idx}") + print(f" Biomass flux — mean: {np.mean(biomass_fluxes):.6f}, " + f"min: {np.min(biomass_fluxes):.6f}, max: {np.max(biomass_fluxes):.6f}") + print(f" Sample shape: {steady_states.shape}") + + # In the free case, we expect a wider range of biomass values + # The minimum should be considerably lower than the max + biomass_range = np.max(biomass_fluxes) - np.min(biomass_fluxes) + print(f" Biomass range: {biomass_range:.6f}") + + # Save results + np.save(os.path.join(os.getcwd(), "tests", "results_biomass_free.npy"), + steady_states) + + print(" [OK] Scenario iii PASSED") + + # ---------------------------------------------------------------- + # Summary comparison across all three scenarios + # ---------------------------------------------------------------- + def test_scenario_comparison(self): + """ + Compare basic statistics across the three sampling scenarios. + This test runs after the individual scenario tests. + """ + print("\n" + "=" * 70) + print("SUMMARY: Comparing all three scenarios") + print("=" * 70) + + model, max_biomass = self._load_model() + biomass_idx = model.biomass_index + reactions = model.reactions + + results = {} + for name, opt_pct in [("optimal", 100), ("half_optimal", 50), ("biomass_free", 0)]: + m, _ = self._load_model() + m.set_opt_percentage(opt_pct) + ss, _ = self._sample_steady_states(m, n_samples=300) + results[name] = ss + + # Print comparison table + print(f"\n {'Scenario':<20} {'Mean Biomass':>14} {'Min Biomass':>14} " + f"{'Max Biomass':>14} {'Flux Std Mean':>14}") + print(" " + "-" * 78) + + for name, ss in results.items(): + bm = ss[biomass_idx, :] + flux_std = np.mean(np.std(ss, axis=1)) + print(f" {name:<20} {np.mean(bm):>14.6f} {np.min(bm):>14.6f} " + f"{np.max(bm):>14.6f} {flux_std:>14.6f}") + + # Verify that the biomass-free scenario has wider flux variability + std_optimal = np.mean(np.std(results["optimal"], axis=1)) + std_free = np.mean(np.std(results["biomass_free"], axis=1)) + + print(f"\n Flux variability ratio (free/optimal): {std_free/std_optimal:.2f}") + print(" [OK] Comparison complete") + + +if __name__ == "__main__": + # Only treat positional (non-flag) arguments as solver name + if len(sys.argv) > 1 and not sys.argv[1].startswith("-"): + set_default_solver(sys.argv[1]) + sys.argv.pop(1) + unittest.main(verbosity=2) diff --git a/tests/nonconvex_sampling_pipeline.py b/tests/nonconvex_sampling_pipeline.py new file mode 100644 index 0000000..b1c0aea --- /dev/null +++ b/tests/nonconvex_sampling_pipeline.py @@ -0,0 +1,258 @@ +# dingo : a python library for metabolic networks sampling and analysis +# dingo is part of GeomScale project + +# Copyright (c) 2024 + +# Licensed under GNU LGPL.3, see LICENCE file + +# ============================================================================= +# Non-Convex Sampling Pipeline +# ============================================================================= +# End-to-end pipeline demonstrating: +# 1. Flux sampling under three biomass scenarios +# 2. Loop analysis and loopless sampling +# 3. Reaction clustering & biological interpretation +# 4. Comparison of standard vs loopless sample properties +# ============================================================================= + +import unittest +import os +import sys +import numpy as np +from dingo import MetabolicNetwork, PolytopeSampler +from dingo.LooplessFluxSampler import LooplessFluxSampler +from dingo.utils import correlated_reactions, cluster_corr_reactions +from dingo.pyoptinterface_based_impl import set_default_solver + + +class TestNonConvexPipeline(unittest.TestCase): + """End-to-end non-convex sampling pipeline on E. coli core.""" + + MODEL_PATH = os.path.join(os.getcwd(), "ext_data", "e_coli_core.json") + + # ---------------------------------------------------------------- + # Full pipeline test + # ---------------------------------------------------------------- + def test_full_pipeline(self): + """ + Complete pipeline: sample → detect loops → filter → cluster → interpret. + """ + print("\n" + "=" * 70) + print("PIPELINE: Non-Convex Sampling End-to-End") + print("=" * 70) + + model = MetabolicNetwork.from_json(self.MODEL_PATH) + reactions = model.reactions + fba_result = model.fba() + max_biomass = fba_result[1] + biomass_idx = model.biomass_index + + print(f"\n Model: E. coli core") + print(f" Reactions: {model.num_of_reactions()}") + print(f" Metabolites: {model.num_of_metabolites()}") + print(f" FBA optimal biomass: {max_biomass:.6f}") + + # ============================================================== + # STEP 1: Standard sampling under 3 biomass scenarios + # ============================================================== + print("\n" + "-" * 50) + print(" STEP 1: Standard flux sampling (3 scenarios)") + print("-" * 50) + + scenarios = { + "optimal (100%)": 100, + "half-optimal (50%)": 50, + "biomass-free (0%)": 0, + } + + standard_results = {} + for name, opt_pct in scenarios.items(): + m = MetabolicNetwork.from_json(self.MODEL_PATH) + m.set_opt_percentage(opt_pct) + sampler = PolytopeSampler(m) + ss = sampler.generate_steady_states_no_multiphase( + method="billiard_walk", n=300, burn_in=50, thinning=2 + ) + standard_results[name] = ss + bm = ss[biomass_idx, :] + print(f" {name}: {ss.shape[1]} samples, " + f"biomass mean={np.mean(bm):.4f}") + + # ============================================================== + # STEP 2: Loop analysis on standard samples + # ============================================================== + print("\n" + "-" * 50) + print(" STEP 2: Loop analysis on standard samples") + print("-" * 50) + + lfs = LooplessFluxSampler(model, loop_tolerance=1e-3) + + print(f" Internal reactions: {len(lfs.internal_indices)}") + print(f" Exchange reactions: {len(lfs.exchange_indices)}") + print(f" Cycle nullspace dim: " + f"{lfs.cycle_nullspace.shape[1] if lfs.cycle_nullspace.size > 0 else 0}") + + for name, ss in standard_results.items(): + analysis = lfs.analyze_loops(ss) + print(f" {name}: " + f"loops={analysis['fraction_with_loops']*100:.1f}%, " + f"mean_norm={analysis['mean_loop_norm']:.6f}") + + # ============================================================== + # STEP 3: Loopless sampling (rejection-based) + # ============================================================== + print("\n" + "-" * 50) + print(" STEP 3: Loopless sampling (rejection)") + print("-" * 50) + + loopless_ss, rej_stats = lfs.sample_loopless( + n_samples=200, + max_attempts_factor=20, + method="billiard_walk", + opt_percentage=100, + ) + + print(f" Accepted: {rej_stats['total_accepted']}") + print(f" Rejected: {rej_stats['total_rejected']}") + print(f" Acceptance rate: {rej_stats['acceptance_rate']:.4f}") + + # ============================================================== + # STEP 4: Penalty-weighted sampling + # ============================================================== + print("\n" + "-" * 50) + print(" STEP 4: Penalty-weighted sampling") + print("-" * 50) + + penalized_ss, weights, pen_stats = lfs.sample_penalized( + n_samples=300, + penalty_weight=10.0, + method="billiard_walk", + opt_percentage=100, + ) + + print(f" Effective sample size: {pen_stats['effective_sample_size']:.1f}") + print(f" N loopless: {pen_stats['n_effectively_loopless']}") + + # ============================================================== + # STEP 5: Reaction clustering on standard vs loopless + # ============================================================== + print("\n" + "-" * 50) + print(" STEP 5: Reaction clustering comparison") + print("-" * 50) + + for label, ss in [("standard", standard_results["optimal (100%)"]), + ("loopless", loopless_ss)]: + if ss.shape[1] < 10: + print(f" {label}: insufficient samples, skipping") + continue + + corr_matrix = correlated_reactions( + ss, + reactions=reactions, + pearson_cutoff=0.0, + indicator_cutoff=0, + cells=10, + lower_triangle=False, + ) + + _, labels, clusters = cluster_corr_reactions( + corr_matrix, reactions, linkage="ward", t=4.0 + ) + + multi_clusters = [c for c in clusters if len(c) > 1] + print(f" {label}: {len(clusters)} total clusters, " + f"{len(multi_clusters)} non-trivial") + for i, c in enumerate(multi_clusters[:3]): + print(f" Cluster {i+1}: {c[:4]}{'...' if len(c)>4 else ''} " + f"(n={len(c)})") + + # ============================================================== + # STEP 6: Biological interpretation summary + # ============================================================== + print("\n" + "-" * 50) + print(" STEP 6: Biological interpretation") + print("-" * 50) + + # Compare key metabolic pathway statistics + pathways = { + "Glycolysis": ["PFK", "PYK", "GAPD", "PGK", "ENO", "PGM"], + "TCA Cycle": ["CS", "ACONTa", "ACONTb", "AKGDH", "SUCOAS", + "FUM", "MDH", "ICDHyr"], + "Pentose Phosphate": ["G6PDH2r", "GND", "RPI", "TKT1", "TKT2", + "TALA"], + "Respiration": ["NADH16", "CYTBD", "ATPS4r"], + } + + standard_opt = standard_results["optimal (100%)"] + + print(f"\n {'Pathway':<22} {'Mean±Std (Standard)':>24} " + f"{'Mean±Std (Loopless)':>24}") + print(" " + "-" * 72) + + for pathway_name, rxn_list in pathways.items(): + std_fluxes = [] + ll_fluxes = [] + for rxn in rxn_list: + if rxn in reactions: + idx = reactions.index(rxn) + std_fluxes.append(np.mean(np.abs(standard_opt[idx, :]))) + if loopless_ss.shape[1] > 0: + ll_fluxes.append( + np.mean(np.abs(loopless_ss[idx, :])) + ) + + if std_fluxes: + std_mean = np.mean(std_fluxes) + std_std = np.std(std_fluxes) + if ll_fluxes: + ll_mean = np.mean(ll_fluxes) + ll_std = np.std(ll_fluxes) + print(f" {pathway_name:<22} " + f"{std_mean:>10.4f} ± {std_std:<10.4f} " + f"{ll_mean:>10.4f} ± {ll_std:<10.4f}") + else: + print(f" {pathway_name:<22} " + f"{std_mean:>10.4f} ± {std_std:<10.4f} " + f"{'N/A':>24}") + + # ============================================================== + # STEP 7: Assessment of non-convex sampling necessity + # ============================================================== + print("\n" + "-" * 50) + print(" STEP 7: Assessment — Are non-convex methods needed?") + print("-" * 50) + + opt_analysis = lfs.analyze_loops(standard_opt) + loop_fraction = opt_analysis["fraction_with_loops"] + mean_norm = opt_analysis["mean_loop_norm"] + + print(f" Fraction of standard samples with loops: " + f"{loop_fraction:.4f}") + print(f" Mean loop norm: {mean_norm:.6f}") + + if loop_fraction > 0.1: + print(f"\n CONCLUSION: {loop_fraction*100:.1f}% of standard " + f"samples contain internal loops.") + print(" Non-convex sampling IS recommended for this model.") + print(" Loopless filtering removes thermodynamically " + "infeasible solutions.") + elif loop_fraction > 0.01: + print(f"\n CONCLUSION: {loop_fraction*100:.1f}% of samples " + f"have loops — moderate impact.") + print(" Non-convex sampling is BENEFICIAL but not critical.") + else: + print(f"\n CONCLUSION: Only {loop_fraction*100:.2f}% of samples " + f"have loops — negligible impact.") + print(" Standard convex sampling appears sufficient for " + "this model.") + + print("\n" + "=" * 70) + print(" [OK] FULL PIPELINE COMPLETED SUCCESSFULLY") + print("=" * 70) + + +if __name__ == "__main__": + if len(sys.argv) > 1 and not sys.argv[1].startswith("-"): + set_default_solver(sys.argv[1]) + sys.argv.pop(1) + unittest.main(verbosity=2) diff --git a/tests/reaction_clustering_analysis.py b/tests/reaction_clustering_analysis.py new file mode 100644 index 0000000..a3ed861 --- /dev/null +++ b/tests/reaction_clustering_analysis.py @@ -0,0 +1,313 @@ +# dingo : a python library for metabolic networks sampling and analysis +# dingo is part of GeomScale project + +# Copyright (c) 2024 + +# Licensed under GNU LGPL.3, see LICENCE file + +# ============================================================================= +# Medium Test: Reaction Clustering Analysis +# ============================================================================= +# This script exploits dingo's functionalities to discover reaction clusters +# that differentiate the three biomass constraint scenarios: +# i. opt_percentage = 100 (optimal biomass) +# ii. opt_percentage = 50 (half-optimal) +# iii. opt_percentage = 0 (biomass free) +# +# Methods used: +# - Pearson correlation with copula-based filtering +# - Hierarchical clustering +# - Cross-scenario comparison of cluster composition +# ============================================================================= + +import unittest +import os +import sys +import numpy as np +from dingo import MetabolicNetwork, PolytopeSampler +from dingo.utils import correlated_reactions, cluster_corr_reactions +from dingo.pyoptinterface_based_impl import set_default_solver + + +class TestReactionClustering(unittest.TestCase): + """Reaction clustering analysis across three biomass scenarios.""" + + MODEL_PATH = os.path.join(os.getcwd(), "ext_data", "e_coli_core.json") + + # ---------------------------------------------------------------- + # Helpers + # ---------------------------------------------------------------- + def _sample_scenario(self, opt_percentage, n_samples=500): + """Sample steady states for a given opt_percentage.""" + model = MetabolicNetwork.from_json(self.MODEL_PATH) + model.set_opt_percentage(opt_percentage) + sampler = PolytopeSampler(model) + steady_states = sampler.generate_steady_states_no_multiphase( + method="billiard_walk", n=n_samples, burn_in=100, thinning=2 + ) + return model, steady_states + + # ---------------------------------------------------------------- + # Test: Correlation matrices for each scenario + # ---------------------------------------------------------------- + def test_correlation_matrices(self): + """ + Compute and validate correlation matrices for all three scenarios. + """ + print("\n" + "=" * 70) + print("MEDIUM: Computing correlation matrices for each scenario") + print("=" * 70) + + scenarios = { + "optimal (100%)": 100, + "half_optimal (50%)": 50, + "biomass_free (0%)": 0, + } + + for name, opt_pct in scenarios.items(): + print(f"\n --- {name} ---") + model, steady_states = self._sample_scenario(opt_pct) + reactions = model.reactions + + # Compute correlation matrix without copula filtering + corr_matrix = correlated_reactions( + steady_states, + reactions=reactions, + pearson_cutoff=0.0, + indicator_cutoff=0, + cells=10, + cop_coeff=0.3, + lower_triangle=False, + ) + + # Validate the correlation matrix + self.assertEqual(corr_matrix.shape[0], len(reactions)) + self.assertEqual(corr_matrix.shape[1], len(reactions)) + self.assertAlmostEqual(np.trace(corr_matrix), len(reactions), places=5) + + # Count highly correlated pairs (|r| > 0.9) + upper_tri = np.triu(corr_matrix, k=1) + n_high_corr = np.sum(np.abs(upper_tri) > 0.9) + + print(f" Matrix shape: {corr_matrix.shape}") + print(f" Highly correlated pairs (|r|>0.9): {n_high_corr}") + print(f" Mean abs correlation: {np.mean(np.abs(upper_tri)):.4f}") + + print("\n [OK] Correlation matrix computation PASSED for all scenarios") + + # ---------------------------------------------------------------- + # Test: Copula-based filtered correlation + # ---------------------------------------------------------------- + def test_copula_filtered_correlation(self): + """ + Compute correlation matrices with copula indicator filtering for + each scenario and compare the number of truly correlated pairs. + """ + print("\n" + "=" * 70) + print("MEDIUM: Copula-filtered correlation analysis") + print("=" * 70) + + scenarios = { + "optimal (100%)": 100, + "half_optimal (50%)": 50, + "biomass_free (0%)": 0, + } + + corr_results = {} + + for name, opt_pct in scenarios.items(): + print(f"\n --- {name} ---") + model, steady_states = self._sample_scenario(opt_pct) + reactions = model.reactions + + # Compute correlation matrix WITH copula indicator filtering + corr_matrix, indicator_dict = correlated_reactions( + steady_states, + reactions=reactions, + pearson_cutoff=0.90, + indicator_cutoff=5, + cells=10, + cop_coeff=0.3, + lower_triangle=False, + verbose=False, + ) + + corr_results[name] = { + "corr_matrix": corr_matrix, + "indicator_dict": indicator_dict, + "reactions": reactions, + } + + # Count positive and negative correlations + n_positive = sum( + 1 for v in indicator_dict.values() + if v["classification"] == "positive" + ) + n_negative = sum( + 1 for v in indicator_dict.values() + if v["classification"] == "negative" + ) + n_none = sum( + 1 for v in indicator_dict.values() + if v["classification"] == "no correlation" + ) + + print(f" Positively correlated pairs: {n_positive}") + print(f" Negatively correlated pairs: {n_negative}") + print(f" No correlation (filtered): {n_none}") + + self.assertEqual(corr_matrix.shape[0], len(reactions)) + + print("\n [OK] Copula-filtered correlation analysis PASSED") + + # ---------------------------------------------------------------- + # Test: Hierarchical clustering and differentiation + # ---------------------------------------------------------------- + def test_hierarchical_clustering(self): + """ + Perform hierarchical clustering on each scenario and compare + which reaction clusters appear, disappear, or change across scenarios. + """ + print("\n" + "=" * 70) + print("MEDIUM: Hierarchical clustering & cluster differentiation") + print("=" * 70) + + scenarios = { + "optimal": 100, + "half_optimal": 50, + "biomass_free": 0, + } + + all_clusters = {} + + for name, opt_pct in scenarios.items(): + print(f"\n --- {name} (opt={opt_pct}%) ---") + model, steady_states = self._sample_scenario(opt_pct, n_samples=500) + reactions = model.reactions + + # Compute correlation matrix (unfiltered for clustering) + corr_matrix = correlated_reactions( + steady_states, + reactions=reactions, + pearson_cutoff=0.0, + indicator_cutoff=0, + cells=10, + lower_triangle=False, + ) + + # Perform hierarchical clustering + dissimilarity, labels, clusters = cluster_corr_reactions( + corr_matrix, + reactions, + linkage="ward", + t=4.0, + correction=True, + ) + + all_clusters[name] = clusters + + print(f" Number of clusters: {len(clusters)}") + for i, cluster in enumerate(clusters): + if len(cluster) > 1: + print(f" Cluster {i+1} ({len(cluster)} reactions): " + f"{cluster[:5]}{'...' if len(cluster) > 5 else ''}") + + # Validate clustering output + self.assertTrue(len(clusters) > 0, "Should find at least one cluster") + total_reactions_in_clusters = sum(len(c) for c in clusters) + self.assertEqual(total_reactions_in_clusters, len(reactions), + "All reactions should be assigned to clusters") + + # --------------------------------------------------------------- + # Compare clusters across scenarios + # --------------------------------------------------------------- + print("\n --- Cross-scenario cluster comparison ---") + + # Convert clusters to sets for comparison + for name, clusters in all_clusters.items(): + cluster_sets = [frozenset(c) for c in clusters] + all_clusters[name] = cluster_sets + + # Find clusters unique to each scenario + for name in all_clusters: + other_names = [n for n in all_clusters if n != name] + unique_clusters = [] + for c in all_clusters[name]: + if len(c) > 1: + is_unique = True + for other in other_names: + if c in all_clusters[other]: + is_unique = False + break + if is_unique: + unique_clusters.append(c) + + if unique_clusters: + print(f"\n Clusters unique to {name}:") + for c in unique_clusters[:3]: + print(f" {list(c)[:5]}{'...' if len(c) > 5 else ''} " + f"(size={len(c)})") + else: + print(f"\n No strictly unique clusters for {name} (clusters may overlap)") + + print("\n [OK] Hierarchical clustering & differentiation PASSED") + + # ---------------------------------------------------------------- + # Test: Key reaction flux distributions across scenarios + # ---------------------------------------------------------------- + def test_flux_distribution_comparison(self): + """ + Compare the flux distributions of key reactions across scenarios + to identify biologically meaningful differences. + """ + print("\n" + "=" * 70) + print("MEDIUM: Flux distribution comparison for key reactions") + print("=" * 70) + + # Key reactions in E. coli core to examine + key_reaction_names = [ + "PFK", # Phosphofructokinase (glycolysis) + "PYK", # Pyruvate kinase (glycolysis) + "CS", # Citrate synthase (TCA cycle) + "AKGDH", # Alpha-ketoglutarate dehydrogenase (TCA) + "PPC", # Phosphoenolpyruvate carboxylase + "ATPM", # ATP maintenance requirement + ] + + scenarios = {"optimal": 100, "half_optimal": 50, "biomass_free": 0} + scenario_results = {} + + for name, opt_pct in scenarios.items(): + model, steady_states = self._sample_scenario(opt_pct, n_samples=500) + scenario_results[name] = { + "model": model, + "steady_states": steady_states, + } + + # Compare distributions + print(f"\n {'Reaction':<12} {'Scenario':<18} {'Mean':>10} " + f"{'Std':>10} {'Min':>10} {'Max':>10}") + print(" " + "-" * 72) + + reactions = scenario_results["optimal"]["model"].reactions + + for rxn_name in key_reaction_names: + if rxn_name in reactions: + rxn_idx = reactions.index(rxn_name) + + for s_name, data in scenario_results.items(): + fluxes = data["steady_states"][rxn_idx, :] + print(f" {rxn_name:<12} {s_name:<18} {np.mean(fluxes):>10.4f} " + f"{np.std(fluxes):>10.4f} {np.min(fluxes):>10.4f} " + f"{np.max(fluxes):>10.4f}") + + print() # blank line between reactions + + print(" [OK] Flux distribution comparison PASSED") + + +if __name__ == "__main__": + if len(sys.argv) > 1 and not sys.argv[1].startswith("-"): + set_default_solver(sys.argv[1]) + sys.argv.pop(1) + unittest.main(verbosity=2) diff --git a/tests/test_loopless_sampling.py b/tests/test_loopless_sampling.py new file mode 100644 index 0000000..04dc140 --- /dev/null +++ b/tests/test_loopless_sampling.py @@ -0,0 +1,306 @@ +# dingo : a python library for metabolic networks sampling and analysis +# dingo is part of GeomScale project + +# Copyright (c) 2024 + +# Licensed under GNU LGPL.3, see LICENCE file + +# ============================================================================= +# Hard Test: Loopless (Non-Convex) Sampling +# ============================================================================= +# Tests the LooplessFluxSampler module for: +# - Internal reaction classification +# - Cycle nullspace computation +# - Rejection-based loopless sampling +# - Penalty-weighted sampling +# - Comparison of loopless vs standard samples +# ============================================================================= + +import unittest +import os +import sys +import numpy as np +from dingo import MetabolicNetwork, PolytopeSampler +from dingo.LooplessFluxSampler import LooplessFluxSampler +from dingo.pyoptinterface_based_impl import set_default_solver + + +class TestLooplessSampling(unittest.TestCase): + """Tests for the LooplessFluxSampler module.""" + + MODEL_PATH = os.path.join(os.getcwd(), "ext_data", "e_coli_core.json") + + # ---------------------------------------------------------------- + # Test: Reaction classification + # ---------------------------------------------------------------- + def test_reaction_classification(self): + """ + Verify that reactions are correctly classified as internal vs exchange. + """ + print("\n" + "=" * 70) + print("HARD: Reaction classification (internal vs exchange)") + print("=" * 70) + + model = MetabolicNetwork.from_json(self.MODEL_PATH) + lfs = LooplessFluxSampler(model) + + n_internal = len(lfs.internal_indices) + n_exchange = len(lfs.exchange_indices) + n_total = model.num_of_reactions() + + print(f" Total reactions: {n_total}") + print(f" Internal reactions: {n_internal}") + print(f" Exchange reactions: {n_exchange}") + + # All reactions must be classified + self.assertEqual(n_internal + n_exchange, n_total, + "All reactions should be classified") + + # E. coli core has ~20 exchange reactions + self.assertTrue(n_exchange > 10, + "E. coli core should have >10 exchange reactions") + self.assertTrue(n_internal > 50, + "E. coli core should have >50 internal reactions") + + # List some exchange reactions for verification + reactions = model.reactions + print("\n Exchange reactions:") + for idx in lfs.exchange_indices[:10]: + print(f" {reactions[idx]}") + + print("\n [OK] Reaction classification PASSED") + + # ---------------------------------------------------------------- + # Test: Cycle nullspace computation + # ---------------------------------------------------------------- + def test_cycle_nullspace(self): + """ + Verify that the internal cycle nullspace is computed correctly. + """ + print("\n" + "=" * 70) + print("HARD: Cycle nullspace computation") + print("=" * 70) + + model = MetabolicNetwork.from_json(self.MODEL_PATH) + lfs = LooplessFluxSampler(model) + + N = lfs.cycle_nullspace + + if N.size > 0: + print(f" Cycle nullspace shape: {N.shape}") + print(f" Number of potential cycle directions: {N.shape[1]}") + + # Verify orthogonality of nullspace vectors + NtN = N.T @ N + identity_check = np.allclose(NtN, np.eye(N.shape[1]), atol=1e-10) + print(f" Nullspace vectors orthonormal: {identity_check}") + + # Verify that nullspace vectors are in the nullspace of S_int + S_int = model.S[:, lfs.internal_indices] + residual = np.linalg.norm(S_int @ N) + print(f" S_int * N residual norm: {residual:.2e}") + self.assertTrue(residual < 1e-10, + "Nullspace vectors should satisfy S_int * N ≈ 0") + else: + print(" No internal cycles found (nullspace is empty)") + + print("\n [OK] Cycle nullspace computation PASSED") + + # ---------------------------------------------------------------- + # Test: Loop detection on known samples + # ---------------------------------------------------------------- + def test_loop_detection(self): + """ + Generate standard samples and check what fraction contains loops. + """ + print("\n" + "=" * 70) + print("HARD: Loop detection in standard samples") + print("=" * 70) + + model = MetabolicNetwork.from_json(self.MODEL_PATH) + lfs = LooplessFluxSampler(model) + + # Generate standard (possibly loopy) samples + sampler = PolytopeSampler(model) + steady_states = sampler.generate_steady_states_no_multiphase( + method="billiard_walk", n=300, burn_in=50, thinning=2 + ) + + # Analyze loops + analysis = lfs.analyze_loops(steady_states) + + print(f" Samples analyzed: {steady_states.shape[1]}") + print(f" Fraction with loops: {analysis['fraction_with_loops']:.4f}") + print(f" Mean loop norm: {analysis['mean_loop_norm']:.6f}") + print(f" Median loop norm: {analysis['median_loop_norm']:.6f}") + print(f" Max loop norm: {analysis['max_loop_norm']:.6f}") + print(f" Internal reactions: {analysis['n_internal_reactions']}") + print(f" Exchange reactions: {analysis['n_exchange_reactions']}") + print(f" Cycle nullspace dim: {analysis['cycle_nullspace_dim']}") + + # Validate analysis output + self.assertTrue(0 <= analysis["fraction_with_loops"] <= 1) + self.assertTrue(analysis["mean_loop_norm"] >= 0) + + print("\n [OK] Loop detection PASSED") + + # ---------------------------------------------------------------- + # Test: Rejection-based loopless sampling + # ---------------------------------------------------------------- + def test_rejection_sampling(self): + """ + Test the rejection-based loopless sampling approach. + """ + print("\n" + "=" * 70) + print("HARD: Rejection-based loopless sampling") + print("=" * 70) + + model = MetabolicNetwork.from_json(self.MODEL_PATH) + lfs = LooplessFluxSampler(model, loop_tolerance=1e-3) + + loopless_samples, stats = lfs.sample_loopless( + n_samples=100, + max_attempts_factor=20, + method="billiard_walk", + burn_in=50, + thinning=2, + opt_percentage=100, + ) + + print(f" Total sampled: {stats['total_sampled']}") + print(f" Total accepted: {stats['total_accepted']}") + print(f" Total rejected: {stats['total_rejected']}") + print(f" Acceptance rate: {stats['acceptance_rate']:.4f}") + print(f" Mean loop norm: {stats['mean_loop_norm']:.6f}") + print(f" Max loop norm: {stats['max_loop_norm']:.6f}") + print(f" Loop tolerance: {stats['loop_tolerance']}") + print(f" Output shape: {loopless_samples.shape}") + + # Validate that accepted samples have correct shape + self.assertEqual(loopless_samples.shape[0], 95, + "Expected 95 reactions") + + # Verify no loops in accepted samples + if loopless_samples.shape[1] > 0: + for j in range(loopless_samples.shape[1]): + has_loop, norm = lfs._has_active_loop(loopless_samples[:, j]) + self.assertFalse(has_loop, + f"Sample {j} should be loopless (norm={norm:.6f})") + + print(f"\n [OK] All {loopless_samples.shape[1]} accepted samples " + f"verified loopless") + else: + print("\n [!] No samples accepted (acceptance rate may be very low)") + + print(" [OK] Rejection-based loopless sampling PASSED") + + # ---------------------------------------------------------------- + # Test: Penalty-weighted sampling + # ---------------------------------------------------------------- + def test_penalty_sampling(self): + """ + Test the penalty-weighted (importance sampling) approach. + """ + print("\n" + "=" * 70) + print("HARD: Penalty-weighted loopless sampling") + print("=" * 70) + + model = MetabolicNetwork.from_json(self.MODEL_PATH) + lfs = LooplessFluxSampler(model) + + samples, weights, stats = lfs.sample_penalized( + n_samples=300, + penalty_weight=10.0, + method="billiard_walk", + burn_in=50, + thinning=2, + opt_percentage=100, + ) + + print(f" Samples shape: {samples.shape}") + print(f" Weights shape: {weights.shape}") + print(f" Mean loop norm: {stats['mean_loop_norm']:.6f}") + print(f" Max loop norm: {stats['max_loop_norm']:.6f}") + print(f" N effectively loopless: {stats['n_effectively_loopless']}") + print(f" Effective sample size: {stats['effective_sample_size']:.1f}") + print(f" Penalty weight: {stats['penalty_weight']}") + + # Weights must sum to 1 + self.assertAlmostEqual(np.sum(weights), 1.0, places=10) + + # All weights should be non-negative + self.assertTrue(np.all(weights >= 0)) + + # Effective sample size should be > 0 + self.assertTrue(stats["effective_sample_size"] > 0) + + # Weighted mean of biomass should be reasonable + biomass_idx = model.biomass_index + weighted_mean_biomass = np.average( + samples[biomass_idx, :], weights=weights + ) + print(f" Weighted mean biomass flux: {weighted_mean_biomass:.6f}") + + print("\n [OK] Penalty-weighted sampling PASSED") + + # ---------------------------------------------------------------- + # Test: Standard vs loopless comparison + # ---------------------------------------------------------------- + def test_standard_vs_loopless(self): + """ + Compare standard and loopless samples on key flux statistics. + """ + print("\n" + "=" * 70) + print("HARD: Standard vs Loopless sampling comparison") + print("=" * 70) + + model = MetabolicNetwork.from_json(self.MODEL_PATH) + lfs = LooplessFluxSampler(model, loop_tolerance=1e-3) + + # Standard sampling + sampler = PolytopeSampler(model) + standard_ss = sampler.generate_steady_states_no_multiphase( + method="billiard_walk", n=300, burn_in=50, thinning=2 + ) + + # Loopless sampling + loopless_ss, stats = lfs.sample_loopless( + n_samples=200, + max_attempts_factor=20, + method="billiard_walk", + burn_in=50, + thinning=2, + opt_percentage=100, + ) + + # Analyze loops in both + standard_analysis = lfs.analyze_loops(standard_ss) + + print(f"\n {'Metric':<30} {'Standard':>15} {'Loopless':>15}") + print(" " + "-" * 62) + print(f" {'N samples':<30} {standard_ss.shape[1]:>15} " + f"{loopless_ss.shape[1]:>15}") + print(f" {'Fraction with loops':<30} " + f"{standard_analysis['fraction_with_loops']:>15.4f} {'0.0000':>15}") + + if loopless_ss.shape[1] > 0: + reactions = model.reactions + key_rxns = ["PFK", "CS", "ATPM"] + + for rxn in key_rxns: + if rxn in reactions: + idx = reactions.index(rxn) + std_mean = np.mean(standard_ss[idx, :]) + ll_mean = np.mean(loopless_ss[idx, :]) + print(f" {'Mean flux ' + rxn:<30} " + f"{std_mean:>15.4f} {ll_mean:>15.4f}") + + print(f"\n Acceptance rate: {stats['acceptance_rate']:.4f}") + print("\n [OK] Standard vs Loopless comparison PASSED") + + +if __name__ == "__main__": + if len(sys.argv) > 1 and not sys.argv[1].startswith("-"): + set_default_solver(sys.argv[1]) + sys.argv.pop(1) + unittest.main(verbosity=2) diff --git a/volestipy.py b/volestipy.py new file mode 100644 index 0000000..c3c2353 --- /dev/null +++ b/volestipy.py @@ -0,0 +1,133 @@ +# volestipy mock/stub module +# ============================================================================= +# This provides a minimal mock of the volestipy C++ extension when the +# native build is not available (e.g. on Windows without Boost/lp_solve). +# The mock HPolytope uses cobra's sampling when available, or falls back +# to a simple hit-and-run sampler implemented in pure Python. +# ============================================================================= + +import numpy as np +from scipy.optimize import linprog + + +class HPolytope: + """Mock HPolytope when volestipy C++ extension is not available.""" + + def __init__(self, A, b): + self.A = np.array(A, dtype=np.float64) + self.b = np.array(b, dtype=np.float64) + self.dim = A.shape[1] + + def _find_interior_point(self): + """Find an interior point using multiple strategies.""" + m, n = self.A.shape + + # Strategy 1: Chebyshev center (max inscribed ball) + try: + norms = np.linalg.norm(self.A, axis=1) + c_lp = np.zeros(n + 1) + c_lp[-1] = -1 # maximize r => minimize -r + + A_lp = np.hstack([self.A, norms.reshape(-1, 1)]) + bounds = [(None, None)] * n + [(0, None)] + + res = linprog(c_lp, A_ub=A_lp, b_ub=self.b, bounds=bounds, method='highs') + if res.success and res.x[-1] > 1e-12: + return res.x[:n] + except Exception: + pass + + # Strategy 2: Find a feasible point via linprog + try: + c_feas = np.zeros(n) + res = linprog(c_feas, A_ub=self.A, b_ub=self.b, method='highs') + if res.success: + x = res.x + # Check strict feasibility + if np.all(self.A @ x <= self.b + 1e-10): + return x + except Exception: + pass + + # Strategy 3: Try origin or midpoint of bounds + if np.all(self.A @ np.zeros(n) <= self.b + 1e-10): + return np.zeros(n) + + return None + + def _hit_and_run_step(self, x): + """One step of coordinate hit-and-run from point x.""" + n = self.dim + # Random direction (coordinate direction for CDHR) + coord = np.random.randint(n) + d = np.zeros(n) + d[coord] = 1.0 + + # Compute step range: A*(x + t*d) <= b => t*(A*d) <= b - A*x + Ad = self.A @ d + residual = self.b - self.A @ x + + t_min = -np.inf + t_max = np.inf + + for i in range(len(Ad)): + if Ad[i] > 1e-12: + t_max = min(t_max, residual[i] / Ad[i]) + elif Ad[i] < -1e-12: + t_min = max(t_min, residual[i] / Ad[i]) + + if t_min >= t_max: + return x # infeasible direction, stay + + t = np.random.uniform(t_min, t_max) + return x + t * d + + def generate_samples(self, method, n, burn_in, thinning, variance, bias_vector, solver=None): + """Generate samples using pure-Python CDHR (coordinate hit-and-run). + Returns shape (n_samples, dim) to match real volestipy API. + """ + x = self._find_interior_point() + if x is None: + raise RuntimeError("Could not find interior point for polytope") + + # Burn-in + for _ in range(max(burn_in, 100)): + x = self._hit_and_run_step(x) + + # Sample — real volestipy returns (n_samples, dim) + samples = np.zeros((n, self.dim)) + for i in range(n): + for _ in range(max(thinning, 1)): + x = self._hit_and_run_step(x) + samples[i, :] = x + + return samples + + def mmcs(self, ess, psrf, parallel_mmcs, num_threads, solver=None): + """Mock MMCS: falls back to hit-and-run with identity rounding.""" + n_samples = max(ess * 2, 1000) + x = self._find_interior_point() + if x is None: + raise RuntimeError("Could not find interior point for polytope") + + # Burn-in + for _ in range(200): + x = self._hit_and_run_step(x) + + # Sample + samples = np.zeros((self.dim, n_samples)) + for i in range(n_samples): + for _ in range(2): + x = self._hit_and_run_step(x) + samples[:, i] = x + + Tr = np.eye(self.dim) + Tr_shift = np.zeros(self.dim) + + return self.A, self.b, Tr, Tr_shift, samples + + def rounding(self, method="john_position", solver=None): + """Mock rounding: returns identity transformation.""" + Tr = np.eye(self.dim) + Tr_shift = np.zeros(self.dim) + return self.A, self.b, Tr, Tr_shift, 1.0