dingo implementations Shake and Bake / Billiard Shake and Bake

.gitignore

@@ -12,4 +12,4 @@ volestipy.egg-info
 venv
 lp_solve_5.5/
 .devcontainer/
-.github/dependabot.yml
+.github/dependabot.yml.venv/

.python-version

@@ -0,0 +1 @@
+3.8.19

dingo/PolytopeSampler.py

@@ -6,10 +6,13 @@
 
 # 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
+from typing import Optional
 import math
+import time
 from dingo.MetabolicNetwork import MetabolicNetwork
 from dingo.utils import (
     map_samples_to_steady_states,
@@ -29,12 +32,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[
@@ -46,21 +49,23 @@ def __init__(self, metabol_net):
 
         self._parameters["tol"] = 1e-06
         self._parameters["solver"] = None
+        self._last_hpoly = None
 
     def get_polytope(self):
         """A member function to derive the corresponding full dimensional polytope
         and a isometric linear transformation that maps the latter to the initial space.
         """
 
         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 +171,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}
         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 +195,50 @@ def generate_steady_states_no_multiphase(
 
         return steady_states
 
+    def generate_steady_states_sb_once(self, n=1000, burn_in=0, sampler="sb", nreflections=None, walk_len=None):
+        """
+        Single-phase boundary sampler with clear separation of concerns.
+        - sampler: "sb"/"shake_and_bake" or "bsb"/"billiard_shake_and_bake"
+        - walk_len: default ceil(sqrt(d)) with a minimum of 5
+        - nreflections: only used for BSB; defaults to ceil(sqrt(d))
+        """
+
+        self.get_polytope()
+        P = HPolytope(self._A, self._b)
+        d = int(self._A.shape[1])
+
+        wl = int(walk_len) if walk_len is not None else max(5, int(np.sqrt(d)))
+        use_bsb = str(sampler).lower() in ("bsb", "billiard_shake_and_bake")
+        nref = (
+            int(nreflections)
+            if nreflections is not None
+            else (int(np.ceil(np.sqrt(d))) if use_bsb else 0)
+        )
+
+        S = P.boundary_sample(
+            sampler="bsb" if use_bsb else "sb",
+            number_of_points=int(n),
+            number_of_points_to_burn=int(burn_in),
+            walk_len=int(wl),
+            nreflections=int(nref),
+        )
+
+        finite_cols = np.isfinite(S).all(axis=0)
+        if not finite_cols.all():
+            S = S[:, finite_cols]
+        if S.shape[1] == 0:
+            raise RuntimeError(
+                "All sample columns contain NaN/Inf; check walk_len, nreflections, or polytope constraints."
+            )
+
+        steady_states = map_samples_to_steady_states(S, self._N, self._N_shift)
+
+        self._last_hpoly = P  # cache if needed later
+        return steady_states
+
+
+
+
     @staticmethod
     def sample_from_polytope(
         A, b, ess=1000, psrf=False, parallel_mmcs=False, num_threads=1, solver=None
@@ -223,7 +272,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 +288,65 @@ def sample_from_polytope_no_multiphase(
 
         samples_T = samples.T
         return samples_T
+
+    @staticmethod
+    def sample_from_polytope_sb_once(
+        A,
+        b,
+        n: int = 1000,
+        burn_in: int = 0,
+        sampler: str = "sb",
+        walk_len: Optional[int] = None,
+        nreflections: Optional[int] = None,
+    ):
+        """
+        Perform a single boundary-sampling phase on Ax <= b using SB or BSB.
+        Returns only the sample matrix (d x N). No timing or diagnostics inside.
+        """
+
+        import numpy as np
+
+        A = np.ascontiguousarray(A, dtype=np.float64)
+        b = np.ascontiguousarray(b, dtype=np.float64)
+        P = HPolytope(A, b)
+        d = int(A.shape[1])
+        wl = int(walk_len) if walk_len is not None else max(5, int(np.sqrt(d)))
+        use_bsb = str(sampler).lower() in ("bsb", "billiard_shake_and_bake")
+        nref = (
+            int(nreflections)
+            if nreflections is not None
+            else (int(np.ceil(np.sqrt(d))) if use_bsb else 0)
+        )
+
+        S = P.boundary_sample(
+            sampler="bsb" if use_bsb else "sb",
+            number_of_points=int(n),
+            number_of_points_to_burn=int(burn_in),
+            walk_len=int(wl),
+            nreflections=int(nref),
+        )
+
+        finite_cols = np.isfinite(S).all(axis=0)
+        if not finite_cols.all():
+            S = S[:, finite_cols]
+        if S.shape[1] == 0:
+            raise RuntimeError(
+                "All sample columns contain NaN/Inf; adjust walk_len, nreflections, or check constraints."
+            )
+
+        return S
+
+    def boundary_diagnostics(self, S):
+        """
+        Compute diagnostics (minESS, maxPSRF, N) for a given sample matrix S
+        using the current polytope stored in this sampler instance.
+        """
+        S = np.ascontiguousarray(S, dtype=np.float64)
+        P = HPolytope(self._A, self._b)
+        diag = P.boundary_diag(S)
+        return diag
+
+
 
     @staticmethod
     def round_polytope(
@@ -352,4 +460,4 @@ def set_tol(self, value):
 
     def set_opt_percentage(self, value):
 
-        self._parameters["opt_percentage"] = value
+        self._parameters["opt_percentage"] = value

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,9 @@
 #include <stdexcept>
 #include "bindings.h"
 #include "hmc_sampling.h"
-
+#include "random_walks/shake_and_bake_walk.hpp"
+#include "random_walks/billiard_shake_and_bake_walk.hpp"
+#include "diagnostics/scaling_ratio.hpp" 
 
 using namespace std;
 
@@ -92,7 +95,8 @@ double HPolytopeCPP::apply_sampling(int walk_len,
                                     double* samples,
                                     double variance_value,
                                     double* bias_vector_,
-                                    int ess){
+                                    int ess,
+                                    int nreflections) {
 
    RNGType rng(HP.dimension());
    HP.normalize();
@@ -167,7 +171,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) {
@@ -181,7 +185,67 @@ double HPolytopeCPP::apply_sampling(int walk_len,
    }
    return 0.0;
 }
+
+// Boundary sampling: shakre-and-bake and billiard shake-and-bake
+int HPolytopeCPP::apply_boundary_sampling(int walk_len,
+                                          int number_of_points,
+                                          int number_of_points_to_burn,
+                                          const char* sampler,
+                                          int nreflections,
+                                          double* samples) {
+    RNGType rng(HP.dimension());
+    HP.normalize();
+
+    auto [boundary_pt, facet_idx] = compute_boundary_point<Point>(HP, rng, static_cast<NT>(1e-8));
+
+    const int d = HP.dimension();
+    std::list<Point> batch;
+
+    if (strcmp(sampler, "shake_and_bake") == 0 || strcmp(sampler, "sb") == 0) {
+        shakeandbake_sampling<ShakeAndBakeWalk>(batch, HP, rng, walk_len, number_of_points, boundary_pt, number_of_points_to_burn, facet_idx);
+
+    } else if (strcmp(sampler, "billiard_shake_and_bake") == 0 || strcmp(sampler, "bsb") == 0) {
+        billiard_shakeandbake_sampling<BilliardShakeAndBakeWalk>(batch, HP, rng, walk_len,nreflections, number_of_points,boundary_pt, number_of_points_to_burn,facet_idx);
+
+    } else {
+        throw std::runtime_error(std::string(sampler) + " is not a boundary sampler.");
+    }
+
+    int n_si = 0;
+    for (auto it = batch.cbegin(); it != batch.cend(); ++it)
+        for (int j = 0; j < d; ++j)
+            samples[n_si++] = (*it)[j];
+
+    return static_cast<int>(batch.size());
+}
+
 //////////         End of "generate_samples()"          //////////
+SBDiagnostics HPolytopeCPP::sb_diagnostics(int d, int N, const double* samples) {
+    MT S(d, N);
+    for (int c = 0; c < N; ++c)
+        for (int r = 0; r < d; ++r)
+            S(r, c) = samples[c*d + r];
+
+    unsigned int min_ess_u = 0;
+    VT ess_vec  = effective_sample_size<NT, VT, MT>(S, min_ess_u);
+    VT rhat_vec = univariate_psrf<NT, VT, MT>(S);
+
+    SBDiagnostics out;
+    out.minESS  = static_cast<double>(ess_vec.minCoeff());
+    out.maxPSRF = static_cast<double>(rhat_vec.maxCoeff());
+    out.N       = N;
+    return out;
+}
+
+void HPolytopeCPP::set_sb_state_from_buffer(int d, int N, const double* samples, const SBDiagnostics& diag) {
+    sb_samples_.resize(d, N);
+    for (int j = 0; j < N; ++j)
+        for (int i = 0; i < d; ++i)
+            sb_samples_(i, j) = samples[i + j * d];
+
+    sb_diag_ = diag; 
+}
+
 
 
 void HPolytopeCPP::get_polytope_as_matrices(double* new_A, double* new_b) const {
@@ -428,6 +492,21 @@ 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* out) const {
+    const int d = static_cast<int>(sb_samples_.rows());
+    const int N = static_cast<int>(sb_samples_.cols());
+    for (int j = 0; j < N; ++j)
+        for (int i = 0; i < d; ++i)
+            out[i + j * d] = sb_samples_(i, j);
+}
+
+void HPolytopeCPP::get_sb_diagnostics(double* out3) const {
+    out3[0] = sb_diag_.minESS;
+    out3[1] = sb_diag_.maxPSRF;
+    out3[2] = static_cast<double>(sb_diag_.N);
+
+}
+
 
 //////////         Start of "rounding()"          //////////
 void HPolytopeCPP::apply_rounding(int rounding_method, double* new_A, double* new_b,