diff --git a/dingo/dynamic_volume.py b/dingo/dynamic_volume.py new file mode 100644 index 0000000..7ecc33a --- /dev/null +++ b/dingo/dynamic_volume.py @@ -0,0 +1,443 @@ +"""Wrapper for the supplied static relative-volume routine. + +This module does not implement a second volume algorithm. For every incoming +linear cut it calls MATLAB ``relative_volume`` from the supplied +MetabolicPolytopesCuttingPlanes source, adds that cut to the parent polytope, +and repeats. No FBA objective or biomass floor is introduced. +""" +from __future__ import annotations + +from dataclasses import dataclass +import os +from pathlib import Path +import shutil +import subprocess +import tempfile +from typing import Iterable, Mapping + +import numpy as np +from scipy.io import loadmat, savemat + + +_STATIC_OPTION_NAMES = { + "r", + "delta", + "alpha", + "nu", + "N_utest", + "window_size", + "epsilon", + "max_annealing_iter", + "max_bisect_iter", + "simdLen", + "t_max_safety", + "verb", +} +_WRAPPER_OPTION_NAMES = {"algorithm", "wrapper_timeout_sec"} +_SOURCE_DEFAULTS = { + "r": 0.10, + "delta": 0.05, + "alpha": 0.20, + "nu": 10, + "N_utest": 320, + "window_size": 5000, + "epsilon": 0.10, + "max_annealing_iter": 200, + "max_bisect_iter": 100, + "simdLen": 8, + "t_max_safety": 1.01, + "verb": 1, +} +_POSITIVE_INTEGER_OPTIONS = { + "nu", + "N_utest", + "window_size", + "max_annealing_iter", + "max_bisect_iter", + "simdLen", +} + + +@dataclass(frozen=True) +class HalfspaceCut: + normal: np.ndarray + threshold: float + label: str = "" + + def __post_init__(self) -> None: + normal = np.asarray(self.normal, dtype=float).reshape(-1) + if normal.size == 0 or not np.all(np.isfinite(normal)): + raise ValueError("cut normal must be a nonempty finite vector") + if np.linalg.norm(normal) == 0: + raise ValueError("cut normal must be nonzero") + if not np.isfinite(self.threshold): + raise ValueError("cut threshold must be finite") + object.__setattr__(self, "normal", normal) + + +@dataclass(frozen=True) +class StaticVolumeResult: + rho: float + log_rho: float + cut_ratios: np.ndarray + cut_log_ratios: np.ndarray + cut_n_phases: np.ndarray + cut_total_steps: np.ndarray + cut_converged: np.ndarray + cut_labels: tuple[str, ...] + algorithm: str = "static_relative_volume_sequential" + objective_floor: bool = False + + +def cuts_from_bounds( + parent_lb, + parent_ub, + child_lb, + child_ub, + *, + atol: float = 1e-12, +) -> list[HalfspaceCut]: + """Convert a nested bound update into sequential scalar halfspace cuts.""" + parent_lb = _vector(parent_lb, "parent_lb") + parent_ub = _vector(parent_ub, "parent_ub") + child_lb = _vector(child_lb, "child_lb") + child_ub = _vector(child_ub, "child_ub") + n = parent_lb.size + if any(x.size != n for x in (parent_ub, child_lb, child_ub)): + raise ValueError("parent and child bound vectors have different lengths") + if atol < 0 or not np.isfinite(atol): + raise ValueError("atol must be nonnegative and finite") + if np.any(parent_lb > parent_ub) or np.any(child_lb > child_ub): + raise ValueError("invalid parent or child bounds") + if np.any(child_lb < parent_lb - atol) or np.any(child_ub > parent_ub + atol): + raise ValueError("child bounds must define a subset of the parent bounds") + cuts: list[HalfspaceCut] = [] + for j in range(n): + if child_ub[j] < parent_ub[j] - atol: + normal = np.zeros(n) + normal[j] = 1.0 + cuts.append(HalfspaceCut(normal, float(child_ub[j]), f"ub[{j}]")) + if child_lb[j] > parent_lb[j] + atol: + normal = np.zeros(n) + normal[j] = -1.0 + cuts.append(HalfspaceCut(normal, float(-child_lb[j]), f"lb[{j}]")) + return cuts + + +def _vector(value, name: str) -> np.ndarray: + result = np.asarray(value, dtype=float).reshape(-1) + if result.size == 0 or not np.all(np.isfinite(result)): + raise ValueError(f"{name} must be a nonempty finite vector") + return result + + +def _matlab_quote(path: str | Path) -> str: + return str(path).replace("'", "''") + + +def _string_from_mat(value) -> str: + if isinstance(value, str): + return value + array = np.asarray(value) + if array.dtype.kind in {"U", "S"}: + return "".join(array.reshape(-1).astype(str).tolist()).strip() + return str(value).strip() + + +class StaticVolumeUpdater: + """Sequential wrapper around the unchanged supplied static estimator.""" + + def __init__( + self, + backend: str = "batch", + *, + matlab_executable: str | None = None, + matlab_source_dir: str | Path | None = None, + polytope_sampler_matlab: str | Path | None = None, + ) -> None: + if backend != "batch": + raise ValueError("only backend='batch' is supported") + self.backend = backend + self.matlab_executable = matlab_executable or os.environ.get( + "DINGO_MATLAB_EXECUTABLE", "matlab" + ) + self.matlab_source_dir = Path( + matlab_source_dir + or os.environ.get("DINGO_VOLUME_UPDATING_MATLAB_DIR", "") + or Path(__file__).resolve().parent / "matlab" / "volume_updating" + ).resolve() + sampler = polytope_sampler_matlab or os.environ.get( + "DINGO_POLYTOPE_SAMPLER_MATLAB" + ) + self.polytope_sampler_matlab = Path(sampler).resolve() if sampler else None + self._validate_installation() + + def _validate_installation(self) -> None: + required = [ + "volume_update_batch.m", + "relative_volume.m", + "default_params.m", + "ul_test.m", + "new_window.m", + "update_window.m", + "init.m", + ] + missing = [name for name in required if not (self.matlab_source_dir / name).is_file()] + if missing: + raise FileNotFoundError( + f"incomplete volume-updating MATLAB folder {self.matlab_source_dir}: {missing}" + ) + executable = shutil.which(self.matlab_executable) + if executable is None and not Path(self.matlab_executable).is_file(): + raise FileNotFoundError( + f"MATLAB executable not found: {self.matlab_executable}. Set " + "DINGO_MATLAB_EXECUTABLE if needed." + ) + if self.polytope_sampler_matlab is not None: + missing = [ + name + for name in ("code", "bin") + if not (self.polytope_sampler_matlab / name).is_dir() + ] + if missing: + raise FileNotFoundError( + "PolytopeSamplerMatlab must contain code/ and bin/: " + f"{self.polytope_sampler_matlab} (missing {missing})" + ) + + @property + def _uses_windows_matlab_from_wsl(self) -> bool: + return os.name == "posix" and self.matlab_executable.lower().endswith(".exe") + + def _path_for_matlab(self, path: str | Path) -> str: + path = str(Path(path).resolve()) + if not self._uses_windows_matlab_from_wsl: + return path + result = subprocess.run( + ["wslpath", "-w", path], + check=True, + capture_output=True, + text=True, + ) + return result.stdout.strip() + + @staticmethod + def _normalize_options(options: Mapping | None) -> tuple[dict, float]: + options = dict(options or {}) + algorithm = options.pop("algorithm", "static") + if algorithm not in {"static", "baseline", "relative_volume"}: + raise ValueError( + "algorithm must be 'static'/'baseline'. The wrapper does not " + "provide a separate reuse estimator." + ) + timeout = float(options.pop("wrapper_timeout_sec", 7200)) + unknown = set(options) - _STATIC_OPTION_NAMES + if unknown: + raise ValueError( + "options are not part of the supplied static API: " + + ", ".join(sorted(unknown)) + ) + for name, value in options.items(): + if isinstance(value, (bool, np.bool_)): + raise ValueError(f"{name} must be numeric, not boolean") + array = np.asarray(value) + if array.size != 1: + raise ValueError(f"{name} must be a numeric scalar") + try: + numeric = float(array.reshape(-1)[0]) + except (TypeError, ValueError) as error: + raise ValueError(f"{name} must be a numeric scalar") from error + if not np.isfinite(numeric): + raise ValueError(f"{name} must be finite") + if name in _POSITIVE_INTEGER_OPTIONS: + if numeric <= 0 or not numeric.is_integer(): + raise ValueError(f"{name} must be a positive integer") + options[name] = int(numeric) + else: + options[name] = numeric + if "r" in options and not 0 < options["r"] < 1: + raise ValueError("r must be strictly between 0 and 1") + if "delta" in options and options["delta"] <= 0: + raise ValueError("delta must be positive") + effective_r = options.get("r", _SOURCE_DEFAULTS["r"]) + effective_delta = options.get("delta", _SOURCE_DEFAULTS["delta"]) + if effective_r + effective_delta >= 1: + raise ValueError("r + delta must be less than 1") + if "alpha" in options and not 0 < options["alpha"] < 1: + raise ValueError("alpha must be strictly between 0 and 1") + if "epsilon" in options and options["epsilon"] <= 0: + raise ValueError("epsilon must be positive") + if "t_max_safety" in options and options["t_max_safety"] < 1: + raise ValueError("t_max_safety must be at least 1") + if "verb" in options and options["verb"] not in {0, 1, 2}: + raise ValueError("verb must be 0, 1, or 2") + if timeout <= 0 or not np.isfinite(timeout): + raise ValueError("wrapper_timeout_sec must be positive and finite") + return options, timeout + + @staticmethod + def _normalize_cuts(cuts: Iterable[HalfspaceCut], dimension: int) -> list[HalfspaceCut]: + result = [] + for cut in cuts: + if not isinstance(cut, HalfspaceCut): + raise TypeError("every cut must be a HalfspaceCut") + if cut.normal.size != dimension: + raise ValueError( + f"cut {cut.label!r} has dimension {cut.normal.size}, expected {dimension}" + ) + result.append(cut) + return result + + def _run_matlab( + self, + request_path: Path, + output_path: Path, + timeout: float, + ) -> None: + source = self._path_for_matlab(self.matlab_source_dir) + request = self._path_for_matlab(request_path) + output = self._path_for_matlab(output_path) + sampler_setup = "" + if self.polytope_sampler_matlab is not None: + sampler = self._path_for_matlab(self.polytope_sampler_matlab) + sampler_setup = ( + "setenv('DINGO_POLYTOPE_SAMPLER_MATLAB', " + f"'{_matlab_quote(sampler)}'); " + ) + command = ( + sampler_setup + + f"addpath(genpath('{_matlab_quote(source)}')); " + f"volume_update_batch('{_matlab_quote(request)}', " + f"'{_matlab_quote(output)}');" + ) + env = os.environ.copy() + completed = subprocess.run( + [self.matlab_executable, "-batch", command], + check=False, + capture_output=True, + text=True, + timeout=timeout, + env=env, + ) + if completed.returncode != 0: + transcript = (completed.stdout + "\n" + completed.stderr).strip() + raise RuntimeError( + f"MATLAB volume update failed with code {completed.returncode}:\n" + + transcript[-12000:] + ) + if not output_path.is_file(): + raise RuntimeError("MATLAB exited successfully but produced no output MAT file") + + def estimate( + self, + S, + lb, + ub, + cuts: Iterable[HalfspaceCut], + *, + options: Mapping | None = None, + ) -> StaticVolumeResult: + S = np.asarray(S, dtype=float) + lb = _vector(lb, "lb") + ub = _vector(ub, "ub") + if S.ndim != 2 or S.shape[1] != lb.size or ub.size != lb.size: + raise ValueError("S, lb and ub dimensions are inconsistent") + if not np.all(np.isfinite(S)): + raise ValueError("S must be finite") + if np.any(lb > ub): + raise ValueError("a lower bound exceeds its upper bound") + cuts = self._normalize_cuts(cuts, lb.size) + static_options, timeout = self._normalize_options(options) + if not cuts: + empty_float = np.empty(0, dtype=float) + return StaticVolumeResult( + 1.0, + 0.0, + empty_float, + empty_float, + np.empty(0, dtype=int), + np.empty(0, dtype=int), + np.empty(0, dtype=bool), + (), + ) + + normals = np.vstack([cut.normal for cut in cuts]) + thresholds = np.asarray([cut.threshold for cut in cuts], dtype=float) + with tempfile.TemporaryDirectory(prefix="volume_update_") as temp: + temp_dir = Path(temp) + request_path = temp_dir / "request.mat" + output_path = temp_dir / "result.mat" + savemat( + request_path, + { + "S": S, + "lb": lb.reshape(-1, 1), + "ub": ub.reshape(-1, 1), + "cut_normals": normals, + "cut_thresholds": thresholds.reshape(-1, 1), + "options": static_options, + }, + do_compression=True, + ) + self._run_matlab(request_path, output_path, timeout) + result = loadmat(output_path, squeeze_me=True, struct_as_record=False) + + cut_ratios = np.atleast_1d(result["cut_ratios"]).astype(float) + cut_logs = np.atleast_1d(result["cut_log_ratios"]).astype(float) + cut_phases = np.atleast_1d(result["cut_n_phases"]).astype(int) + cut_steps = np.atleast_1d(result["cut_total_steps"]).astype(int) + cut_converged = np.atleast_1d(result["cut_converged"]).astype(bool) + expected = len(cuts) + arrays = [cut_ratios, cut_logs, cut_phases, cut_steps, cut_converged] + if any(array.size != expected for array in arrays): + raise RuntimeError("MATLAB result arrays do not match the number of cuts") + if ( + not np.all(np.isfinite(cut_ratios)) + or np.any(cut_ratios < 0) + or np.any(cut_ratios > 1 + 1e-6) + ): + raise RuntimeError("MATLAB returned an invalid per-cut ratio") + expected_logs = np.log(np.maximum(cut_ratios, 1e-300)) + if not np.allclose(cut_logs, expected_logs, atol=1e-12, rtol=1e-10): + raise RuntimeError("MATLAB per-cut ratios and log ratios are inconsistent") + if np.any(cut_phases < 0) or np.any(cut_steps < 0): + raise RuntimeError("MATLAB returned a negative phase or step count") + rho = float(np.asarray(result["rho"]).squeeze()) + log_rho = float(np.asarray(result["log_rho"]).squeeze()) + if not np.isfinite(rho) or rho < 0 or rho > 1 + 1e-6: + raise RuntimeError(f"invalid cumulative ratio returned by MATLAB: {rho}") + if not np.isclose(log_rho, np.sum(cut_logs), atol=1e-10, rtol=1e-10): + raise RuntimeError("MATLAB cumulative log ratio is inconsistent with cut ratios") + algorithm = _string_from_mat(result.get("algorithm", "static_relative_volume_sequential")) + if algorithm != "static_relative_volume_sequential": + raise RuntimeError(f"unexpected MATLAB algorithm identifier: {algorithm!r}") + objective_floor = bool(np.asarray(result.get("objective_floor", False)).squeeze()) + if objective_floor: + raise RuntimeError("MATLAB result unexpectedly reports an objective floor") + return StaticVolumeResult( + rho=rho, + log_rho=log_rho, + cut_ratios=cut_ratios, + cut_log_ratios=cut_logs, + cut_n_phases=cut_phases, + cut_total_steps=cut_steps, + cut_converged=cut_converged, + cut_labels=tuple(cut.label for cut in cuts), + algorithm=algorithm, + objective_floor=False, + ) + + +# Backward-compatible class name used by the earlier wrapper code. The alias +# now points to the supplied static implementation; it does not imply a second +# dynamic or reuse algorithm. +MatlabVolumeUpdater = StaticVolumeUpdater + + +__all__ = [ + "StaticVolumeResult", + "StaticVolumeUpdater", + "HalfspaceCut", + "MatlabVolumeUpdater", + "cuts_from_bounds", +] diff --git a/dingo/matlab/volume_updating/LICENSE b/dingo/matlab/volume_updating/LICENSE new file mode 100644 index 0000000..1ef6bf4 --- /dev/null +++ b/dingo/matlab/volume_updating/LICENSE @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2026 GeomScale + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/dingo/matlab/volume_updating/default_params.m b/dingo/matlab/volume_updating/default_params.m new file mode 100644 index 0000000..fcfb456 --- /dev/null +++ b/dingo/matlab/volume_updating/default_params.m @@ -0,0 +1,29 @@ +function params = default_params() + %DEFAULT_PARAMS Return the default parameter struct for relative_volume. + % + % Usage: + % params = default_params(); + % params = Setfield(default_params(), user_overrides); + + % ---- Annealing schedule parameters ---- + params.r = 0.1; % target lower bound for each volume ratio + params.delta = 0.05; % half-width of target ratio interval [r, r+delta] + params.alpha = 0.20; % significance level for U/L t-tests + params.nu = 10; % number of sublists for U/L tests + params.N_utest = 320; % points per sublist (total = nu * N_utest) + params.window_size = 5000; % size of the sliding window (We now set it in the new_window()!!) + + % ---- Convergence parameters ---- + params.epsilon = 0.10; % target relative error for final estimate + params.max_annealing_iter = 200; % max Phase 1 construction iterations + params.max_bisect_iter = 100; % max binary search iterations per body + + % ---- Sampler parameters ---- + params.simdLen = 8; % number of parallel chains per step() + + % ---- t_max estimation ---- + params.t_max_safety = 1.01; % inflate sampled t_max by this factor + + % ---- Verbosity ---- + params.verb = 1; % 0 = quiet, 1 = normal, 2 = detailed +end diff --git a/dingo/matlab/volume_updating/init.m b/dingo/matlab/volume_updating/init.m new file mode 100644 index 0000000..1a59eff --- /dev/null +++ b/dingo/matlab/volume_updating/init.m @@ -0,0 +1,81 @@ +function simdLen = init() + %INIT Add PolytopeSamplerMatlab and src/ to the MATLAB path and compile solvers. + % + % simdLen = init() + % + % Returns the maximum usable simdLen (number of parallel SIMD chains). + % Tries simdLen=8 first; falls back to simdLen=4 if the MEX solver for 8 + % cannot be compiled (e.g. no C++ compiler is available). simdLen=4 has + % pre-compiled binaries in PolytopeSamplerMatlab/bin/ for Linux, macOS, + % and Windows on x86-64. + % + % This function is idempotent — subsequent calls return the cached value + % without recompiling. Called automatically by relative_volume() and + % run_relative_volume(). + + persistent cached_simdLen; + if ~isempty(cached_simdLen) + simdLen = cached_simdLen; + return; + end + + src_dir = fileparts(mfilename('fullpath')); + root = fileparts(src_dir); + + % Add src/ so default_params.m, sliding_window.m, ul_test.m are findable + addpath(src_dir); + + % The original repository keeps PolytopeSamplerMatlab beside src/. The + % Dingo overlay can be installed in a different layout, so first honor an + % explicit environment variable and then check deterministic locations. + sampler_root = getenv('DINGO_POLYTOPE_SAMPLER_MATLAB'); + candidates = { + sampler_root, ... + fullfile(root, 'PolytopeSamplerMatlab'), ... + fullfile(fileparts(fileparts(fileparts(src_dir))), 'PolytopeSamplerMatlab'), ... + fullfile(fileparts(fileparts(fileparts(fileparts(src_dir)))), 'PolytopeSamplerMatlab') ... + }; + sampler_root = ''; + for candidate_idx = 1:numel(candidates) + candidate = candidates{candidate_idx}; + if ~isempty(candidate) && exist(fullfile(candidate, 'code'), 'dir') + sampler_root = candidate; + break; + end + end + if isempty(sampler_root) + error(['PolytopeSamplerMatlab was not found. Set ', ... + 'DINGO_POLYTOPE_SAMPLER_MATLAB to its repository root.']); + end + sampler_code = fullfile(sampler_root, 'code'); + sampler_bin = fullfile(sampler_root, 'bin'); + + addpath(genpath(sampler_code)); + addpath(sampler_bin); + + % Scalar solver — pre-compiled binaries exist for all platforms + compile_solver(0); + + % SIMD solver. Try simdLen=8 first (best throughput); fall back to + % simdLen=4 if the MEX binary cannot be found or compiled. simdLen=4 + % has pre-compiled binaries shipped with PolytopeSamplerMatlab. + for trial = [8, 4] + try + compile_solver(trial); + cached_simdLen = trial; + if trial == 4 + warning('init:fallback', ... + ['Cannot compile MEX solver for simdLen=8 ' ... + '(no C++ compiler?). Falling back to simdLen=4 ' ... + 'which uses pre-compiled binaries.']); + end + simdLen = cached_simdLen; + return; + catch ME + if trial == 4 + rethrow(ME); + end + % simdLen=8 failed — try simdLen=4 + end + end +end diff --git a/dingo/matlab/volume_updating/new_window.m b/dingo/matlab/volume_updating/new_window.m new file mode 100644 index 0000000..afdfe65 --- /dev/null +++ b/dingo/matlab/volume_updating/new_window.m @@ -0,0 +1,82 @@ +function window = new_window(params, dim, eps_i) + %NEW_WINDOW Create a sliding window convergence tracker. + % + % Input: + % dim - effective dimension of the polytope + % eps_i - per-ratio error budget + % + % Output: + % window - struct with fields: + % .W window size = ceil(4*dim^2 + 500) + % .min_val current minimum in the window + % .max_val current maximum in the window + % .min_index index of min in circular buffer + % .max_index index of max in circular buffer + % .last_W circular buffer of length W + % .index write pointer (1-based) + % .converged logical flag + % .eps_i per-ratio error budget + % + % Reference: Cousins & Vempala (2016), Section 3.3. + % Adapted from Volume-and-Sampling/Volume.m lines 246-266. + + %window.W = ceil(4 * dim^2 + 500); + window.W = ceil(params.simdLen * dim^0.52)+400; + %window.W = params.window_size; + window.min_val = -realmax; + window.max_val = realmax; + window.min_index = window.W; + window.max_index = window.W; + window.last_W = zeros(window.W, 1); + window.index = 1; + window.converged = false; + window.eps_i = eps_i; + window.filled = false; +end + +function window = update_window(window, val) + %UPDATE_WINDOW Push a new cumulative running average into the window. + % + % Input: + % window - window struct from new_window + % val - current cumulative running average (scalar) + % + % Output: + % window - updated struct with .converged set to true if converged + % + % Maintains O(1) amortized min/max tracking via lazy eviction: + % when the outgoing value was NOT the min or max, updates are O(1). + % When the outgoing value WAS the min or max, a full scan recomputes. + % + % Adapted from Volume-and-Sampling/Volume.m lines 268-301. + + window.last_W(window.index) = val; + + % Update minimum + if val <= window.min_val + window.min_val = val; + window.min_index = window.index; + elseif window.min_index == window.index + % The outgoing element was the min -- recompute + window.filled = true; + [window.min_val, window.min_index] = min(window.last_W); + end + + % Update maximum + if val >= window.max_val + window.max_val = val; + window.max_index = window.index; + elseif window.max_index == window.index + % The outgoing element was the max -- recompute + window.filled = true; + [window.max_val, window.max_index] = max(window.last_W); + end + + % Check convergence: relative spread within eps_i/2 + if (window.max_val - window.min_val) / window.max_val <= window.eps_i / 2 && window.filled + %if (max(window.last_W) - min(window.last_W)) / max(window.last_W) <= window.eps_i / 2 + window.converged = true; + end + % Advance circular buffer pointer + window.index = mod(window.index, window.W) + 1; +end diff --git a/dingo/matlab/volume_updating/relative_volume.m b/dingo/matlab/volume_updating/relative_volume.m new file mode 100644 index 0000000..9f110fc --- /dev/null +++ b/dingo/matlab/volume_updating/relative_volume.m @@ -0,0 +1,706 @@ +function [rho, result] = relative_volume(problem, c, c0, user_opts) + %RELATIVE_VOLUME Estimate vol(P ∩ H(c0)) / vol(P) via halfspace simulated annealing. + % + % Input: + % problem - struct with fields .Aeq (=S), .beq (=0), .lb, .ub + % (and optionally .Aineq, .bineq) defining the metabolic polytope + % P = {v | Sv=0, lb <= v <= ub}. + % c - cost vector (d x 1 or 1 x d) + % c0 - scalar target threshold for halfspace c'*v <= c0 + % user_opts - optional struct overriding default parameters (see default_params.m) + % + % Output: + % rho - estimate of vol(P ∩ H(c0)) / vol(P) + % result - struct with diagnostic fields: + % .thresholds cell array of thresholds [t_max_or_Inf; t_k; ...; t_1] + % .ratios individual ratio estimates [r_k, ..., r_0] + % .log_ratios log of each ratio estimate + % .steps_per_phase HMC step() counts per ratio + % .total_steps total HMC step() calls + % .n_phases number of intermediate bodies (= length(thresholds)-1) + % .log_rho log of final estimate + % .converged whether sliding window converged for each ratio + % + % Algorithm: Chalkis et al. (2023) annealing schedule with + % Cousins & Vempala (2016) sliding window convergence, + % using Kook et al. (2022) CRHMC via PolytopeSamplerMatlab. + + % ===================================================================== + % 0. Input validation and setup + % ===================================================================== + c = c(:); % ensure column vector + + % Shuffle the global RNG so each run produces different seeds / chains + rng('shuffle'); + + % Must initialize paths BEFORE using Setfield or default_options. + % init() returns the maximum usable simdLen (8 if MEX compiled, 4 if + % only pre-compiled binaries are available). + avail_simdLen = init(); + + if nargin < 4 || isempty(user_opts) + params = default_params(); + else + params = Setfield(default_params(), user_opts); + end + params + + % Use the best simdLen that is actually available + if params.simdLen > avail_simdLen + if params.verb >= 1 + fprintf('Note: simdLen reduced from %d to %d (MEX solver unavailable).\n', ... + params.simdLen, avail_simdLen); + end + params.simdLen = avail_simdLen; + end + + [polytope_opts, sampler_opts] = build_opts(params); + + % ===================================================================== + % Edge case: empty polytope or invalid problem + % ===================================================================== + n_vars = max([size(problem.Aeq,2), length(problem.lb), length(problem.ub)]); + if n_vars == 0 + error('Problem has zero variables.'); + end + + % ===================================================================== + % Phase 1: Build the annealing schedule + % ===================================================================== + if params.verb >= 1 + fprintf('=== Phase 1: Building Annealing Schedule ===\n'); + end + + % Create full polytope P (no extra halfspace) + registry = struct(); + registry.thresholds = {}; + registry.polytopes = {}; + registry.samplers = {}; + registry.initial_ratios = {}; % warm-start values for Phase 2 + registry.dim = 0; + registry.c = c; + + registry = registry_store(registry, Inf, c, problem, polytope_opts, sampler_opts); + sampler_full = registry.samplers{1}; + registry.dim = size(sampler_full.x, 2); + + % ---- Phase 1A+B: Collect first batch, determine t_max, find t_star ---- + % Draw nu*N points from the full polytope P. These serve double duty: + % (a) determine t_max = max_i (c'*x_i) * safety + % (b) reuse as the first U/L-test batch in the binary search for t_star + if params.verb >= 1 + fprintf(' Sampling %d points from P to estimate t_max...\n', ... + params.nu * params.N_utest); + end + + n_needed = params.nu * params.N_utest; + [points, sampler_full] = advance_sampler(sampler_full, n_needed); + c_vals = points * c; + t_max = max(c_vals) * params.t_max_safety; + + if params.verb >= 1 + fprintf(' t_max = %.4f (safety factor %.2f)\n', t_max, params.t_max_safety); + fprintf(' target c0 = %.4f\n', c0); + end + + % Edge case: halfspace contains all of P + if c0 >= t_max + if params.verb >= 1 + fprintf(' c0 >= t_max: halfspace contains entire polytope. rho = 1.\n'); + end + rho = 1.0; + result = build_result(registry, [], [], [], 0, 0, true(0,1)); + return; + end + + % ---- Check: is the ratio to P0 already large enough? ---- + % Run the L-test on the initial points: H0: theta <= r vs H1: theta > r. + % If it passes, vol(P0)/vol(P) >= r with confidence 1-alpha, i.e. the + % ratio is already large enough for direct rejection sampling. Skip the + % binary search and all intermediate body construction. + [~, L_pass_direct, theta_bar_direct, ~] = ul_test(points, c0, c, params); + + if params.verb >= 1 + fprintf(' L-test to P0 on initial points: theta_bar=%.4f, L_pass=%d\n', ... + theta_bar_direct, L_pass_direct); + end + + if L_pass_direct + if params.verb >= 1 + fprintf(' Ratio is already >= r — direct estimation (1 phase).\n'); + end + % sampler_full is the only sampler we need; Phase 2 estimates + % vol(P0) / vol(P) directly via rejection sampling. + direct_estimation = true; + % Store the L-test theta_bar as a warm-start for Phase 2 + registry.initial_ratios{1} = theta_bar_direct; + else + direct_estimation = false; + + % ---- Phase 1B: Binary search for t_star ---- + if params.verb >= 1 + fprintf(' Binary search for t* in [%.4f, %.4f]...\n', c0, t_max); + end + + [t_star, ~, r_init_star] = binary_search_t(sampler_full, c, c0, t_max, params, points); + r_init_star + if params.verb >= 1 + fprintf(' t* = %.4f\n', t_star); + end + + % Create body at t_star + if t_star < t_max - 1e-8 * abs(t_max) + registry = registry_store(registry, t_star, c, problem, ... + polytope_opts, sampler_opts, sampler_full); + % Store warm-start ratio for body P → t_star + r_init_star + if ~isempty(r_init_star) + registry.initial_ratios{1} = r_init_star; + end + sampler_curr = registry.samplers{end}; + t_curr = t_star; + if params.verb >= 1 + fprintf(' Created Polytope + Sampler for P ∩ H(t*).\n'); + end + else + % t_star ≈ t_max: use the full P sampler, skip intermediate bodies + sampler_curr = sampler_full; + t_curr = t_max; + % Store the r_init for the full-P → P0 ratio (if available) + if ~isempty(r_init_star) + registry.initial_ratios{1} = r_init_star; + end + if params.verb >= 1 + fprintf(' t* ≈ t_max, using full P sampler.\n'); + end + end + end + + % ---- Phase 1C-D: Iteratively build intermediate bodies ---- + % Only run when we have an intermediate body and the L-test to P0 + % did not already pass on the initial batch. + if ~direct_estimation && exist('sampler_curr', 'var') + if params.verb >= 1 + fprintf(' Iterative construction (interval [c0=%.4f, t_curr=%.4f])...\n', ... + c0, t_curr); + end + + for iter = 1:params.max_annealing_iter + % Step D: check if we can reach P_0 directly + [pts_check, sampler_curr] = advance_sampler(sampler_curr, n_needed); + [~, L_pass, theta_bar, ~] = ul_test(pts_check, c0, c, params); + + if params.verb >= 2 + fprintf(' Iter %d: L-test to P0: theta_bar=%.4f, L_pass=%d\n', ... + iter, theta_bar, L_pass); + end + + if L_pass + if params.verb >= 1 + fprintf(' L-test to target body passed. Schedule complete.\n'); + end + registry.initial_ratios{end+1} = theta_bar; + break; + end + + % Binary search for next threshold in [c0, t_curr] + [t_next, last_body, r_init_next] = binary_search_t(sampler_curr, c, c0, t_curr, params, pts_check); + + if params.verb >= 1 + fprintf(' Iter %d: t_next = %.4f\n', iter, t_next); + end + + if last_body + if params.verb >= 1 + fprintf(' Last body reached; stopping phase construction.\n'); + end + if ~isempty(r_init_next) + % registry.initial_ratios{k} is the ratio from body k to body k+1 + r_init_next + registry.initial_ratios{end+1} = r_init_next; + end + break; + end + + % Stop if t_next converged to an interval boundary + %if t_next >= t_curr - 1e-8 * abs(t_curr) || t_next <= c0 + 1e-8 + % if params.verb >= 1 + % fprintf(' t_next converged to boundary; stopping iteration.\n'); + % end + % break; + %end + + % Create body at t_next + registry = registry_store(registry, t_next, c, problem, ... + polytope_opts, sampler_opts, sampler_curr); + % Store warm-start ratio for this body pair + if ~isempty(r_init_next) + % registry.initial_ratios{k} is the ratio from body k to body k+1 + r_init_next + registry.initial_ratios{end+1} = r_init_next; + end + sampler_curr = registry.samplers{end}; + t_curr = t_next; + end + + if iter >= params.max_annealing_iter && ~L_pass + warning('Annealing reached max iterations (%d) without L-test passing.', ... + params.max_annealing_iter); + end + elseif ~direct_estimation + % t_star ≈ t_max and no sampler_curr: this is a degenerate case. + % Fall back to direct estimation from full P. + if params.verb >= 1 + fprintf(' No intermediate bodies; estimating directly from P.\n'); + end + end + + n_bodies = length(registry.samplers); + n_ratios = n_bodies; % one ratio per sampler (outer body) + + if params.verb >= 1 + fprintf(' Annealing complete: %d samplers, %d ratios to estimate.\n', ... + n_bodies, n_ratios); + fprintf(' Thresholds (descending): '); + for i = 1:length(registry.thresholds) + fprintf('%.4f ', registry.thresholds{i}); + end + fprintf('\n'); + end + + % ===================================================================== + % Phase 2: Estimate each volume ratio with sliding window + % ===================================================================== + if params.verb >= 1 + fprintf('\n=== Phase 2: Estimating Volume Ratios ===\n'); + end + + eps_i = params.epsilon / sqrt(n_ratios + 1); + log_ratios = zeros(n_ratios, 1); + ratios = zeros(n_ratios, 1); + steps_per_phase = zeros(n_ratios, 1); + converged_flags = false(n_ratios, 1); + + for pair_idx = 1:n_ratios + % Outer body sampler + sampler = registry.samplers{pair_idx}; + + % Inner threshold + if pair_idx < n_ratios + t_inner = registry.thresholds{pair_idx + 1}; + else + t_inner = c0; % target body P_0 + end + + if params.verb >= 1 + fprintf(' Ratio %d/%d: outer threshold=%s, inner threshold=%.4f\n', ... + pair_idx, n_ratios, ... + format_threshold(registry.thresholds{pair_idx}), t_inner); + end + + % Retrieve warm-start ratio from Phase 1 (if any) + %registry.initial_ratios + if pair_idx <= length(registry.initial_ratios) + r_init = registry.initial_ratios{pair_idx}; + else + r_init = []; + end + + % Estimate ratio via rejection sampling with sliding window + r_init + [r_hat, converged, n_steps] = estimate_ratio_sliding(... + sampler, c, t_inner, eps_i, params, r_init); + + ratios(pair_idx) = r_hat; + log_ratios(pair_idx) = log(max(r_hat, 1e-300)); + steps_per_phase(pair_idx) = n_steps; + converged_flags(pair_idx) = converged; + + if params.verb >= 1 + fprintf(' r_hat = %.4f, converged = %d, steps = %d\n', ... + r_hat, converged, n_steps); + end + + % Early exit if ratio is numerically zero + if r_hat < 1e-300 + if params.verb >= 1 + fprintf(' Ratio numerically zero; stopping Phase 2.\n'); + end + % Zero out remaining ratios + log_ratios(pair_idx+1:end) = 0; + ratios(pair_idx+1:end) = 0; + converged_flags(pair_idx+1:end) = false; + break; + end + end + + % ===================================================================== + % Phase 3: Assemble result + % ===================================================================== + log_rho = sum(log_ratios); + + if log_rho < log(realmin) + 10 + warning('log(rho) = %.2f is near underflow; returning rho = 0.', log_rho); + rho = 0.0; + else + rho = exp(log_rho); + end + + result = build_result(registry, ratios, log_ratios, steps_per_phase, ... + n_ratios, log_rho, converged_flags); + + if params.verb >= 1 + fprintf('\n=== Complete ===\n'); + fprintf(' rho_hat = %.16f\n', rho); + fprintf(' log_rho = %.16f\n', log_rho); + fprintf(' n_phases = %d\n', n_ratios); + fprintf(' total HMC steps = %d\n', sum(steps_per_phase)); + end +end + +% ========================================================================= +% Subfunctions +% ========================================================================= + +function [polytope_opts, sampler_opts] = build_opts(params) + %BUILD_OPTS Construct Polytope and Sampler options structs. + % + % Both Polytope() and Sampler() constructors accept the full options + % struct and silently ignore fields they don't consume. We build one + % base struct from default_options() and fork it for each constructor. + + base_opts = default_options(); + + % ---- Polytope options ---- + polytope_opts = base_opts; + polytope_opts.presolve.runSimplify = true; + polytope_opts.presolve.logFunc = @(tag, msg) 0; % suppress constructor output + + % ---- Sampler options ---- + sampler_opts = base_opts; + sampler_opts.maxTime = 86400 * 7; % 1 week (prevent timeout) + sampler_opts.maxStep = Inf; % never stop on step count + sampler_opts.N = Inf; % never stop on ESS + sampler_opts.simdLen = params.simdLen; + sampler_opts.initalStepSize = 0.2; + sampler_opts.freezeMCMCAfterSamples = Inf; + + % Memory limit for chain storage (Module: MemoryStorage) + sampler_opts.MemoryStorage.memoryLimit = 32 * 1024^3; % 32 GB + + % Remove ProgressBar (cosmetic) but keep MixingTimeEstimator and MemoryStorage + sampler_opts.module = setdiff(sampler_opts.module, {'ProgressBar'}, 'stable'); + + % Suppress all sampler logging. We control convergence ourselves via + % the sliding window — the built-in ESS/step logging is not needed. + % Note: this also suppresses diagnostic messages from DynamicStepSize + % (e.g. step-shrink events). Set params.verb >= 3 in the future to + % route those through a function-handle logger. + sampler_opts.logging = []; +end + +function registry = registry_store(registry, t, c, problem, polytope_opts, sampler_opts, ... + prev_sampler) + %REGISTRY_STORE Create Polytope + Sampler for P ∩ H(t) and store in registry. + % + % If t == Inf, creates the full polytope P without adding a halfspace. + % Polytope and Sampler are handle objects — mutations persist through the registry. + % + % Optional prev_sampler — if provided, its adapted stepSize is used as the + % initial step size for the new sampler (warm-start across nested bodies). + + if isinf(t) && t > 0 + % Full polytope P (no halfspace constraint) + problem_new = problem; + else + % Add halfspace c'*x <= t + problem_new = problem; + c_row = c(:)'; % ensure row vector for Aineq + if ~isfield(problem_new, 'Aineq') || isempty(problem_new.Aineq) + problem_new.Aineq = c_row; + problem_new.bineq = t; + else + problem_new.Aineq = [problem_new.Aineq; c_row]; + problem_new.bineq = [problem_new.bineq(:); t]; + end + end + + % Set a random seed for each new body (avoid identical chains across bodies) + % Use uint32 to match the SimdTwister RNG expected by the Sampler constructor + sampler_opts.seed = randi(intmax('uint32'), 'uint32'); + + % Carry forward adapted step size from previous (outer) body. + % Must be set before Sampler() reads opts.initalStepSize. + if nargin >= 8 && ~isempty(prev_sampler) && isvalid(prev_sampler) + sampler_opts.initalStepSize = prev_sampler.stepSize; + end + + % PolytopeSamplerMatlab's top-level sample() initializes a shared timeout + % clock before constructing Polytope and Sampler. This routine constructs + % those objects directly, so reproduce that required interface setup here. + start_time = tic; + polytope_opts.startTime = start_time; + sampler_opts.startTime = start_time; + + polytope = Polytope(problem_new, polytope_opts); + % sample() also initializes the SIMD random stream before Sampler(). + rng(sampler_opts.seed, 'simdTwister'); + sampler = Sampler(polytope, sampler_opts); + + % ---- Burn-in: settle step size and approach stationarity, then freeze ---- + % Kook et al. (2022) report empirically O(d^0.52) steps per effective + % sample for CRHMC on metabolic polytopes. 15 effective samples gives + % the chain time to mix and DynamicStepSize time to stabilise. + tic + d = size(sampler.x, 2); + burn_in_steps = ceil(15 * d^0.52); + for k = 1:burn_in_steps + sampler.step(); + if sampler.terminate == 3 + error('Sampler terminated with error during burn-in at step %d.', k); + elseif sampler.terminate ~= 0 + warning('Sampler terminated with code %d during burn-in at step %d.', ... + sampler.terminate, k); + break; + end + end + time_to_burn_in = toc + sampler.freezed = true; + + registry.thresholds{end+1} = t; + registry.polytopes{end+1} = polytope; + registry.samplers{end+1} = sampler; + + if registry.dim == 0 + registry.dim = size(sampler.x, 2); + end +end + +function [points, sampler] = advance_sampler(sampler, n_points) + %ADVANCE_SAMPLER Advance the sampler and collect transformed points. + % + % Input: + % sampler - Sampler handle object + % n_points - number of points to collect + % + % Output: + % points - M x d_orig matrix of points in ORIGINAL space + % where M = ceil(n_points / simdLen) * simdLen + % sampler - the same handle object (mutated in place) + + simdLen = size(sampler.x, 1); + if isempty(simdLen) || simdLen == 0 + simdLen = 1; + end + n_steps = ceil(n_points / simdLen); + n_total = n_steps * simdLen; + + % Pre-allocate in original space + dim_orig = size(sampler.problem.T, 1); + points = zeros(n_total, dim_orig); + + idx = 0; % last index written + + for step = 1:n_steps + sampler.step(); + + % Check for sampler errors + if sampler.terminate == 3 + error('Sampler terminated with error at iteration %d.', sampler.i); + elseif sampler.terminate ~= 0 + warning('Sampler terminated with code %d at iteration %d.', ... + sampler.terminate, sampler.i); + break; + end + + % Transform each chain's point to original space + % x_orig = T * x_internal + y + for chain = 1:simdLen + idx = (step-1)*simdLen + chain; + x_int = sampler.x(chain, :)'; + points(idx, :) = (sampler.problem.T * x_int + sampler.problem.y)'; + end + end + + % Trim if we collected fewer than allocated + if idx < n_total + points = points(1:idx, :); + end +end + +function [t_found, last_body, r_init] = binary_search_t(sampler, c, t_lo, t_hi, params, initial_points) + %BINARY_SEARCH_T Find threshold t such that vol ratio ∈ [r, r+δ]. + % + % Searches in [t_lo, t_hi] using U/L tests on points from the sampler. + % The sampler is already sampling from the OUTER body. + % + % Optional: initial_points (M×d matrix) — if provided, used for the + % first bisection iteration instead of collecting fresh points. + % + % Output: + % t_found - threshold satisfying U-test and L-test + % last_body - true if the L-test to t_lo passes (no more bodies needed) + % r_init - theta_bar at t_found (warm-start for Phase 2) + + last_body = false; + r_init = []; % warm-start ratio for Phase 2 + lo = t_lo; + hi = t_hi; + t_found = (lo + hi) / 2; % fallback + have_initial = (nargin >= 6 && ~isempty(initial_points)); + % Use pre-collected points for first iteration, otherwise collect fresh + if have_initial + pts = initial_points; + else + n_needed = params.nu * params.N_utest; + [pts, ~] = advance_sampler(sampler, n_needed); + end + %[~, L_pass, theta_bar_lo, ~] = ul_test(pts, t_lo, c, params); + %if L_pass + % t_found = t_lo; + % last_body = true; + % r_init = theta_bar_lo; + % return; + %end + theta_bar = 0; + for i = 1:params.max_bisect_iter + if i==params.max_bisect_iter + fprintf(' hello world max_bisect_iter reached\n'); + end + if hi - lo < 1e-8 * max(abs(hi), 1) + t_found = (lo + hi) / 2; + fprintf(' hello world t_found = %.4f\n', t_found); + r_init = theta_bar; + break; + end + + t_mid = (lo + hi) / 2; + + [U_pass, L_pass, theta_bar, s_theta] = ul_test(pts, t_mid, c, params); + + if params.verb >= 2 + fprintf(' bisect %d: t=%.4f, theta_bar=%.4f, s=%.4f, U=%d, L=%d\n', ... + i, t_mid, theta_bar, s_theta, U_pass, L_pass); + end + + if U_pass && L_pass + t_found = t_mid; + r_init = theta_bar; + return; + elseif ~U_pass + % Ratio >= r+delta → inner body too large → tighten (decrease t) + hi = t_mid; + else % ~L_pass + % Ratio <= r → inner body too small → loosen (increase t) + lo = t_mid; + end + end + + % If we exit without exact convergence, use midpoint + if ~exist('U_pass', 'var') || ~(U_pass && L_pass) + t_found = (lo + hi) / 2; + end +end + +function [r_hat, converged, n_steps] = estimate_ratio_sliding(sampler, c, t_inner, ... + eps_i, params, r_init) + %ESTIMATE_RATIO_SLIDING Estimate volume ratio via rejection sampling with sliding window. + % + % Samples from the outer body (sampler) and tests membership in the + % inner body (c'*x <= t_inner). Convergence via sliding window. + % + % Optional r_init (scalar in (0,1)) warm-starts the counters and + % pre-fills the sliding window from a Phase 1 U/L-test ratio estimate. + + % Use this sampler's internal dimension for window sizing + dim = size(sampler.x, 2); + window = new_window(params, dim, eps_i); + simdLen = size(sampler.x, 1); + if isempty(simdLen) || simdLen == 0 + simdLen = 1; + end + c_vec = c(:); + + % ---- Warm-start from Phase 1 ratio estimate ---- + if nargin >= 7 && ~isempty(r_init) && r_init > 0 && r_init < 1 + % Use the same sample size as Phase 1 for consistent weighting + n_warm = params.nu * params.N_utest; + count_in = round(r_init * n_warm); + total = n_warm; + % Pre-fill the sliding window with the initial estimate + %window.last_W(:) = r_init; + %window.min_val = r_init; + %window.max_val = r_init; + %window.min_index = 1; + %window.max_index = 1; + %window.index = 1; + if params.verb >= 2 + fprintf(' (warm-start: r_init=%.4f, n_warm=%d)\n', r_init, n_warm); + end + else + count_in = 0; + total = 0; + end + n_steps = 0; + converged = false; + + while ~converged + sampler.step(); + n_steps = n_steps + simdLen; + + % Check for errors + if sampler.terminate == 3 + error('Sampler terminated with error at iteration %d.', sampler.i); + elseif sampler.terminate ~= 0 + warning('Sampler terminated with code %d at iteration %d.', ... + sampler.terminate, sampler.i); + break; + end + + for chain = 1:simdLen + % Transform from internal to original space + x_int = sampler.x(chain, :)'; + x_orig = sampler.problem.T * x_int + sampler.problem.y; + + % Membership test + in_smaller = (c_vec' * x_orig <= t_inner); + + count_in = count_in + in_smaller; + total = total + 1; + + % Update running average and sliding window + val = count_in / total; + window = update_window(window, val); + + if window.converged + converged = true; + break; + end + end + end + + r_hat = count_in / total; +end + +function result = build_result(registry, ratios, log_ratios, steps_per_phase, ... + n_ratios, log_rho, converged) + %BUILD_RESULT Assemble the output result struct. + result.thresholds = registry.thresholds; + result.ratios = ratios; + result.log_ratios = log_ratios; + result.steps_per_phase = steps_per_phase; + result.total_steps = sum(steps_per_phase); + result.n_phases = n_ratios; + result.log_rho = log_rho; + result.converged = converged; + result.samplers = registry.samplers; +end + +function s = format_threshold(t) + %FORMAT_THRESHOLD String representation of a threshold value. + if isinf(t) && t > 0 + s = '+Inf'; + else + s = sprintf('%.4f', t); + end +end diff --git a/dingo/matlab/volume_updating/run_relative_volume.m b/dingo/matlab/volume_updating/run_relative_volume.m new file mode 100644 index 0000000..a5a49b0 --- /dev/null +++ b/dingo/matlab/volume_updating/run_relative_volume.m @@ -0,0 +1,247 @@ +function result = run_relative_volume(model_path, lambda, user_opts) + %RUN_RELATIVE_VOLUME Load a metabolic model, generate a random halfspace, + % and compute the relative volume vol(P ∩ H(c0)) / vol(P). + % + % Input: + % model_path - path to a .mat file containing a 'problem' struct + % with fields .Aeq (=S), .beq, .lb, .ub. + % Relative to repo root, e.g. 'data/e_coli_core' + % (the .mat extension is optional). + % lambda - scalar in (0,1) determining the halfspace threshold: + % c0 = (1-lambda)*c_max + lambda*c_min + % lambda=0 → c0 = c_max (halfspace contains all of P) + % lambda=1 → c0 = c_min (halfspace barely touches P) + % user_opts - optional struct overriding default parameters for + % relative_volume (see default_params.m). + % + % Output: + % result - struct containing all fields from relative_volume plus: + % .rho the relative volume estimate + % .c random cost vector used (d x 1, unit norm) + % .c_min min c'*v over P (from LP) + % .c_max max c'*v over P (from LP) + % .c0 halfspace threshold c'*v <= c0 + % .lambda the lambda value used + + % ----------------------------------------------------------------- + % 0. Path setup & load the model + % ----------------------------------------------------------------- + avail_simdLen = init(); %#ok % sets up paths, compiles solvers + + root = fileparts(fileparts(mfilename('fullpath'))); + + % Accept path with or without .mat extension + if ~endsWith(model_path, '.mat') + full_path = fullfile(root, [model_path '.mat']); + else + full_path = fullfile(root, model_path); + end + + if ~exist(full_path, 'file') + error('Model file not found: %s', full_path); + end + + data = load(full_path); + + if isfield(data, 'problem') + % --- Existing format: file contains a 'problem' struct --- + problem = data.problem; + else + % --- Try to auto-detect a COBRA-style metabolic model --- + model = detect_metabolic_model(data); + if isempty(model) + error('File does not contain a ''problem'' struct or a recognized metabolic model: %s', full_path); + end + fprintf(' Detected metabolic model: %d rxns, %d mets\n', ... + length(model.rxns), length(model.mets)); + problem = metabolic_to_problem(model); + end + problem = standardize_problem(problem); + + % ----------------------------------------------------------------- + % Detect and handle netlib-style (unbounded) problems. + % Netlib .mat files ship with a non-empty df vector (the LP cost). + % Metabolic models have df = [] (uniform). Netlib polytopes are + % typically unbounded, so we apply the same bounding logic used by + % PolytopeSampler's loadProblem.m. + % ----------------------------------------------------------------- + is_netlib = ~isempty(problem.df) && isnumeric(problem.df); + if is_netlib + fprintf(' Detected netlib-style problem (non-empty df). Applying bounds...\n'); + + Aeq_net = problem.Aeq; + beq_net = problem.beq(:); + lb_net = problem.lb(:); + ub_net = problem.ub(:); + if isfield(problem, 'Aineq') && ~isempty(problem.Aineq) + Aineq_net = problem.Aineq; + bineq_net = problem.bineq(:); + else + Aineq_net = []; + bineq_net = []; + end + + lp_opts = optimoptions('linprog', 'Display', 'none'); + + % Solve the netlib LP: min df'*x over the original constraints + df_net = problem.df(:); + [x_opt, ~, exitflag] = linprog(df_net, Aineq_net, bineq_net, ... + Aeq_net, beq_net, lb_net, ub_net, lp_opts); + + if exitflag <= 0 + warning('Netlib bounding LP did not converge (exitflag=%d). Using problem as-is.', exitflag); + else + % Add a bounding inequality: df'*x <= df'*x_opt + |df|'*|x_opt| + threshold = df_net' * x_opt + abs(df_net)' * abs(x_opt); + problem.Aineq = [Aineq_net; df_net']; + problem.bineq = [bineq_net; threshold]; + + % Clamp variable bounds around the optimal vertex + bound = 2 * max(abs(x_opt)); + problem.lb = max(lb_net, -bound); + problem.ub = min(ub_net, bound); + + % Switch to uniform sampling + problem.f = []; + problem.df = []; + problem.ddf = []; + + fprintf(' Bounding LP solved. Added df''*x <= %.4f, bounds clamped to ±%.4f.\n', ... + threshold, bound); + end + end + + % Effective dimension (before presolve — number of original variables) + d = max([size(problem.Aeq, 2), length(problem.lb), length(problem.ub)]); + + % ----------------------------------------------------------------- + % 1. Generate random cost vector + % ----------------------------------------------------------------- + c = randn(d, 1); + c = c / norm(c); + + % ----------------------------------------------------------------- + % 2. Compute c_min and c_max via LP + % ----------------------------------------------------------------- + Aeq = problem.Aeq; + beq = problem.beq(:); + lb = problem.lb(:); + ub = problem.ub(:); + + if isfield(problem, 'Aineq') && ~isempty(problem.Aineq) + Aineq = problem.Aineq; + bineq = problem.bineq(:); + else + Aineq = []; + bineq = []; + end + + lp_opts = optimoptions('linprog', 'Display', 'none'); + + % min c'*x over P + [~, c_min] = linprog(c, Aineq, bineq, Aeq, beq, lb, ub, lp_opts); + + % max c'*x = -min -c'*x over P + [~, c_max_neg] = linprog(-c, Aineq, bineq, Aeq, beq, lb, ub, lp_opts); + c_max = -c_max_neg; + + % ----------------------------------------------------------------- + % 3. Set threshold: interpolate between c_min and c_max + % ----------------------------------------------------------------- + c0 = (1 - lambda) * c_max + lambda * c_min; + + fprintf('=== run_relative_volume ===\n'); + fprintf(' model : %s\n', model_path); + fprintf(' dimension : %d\n', d); + fprintf(' lambda : %.4f\n', lambda); + fprintf(' c_min : %.4f\n', c_min); + fprintf(' c_max : %.4f\n', c_max); + fprintf(' c0 : %.4f\n', c0); + fprintf(' (lambda=0 => c0=c_max=%.4f, lambda=1 => c0=c_min=%.4f)\n', c_max, c_min); + + % ----------------------------------------------------------------- + % 4. Call relative_volume + % ----------------------------------------------------------------- + if nargin < 3 || isempty(user_opts) + [rho, vol_result] = relative_volume(problem, c, c0); + else + [rho, vol_result] = relative_volume(problem, c, c0, user_opts); + end + + % ----------------------------------------------------------------- + % 5. Assemble output + % ----------------------------------------------------------------- + result = vol_result; + result.rho = rho; + result.c = c; + result.c_min = c_min; + result.c_max = c_max; + result.c0 = c0; + result.lambda = lambda; +end + +% ========================================================================= +% Helper: detect a COBRA-style metabolic model in a loaded .mat struct +% ========================================================================= +function model = detect_metabolic_model(data) + %DETECT_METABOLIC_MODEL Find the first COBRA-style model in loaded data. + % + % A metabolic model is identified by having fields: S, lb, ub, and rxns. + + fields = fieldnames(data); + for i = 1:numel(fields) + var = data.(fields{i}); + if isstruct(var) && all(isfield(var, {'S', 'lb', 'ub', 'rxns'})) + % Found a COBRA model — store its variable name for reporting + model = var; + model.model_var_name = fields{i}; %#ok + return; + end + end + + % No metabolic model found + model = []; +end + +% ========================================================================= +% Helper: convert a COBRA metabolic model to the PolytopeSampler problem +% struct (Aeq·x = beq, lb ≤ x ≤ ub). +% ========================================================================= +function problem = metabolic_to_problem(model) + %METABOLIC_TO_PROBLEM Convert COBRA model to PolytopeSampler format. + % + % COBRA model fields used: + % .S - stoichiometric matrix (mets × rxns), S·v = b + % .b - right-hand side (default: zeros) + % .lb - lower bounds on reaction fluxes + % .ub - upper bounds on reaction fluxes + % + % Output problem fields: + % .Aeq = S (after removing all-zero metabolite rows) + % .beq = b (after removing the same rows) + % .lb = lb(:) + % .ub = ub(:) + % .f, .df, .ddf = [] (uniform sampling) + + % --- Stoichiometric matrix & right-hand side --- + S = model.S; + if isfield(model, 'b') && ~isempty(model.b) + b_vec = model.b(:); + else + b_vec = zeros(size(S, 1), 1); + end + + % Remove metabolite rows that impose no constraint (all-zero rows) + nonzero_rows = any(S, 2); + S = S(nonzero_rows, :); + b_vec = b_vec(nonzero_rows); + + % --- Assemble problem struct --- + problem.Aeq = S; + problem.beq = b_vec; + problem.lb = model.lb(:); + problem.ub = model.ub(:); + problem.f = []; + problem.df = []; + problem.ddf = []; +end diff --git a/dingo/matlab/volume_updating/sliding_window.m b/dingo/matlab/volume_updating/sliding_window.m new file mode 100644 index 0000000..d658eb0 --- /dev/null +++ b/dingo/matlab/volume_updating/sliding_window.m @@ -0,0 +1,77 @@ +function window = new_window(dim, eps_i) + %NEW_WINDOW Create a sliding window convergence tracker. + % + % Input: + % dim - effective dimension of the polytope + % eps_i - per-ratio error budget + % + % Output: + % window - struct with fields: + % .W window size = ceil(4*dim^2 + 500) + % .min_val current minimum in the window + % .max_val current maximum in the window + % .min_index index of min in circular buffer + % .max_index index of max in circular buffer + % .last_W circular buffer of length W + % .index write pointer (1-based) + % .converged logical flag + % .eps_i per-ratio error budget + % + % Reference: Cousins & Vempala (2016), Section 3.3. + % Adapted from Volume-and-Sampling/Volume.m lines 246-266. + + window.W = ceil(4 * dim^2 + 500); + window.min_val = -realmax; + window.max_val = realmax; + window.min_index = window.W; + window.max_index = window.W; + window.last_W = zeros(window.W, 1); + window.index = 1; + window.converged = false; + window.eps_i = eps_i; +end + +function window = update_window(window, val) + %UPDATE_WINDOW Push a new cumulative running average into the window. + % + % Input: + % window - window struct from new_window + % val - current cumulative running average (scalar) + % + % Output: + % window - updated struct with .converged set to true if converged + % + % Maintains O(1) amortized min/max tracking via lazy eviction: + % when the outgoing value was NOT the min or max, updates are O(1). + % When the outgoing value WAS the min or max, a full scan recomputes. + % + % Adapted from Volume-and-Sampling/Volume.m lines 268-301. + + window.last_W(window.index) = val; + + % Update minimum + if val <= window.min_val + window.min_val = val; + window.min_index = window.index; + elseif window.min_index == window.index + % The outgoing element was the min -- recompute + [window.min_val, window.min_index] = min(window.last_W); + end + + % Update maximum + if val >= window.max_val + window.max_val = val; + window.max_index = window.index; + elseif window.max_index == window.index + % The outgoing element was the max -- recompute + [window.max_val, window.max_index] = max(window.last_W); + end + + % Check convergence: relative spread within eps_i/2 + if (window.max_val - window.min_val) / window.max_val <= window.eps_i / 2 + window.converged = true; + end + + % Advance circular buffer pointer + window.index = mod(window.index, window.W) + 1; +end diff --git a/dingo/matlab/volume_updating/ul_test.m b/dingo/matlab/volume_updating/ul_test.m new file mode 100644 index 0000000..ef4232e --- /dev/null +++ b/dingo/matlab/volume_updating/ul_test.m @@ -0,0 +1,63 @@ +function [U_pass, L_pass, theta_bar, s_theta] = ul_test(points, t_inner, c, params) + %UL_TEST U-test and L-test for volume ratio interval verification. + % + % Input: + % points - M x d matrix of points in ORIGINAL space from the OUTER body + % t_inner - threshold defining the INNER body (c'*x <= t_inner) + % c - cost vector (d x 1) + % params - struct with fields: .r, .delta, .alpha, .nu, .N_utest + % + % Output: + % U_pass - true if U-test rejects H0: theta >= r+delta + % L_pass - true if L-test rejects H0: theta <= r + % theta_bar - sample mean of per-sublist empirical ratios + % s_theta - sample standard deviation of per-sublist ratios + % + % Reference: Chalkis et al. (2023), Section 4.1, Eqs. 9-10. + + nu = params.nu; + N = params.N_utest; + r = params.r; + delta = params.delta; + alpha = params.alpha; + + n_available = size(points, 1); + n_needed = nu * N; + + if n_available < n_needed + warning('ul_test: only %d points available, need %d. Returning not-passed.', ... + n_available, n_needed); + U_pass = false; + L_pass = false; + theta_bar = NaN; + s_theta = NaN; + return; + end + + % Compute c'*x for the first nu*N points + c_vals = points(1:n_needed, :) * c(:); + + % Partition into nu sublists and compute per-sublist ratios + theta_hat = zeros(nu, 1); + for j = 1:nu + idx_start = (j-1)*N + 1; + idx_end = j*N; + n_in = sum(c_vals(idx_start:idx_end) <= t_inner); + theta_hat(j) = n_in / N; + end + + % Sample statistics + theta_bar = mean(theta_hat); + s_theta = std(theta_hat); + + % Critical value: upper alpha-quantile of t-distribution with nu-1 df + t_crit = tinv(1 - alpha, nu - 1); + + se = s_theta / sqrt(nu); % standard error + + % U-test: reject H0: theta >= r+delta if theta_bar <= r+delta - t_crit*se + U_pass = (theta_bar <= r + delta - t_crit * se); + + % L-test: reject H0: theta <= r if theta_bar >= r + t_crit*se + L_pass = (theta_bar >= r + t_crit * se); +end diff --git a/dingo/matlab/volume_updating/update_window.m b/dingo/matlab/volume_updating/update_window.m new file mode 100644 index 0000000..16d0ab2 --- /dev/null +++ b/dingo/matlab/volume_updating/update_window.m @@ -0,0 +1,44 @@ +function window = update_window(window, val) + %UPDATE_WINDOW Push a new cumulative running average into the window. + % + % Input: + % window - window struct from new_window + % val - current cumulative running average (scalar) + % + % Output: + % window - updated struct with .converged set to true if converged + % + % Maintains O(1) amortized min/max tracking via lazy eviction: + % when the outgoing value was NOT the min or max, updates are O(1). + % When the outgoing value WAS the min or max, a full scan recomputes. + % + % Adapted from Volume-and-Sampling/Volume.m lines 268-301. + + window.last_W(window.index) = val; + + % Update minimum + if val <= window.min_val + window.min_val = val; + window.min_index = window.index; + elseif window.min_index == window.index + % The outgoing element was the min -- recompute + [window.min_val, window.min_index] = min(window.last_W); + end + + % Update maximum + if val >= window.max_val + window.max_val = val; + window.max_index = window.index; + elseif window.max_index == window.index + % The outgoing element was the max -- recompute + [window.max_val, window.max_index] = max(window.last_W); + end + + % Check convergence: relative spread within eps_i/2 + if (window.max_val - window.min_val) / window.max_val <= window.eps_i / 2 + window.converged = true; + end + + % Advance circular buffer pointer + window.index = mod(window.index, window.W) + 1; +end diff --git a/dingo/matlab/volume_updating/volume_update_batch.m b/dingo/matlab/volume_updating/volume_update_batch.m new file mode 100644 index 0000000..3aac124 --- /dev/null +++ b/dingo/matlab/volume_updating/volume_update_batch.m @@ -0,0 +1,102 @@ +function volume_update_batch(request_path, output_path) +%VOLUME_UPDATE_BATCH Apply sequential cuts using relative_volume unchanged. +% +% The request MAT file must contain S, lb, ub, cut_normals and +% cut_thresholds. Each new cut becomes part of the parent polytope before the +% next call. No objective, biomass floor, FBA or alternative volume logic is +% introduced here. + + req = load(request_path); + required = {'S', 'lb', 'ub', 'cut_normals', 'cut_thresholds'}; + for k = 1:numel(required) + if ~isfield(req, required{k}) + error('Missing request field: %s', required{k}); + end + end + + S = double(req.S); + lb = double(req.lb(:)); + ub = double(req.ub(:)); + cut_normals = double(req.cut_normals); + cut_thresholds = double(req.cut_thresholds(:)); + n = numel(lb); + + if size(S, 2) ~= n || numel(ub) ~= n + error('Dimension mismatch among S, lb and ub.'); + end + if size(cut_normals, 2) ~= n + error('cut_normals must have one column per reaction.'); + end + if size(cut_normals, 1) ~= numel(cut_thresholds) + error('cut_normals rows and cut_thresholds length differ.'); + end + if any(~isfinite(S(:))) || any(~isfinite(lb)) || any(~isfinite(ub)) || ... + any(~isfinite(cut_normals(:))) || any(~isfinite(cut_thresholds)) + error('All wrapper inputs must be finite.'); + end + if any(lb > ub) + error('At least one lower bound exceeds its upper bound.'); + end + + problem.Aeq = S; + problem.beq = zeros(size(S, 1), 1); + problem.Aineq = zeros(0, n); + problem.bineq = zeros(0, 1); + problem.lb = lb; + problem.ub = ub; + problem.f = []; + problem.df = []; + problem.ddf = []; + + if isfield(req, 'options') + options = req.options; + else + options = struct(); + end + + n_cuts = size(cut_normals, 1); + cut_ratios = ones(n_cuts, 1); + cut_log_ratios = zeros(n_cuts, 1); + cut_n_phases = zeros(n_cuts, 1); + cut_total_steps = zeros(n_cuts, 1); + cut_converged = false(n_cuts, 1); + cut_diagnostics = cell(n_cuts, 1); + + for j = 1:n_cuts + c = cut_normals(j, :)'; + c0 = cut_thresholds(j); + + [rho_j, diag_j] = relative_volume(problem, c, c0, options); + if ~isscalar(rho_j) || ~isfinite(rho_j) || rho_j < 0 || rho_j > 1 + 1e-6 + error('Invalid ratio returned for cut %d: %.17g', j, rho_j); + end + + cut_ratios(j) = rho_j; + cut_log_ratios(j) = log(max(rho_j, 1e-300)); + if isfield(diag_j, 'n_phases'); cut_n_phases(j) = diag_j.n_phases; end + if isfield(diag_j, 'total_steps'); cut_total_steps(j) = diag_j.total_steps; end + if isfield(diag_j, 'converged') + cut_converged(j) = all(diag_j.converged(:)); + end + % MATLAB handle objects cannot be serialized reliably to a portable + % v7 MAT file. They are internal sampler state, not wrapper output. + if isfield(diag_j, 'samplers') + diag_j = rmfield(diag_j, 'samplers'); + end + cut_diagnostics{j} = diag_j; + + % This is the entire sequential-update rule: the accepted cut becomes + % part of the next parent, and the same static Apostolos routine is + % called again for the next cut. + problem.Aineq(end + 1, :) = c'; %#ok + problem.bineq(end + 1, 1) = c0; %#ok + end + + log_rho = sum(cut_log_ratios); + rho = exp(log_rho); + algorithm = 'static_relative_volume_sequential'; + objective_floor = false; + save(output_path, 'rho', 'log_rho', 'cut_ratios', 'cut_log_ratios', ... + 'cut_n_phases', 'cut_total_steps', 'cut_converged', ... + 'cut_diagnostics', 'algorithm', 'objective_floor', '-v7'); +end diff --git a/dingo/volume_updating.py b/dingo/volume_updating.py new file mode 100644 index 0000000..c4c9127 --- /dev/null +++ b/dingo/volume_updating.py @@ -0,0 +1,21 @@ +"""Public imports for sequential relative-volume updating. + +Existing scripts imported ``dingo.volume_updating``. Keep that import path, +but expose only the sequential wrapper around the supplied static routine. +""" + +from dingo.dynamic_volume import ( + StaticVolumeResult, + StaticVolumeUpdater, + HalfspaceCut, + MatlabVolumeUpdater, + cuts_from_bounds, +) + +__all__ = [ + "StaticVolumeResult", + "StaticVolumeUpdater", + "HalfspaceCut", + "MatlabVolumeUpdater", + "cuts_from_bounds", +] diff --git a/doc/dynamic_volume_updating.md b/doc/dynamic_volume_updating.md new file mode 100644 index 0000000..7c9e389 --- /dev/null +++ b/doc/dynamic_volume_updating.md @@ -0,0 +1,213 @@ +# Corrected Summer Project Report + +## Dynamic Volume Updating of Convex Polytopes via Random-Walk Sampling + +## 1. Correct project scope + +The project concerns a sequence of nested polytopes formed by adding linear +constraints. In a metabolic model the initial feasible body is + +`P = {v : S v = 0, lb <= v <= ub}`. + +No biomass objective floor or FBA-optimality restriction is inserted. Boundary +sampling is a separate project and has no role here. + +The numerical foundation is the supplied static `relative_volume.m` method: + +`rho = vol(P intersect {c^T v <= c0}) / vol(P)`. + +That method builds a simulated-annealing schedule of parallel halfspaces, +estimates a telescoping product of volume ratios, and obtains the necessary +samples with PolytopeSamplerMatlab CRHMC. + +This integration does not replace that logic. For a stream of cuts it invokes +the same static routine repeatedly, adding each accepted cut to the parent +before the next call. + +## 2. What was wrong in the previous Dingo wrapper + +The earlier folder contained three incompatible ideas: + +1. `volume_reference.py` manually projected `S v = 0` into a SciPy nullspace + and divided two independent volesti absolute-volume estimates. This bypassed + the supplied method completely. It also did not robustly cover fixed-width reactions, + infinite bounds, additional affine equalities, or dimension-changing cuts. +2. `MatlabVolumeEngine` attempted to import `dingo.volume_updating`, but that + module was absent from the delivered bundle. +3. The engine expected diagnostics such as `cut_reused_pilot` and + `cut_fresh_samples` that the supplied MATLAB code does not return. + +Consequently, the wrapper was neither self-contained nor a faithful wrapper of +the static method. + +## 3. Corrected Dingo architecture + +The corrected integration contains two Python modules: + +- `dingo.dynamic_volume`: the real implementation; +- `dingo.volume_updating`: a backward-compatible import path. + +The public objects are: + +- `HalfspaceCut`; +- `cuts_from_bounds`; +- `StaticVolumeUpdater`; +- backward-compatible alias `MatlabVolumeUpdater`; +- `StaticVolumeResult`. + +The wrapper validates dimensions, nested bounds, options, executable paths and +all MATLAB outputs. It supports Linux MATLAB and Windows MATLAB invoked from +WSL. Temporary request/result files are isolated for each call. + +## 4. Sequential static algorithm + +Let `P^(0)` be the supplied parent and let cut `j` be + +`C_j = {v : c_j^T v <= t_j}`. + +For each `j`: + +1. call `relative_volume(P^(j-1), c_j, t_j)`; +2. obtain `rho_j = vol(P^(j-1) intersect C_j)/vol(P^(j-1))`; +3. set `P^j = P^(j-1) intersect C_j`; +4. continue with the next cut. + +The final estimate is + +`rho = product_j rho_j`, + +computed in log space. This is only the sequential composition of the static +volume estimator. + +## 5. Relationship to the supplied cutting-planes ZIP + +The installed MATLAB directory contains the eight original source files. + +Byte-for-byte unchanged: + +- `default_params.m`; +- `new_window.m`; +- `run_relative_volume.m`; +- `sliding_window.m`; +- `ul_test.m`; +- `update_window.m`. + +Integration-only modification: + +- `init.m` retains its compilation/setup behavior but additionally searches + `DINGO_POLYTOPE_SAMPLER_MATLAB` and common Dingo-relative locations. +- `relative_volume.m` retains the supplied estimator and adds the two interface + initializations normally performed by PolytopeSamplerMatlab `sample.m`: + the shared `startTime` timeout clock and the seeded `simdTwister` stream. + +New bridge: + +- `volume_update_batch.m` reads a Python-generated MAT request, invokes + `relative_volume` once per cut, appends the cut to the next parent, strips + nonportable MATLAB handle objects from saved diagnostics, and writes a v7 MAT + result readable by SciPy. + +No annealing, ratio-estimation, convergence, CRHMC-transition, or sequential +update mathematics in `relative_volume.m` was rewritten. + +The ZIP's model data, standalone random-cut runner documentation, HTML/Markdown +exposition, and Git metadata are not installed into the Dingo Python package +because they are not runtime dependencies. The supplied TeX theory is +preserved under `docs/volume_updating/theory.tex`, the original MIT license is +preserved under `docs/volume_updating/MetabolicPolytopesCuttingPlanes-LICENSE`, +and exact source hashes and both complete integration patches are under `docs/`. + +## 6. Bound and GPR handling + +The volume wrapper accepts linear cuts. A biological application may first use +a GPR rule and severity model to produce child reaction bounds. The helper +`cuts_from_bounds` converts every nested bound change into the correct scalar +halfspace: + +- upper bound `v_k <= u'_k` becomes normal `+e_k` and threshold `u'_k`; +- lower bound `v_k >= l'_k` becomes normal `-e_k` and threshold `-l'_k`. + +Both sides are emitted for a tightened reversible reaction. Any attempted bound +relaxation is rejected because the static telescoping update requires a nested +child. + +GPR semantics, severity sweeps, angles and RL outputs remain application +layers. They do not change the volume estimator. + +## 7. Parameter policy + +The wrapper accepts only parameters present in the supplied static API: + +`r`, `delta`, `alpha`, `nu`, `N_utest`, `window_size`, `epsilon`, +`max_annealing_iter`, `max_bisect_iter`, `simdLen`, `t_max_safety`, and `verb`. + +Wrapper-only options are `algorithm` and `wrapper_timeout_sec`. +`algorithm` must be `static`, `baseline`, or `relative_volume`. A requested +`reuse` algorithm raises an error rather than silently running a different +method. Unsupported fields such as `ratio_ess_min` are rejected. + +This means parameter changes are explicit while the underlying static logic +remains the supplied logic. + +When no override is provided, the wrapper uses the ZIP's exact defaults: +`r=0.10`, `delta=0.05`, `alpha=0.20`, `nu=10`, `N_utest=320`, +`window_size=5000`, `epsilon=0.10`, `max_annealing_iter=200`, +`max_bisect_iter=100`, `simdLen=8`, `t_max_safety=1.01`, and `verb=1`. +The prior `alpha=0.10`, `N_utest=125`, `epsilon=0.05` profile is therefore not +silently imposed by this wrapper; it must be requested as a named experiment. + +## 8. Result contract + +The wrapper returns: + +- cumulative `rho` and `log_rho`; +- per-cut ratios and log ratios; +- number of annealing phases per cut; +- total CRHMC step counts per cut; +- per-cut convergence flags; +- cut labels; +- `algorithm='static_relative_volume_sequential'`; +- `objective_floor=False`. + +It does not invent fresh-sample, reused-pilot or cross-update-state statistics +that the static implementation does not expose. + +## 9. Testing completed + +Nine Python wrapper tests pass: + +- two-cut sequential result parsing; +- identity behavior for an empty cut stream; +- bound-to-halfspace conversion including reversible two-sided tightening; +- rejection of bound relaxation; +- rejection of dimension mismatch; +- rejection of a fake `reuse` mode; +- rejection of unsupported options; +- rejection of malformed source numeric options. +- verification that the required sampler clock and random stream are initialized. + +The Python modules compile successfully. A real sequential integration test +was executed through WSL using Windows MATLAB and PolytopeSamplerMatlab. For +the exact two-cut toy problem, the estimator returned rho=0.0721726401 versus +11/150=0.0733333333, corresponding to 1.58% relative error. Both cuts +converged, the phase counts were [1, 2], the CRHMC step counts were +[2488, 4072], and objective_floor=False. Total wall-clock time was 15.25 s. + +## 10. Correct experimental comparison + +The three scientific columns should remain distinct: + +1. objective-free Dingo/volesti reference; +2. supplied static CRHMC, freshly applied to each sequential cut; +3. any future optimized reuse implementation, if built, tested against column + 2 and clearly separated from the static wrapper. + +The wrapper delivered here implements column 2 only. This is intentional and +matches the requested architecture. + +## 11. Deliverable + +The `dynamic_volume_updating.zip` archive can be unzipped directly into the +repository. It places the wrapper files into the existing Python package and +includes the report, tests, doctor, theory, source hashes, and exact integration +diff. diff --git a/pyproject.toml b/pyproject.toml index d0cdf08..c4f1f8e 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,6 +4,10 @@ version = "0.1.0" description = "A python library for metabolic networks sampling and analysis. dingo is part of GeomScale project" authors = ["Apostolos Chalkis "] build = 'build.py' +include = [ + { path = "dingo/matlab/volume_updating/*.m", format = ["sdist", "wheel"] }, + { path = "dingo/matlab/volume_updating/LICENSE", format = ["sdist", "wheel"] }, +] [tool.poetry.dependencies] python = "^3.8" diff --git a/setup.py b/setup.py index 77ec7ec..f019a5d 100755 --- a/setup.py +++ b/setup.py @@ -109,6 +109,12 @@ author_email=author_email, name=name, packages=packages, + package_data={ + "dingo": [ + "matlab/volume_updating/*.m", + "matlab/volume_updating/LICENSE", + ], + }, ext_modules=ext_modules, ) diff --git a/tests/test_volume_updating.py b/tests/test_volume_updating.py new file mode 100644 index 0000000..cbc9710 --- /dev/null +++ b/tests/test_volume_updating.py @@ -0,0 +1,174 @@ +from __future__ import annotations + +import importlib.util +from pathlib import Path +import sys +import unittest + +import numpy as np +from scipy.io import loadmat, savemat + + +ROOT = Path(__file__).resolve().parents[1] +MODULE_PATH = ROOT / "dingo" / "dynamic_volume.py" +if not MODULE_PATH.is_file(): + MODULE_PATH = ROOT / "dingo" / "dynamic_volume.py" +SPEC = importlib.util.spec_from_file_location("volume_updating_test", MODULE_PATH) +MODULE = importlib.util.module_from_spec(SPEC) +sys.modules[SPEC.name] = MODULE +SPEC.loader.exec_module(MODULE) + +HalfspaceCut = MODULE.HalfspaceCut +StaticVolumeUpdater = MODULE.StaticVolumeUpdater +cuts_from_bounds = MODULE.cuts_from_bounds + + +class FakeUpdater(StaticVolumeUpdater): + last_request = None + + def _run_matlab(self, request_path, output_path, timeout): + self.last_request = loadmat(request_path, squeeze_me=True, struct_as_record=False) + n_cuts = np.atleast_2d(self.last_request["cut_normals"]).shape[0] + ratios = np.array([0.5, 0.25], dtype=float)[:n_cuts] + logs = np.log(ratios) + savemat( + output_path, + { + "rho": float(np.prod(ratios)), + "log_rho": float(np.sum(logs)), + "cut_ratios": ratios, + "cut_log_ratios": logs, + "cut_n_phases": np.arange(1, n_cuts + 1), + "cut_total_steps": 10 * np.arange(1, n_cuts + 1), + "cut_converged": np.ones(n_cuts, dtype=np.uint8), + "algorithm": "static_relative_volume_sequential", + "objective_floor": np.uint8(0), + }, + ) + + +class VolumeUpdatingWrapperTests(unittest.TestCase): + @classmethod + def setUpClass(cls): + cls.matlab_dir = MODULE_PATH.parent / "matlab" / "volume_updating" + + def updater(self): + return FakeUpdater( + matlab_executable="/bin/true", + matlab_source_dir=self.matlab_dir, + ) + + def test_sequential_static_result_contract(self): + updater = self.updater() + cuts = [ + HalfspaceCut(np.array([1.0, 0.0]), 0.5, "x"), + HalfspaceCut(np.array([0.0, 1.0]), 0.25, "y"), + ] + result = updater.estimate( + np.zeros((0, 2)), + np.zeros(2), + np.ones(2), + cuts, + options={"epsilon": 0.1, "algorithm": "static"}, + ) + self.assertAlmostEqual(result.rho, 0.125) + self.assertAlmostEqual(result.log_rho, np.log(0.125)) + np.testing.assert_allclose(result.cut_ratios, [0.5, 0.25]) + self.assertEqual(result.cut_labels, ("x", "y")) + self.assertFalse(result.objective_floor) + np.testing.assert_allclose( + np.atleast_2d(updater.last_request["cut_normals"]), + [[1.0, 0.0], [0.0, 1.0]], + ) + + def test_empty_cut_stream_is_identity(self): + result = self.updater().estimate( + np.zeros((0, 2)), np.zeros(2), np.ones(2), [] + ) + self.assertEqual(result.rho, 1.0) + self.assertEqual(result.cut_ratios.size, 0) + + def test_reuse_algorithm_is_not_faked(self): + with self.assertRaisesRegex(ValueError, "does not provide a separate reuse"): + self.updater().estimate( + np.zeros((0, 2)), + np.zeros(2), + np.ones(2), + [HalfspaceCut(np.ones(2), 1.0)], + options={"algorithm": "reuse"}, + ) + + def test_unknown_non_source_option_is_rejected(self): + with self.assertRaisesRegex(ValueError, "not part of the supplied static API"): + self.updater().estimate( + np.zeros((0, 2)), + np.zeros(2), + np.ones(2), + [HalfspaceCut(np.ones(2), 1.0)], + options={"ratio_ess_min": 3000}, + ) + + def test_invalid_static_option_is_rejected(self): + with self.assertRaisesRegex(ValueError, "positive integer"): + self.updater().estimate( + np.zeros((0, 2)), + np.zeros(2), + np.ones(2), + [HalfspaceCut(np.ones(2), 1.0)], + options={"N_utest": 12.5}, + ) + with self.assertRaisesRegex(ValueError, r"r \+ delta"): + self.updater().estimate( + np.zeros((0, 2)), + np.zeros(2), + np.ones(2), + [HalfspaceCut(np.ones(2), 1.0)], + options={"r": 0.8, "delta": 0.3}, + ) + with self.assertRaisesRegex(ValueError, r"r \+ delta"): + self.updater().estimate( + np.zeros((0, 2)), + np.zeros(2), + np.ones(2), + [HalfspaceCut(np.ones(2), 1.0)], + options={"r": 0.98}, + ) + + def test_cut_dimension_is_checked(self): + with self.assertRaisesRegex(ValueError, "expected 2"): + self.updater().estimate( + np.zeros((0, 2)), + np.zeros(2), + np.ones(2), + [HalfspaceCut(np.ones(3), 1.0)], + ) + + def test_bound_update_compiles_to_two_sided_cuts(self): + cuts = cuts_from_bounds( + np.array([-2.0, 0.0]), + np.array([2.0, 4.0]), + np.array([-1.0, 0.0]), + np.array([1.5, 3.0]), + ) + self.assertEqual([cut.label for cut in cuts], ["ub[0]", "lb[0]", "ub[1]"]) + np.testing.assert_allclose(cuts[0].normal, [1.0, 0.0]) + np.testing.assert_allclose(cuts[1].normal, [-1.0, 0.0]) + self.assertEqual(cuts[1].threshold, 1.0) + + def test_bound_relaxation_is_rejected(self): + with self.assertRaisesRegex(ValueError, "subset"): + cuts_from_bounds([0.0], [1.0], [-0.1], [1.0]) + + def test_sampler_interface_initialization_is_present(self): + source = ( + MODULE_PATH.parent + / "matlab" + / "volume_updating" + / "relative_volume.m" + ).read_text(encoding="utf-8") + self.assertIn("sampler_opts.startTime = start_time;", source) + self.assertIn("rng(sampler_opts.seed, 'simdTwister');", source) + + +if __name__ == "__main__": + unittest.main()