diff --git a/.gitignore b/.gitignore index 23c482eb..d8e930b6 100644 --- a/.gitignore +++ b/.gitignore @@ -12,4 +12,4 @@ volestipy.egg-info venv lp_solve_5.5/ .devcontainer/ -.github/dependabot.yml \ No newline at end of file +.github/dependabot.yml.venv/ diff --git a/.python-version b/.python-version new file mode 100644 index 00000000..143c2f5d --- /dev/null +++ b/.python-version @@ -0,0 +1 @@ +3.8.19 diff --git a/dingo/PolytopeSampler.py b/dingo/PolytopeSampler.py index 27f37853..a683ef92 100644 --- a/dingo/PolytopeSampler.py +++ b/dingo/PolytopeSampler.py @@ -6,6 +6,7 @@ # Licensed under GNU LGPL.3, see LICENCE file +# Contributed and/or modified by Iva Janković, as part of Google Summer of Code 2025 program. import numpy as np import warnings @@ -29,12 +30,12 @@ def __init__(self, metabol_net): raise Exception("An unknown input object given for initialization.") self._metabolic_network = metabol_net - self._A = [] - self._b = [] - self._N = [] - self._N_shift = [] - self._T = [] - self._T_shift = [] + self._A = None + self._b = None + self._N = None + self._N_shift = None + self._T = None + self._T_shift = None self._parameters = {} self._parameters["nullspace_method"] = "sparseQR" self._parameters["opt_percentage"] = self.metabolic_network.parameters[ @@ -53,14 +54,15 @@ def get_polytope(self): """ if ( - self._A == [] - or self._b == [] - or self._N == [] - or self._N_shift == [] - or self._T == [] - or self._T_shift == [] + self._A is None + or self._b is None + or self._N is None + or self._N_shift is None + or self._T is None + or self._T_shift is None ): + ( max_flux_vector, max_objective, @@ -166,7 +168,7 @@ def generate_steady_states_no_multiphase( """A member function to sample steady states. Keyword arguments: - method -- An MCMC method to sample, i.e. {'billiard_walk', 'cdhr', 'rdhr', 'ball_walk', 'dikin_walk', 'john_walk', 'vaidya_walk', 'gaussian_hmc_walk', 'exponential_hmc_walk', 'hmc_leapfrog_gaussian', 'hmc_leapfrog_exponential'} + method -- An MCMC method to sample, i.e. {'billiard_walk', 'cdhr', 'rdhr', 'ball_walk', 'dikin_walk', 'john_walk', 'vaidya_walk', 'gaussian_hmc_walk', 'exponential_hmc_walk', 'hmc_leapfrog_gaussian', 'hmc_leapfrog_exponential', 'shake_and_bake', 'billiard_shake_and_bake'} n -- the number of steady states to sample burn_in -- the number of points to burn before sampling thinning -- the walk length of the chain @@ -190,6 +192,40 @@ def generate_steady_states_no_multiphase( return steady_states + def generate_steady_states_sb_once(self,n=1000,burn_in=0,thinning=1,variance=1.0,bias_vector=None,ess=0): + """ + One Shake and Bake phase: samples n points, returns (steady_states, diagnostics) + diagnostics = {'minESS':..., 'maxPSRF':..., 'N':..., 'phases': 1, 'seconds': ...}. + """ + self.get_polytope() + P = HPolytope(self._A, self._b) + + if bias_vector is None: + bias_vector = np.ones(self._A.shape[1], dtype=np.float64) + else: + bias_vector = np.asarray(bias_vector, dtype=np.float64) + if bias_vector.shape[0] != self._A.shape[1]: + raise ValueError(f"bias_vector length {bias_vector.shape[0]} != {self._A.shape[1]}") + + samples = P.generate_samples(b"shake_and_bake",int(n),int(burn_in),int(thinning),float(variance), bias_vector,self._parameters["solver"],int(ess),) + diag_buf = np.empty(5, dtype=np.float64) + minESS = maxPSRF = seconds = np.nan + Ncpp = phases = 0 + + P.get_sb_diagnostics(np.ascontiguousarray(diag_buf)) + minESS, maxPSRF, Ncpp, phases, seconds = diag_buf + + steady_states = map_samples_to_steady_states(samples.T, self._N, self._N_shift) + + diagnostics = { + "minESS": float(minESS), + "maxPSRF": float(maxPSRF), + "N": int(Ncpp), + "phases": int(phases), + "seconds": float(seconds) if np.isfinite(seconds) else None, + } + return steady_states, diagnostics + @staticmethod def sample_from_polytope( A, b, ess=1000, psrf=False, parallel_mmcs=False, num_threads=1, solver=None @@ -223,7 +259,7 @@ def sample_from_polytope_no_multiphase( Keyword arguments: A -- an mxn matrix that contains the normal vectors of the facets of the polytope row-wise b -- a m-dimensional vector, s.t. A*x <= b - method -- An MCMC method to sample, i.e. {'billiard_walk', 'cdhr', 'rdhr', 'ball_walk', 'dikin_walk', 'john_walk', 'vaidya_walk', 'gaussian_hmc_walk', 'exponential_hmc_walk', 'hmc_leapfrog_gaussian', 'hmc_leapfrog_exponential'} + method -- An MCMC method to sample, i.e. {'billiard_walk', 'cdhr', 'rdhr', 'ball_walk', 'dikin_walk', 'john_walk', 'vaidya_walk', 'gaussian_hmc_walk', 'exponential_hmc_walk', 'hmc_leapfrog_gaussian', 'hmc_leapfrog_exponential', 'shake_and_bake', 'billiard_shake_and_bake'} n -- the number of steady states to sample burn_in -- the number of points to burn before sampling thinning -- the walk length of the chain @@ -239,6 +275,36 @@ def sample_from_polytope_no_multiphase( samples_T = samples.T return samples_T + + @staticmethod + def sample_from_polytope_sb_once( + A, b, n=1000, burn_in=0, thinning=1, variance=1.0, bias_vector=None, solver=None, ess=0 + ): + """ + One Shake and Bake phase for polytope defined by A,b. + Returns (samples_T_d_by_N, diagnostics_dict). + """ + if bias_vector is None: + bias_vector = np.ones(A.shape[1], dtype=np.float64) + else: + bias_vector = bias_vector.astype("float64") + + P = HPolytope(A, b) + samples = P.generate_samples(b"shake_and_bake",n,burn_in,thinning,variance,bias_vector,solver,ess) + + diag = np.zeros(5, dtype=np.float64) + P.get_sb_diagnostics(diag) + minESS, maxPSRF, Ncpp, phases, seconds = diag + + diagnostics = { + "minESS": float(minESS), + "maxPSRF": float(maxPSRF), + "N": int(Ncpp), + "phases": int(phases), + "seconds": float(seconds) if not np.isnan(seconds) else None, + } + return samples.T, diagnostics + @staticmethod def round_polytope( @@ -352,4 +418,4 @@ def set_tol(self, value): def set_opt_percentage(self, value): - self._parameters["opt_percentage"] = value + self._parameters["opt_percentage"] = value \ No newline at end of file diff --git a/dingo/bindings/bindings.cpp b/dingo/bindings/bindings.cpp index cbdb1246..0fd19cf4 100644 --- a/dingo/bindings/bindings.cpp +++ b/dingo/bindings/bindings.cpp @@ -6,6 +6,7 @@ // Contributed and/or modified by Haris Zafeiropoulos // Contributed and/or modified by Pedro Zuidberg Dos Martires +// Contributed and/or modified by Iva Janković, as part of Google Summer of Code 2025 program. // Licensed under GNU LGPL.3, see LICENCE file @@ -14,7 +15,7 @@ #include #include "bindings.h" #include "hmc_sampling.h" - +#include "random_walks/shake_and_bake_walk.hpp" using namespace std; @@ -134,6 +135,43 @@ double HPolytopeCPP::apply_sampling(int walk_len, } else if (strcmp(method, "vaidya_walk") == 0) { // vaidya walk uniform_sampling(rand_points, HP, rng, walk_len, number_of_points, starting_point, number_of_points_to_burn); + } else if (strcmp(method, "shake_and_bake") == 0) { + using clock = std::chrono::high_resolution_clock; + auto t0 = clock::now(); + + auto [boundary_pt, facet_idx] = compute_boundary_point(HP, rng, static_cast(1e-8)); + + const int d = HP.dimension(); + const int n_phase = number_of_points; + int nburn = number_of_points_to_burn; + + std::list batch; + shakeandbake_sampling(batch, HP, rng, walk_len, n_phase, boundary_pt, nburn, facet_idx); + rand_points.swap(batch); + const int N = static_cast(rand_points.size()); + MT S(d, N); + int c = 0; + for (auto it = rand_points.cbegin(); it != rand_points.cend(); ++it, ++c) + for (int j = 0; j < d; ++j) + S(j, c) = (*it)[j]; + + unsigned int min_ess_u = 0; + VT ess_vec = effective_sample_size(S, min_ess_u); + NT min_ess = ess_vec.minCoeff(); + + VT rhat_vec = univariate_psrf(S); + NT max_psrf = rhat_vec.maxCoeff(); + + auto t1 = clock::now(); + double secs = std::chrono::duration(t1 - t0).count(); + sb_samples_ = S; // d x N + sb_diag_.minESS = static_cast(min_ess); + sb_diag_.maxPSRF = static_cast(max_psrf); + sb_diag_.N = N; + sb_diag_.phases = 1; + sb_diag_.seconds = secs; + + } else if (strcmp(method, "mmcs") == 0) { // vaidya walk MT S; int total_ess; @@ -167,7 +205,7 @@ double HPolytopeCPP::apply_sampling(int walk_len, } else { - throw std::runtime_error("This function must not be called."); + throw std::runtime_error(method + std::string(" is not recognized as a valid sampling method.")); } if (strcmp(method, "mmcs") != 0) { @@ -428,6 +466,23 @@ void HPolytopeCPP::get_mmcs_samples(double* T_matrix, double* T_shift, double* s mmcs_set_of_parameters.samples.resize(0,0); } +void HPolytopeCPP::get_sb_samples(double* samples) const { + const int d = sb_samples_.rows(); + const int N = sb_samples_.cols(); + int t = 0; + for (int i = 0; i < d; ++i) + for (int j = 0; j < N; ++j) + samples[t++] = sb_samples_(i, j); +} + +void HPolytopeCPP::get_sb_diagnostics(double* out5) const { + out5[0] = sb_diag_.minESS; + out5[1] = sb_diag_.maxPSRF; + out5[2] = static_cast(sb_diag_.N); + out5[3] = static_cast(sb_diag_.phases); + out5[4] = sb_diag_.seconds; +} + ////////// Start of "rounding()" ////////// void HPolytopeCPP::apply_rounding(int rounding_method, double* new_A, double* new_b, diff --git a/dingo/bindings/bindings.h b/dingo/bindings/bindings.h index 553ea99e..5bfe3ac5 100644 --- a/dingo/bindings/bindings.h +++ b/dingo/bindings/bindings.h @@ -5,7 +5,8 @@ // Copyright (c) 2018-2021 Apostolos Chalkis // Contributed and/or modified by Haris Zafeiropoulos -// Contributed and/or modified by Pedro Zuidberg Dos Martires +// Contributed and/or modified by Pedro Zuidberg Dos +// Contributed and/or modified by Iva Janković, as part of Google Summer of Code 2025 program. // Licensed under GNU LGPL.3, see LICENCE file @@ -29,16 +30,20 @@ #include "sampling/mmcs.hpp" #include "sampling/parallel_mmcs.hpp" #include "diagnostics/univariate_psrf.hpp" +#include "diagnostics/effective_sample_size.hpp" //from generate_samples, some extra headers not already included #include +#include #include "sampling/sampling.hpp" #include "ode_solvers/ode_solvers.hpp" +#include "preprocess/feasible_point.hpp" // for rounding #include "preprocess/min_sampling_covering_ellipsoid_rounding.hpp" #include "preprocess/svd_rounding.hpp" #include "preprocess/inscribed_ellipsoid_rounding.hpp" +#include "preprocess/feasible_point.hpp" typedef double NT; typedef Cartesian Kernel; @@ -48,6 +53,13 @@ typedef typename Hpolytope::MT MT; typedef typename Hpolytope::VT VT; typedef BoostRandomNumberGenerator RNGType; +struct Diagnostics { + double minESS = 0.0; + double maxPSRF = 0.0; + long long N = 0; + int phases = 0; + double seconds = 0.0; +}; template struct mmcs_parameters @@ -155,6 +167,16 @@ class HPolytopeCPP{ void apply_rounding(int rounding_method, double* new_A, double* new_b, double* T_matrix, double* shift, double &round_value, double* inner_point, double radius); + // Shake and Bake one-phase sampling methods + inline const Diagnostics& sb_diagnostics() const { return sb_diag_; } + inline const MT& sb_samples() const { return sb_samples_; } + void get_sb_samples(double* samples) const; + void get_sb_diagnostics(double* out5) const; + + private: + Diagnostics sb_diag_; + MT sb_samples_; + }; diff --git a/dingo/volestipy.pyx b/dingo/volestipy.pyx index 6e8d3a98..ddade0d0 100644 --- a/dingo/volestipy.pyx +++ b/dingo/volestipy.pyx @@ -9,6 +9,8 @@ # Licensed under GNU LGPL.3, see LICENCE file +# Contributed and/or modified by Iva Janković, as part of Google Summer of Code 2025 program. + #!python #cython: language_level=3 #cython: boundscheck=False @@ -72,6 +74,9 @@ cdef extern from "bindings.h": void apply_rounding(int rounding_method, double* new_A, double* new_b, double* T_matrix, \ double* shift, double &round_value, double* inner_point, double radius); + void get_sb_diagnostics(double* out5) const + void get_sb_samples(double* samples) const + # The lowDimPolytopeCPP class along with its functions cdef cppclass lowDimHPolytopeCPP: @@ -85,7 +90,7 @@ cdef extern from "bindings.h": # Lists with the methods supported by volesti for volume approximation and random walk volume_methods = ["sequence_of_balls".encode("UTF-8"), "cooling_gaussian".encode("UTF-8"), "cooling_balls".encode("UTF-8")] walk_methods = ["uniform_ball".encode("UTF-8"), "CDHR".encode("UTF-8"), "RDHR".encode("UTF-8"), "gaussian_ball".encode("UTF-8"), \ - "gaussian_CDHR".encode("UTF-8"), "gaussian_RDHR".encode("UTF-8"), "uniform_ball".encode("UTF-8"), "billiard".encode("UTF-8")] + "gaussian_CDHR".encode("UTF-8"), "gaussian_RDHR".encode("UTF-8"), "uniform_ball".encode("UTF-8"), "billiard".encode("UTF-8"),"shake_and_bake".encode("UTF-8"),"billiard_shake_and_bake".encode("UTF-8") ] rounding_methods = ["min_ellipsoid".encode("UTF-8"), "svd".encode("UTF-8"), "max_ellipsoid".encode("UTF-8")] # Build the HPolytope class @@ -136,6 +141,33 @@ cdef class HPolytope: variance_value, &bias_vector_[0], ess) return np.asarray(samples) + def get_sb_diagnostics(self, out): + """ + out: np.ndarray float64, shape (5,), redosled: + [minESS, maxPSRF, N, phases, seconds] + """ + cdef double[::1] buf = np.ascontiguousarray(out, dtype=np.float64) + if buf.shape[0] < 5: + raise ValueError("out must have length >= 5") + self.polytope_cpp.get_sb_diagnostics(&buf[0]) + return np.asarray(buf) + + def get_sb_samples(self): + """ + Return d x N numpy matrix of samples from Shake and Bake sampling + """ + cdef int d = self._A.shape[1] + + cdef double[::1] tmp = np.zeros(5, dtype=np.float64, order="C") + self.polytope_cpp.get_sb_diagnostics(&tmp[0]) + cdef Py_ssize_t N = tmp[2] + if N <= 0: + return np.zeros((d, 0), dtype=np.float64) + + cdef double[:,::1] S = np.zeros((d, N), dtype=np.float64, order="C") + self.polytope_cpp.get_sb_samples(&S[0,0]) + return np.asarray(S) + # The rounding() function; as in compute_volume, more than one method is available for this step def rounding(self, rounding_method = 'john_position', solver = None): diff --git a/eigen b/eigen index 02f42001..5e8edd21 160000 --- a/eigen +++ b/eigen @@ -1 +1 @@ -Subproject commit 02f420012a169ed9267a8a78083aaa588e713353 +Subproject commit 5e8edd21863b8321fc6b9c82322e6cc8cfc47c14 diff --git a/pyproject.toml b/pyproject.toml index d0cdf087..f0c3cb74 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -24,4 +24,4 @@ networkx = "3.1" [build-system] requires = ["poetry-core>=1.0.0", "cython", "numpy==1.20.1"] -build-backend = "poetry.core.masonry.api" +build-backend = "poetry.core.masonry.api" \ No newline at end of file diff --git a/setup.py b/setup.py old mode 100755 new mode 100644 diff --git a/tests/sampling_no_multiphase.py b/tests/sampling_no_multiphase.py index 86e304f3..4bc54e52 100644 --- a/tests/sampling_no_multiphase.py +++ b/tests/sampling_no_multiphase.py @@ -76,4 +76,4 @@ def test_sample_sbml(self): if len(sys.argv) > 1: set_default_solver(sys.argv[1]) sys.argv.pop(1) - unittest.main() + unittest.main() \ No newline at end of file diff --git a/volesti b/volesti index c3109bba..a33e5f45 160000 --- a/volesti +++ b/volesti @@ -1 +1 @@ -Subproject commit c3109bba06a9b623446bdde4c5fadb02722de876 +Subproject commit a33e5f450a8badb506384d7a619344dca855fe75