From eeb6ff8b21f67ebedf0e05cfb98f0825a8dad1b8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Iva=20Jankovi=C4=87?= Date: Mon, 18 Aug 2025 20:00:53 +0200 Subject: [PATCH 01/14] Dingo SB implementation --- .python-version | 1 + dingo/bindings/bindings.cpp | 4 ++++ dingo/bindings/bindings.h | 1 + dingo/volestipy.pyx | 2 +- eigen | 2 +- tests/shake_and_bake_test.py | 39 ++++++++++++++++++++++++++++++++++++ volesti | 2 +- 7 files changed, 48 insertions(+), 3 deletions(-) create mode 100644 .python-version create mode 100644 tests/shake_and_bake_test.py 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/bindings/bindings.cpp b/dingo/bindings/bindings.cpp index 0e3a62f4..a45975ef 100644 --- a/dingo/bindings/bindings.cpp +++ b/dingo/bindings/bindings.cpp @@ -134,6 +134,10 @@ double HPolytopeCPP::apply_sampling(int walk_len, } else if (strcmp(method, "vaidya_walk")) { // 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_walk") == 0) { // shake and bake walk + auto [boundary_pt, facet_idx] = compute_boundary_point(HP, rng, static_cast(1e-4)); + shakeandbake_sampling(rand_points,HP, rng, walk_len,number_of_points, + boundary_pt,number_of_points_to_burn, facet_idx); } else if (strcmp(method, "mmcs")) { // vaidya walk MT S; int total_ess; diff --git a/dingo/bindings/bindings.h b/dingo/bindings/bindings.h index 553ea99e..43f38d96 100644 --- a/dingo/bindings/bindings.h +++ b/dingo/bindings/bindings.h @@ -34,6 +34,7 @@ #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" diff --git a/dingo/volestipy.pyx b/dingo/volestipy.pyx index 6e8d3a98..9c64d003 100644 --- a/dingo/volestipy.pyx +++ b/dingo/volestipy.pyx @@ -85,7 +85,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") ] rounding_methods = ["min_ellipsoid".encode("UTF-8"), "svd".encode("UTF-8"), "max_ellipsoid".encode("UTF-8")] # Build the HPolytope class 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/tests/shake_and_bake_test.py b/tests/shake_and_bake_test.py new file mode 100644 index 00000000..5fc7251d --- /dev/null +++ b/tests/shake_and_bake_test.py @@ -0,0 +1,39 @@ +import numpy as np +from dingo import PolytopeSampler + +# Definišemo kocku [-1,1]^3: +A = np.array([ + [ 1, 0, 0], + [-1, 0, 0], + [ 0, 1, 0], + [ 0,-1, 0], + [ 0, 0, 1], + [ 0, 0,-1] +], dtype=float) +b = np.array([1, 1, 1, 1, 1, 1], dtype=float) + +# Napravi sampler direktno iz (A, b) + +# Generiši uzorke (na primer billiard_walk) +samples = PolytopeSampler.sample_from_polytope_no_multiphase( + A, b, method = 'cdhr', n=1000, burn_in=100, thinning=1, variance=1.0, bias_vector=None, solver=None, ess=0 + ) + +print("Oblik matrice uzoraka:", samples.shape) +print("Min po osi:", samples.min(axis=1)) +print("Max po osi:", samples.max(axis=1)) + +# Ako hoćeš plot +try: + import matplotlib.pyplot as plt + from mpl_toolkits.mplot3d import Axes3D # noqa: F401 + fig = plt.figure() + ax = fig.add_subplot(111, projection="3d") + ax.scatter(samples[0,:], samples[1,:], samples[2,:], s=4, alpha=0.3) + ax.set_xlabel("x") + ax.set_ylabel("y") + ax.set_zlabel("z") + ax.set_title("Samples iz [-1,1]^3 kocke") + plt.show() +except ImportError: + pass diff --git a/volesti b/volesti index c3109bba..a33e5f45 160000 --- a/volesti +++ b/volesti @@ -1 +1 @@ -Subproject commit c3109bba06a9b623446bdde4c5fadb02722de876 +Subproject commit a33e5f450a8badb506384d7a619344dca855fe75 From 3b47c5ef31775eb198665a8aab2eb4d88896245d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Iva=20Jankovi=C4=87?= Date: Mon, 18 Aug 2025 20:01:45 +0200 Subject: [PATCH 02/14] Ignore virtual environment folder --- .gitignore | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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/ From 7c4f4ddd16ec50fd25ee01ed34a7b9a5ab301d15 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Iva=20Jankovi=C4=87?= Date: Mon, 18 Aug 2025 20:02:54 +0200 Subject: [PATCH 03/14] Added .venv in gitignore --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index d8e930b6..cf0beec2 100644 --- a/.gitignore +++ b/.gitignore @@ -13,3 +13,4 @@ venv lp_solve_5.5/ .devcontainer/ .github/dependabot.yml.venv/ +.venv/ From 1712a7a2c6d7f4d4d12d7c0af0883e70f18ee665 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Iva=20Jankovi=C4=87?= Date: Sat, 23 Aug 2025 16:08:43 +0200 Subject: [PATCH 04/14] First version --- dingo/PolytopeSampler.py | 16 ++++++++-------- dingo/bindings/bindings.cpp | 8 ++++++-- dingo/volestipy.pyx | 2 +- 3 files changed, 15 insertions(+), 11 deletions(-) diff --git a/dingo/PolytopeSampler.py b/dingo/PolytopeSampler.py index 27f37853..1d5b21b3 100644 --- a/dingo/PolytopeSampler.py +++ b/dingo/PolytopeSampler.py @@ -53,12 +53,12 @@ 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 ): ( @@ -166,7 +166,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 @@ -223,7 +223,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 diff --git a/dingo/bindings/bindings.cpp b/dingo/bindings/bindings.cpp index a45975ef..072b6553 100644 --- a/dingo/bindings/bindings.cpp +++ b/dingo/bindings/bindings.cpp @@ -134,10 +134,14 @@ double HPolytopeCPP::apply_sampling(int walk_len, } else if (strcmp(method, "vaidya_walk")) { // 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_walk") == 0) { // shake and bake walk - auto [boundary_pt, facet_idx] = compute_boundary_point(HP, rng, static_cast(1e-4)); + } else if (strcmp(method, "shake_and_bake") == 0) { // shake and bake walk + auto [boundary_pt, facet_idx] = compute_boundary_point(HP, rng, static_cast(1e-10)); shakeandbake_sampling(rand_points,HP, rng, walk_len,number_of_points, boundary_pt,number_of_points_to_burn, facet_idx); + } else if (strcmp(method, "billiard_shake_and_bake") == 0) { // billiard shake and bake walk + auto [boundary_pt, facet_idx] = compute_boundary_point(HP, rng, static_cast(1e-10)); + billiard_shakeandbake_sampling(randPoints, P, rng, walkL,nreflections, numpoints, + StartingPoint, nburns, facet_index); } else if (strcmp(method, "mmcs")) { // vaidya walk MT S; int total_ess; diff --git a/dingo/volestipy.pyx b/dingo/volestipy.pyx index 9c64d003..02741956 100644 --- a/dingo/volestipy.pyx +++ b/dingo/volestipy.pyx @@ -85,7 +85,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"),"shake_and_bake".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 From aaa228bbb14e11dd535fa512e10d90c29857af19 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Iva=20Jankovi=C4=87?= Date: Sat, 23 Aug 2025 16:20:12 +0200 Subject: [PATCH 05/14] Remove tests/shake_and_bake_test.py --- tests/shake_and_bake_test.py | 39 ------------------------------------ 1 file changed, 39 deletions(-) delete mode 100644 tests/shake_and_bake_test.py diff --git a/tests/shake_and_bake_test.py b/tests/shake_and_bake_test.py deleted file mode 100644 index 5fc7251d..00000000 --- a/tests/shake_and_bake_test.py +++ /dev/null @@ -1,39 +0,0 @@ -import numpy as np -from dingo import PolytopeSampler - -# Definišemo kocku [-1,1]^3: -A = np.array([ - [ 1, 0, 0], - [-1, 0, 0], - [ 0, 1, 0], - [ 0,-1, 0], - [ 0, 0, 1], - [ 0, 0,-1] -], dtype=float) -b = np.array([1, 1, 1, 1, 1, 1], dtype=float) - -# Napravi sampler direktno iz (A, b) - -# Generiši uzorke (na primer billiard_walk) -samples = PolytopeSampler.sample_from_polytope_no_multiphase( - A, b, method = 'cdhr', n=1000, burn_in=100, thinning=1, variance=1.0, bias_vector=None, solver=None, ess=0 - ) - -print("Oblik matrice uzoraka:", samples.shape) -print("Min po osi:", samples.min(axis=1)) -print("Max po osi:", samples.max(axis=1)) - -# Ako hoćeš plot -try: - import matplotlib.pyplot as plt - from mpl_toolkits.mplot3d import Axes3D # noqa: F401 - fig = plt.figure() - ax = fig.add_subplot(111, projection="3d") - ax.scatter(samples[0,:], samples[1,:], samples[2,:], s=4, alpha=0.3) - ax.set_xlabel("x") - ax.set_ylabel("y") - ax.set_zlabel("z") - ax.set_title("Samples iz [-1,1]^3 kocke") - plt.show() -except ImportError: - pass From d881edab145d02e2f54185a1f458837ccdc9cdb3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Iva=20Jankovi=C4=87?= Date: Sat, 23 Aug 2025 17:18:53 +0200 Subject: [PATCH 06/14] Small fix --- .gitignore | 1 - dingo/PolytopeSampler.py | 12 ++++++------ 2 files changed, 6 insertions(+), 7 deletions(-) diff --git a/.gitignore b/.gitignore index cf0beec2..d8e930b6 100644 --- a/.gitignore +++ b/.gitignore @@ -13,4 +13,3 @@ venv lp_solve_5.5/ .devcontainer/ .github/dependabot.yml.venv/ -.venv/ diff --git a/dingo/PolytopeSampler.py b/dingo/PolytopeSampler.py index 1d5b21b3..012355fb 100644 --- a/dingo/PolytopeSampler.py +++ b/dingo/PolytopeSampler.py @@ -53,12 +53,12 @@ def get_polytope(self): """ if ( - 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 + self._A == [] + or self._b == [] + or self._N == [] + or self._N_shift == [] + or self._T == [] + or self._T_shift == [] ): ( From 6ea096d054c953474f4a1bfe870a4966a6e4653c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Iva=20Jankovi=C4=87?= Date: Sat, 23 Aug 2025 17:27:45 +0200 Subject: [PATCH 07/14] . --- pyproject.toml | 2 +- setup.py | 91 +++++++++++++++++++++----------------------------- 2 files changed, 39 insertions(+), 54 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index d0cdf087..ca7f4735 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -23,5 +23,5 @@ networkx = "3.1" [tool.poetry.dev-dependencies] [build-system] -requires = ["poetry-core>=1.0.0", "cython", "numpy==1.20.1"] +requires = ["poetry-core>=1.0.0", "setuptools", "wheel", "cython", "numpy==1.20.1"] build-backend = "poetry.core.masonry.api" diff --git a/setup.py b/setup.py index 77ec7ec7..77bb0e3f 100755 --- a/setup.py +++ b/setup.py @@ -6,38 +6,29 @@ # Licensed under GNU LGPL.3, see LICENCE file -# This is the setup Python script for building the dingo library - -from distutils.core import setup -from distutils.core import Extension +from setuptools import setup, Extension from Cython.Build import cythonize from os.path import join import numpy -import os -# information about the dingo library version = "0.1.0" -license = ("LGPL3",) packages = ["dingo"] description = "A python library for metabolic networks sampling and analysis" author = "Apostolos Chalkis" author_email = "tolis.chal@gmail.com" name = "dingo" - source_directory_list = ["dingo", join("dingo", "bindings")] compiler_args = ["-std=c++17", "-O3", "-DBOOST_NO_AUTO_PTR", "-ldl", "-lm", "-fopenmp"] -lp_solve_compiler_args = ["-DYY_NEVER_INTERACTIVE", "-DLoadInverseLib=0", "-DLoadLanguageLib=0", -"-DRoleIsExternalInvEngine", "-DINVERSE_ACTIVE=3", "-DLoadableBlasLib=0"] - +lp_solve_compiler_args = [ + "-DYY_NEVER_INTERACTIVE", "-DLoadInverseLib=0", "-DLoadLanguageLib=0", + "-DRoleIsExternalInvEngine", "-DINVERSE_ACTIVE=3", "-DLoadableBlasLib=0" +] link_args = ["-O3", "-fopenmp"] extra_volesti_include_dirs = [ - # include binding files join("dingo", "bindings"), - # the volesti code uses some external classes. - # external directories we need to add join("eigen"), join("boost_1_76_0"), join("boost_1_76_0", "boost"), @@ -49,7 +40,6 @@ join("lp_solve_5.5", "shared"), join("volesti", "external"), join("volesti", "external", "minimum_ellipsoid"), - # include and add the directories on the "include" directory join("volesti", "include"), join("volesti", "include", "convex_bodies"), join("volesti", "include", "random_walks"), @@ -58,50 +48,45 @@ join("volesti", "include", "cartesian_geom"), ] -src_files = ["lp_solve_5.5/bfp/bfp_LUSOL/lp_LUSOL.c" - , "lp_solve_5.5/bfp/bfp_LUSOL/LUSOL/lusol.c" - , "lp_solve_5.5/colamd/colamd.c" - , "lp_solve_5.5/ini.c" - , "lp_solve_5.5/shared/commonlib.c" - , "lp_solve_5.5/shared/mmio.c" - , "lp_solve_5.5/shared/myblas.c" - , "lp_solve_5.5/lp_crash.c" - , "lp_solve_5.5/lp_Hash.c" - , "lp_solve_5.5/lp_lib.c" - , "lp_solve_5.5/lp_matrix.c" - , "lp_solve_5.5/lp_MDO.c" - , "lp_solve_5.5/lp_mipbb.c" - , "lp_solve_5.5/lp_MPS.c" - , "lp_solve_5.5/lp_params.c" - , "lp_solve_5.5/lp_presolve.c" - , "lp_solve_5.5/lp_price.c" - , "lp_solve_5.5/lp_pricePSE.c" - , "lp_solve_5.5/lp_report.c" - , "lp_solve_5.5/lp_scale.c" - , "lp_solve_5.5/lp_simplex.c" - , "lp_solve_5.5/lp_SOS.c" - , "lp_solve_5.5/lp_utils.c" - , "lp_solve_5.5/lp_wlp.c" - , "dingo/volestipy.pyx" - , "dingo/bindings/bindings.cpp"] +src_files = [ + "lp_solve_5.5/bfp/bfp_LUSOL/lp_LUSOL.c", + "lp_solve_5.5/bfp/bfp_LUSOL/LUSOL/lusol.c", + "lp_solve_5.5/colamd/colamd.c", + "lp_solve_5.5/ini.c", + "lp_solve_5.5/shared/commonlib.c", + "lp_solve_5.5/shared/mmio.c", + "lp_solve_5.5/shared/myblas.c", + "lp_solve_5.5/lp_crash.c", + "lp_solve_5.5/lp_Hash.c", + "lp_solve_5.5/lp_lib.c", + "lp_solve_5.5/lp_matrix.c", + "lp_solve_5.5/lp_MDO.c", + "lp_solve_5.5/lp_mipbb.c", + "lp_solve_5.5/lp_MPS.c", + "lp_solve_5.5/lp_params.c", + "lp_solve_5.5/lp_presolve.c", + "lp_solve_5.5/lp_price.c", + "lp_solve_5.5/lp_pricePSE.c", + "lp_solve_5.5/lp_report.c", + "lp_solve_5.5/lp_scale.c", + "lp_solve_5.5/lp_simplex.c", + "lp_solve_5.5/lp_SOS.c", + "lp_solve_5.5/lp_utils.c", + "lp_solve_5.5/lp_wlp.c", + "dingo/volestipy.pyx", + "dingo/bindings/bindings.cpp", +] -# Return the directory that contains the NumPy *.h header files. -# Extension modules that need to compile against NumPy should use this -# function to locate the appropriate include directory. extra_include_dirs = [numpy.get_include()] ext_module = Extension( - "volestipy", + "dingo.volestipy", # Use package-relative name language="c++", sources=src_files, include_dirs=extra_include_dirs + extra_volesti_include_dirs, extra_compile_args=compiler_args + lp_solve_compiler_args, extra_link_args=link_args, ) -print("The Extension function is OK.") - -ext_modules = cythonize([ext_module], gdb_debug=False) -print("The cythonize function ran fine!") setup( version=version, @@ -109,7 +94,7 @@ author_email=author_email, name=name, packages=packages, - ext_modules=ext_modules, -) - -print("Installation of dingo completed.") + ext_modules=cythonize([ext_module], gdb_debug=False), + description=description, + install_requires=["numpy", "cython"], +) \ No newline at end of file From aec145148100ea652360ad5857650a87bc12d16a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Iva=20Jankovi=C4=87?= Date: Sat, 23 Aug 2025 17:43:49 +0200 Subject: [PATCH 08/14] Revert pyproject.toml & setup.py to upstream/develop versions --- pyproject.toml | 2 +- setup.py | 91 +++++++++++++++++++++++++++++--------------------- 2 files changed, 54 insertions(+), 39 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index ca7f4735..d0cdf087 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -23,5 +23,5 @@ networkx = "3.1" [tool.poetry.dev-dependencies] [build-system] -requires = ["poetry-core>=1.0.0", "setuptools", "wheel", "cython", "numpy==1.20.1"] +requires = ["poetry-core>=1.0.0", "cython", "numpy==1.20.1"] build-backend = "poetry.core.masonry.api" diff --git a/setup.py b/setup.py index 77bb0e3f..77ec7ec7 100755 --- a/setup.py +++ b/setup.py @@ -6,29 +6,38 @@ # Licensed under GNU LGPL.3, see LICENCE file -from setuptools import setup, Extension +# This is the setup Python script for building the dingo library + +from distutils.core import setup +from distutils.core import Extension from Cython.Build import cythonize from os.path import join import numpy +import os +# information about the dingo library version = "0.1.0" +license = ("LGPL3",) packages = ["dingo"] description = "A python library for metabolic networks sampling and analysis" author = "Apostolos Chalkis" author_email = "tolis.chal@gmail.com" name = "dingo" + source_directory_list = ["dingo", join("dingo", "bindings")] compiler_args = ["-std=c++17", "-O3", "-DBOOST_NO_AUTO_PTR", "-ldl", "-lm", "-fopenmp"] -lp_solve_compiler_args = [ - "-DYY_NEVER_INTERACTIVE", "-DLoadInverseLib=0", "-DLoadLanguageLib=0", - "-DRoleIsExternalInvEngine", "-DINVERSE_ACTIVE=3", "-DLoadableBlasLib=0" -] +lp_solve_compiler_args = ["-DYY_NEVER_INTERACTIVE", "-DLoadInverseLib=0", "-DLoadLanguageLib=0", +"-DRoleIsExternalInvEngine", "-DINVERSE_ACTIVE=3", "-DLoadableBlasLib=0"] + link_args = ["-O3", "-fopenmp"] extra_volesti_include_dirs = [ + # include binding files join("dingo", "bindings"), + # the volesti code uses some external classes. + # external directories we need to add join("eigen"), join("boost_1_76_0"), join("boost_1_76_0", "boost"), @@ -40,6 +49,7 @@ join("lp_solve_5.5", "shared"), join("volesti", "external"), join("volesti", "external", "minimum_ellipsoid"), + # include and add the directories on the "include" directory join("volesti", "include"), join("volesti", "include", "convex_bodies"), join("volesti", "include", "random_walks"), @@ -48,45 +58,50 @@ join("volesti", "include", "cartesian_geom"), ] -src_files = [ - "lp_solve_5.5/bfp/bfp_LUSOL/lp_LUSOL.c", - "lp_solve_5.5/bfp/bfp_LUSOL/LUSOL/lusol.c", - "lp_solve_5.5/colamd/colamd.c", - "lp_solve_5.5/ini.c", - "lp_solve_5.5/shared/commonlib.c", - "lp_solve_5.5/shared/mmio.c", - "lp_solve_5.5/shared/myblas.c", - "lp_solve_5.5/lp_crash.c", - "lp_solve_5.5/lp_Hash.c", - "lp_solve_5.5/lp_lib.c", - "lp_solve_5.5/lp_matrix.c", - "lp_solve_5.5/lp_MDO.c", - "lp_solve_5.5/lp_mipbb.c", - "lp_solve_5.5/lp_MPS.c", - "lp_solve_5.5/lp_params.c", - "lp_solve_5.5/lp_presolve.c", - "lp_solve_5.5/lp_price.c", - "lp_solve_5.5/lp_pricePSE.c", - "lp_solve_5.5/lp_report.c", - "lp_solve_5.5/lp_scale.c", - "lp_solve_5.5/lp_simplex.c", - "lp_solve_5.5/lp_SOS.c", - "lp_solve_5.5/lp_utils.c", - "lp_solve_5.5/lp_wlp.c", - "dingo/volestipy.pyx", - "dingo/bindings/bindings.cpp", -] +src_files = ["lp_solve_5.5/bfp/bfp_LUSOL/lp_LUSOL.c" + , "lp_solve_5.5/bfp/bfp_LUSOL/LUSOL/lusol.c" + , "lp_solve_5.5/colamd/colamd.c" + , "lp_solve_5.5/ini.c" + , "lp_solve_5.5/shared/commonlib.c" + , "lp_solve_5.5/shared/mmio.c" + , "lp_solve_5.5/shared/myblas.c" + , "lp_solve_5.5/lp_crash.c" + , "lp_solve_5.5/lp_Hash.c" + , "lp_solve_5.5/lp_lib.c" + , "lp_solve_5.5/lp_matrix.c" + , "lp_solve_5.5/lp_MDO.c" + , "lp_solve_5.5/lp_mipbb.c" + , "lp_solve_5.5/lp_MPS.c" + , "lp_solve_5.5/lp_params.c" + , "lp_solve_5.5/lp_presolve.c" + , "lp_solve_5.5/lp_price.c" + , "lp_solve_5.5/lp_pricePSE.c" + , "lp_solve_5.5/lp_report.c" + , "lp_solve_5.5/lp_scale.c" + , "lp_solve_5.5/lp_simplex.c" + , "lp_solve_5.5/lp_SOS.c" + , "lp_solve_5.5/lp_utils.c" + , "lp_solve_5.5/lp_wlp.c" + , "dingo/volestipy.pyx" + , "dingo/bindings/bindings.cpp"] +# Return the directory that contains the NumPy *.h header files. +# Extension modules that need to compile against NumPy should use this +# function to locate the appropriate include directory. extra_include_dirs = [numpy.get_include()] ext_module = Extension( - "dingo.volestipy", # Use package-relative name + "volestipy", language="c++", sources=src_files, include_dirs=extra_include_dirs + extra_volesti_include_dirs, extra_compile_args=compiler_args + lp_solve_compiler_args, extra_link_args=link_args, ) +print("The Extension function is OK.") + +ext_modules = cythonize([ext_module], gdb_debug=False) +print("The cythonize function ran fine!") setup( version=version, @@ -94,7 +109,7 @@ author_email=author_email, name=name, packages=packages, - ext_modules=cythonize([ext_module], gdb_debug=False), - description=description, - install_requires=["numpy", "cython"], -) \ No newline at end of file + ext_modules=ext_modules, +) + +print("Installation of dingo completed.") From f7f5813df9f607107262e5dea378591e80d61340 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Iva=20Jankovi=C4=87?= Date: Fri, 31 Oct 2025 01:53:05 +0000 Subject: [PATCH 09/14] Pipeline for Shake and Bake --- dingo/PolytopeSampler.py | 92 ++++++++++++++++++++++++++++----- dingo/bindings/bindings.cpp | 67 ++++++++++++++++++++---- dingo/bindings/bindings.h | 23 ++++++++- dingo/volestipy.pyx | 32 ++++++++++++ pyproject.toml | 2 +- setup.py | 0 tests/sampling_no_multiphase.py | 2 +- 7 files changed, 192 insertions(+), 26 deletions(-) mode change 100755 => 100644 setup.py diff --git a/dingo/PolytopeSampler.py b/dingo/PolytopeSampler.py index 012355fb..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, @@ -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 @@ -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 9e6521ae..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,14 +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) { // shake and bake walk - auto [boundary_pt, facet_idx] = compute_boundary_point(HP, rng, static_cast(1e-10)); - shakeandbake_sampling(rand_points,HP, rng, walk_len,number_of_points, - boundary_pt,number_of_points_to_burn, facet_idx); - } else if (strcmp(method, "billiard_shake_and_bake") == 0) { // billiard shake and bake walk - auto [boundary_pt, facet_idx] = compute_boundary_point(HP, rng, static_cast(1e-10)); - billiard_shakeandbake_sampling(randPoints, P, rng, walkL,nreflections, numpoints, - StartingPoint, nburns, facet_index); + } 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; @@ -175,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) { @@ -436,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 43f38d96..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,9 +30,11 @@ #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" @@ -40,6 +43,7 @@ #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; @@ -49,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 @@ -156,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 02741956..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: @@ -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/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 From c7fff84de909e48431133a0a3bb38b7ebdd4b8d7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Iva=20Jankovi=C4=87?= Date: Mon, 3 Nov 2025 13:25:38 +0000 Subject: [PATCH 10/14] Billiard Shake and Bake --- dingo/PolytopeSampler.py | 125 ++++++++++++++++++++++++++++-------- dingo/bindings/bindings.cpp | 85 ++++++++++++++++++++++-- dingo/bindings/bindings.h | 3 +- dingo/volestipy.pyx | 55 +++++++++++++++- tests/shake_and_bake.py | 100 +++++++++++++++++++++++++++++ 5 files changed, 332 insertions(+), 36 deletions(-) create mode 100644 tests/shake_and_bake.py diff --git a/dingo/PolytopeSampler.py b/dingo/PolytopeSampler.py index a683ef92..3e32868f 100644 --- a/dingo/PolytopeSampler.py +++ b/dingo/PolytopeSampler.py @@ -47,6 +47,7 @@ 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 @@ -168,7 +169,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', 'shake_and_bake', 'billiard_shake_and_bake'} + 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 @@ -192,40 +193,107 @@ 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): + def generate_steady_states_sb_once(self, n=1000, burn_in=0, sampler="sb", nreflections=None): """ - One Shake and Bake phase: samples n points, returns (steady_states, diagnostics) - diagnostics = {'minESS':..., 'maxPSRF':..., 'N':..., 'phases': 1, 'seconds': ...}. + Single-phase boundary sampler. + - sampler: "sb"/"shake_and_bake" or "bsb"/"billiard_shake_and_bake" + - nreflections: only for BSB; defaults to ceil(sqrt(d)) if None """ + import numpy as np + + # Build H-polytope from current A, b 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) + # Normalize sampler keyword + s = (sampler or "sb").lower() + if s in ("sb", "shake_and_bake"): + method = b"shake_and_bake" + use_bsb = False + elif s in ("bsb", "billiard_shake_and_bake"): + method = b"billiard_shake_and_bake" + use_bsb = True 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]}") + raise ValueError(f"Unknown sampler '{sampler}'") - 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 + d = int(self._A.shape[1]) - P.get_sb_diagnostics(np.ascontiguousarray(diag_buf)) - minESS, maxPSRF, Ncpp, phases, seconds = diag_buf + # IMPORTANT: walk_len must be modest; too large often causes poor mixing or NaNs. + # A robust default is O(sqrt(d)) with a small lower bound. + walk_len = max(5, int(np.sqrt(d))) + + # bias_vector must be a 1D array (never None) + bias_vec = np.ones(d, dtype=np.float64) + + # Reflections for BSB (if not provided, use ceil(sqrt(d))) + if use_bsb: + nref = int(nreflections) if nreflections is not None else int(np.ceil(np.sqrt(d))) + else: + nref = 0 + + # Call into Cython in the exact argument order it expects: + # (method, number_of_points, number_of_points_to_burn, walk_len, + # variance_value, bias_vector, solver=None, ess=1000, nreflections) + samples = P.generate_samples( + method, + int(n), + int(burn_in), + int(walk_len), + 1.0, + bias_vec, + self._parameters.get("solver"), + 0, + int(nref), + ) - steady_states = map_samples_to_steady_states(samples.T, self._N, self._N_shift) + # Pull diagnostics from C++ (could contain NaN if samples had NaNs) + self._last_hpoly = P + diag_buf = np.empty(5, dtype=np.float64) + P.get_sb_diagnostics(diag_buf) + minESS, maxPSRF, Ncpp_reported, phases, seconds = diag_buf + + # Cython returns samples as (N, d); convert to (d, N) + S = samples.T # shape (d, N) + + # Sanitize columns: drop any column with NaN/Inf to avoid propagating NaNs + finite_cols = np.isfinite(S).all(axis=0) + if not finite_cols.all(): + S = S[:, finite_cols] + + N_kept = S.shape[1] + if N_kept == 0: + raise RuntimeError( + "All generated samples are NaN/Inf; check walk_len, inner ball, and constraints." + ) + + # Map to full steady-state space + steady_states = map_samples_to_steady_states(S, self._N, self._N_shift) + + # Diagnostics: report actual kept N; guard PSRF if NaN diagnostics = { "minESS": float(minESS), - "maxPSRF": float(maxPSRF), - "N": int(Ncpp), + "maxPSRF": float(maxPSRF) if np.isfinite(maxPSRF) else 1.0, + "N": int(N_kept), "phases": int(phases), "seconds": float(seconds) if np.isfinite(seconds) else None, } + return steady_states, diagnostics + def sb_scaling_ratio(self, tol=1e-10, min_ratio=0.01): + """ + Prosleđuje na C++ get_sb_scaling_ratio preko POSLEDNJEG HPolytope + koji je generisao SB/BSB uzorke. Moraš prethodno pozvati + generate_steady_states_sb_once(...). + """ + if self._last_hpoly is None: + raise RuntimeError("sb_scaling_ratio: nema keširanog HPolytope sa SB/BSB uzorcima; pozovi sampler prvo.") + # ovo zove Cython metodu koja pakuje numpy bafer-e i zove C++ get_sb_scaling_ratio(...) + return self._last_hpoly.sb_scaling_ratio(tol=tol, min_ratio=min_ratio) + + + @staticmethod def sample_from_polytope( A, b, ess=1000, psrf=False, parallel_mmcs=False, num_threads=1, solver=None @@ -278,20 +346,24 @@ def sample_from_polytope_no_multiphase( @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 + A, b, n=1000, burn_in=0, solver=None, sampler="sb" ): """ - One Shake and Bake phase for polytope defined by A,b. + One boundary-sampling phase for polytope A,b using SB or BSB. 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) + # choose method inline, no extra vars + method = b"billiard_shake_and_bake" if sampler == "bsb" else b"shake_and_bake" + + # call binding with or without solver depending on signature + try: + samples = P.generate_samples(method, int(n), int(burn_in), solver) + except TypeError: + samples = P.generate_samples(method, int(n), int(burn_in)) + + # shared diagnostics for SB/BSB diag = np.zeros(5, dtype=np.float64) P.get_sb_diagnostics(diag) minESS, maxPSRF, Ncpp, phases, seconds = diag @@ -306,6 +378,7 @@ def sample_from_polytope_sb_once( return samples.T, diagnostics + @staticmethod def round_polytope( A, b, method = "john_position", solver = None diff --git a/dingo/bindings/bindings.cpp b/dingo/bindings/bindings.cpp index 0fd19cf4..ae287e59 100644 --- a/dingo/bindings/bindings.cpp +++ b/dingo/bindings/bindings.cpp @@ -16,6 +16,8 @@ #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; @@ -93,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(); @@ -142,11 +145,9 @@ double HPolytopeCPP::apply_sampling(int walk_len, 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); + shakeandbake_sampling(batch, HP, rng, walk_len, number_of_points, boundary_pt, number_of_points_to_burn, facet_idx); rand_points.swap(batch); const int N = static_cast(rand_points.size()); MT S(d, N); @@ -164,13 +165,44 @@ double HPolytopeCPP::apply_sampling(int walk_len, auto t1 = clock::now(); double secs = std::chrono::duration(t1 - t0).count(); - sb_samples_ = S; // d x N + 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_.N = N; sb_diag_.phases = 1; sb_diag_.seconds = secs; + } else if (strcmp(method, "billiard_shake_and_bake") == 0) { // billiard shake and bake walk + 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(); + std::list batch; + billiard_shakeandbake_sampling(batch, HP, rng, walk_len, nreflections, number_of_points, + boundary_pt, number_of_points_to_burn, 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; @@ -439,6 +471,47 @@ double HPolytopeCPP::mmcs_step(double* inner_point, double radius, int &N) { return 0.0; } +void HPolytopeCPP::get_sb_scaling_ratio(double tol, + double min_ratio, + double* scale_out, + double* coverage_out, + double* maxdev_out, + double* avgdev_out) const +{ + const int d = sb_samples_.rows(); + const int N = sb_samples_.cols(); + if (d == 0 || N == 0) { + throw std::runtime_error("get_sb_scaling_ratio: no SB/BSB samples available. " + "Call a SB/BSB sampler first."); + } + + auto result = scaling_ratio_boundary_test(HP, sb_samples_, + static_cast(tol), + static_cast(min_ratio)); + + const auto& scale = std::get<0>(result); // VT (K) + const auto& coverage = std::get<1>(result); // MT (m x K) + const auto& max_dev = std::get<2>(result); // VT (m) + const auto& avg_dev = std::get<3>(result); // VT (m) + + const int K = static_cast(scale.size()); + const int m = coverage.rows(); + const int Kcov = coverage.cols(); + + for (int k = 0; k < K; ++k) scale_out[k] = static_cast(scale[k]); + + for (int i = 0; i < m; ++i) + for (int k = 0; k < Kcov; ++k) + coverage_out[i * Kcov + k] = static_cast(coverage(i, k)); + + for (int i = 0; i < m; ++i) { + maxdev_out[i] = static_cast(max_dev[i]); + avgdev_out[i] = static_cast(avg_dev[i]); + } +} + + + void HPolytopeCPP::get_mmcs_samples(double* T_matrix, double* T_shift, double* samples) { int n_variables = HP.dimension(); diff --git a/dingo/bindings/bindings.h b/dingo/bindings/bindings.h index 5bfe3ac5..6bb4f4b0 100644 --- a/dingo/bindings/bindings.h +++ b/dingo/bindings/bindings.h @@ -153,7 +153,7 @@ class HPolytopeCPP{ // the apply_sampling() function double apply_sampling(int walk_len, int number_of_points, int number_of_points_to_burn, char* method, double* inner_point, double radius, double* samples, - double variance_value, double* bias_vector, int ess); + double variance_value, double* bias_vector, int ess, int nreflections); void mmcs_initialize(int d, int ess, bool psrf_check, bool parallelism, int num_threads); @@ -172,6 +172,7 @@ class HPolytopeCPP{ inline const MT& sb_samples() const { return sb_samples_; } void get_sb_samples(double* samples) const; void get_sb_diagnostics(double* out5) const; + void get_sb_scaling_ratio(double tol, double min_ratio,double* scale_out,double* coverage_out,double* maxdev_out,double* avgdev_out) const; private: Diagnostics sb_diag_; diff --git a/dingo/volestipy.pyx b/dingo/volestipy.pyx index ddade0d0..51f74677 100644 --- a/dingo/volestipy.pyx +++ b/dingo/volestipy.pyx @@ -57,7 +57,7 @@ cdef extern from "bindings.h": # Random sampling double apply_sampling(int walk_len, int number_of_points, int number_of_points_to_burn, \ char* method, double* inner_point, double radius, double* samples, \ - double variance_value, double* bias_vector, int ess) + double variance_value, double* bias_vector, int ess,int nreflections) # Initialize the parameters for the (m)ultiphase (m)onte (c)arlo (s)ampling algorithm void mmcs_initialize(unsigned int d, int ess, int psrf_check, int parallelism, int num_threads); @@ -76,6 +76,12 @@ cdef extern from "bindings.h": void get_sb_diagnostics(double* out5) const void get_sb_samples(double* samples) const + void get_sb_scaling_ratio(double tol, double min_ratio, + double* scale_out, + double* coverage_out, + double* maxdev_out, + double* avgdev_out) const + # The lowDimPolytopeCPP class along with its functions cdef cppclass lowDimHPolytopeCPP: @@ -124,7 +130,7 @@ cdef class HPolytope: # Likewise, the generate_samples() function def generate_samples(self, method, number_of_points, number_of_points_to_burn, walk_len, - variance_value, bias_vector, solver = None, ess = 1000): + variance_value, bias_vector, solver = None, ess = 1000, nreflections=None): n_variables = self._A.shape[1] cdef double[:,::1] samples = np.zeros((number_of_points, n_variables), dtype = np.float64, order = "C") @@ -138,7 +144,7 @@ cdef class HPolytope: self.polytope_cpp.apply_sampling(walk_len, number_of_points, number_of_points_to_burn, \ method, &inner_point_for_c[0], radius, &samples[0,0], \ - variance_value, &bias_vector_[0], ess) + variance_value, &bias_vector_[0], ess, nreflections) return np.asarray(samples) def get_sb_diagnostics(self, out): @@ -168,6 +174,49 @@ cdef class HPolytope: self.polytope_cpp.get_sb_samples(&S[0,0]) return np.asarray(S) + def sb_scaling_ratio(self, tol=1e-10, min_ratio=0.01): + """ + Vraća: (scale, coverage, max_dev, avg_dev) + - scale: (K,) [tipično K=10] + - coverage: (m,K) [m = broj hiper-ravnina = self._A.shape[0]] + - max_dev: (m,) + - avg_dev: (m,) + Potrebno: prethodno pozvati SB/BSB da bi postojali sb_samples_ u C++. + """ + cdef int m = self._A.shape[0] + # Trenutna C++ implementacija pravi tačno 10 skala; ostavi 10 ovde. + cdef int K = 10 + + # Alokacije kao C-kontigvne matrice + cdef np.ndarray[np.float64_t, ndim=1, mode="c"] scale_np = np.empty((K,), dtype=np.float64, order="C") + cdef np.ndarray[np.float64_t, ndim=2, mode="c"] coverage_np = np.empty((m, K), dtype=np.float64, order="C") + cdef np.ndarray[np.float64_t, ndim=1, mode="c"] maxdev_np = np.empty((m,), dtype=np.float64, order="C") + cdef np.ndarray[np.float64_t, ndim=1, mode="c"] avgdev_np = np.empty((m,), dtype=np.float64, order="C") + + # Typed memoryview-ovi nad C-kontigvnim baferima + cdef double[::1] scale_mv = scale_np + cdef double[:,::1] coverage_mv = coverage_np + cdef double[::1] maxdev_mv = maxdev_np + cdef double[::1] avgdev_mv = avgdev_np + + # Poziv C++ wrappera – prosleđuju se sirove adrese bafera + try: + self.polytope_cpp.get_sb_scaling_ratio( + tol, + min_ratio, + &scale_mv[0], + &coverage_mv[0, 0], # (m,K) row-major: i*K + k + &maxdev_mv[0], + &avgdev_mv[0] + ) + except RuntimeError as e: + # C++ baca kada nema SB/BSB uzoraka; prebaci poruku dalje. + raise RuntimeError(str(e)) + + # Vrati NumPy objekte (već su odgovarajućih dimenzija) + return scale_np, coverage_np, maxdev_np, avgdev_np + + # 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/tests/shake_and_bake.py b/tests/shake_and_bake.py new file mode 100644 index 00000000..190bc9e4 --- /dev/null +++ b/tests/shake_and_bake.py @@ -0,0 +1,100 @@ +# dingo : a python library for metabolic networks sampling and analysis +# dingo is part of the GeomScale project + +import unittest +import os +import sys +import time +import warnings +import gc +import numpy as np + +from dingo import MetabolicNetwork, PolytopeSampler +from dingo.pyoptinterface_based_impl import set_default_solver + +# --- environment safety --- +os.environ["OMP_NUM_THREADS"] = "1" +os.environ["OPENBLAS_NUM_THREADS"] = "1" +os.environ["MKL_NUM_THREADS"] = "1" + +# === Global parameters === +N_SAMPLES = 50000 +BURN_IN = 200 +NREFL = None # will default to ceil(sqrt(d)) inside generate_steady_states_sb_once + +def _print_diag(diag: dict, elapsed_s: float, tag: str): + secs_cpp = diag.get("seconds", None) + secs_cpp_str = f"{secs_cpp:.6f}s" if isinstance(secs_cpp, float) else "None" + print( + f"\n[{tag}] Diagnostics:\n" + f" minESS = {diag.get('minESS')}\n" + f" maxPSRF = {diag.get('maxPSRF')}\n" + f" N = {diag.get('N')}\n" + f" phases = {diag.get('phases')}\n" + f" seconds(C++) = {secs_cpp_str}\n" + f" elapsed(py) = {elapsed_s:.6f}s\n" + ) + +def _print_samples_summary(steady_states: np.ndarray, tag: str, n_dims_preview: int = 3): + ss_min = np.nanmin(steady_states) + ss_max = np.nanmax(steady_states) + print(f"[{tag}] Steady states summary:") + print(f" shape = {steady_states.shape}, global min = {ss_min:.6g}, global max = {ss_max:.6g}") + d, N = steady_states.shape + dims = min(n_dims_preview, d) + for i in range(dims): + mean_i = np.nanmean(steady_states[i, :]) + std_i = np.nanstd(steady_states[i, :]) + print(f" dim {i}: mean = {mean_i:.6g}, std = {std_i:.6g}") + print("") + +def _run_bsb_once(model: MetabolicNetwork, tc: unittest.TestCase, tag: str): + """Runs one Billiard Shake-and-Bake phase""" + sampler = PolytopeSampler(model) + warnings.filterwarnings("ignore", category=DeprecationWarning) + + A, b, N, N_shift = sampler.get_polytope() + d_eff = A.shape[1] + n_rxns = N.shape[0] + tc.assertGreater(d_eff, 0, "Effective dimension is zero") + + t0 = time.perf_counter() + steady_states, diag = sampler.generate_steady_states_sb_once( + n=N_SAMPLES, burn_in=BURN_IN, sampler="billiard_shake_and_bake", nreflections=NREFL + ) + elapsed = time.perf_counter() - t0 + + _print_diag(diag, elapsed, f"{tag} :: BILLIARD_SHAKE_AND_BAKE") + _print_samples_summary(steady_states, f"{tag} :: BILLIARD_SHAKE_AND_BAKE") + + + + # sanity + tc.assertEqual(steady_states.shape[0], n_rxns) + tc.assertTrue(np.isfinite(diag["minESS"])) + tc.assertGreater(diag["minESS"], 0) + + del steady_states, sampler, A, b, N, N_shift + gc.collect() + +class TestBilliardShakeAndBake(unittest.TestCase): + def test_ecoli_core_json(self): + input_file = os.path.join(os.getcwd(), "ext_data", "e_coli_core.json") + model = MetabolicNetwork.from_json(input_file) + _run_bsb_once(model, self, tag="e_coli_core.json") + + def test_ecoli_core_mat(self): + input_file = os.path.join(os.getcwd(), "ext_data", "e_coli_core.mat") + model = MetabolicNetwork.from_mat(input_file) + _run_bsb_once(model, self, tag="e_coli_core.mat") + + def test_ecoli_core_sbml(self): + input_file = os.path.join(os.getcwd(), "ext_data", "e_coli_core.xml") + model = MetabolicNetwork.from_sbml(input_file) + _run_bsb_once(model, self, tag="e_coli_core.xml") + +if __name__ == "__main__": + if len(sys.argv) > 1: + set_default_solver(sys.argv[1]) + sys.argv.pop(1) + unittest.main() From 4513966a87cf54d9caec09e2fa21c0bac231599b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Iva=20Jankovi=C4=87?= Date: Tue, 4 Nov 2025 13:00:38 +0000 Subject: [PATCH 11/14] Example of running --- tests/run_sb_exp.py | 50 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 50 insertions(+) create mode 100644 tests/run_sb_exp.py diff --git a/tests/run_sb_exp.py b/tests/run_sb_exp.py new file mode 100644 index 00000000..7463c2a2 --- /dev/null +++ b/tests/run_sb_exp.py @@ -0,0 +1,50 @@ +#!/usr/bin/python3 +import os, sys +import numpy as np +from dingo import PolytopeSampler +from time import process_time +import pickle + + +def sample_on_polyround_processed_polytope(p): + name = os.path.basename(p) + + with open(p, "rb") as f: + obj = pickle.load(f) + polytope = obj[0] + + polyround_A = polytope.A.to_numpy() + polyround_b = polytope.b.to_numpy() + + print(f"Dimensions = {polyround_A.shape[1]}, Ograničenja = {polyround_A.shape[0]}") + start = process_time() + + steady_states, diag = PolytopeSampler.sample_from_polytope_sb_once( + polyround_A, + polyround_b, + n=15000, + burn_in=1000, + sampler="sb", + walk_len=48, + ) + + end = process_time() + + print( + f"[{name}] minESS={diag['minESS']:.3f} maxPSRF={diag['maxPSRF']:.3f} " + f"N={diag['N']} phases={diag['phases']} seconds={diag['seconds']}" + ) + print(f"Total sampling time: {end - start:.2f} s") + + out_path = f"dingo_polyround_no_multiphase_{name}.pckl" + with open(out_path, "wb") as f_out: + pickle.dump({"samples": steady_states, "diagnostics": diag}, f_out) + + print(f"Saved results to {out_path}") + + +if __name__ == "__main__": + if len(sys.argv) < 2: + sys.exit("Usage: python run_bwr_exp.py ") + file_name = sys.argv[1] + sample_on_polyround_processed_polytope(file_name) From 5335e477efa9d624b997ad2bb892cfc7ac7304d8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Iva=20Jankovi=C4=87?= Date: Sun, 9 Nov 2025 17:26:52 +0000 Subject: [PATCH 12/14] Cleared shake and bake wrappers --- dingo/PolytopeSampler.py | 175 ++++++++++++++------------------- dingo/bindings/bindings.cpp | 191 ++++++++++++++---------------------- dingo/bindings/bindings.h | 43 +++++--- dingo/volestipy.pyx | 130 ++++++++++++------------ 4 files changed, 239 insertions(+), 300 deletions(-) diff --git a/dingo/PolytopeSampler.py b/dingo/PolytopeSampler.py index 3e32868f..88f673e4 100644 --- a/dingo/PolytopeSampler.py +++ b/dingo/PolytopeSampler.py @@ -10,7 +10,9 @@ 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, @@ -193,104 +195,47 @@ 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): + def generate_steady_states_sb_once(self, n=1000, burn_in=0, sampler="sb", nreflections=None, walk_len=None): """ - Single-phase boundary sampler. + Single-phase boundary sampler with clear separation of concerns. - sampler: "sb"/"shake_and_bake" or "bsb"/"billiard_shake_and_bake" - - nreflections: only for BSB; defaults to ceil(sqrt(d)) if None + - walk_len: default ceil(sqrt(d)) with a minimum of 5 + - nreflections: only used for BSB; defaults to ceil(sqrt(d)) """ - import numpy as np - # Build H-polytope from current A, b self.get_polytope() P = HPolytope(self._A, self._b) - - # Normalize sampler keyword - s = (sampler or "sb").lower() - if s in ("sb", "shake_and_bake"): - method = b"shake_and_bake" - use_bsb = False - elif s in ("bsb", "billiard_shake_and_bake"): - method = b"billiard_shake_and_bake" - use_bsb = True - else: - raise ValueError(f"Unknown sampler '{sampler}'") - d = int(self._A.shape[1]) - # IMPORTANT: walk_len must be modest; too large often causes poor mixing or NaNs. - # A robust default is O(sqrt(d)) with a small lower bound. - walk_len = max(5, int(np.sqrt(d))) - - # bias_vector must be a 1D array (never None) - bias_vec = np.ones(d, dtype=np.float64) - - # Reflections for BSB (if not provided, use ceil(sqrt(d))) - if use_bsb: - nref = int(nreflections) if nreflections is not None else int(np.ceil(np.sqrt(d))) - else: - nref = 0 - - # Call into Cython in the exact argument order it expects: - # (method, number_of_points, number_of_points_to_burn, walk_len, - # variance_value, bias_vector, solver=None, ess=1000, nreflections) - samples = P.generate_samples( - method, - int(n), - int(burn_in), - int(walk_len), - 1.0, - bias_vec, - self._parameters.get("solver"), - 0, - int(nref), + 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) ) - # Pull diagnostics from C++ (could contain NaN if samples had NaNs) - self._last_hpoly = P - - diag_buf = np.empty(5, dtype=np.float64) - P.get_sb_diagnostics(diag_buf) - minESS, maxPSRF, Ncpp_reported, phases, seconds = diag_buf - - # Cython returns samples as (N, d); convert to (d, N) - S = samples.T # shape (d, N) + 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), + ) - # Sanitize columns: drop any column with NaN/Inf to avoid propagating NaNs finite_cols = np.isfinite(S).all(axis=0) if not finite_cols.all(): S = S[:, finite_cols] - - N_kept = S.shape[1] - if N_kept == 0: + if S.shape[1] == 0: raise RuntimeError( - "All generated samples are NaN/Inf; check walk_len, inner ball, and constraints." + "All sample columns contain NaN/Inf; check walk_len, nreflections, or polytope constraints." ) - # Map to full steady-state space steady_states = map_samples_to_steady_states(S, self._N, self._N_shift) - # Diagnostics: report actual kept N; guard PSRF if NaN - diagnostics = { - "minESS": float(minESS), - "maxPSRF": float(maxPSRF) if np.isfinite(maxPSRF) else 1.0, - "N": int(N_kept), - "phases": int(phases), - "seconds": float(seconds) if np.isfinite(seconds) else None, - } - - return steady_states, diagnostics + self._last_hpoly = P # cache if needed later + return steady_states - def sb_scaling_ratio(self, tol=1e-10, min_ratio=0.01): - """ - Prosleđuje na C++ get_sb_scaling_ratio preko POSLEDNJEG HPolytope - koji je generisao SB/BSB uzorke. Moraš prethodno pozvati - generate_steady_states_sb_once(...). - """ - if self._last_hpoly is None: - raise RuntimeError("sb_scaling_ratio: nema keširanog HPolytope sa SB/BSB uzorcima; pozovi sampler prvo.") - # ovo zove Cython metodu koja pakuje numpy bafer-e i zove C++ get_sb_scaling_ratio(...) - return self._last_hpoly.sb_scaling_ratio(tol=tol, min_ratio=min_ratio) @@ -346,36 +291,60 @@ def sample_from_polytope_no_multiphase( @staticmethod def sample_from_polytope_sb_once( - A, b, n=1000, burn_in=0, solver=None, sampler="sb" + A, + b, + n: int = 1000, + burn_in: int = 0, + sampler: str = "sb", + walk_len: Optional[int] = None, + nreflections: Optional[int] = None, ): """ - One boundary-sampling phase for polytope A,b using SB or BSB. - Returns (samples_T_d_by_N, diagnostics_dict). + 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), + ) - # choose method inline, no extra vars - method = b"billiard_shake_and_bake" if sampler == "bsb" else b"shake_and_bake" - - # call binding with or without solver depending on signature - try: - samples = P.generate_samples(method, int(n), int(burn_in), solver) - except TypeError: - samples = P.generate_samples(method, int(n), int(burn_in)) - - # shared diagnostics for SB/BSB - 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 + 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 diff --git a/dingo/bindings/bindings.cpp b/dingo/bindings/bindings.cpp index ae287e59..af6bf810 100644 --- a/dingo/bindings/bindings.cpp +++ b/dingo/bindings/bindings.cpp @@ -138,72 +138,6 @@ 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(); - - std::list batch; - shakeandbake_sampling(batch, HP, rng, walk_len, number_of_points, boundary_pt, number_of_points_to_burn, 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, "billiard_shake_and_bake") == 0) { // billiard shake and bake walk - 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(); - std::list batch; - billiard_shakeandbake_sampling(batch, HP, rng, walk_len, nreflections, number_of_points, - boundary_pt, number_of_points_to_burn, 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; @@ -251,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(HP, rng, static_cast(1e-8)); + + const int d = HP.dimension(); + std::list batch; + + if (strcmp(sampler, "shake_and_bake") == 0 || strcmp(sampler, "sb") == 0) { + shakeandbake_sampling(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(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(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(S, min_ess_u); + VT rhat_vec = univariate_psrf(S); + + SBDiagnostics out; + out.minESS = static_cast(ess_vec.minCoeff()); + out.maxPSRF = static_cast(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 { @@ -471,47 +465,6 @@ double HPolytopeCPP::mmcs_step(double* inner_point, double radius, int &N) { return 0.0; } -void HPolytopeCPP::get_sb_scaling_ratio(double tol, - double min_ratio, - double* scale_out, - double* coverage_out, - double* maxdev_out, - double* avgdev_out) const -{ - const int d = sb_samples_.rows(); - const int N = sb_samples_.cols(); - if (d == 0 || N == 0) { - throw std::runtime_error("get_sb_scaling_ratio: no SB/BSB samples available. " - "Call a SB/BSB sampler first."); - } - - auto result = scaling_ratio_boundary_test(HP, sb_samples_, - static_cast(tol), - static_cast(min_ratio)); - - const auto& scale = std::get<0>(result); // VT (K) - const auto& coverage = std::get<1>(result); // MT (m x K) - const auto& max_dev = std::get<2>(result); // VT (m) - const auto& avg_dev = std::get<3>(result); // VT (m) - - const int K = static_cast(scale.size()); - const int m = coverage.rows(); - const int Kcov = coverage.cols(); - - for (int k = 0; k < K; ++k) scale_out[k] = static_cast(scale[k]); - - for (int i = 0; i < m; ++i) - for (int k = 0; k < Kcov; ++k) - coverage_out[i * Kcov + k] = static_cast(coverage(i, k)); - - for (int i = 0; i < m; ++i) { - maxdev_out[i] = static_cast(max_dev[i]); - avgdev_out[i] = static_cast(avg_dev[i]); - } -} - - - void HPolytopeCPP::get_mmcs_samples(double* T_matrix, double* T_shift, double* samples) { int n_variables = HP.dimension(); @@ -539,21 +492,19 @@ 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_samples(double* out) const { + const int d = static_cast(sb_samples_.rows()); + const int N = static_cast(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* 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; +void HPolytopeCPP::get_sb_diagnostics(double* out3) const { + out3[0] = sb_diag_.minESS; + out3[1] = sb_diag_.maxPSRF; + out3[2] = static_cast(sb_diag_.N); + } diff --git a/dingo/bindings/bindings.h b/dingo/bindings/bindings.h index 6bb4f4b0..0bdb058a 100644 --- a/dingo/bindings/bindings.h +++ b/dingo/bindings/bindings.h @@ -53,12 +53,10 @@ 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; +struct SBDiagnostics { + double minESS; + double maxPSRF; + long long N; }; template @@ -167,15 +165,34 @@ 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; - void get_sb_scaling_ratio(double tol, double min_ratio,double* scale_out,double* coverage_out,double* maxdev_out,double* avgdev_out) const; + // the boundary sampling function: shake and bake, billiard shake and bake + int apply_boundary_sampling(int walk_len,int number_of_points,int number_of_points_to_burn, const char* sampler, int nreflections, double* samples); + + // Compute diagnostics (minESS, maxPSRF, N, phases, seconds) for a given sample buffer. + // This is a static utility function and does not depend on the internal state. + static SBDiagnostics sb_diagnostics(int d, int N, const double* samples); + + // Return a reference to the last stored diagnostics object (legacy internal storage). + // Used for backward compatibility; modern code should call get_sb_diagnostics(). + inline const SBDiagnostics& sb_diagnostics_legacy() const { return sb_diag_; } + + // Return a reference to the last stored samples matrix (legacy internal storage). + inline const MT& sb_samples_legacy() const { return sb_samples_; } + + // Set the internal Shake-and-Bake state from an external buffer. + // Copies samples (d x N) and diagnostic values into internal members. + void set_sb_state_from_buffer(int d, int N, const double* samples, const SBDiagnostics& diag); + + // Copy the internally stored sample matrix (sb_samples_) into a provided output buffer. + // The output pointer must have enough space for d * N doubles. + void get_sb_samples(double* out) const; + + // Copy the current diagnostics (sb_diag_) into a provided 3-element double buffer. + // Order: [minESS, maxPSRF, N]. + void get_sb_diagnostics(double* out3) const; private: - Diagnostics sb_diag_; + SBDiagnostics sb_diag_; MT sb_samples_; }; diff --git a/dingo/volestipy.pyx b/dingo/volestipy.pyx index 51f74677..e77182be 100644 --- a/dingo/volestipy.pyx +++ b/dingo/volestipy.pyx @@ -44,6 +44,11 @@ def get_time_seed(): # Get classes from the bindings.h file cdef extern from "bindings.h": + cdef struct SBDiagnostics: + double minESS + double maxPSRF + long long N + # The HPolytopeCPP class along with its functions cdef cppclass HPolytopeCPP: @@ -74,13 +79,19 @@ 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 - void get_sb_scaling_ratio(double tol, double min_ratio, - double* scale_out, - double* coverage_out, - double* maxdev_out, - double* avgdev_out) const + int apply_boundary_sampling(int walk_len, + int number_of_points, + int number_of_points_to_burn, + const char* sampler, + int nreflections, + double* samples) + + void set_sb_state_from_buffer(int d, int N, const double* samples, SBDiagnostics diag) + void get_sb_samples(double* out) + void get_sb_diagnostics(double* out5) + + @staticmethod + SBDiagnostics sb_diagnostics(int d, int N, const double* samples) # The lowDimPolytopeCPP class along with its functions @@ -147,75 +158,66 @@ cdef class HPolytope: variance_value, &bias_vector_[0], ess, nreflections) return np.asarray(samples) + def boundary_sample(self, sampler="sb", number_of_points=10000,number_of_points_to_burn=0, walk_len=1, nreflections=0): + cdef int d = self._A.shape[1] + cdef int N = number_of_points + cdef int burn = number_of_points_to_burn + cdef int wl = walk_len + cdef int nref = nreflections + cdef bytes sampler_b = sampler.encode("UTF-8") + + cdef np.ndarray[np.float64_t, ndim=2, mode="fortran"] S = \ + np.zeros((d, N), dtype=np.float64, order="F") + + cdef int made = self.polytope_cpp.apply_boundary_sampling(wl, N, burn, sampler_b, nref, &S[0,0]) + if made < 0: + raise RuntimeError("apply_boundary_sampling failed") + + cdef SBDiagnostics diag = HPolytopeCPP.sb_diagnostics(d, made, &S[0,0]) + self.polytope_cpp.set_sb_state_from_buffer(d, made, &S[0,0], diag) + + return np.asarray(S)[:, :made] + def get_sb_diagnostics(self, out): """ - out: np.ndarray float64, shape (5,), redosled: - [minESS, maxPSRF, N, phases, seconds] + out: np.ndarray float64, shape (>=3,), order: + [minESS, maxPSRF, N] """ - 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) + cdef np.ndarray[np.float64_t, ndim=1] out_arr = np.ascontiguousarray(out, dtype=np.float64) + if out_arr.shape[0] < 3: + raise ValueError("out must have length >= 3") + self.polytope_cpp.get_sb_diagnostics( &out_arr[0]) + return np.asarray(out_arr) + def get_sb_samples(self): """ - Return d x N numpy matrix of samples from Shake and Bake sampling + Returns a d x N matrix from the internal C++ sample buffer. """ - 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") + cdef int d + cdef long long Nll + cdef np.ndarray[np.float64_t, ndim=1] tmp = np.zeros(3, dtype=np.float64) + # Retrieve N from diagnostics (tmp = [minESS, maxPSRF, N]) + self.polytope_cpp.get_sb_diagnostics( &tmp[0]) + Nll = tmp[2] + d = self._A.shape[1] + + cdef np.ndarray[np.float64_t, ndim=2, mode="fortran"] S = \ + np.zeros((d, Nll), dtype=np.float64, order="F") self.polytope_cpp.get_sb_samples(&S[0,0]) return np.asarray(S) - def sb_scaling_ratio(self, tol=1e-10, min_ratio=0.01): + + def boundary_diag(self, S): """ - Vraća: (scale, coverage, max_dev, avg_dev) - - scale: (K,) [tipično K=10] - - coverage: (m,K) [m = broj hiper-ravnina = self._A.shape[0]] - - max_dev: (m,) - - avg_dev: (m,) - Potrebno: prethodno pozvati SB/BSB da bi postojali sb_samples_ u C++. + Compute diagnostics directly from a given sample matrix S. + Returns a Python dictionary with keys: minESS, maxPSRF, and N. """ - cdef int m = self._A.shape[0] - # Trenutna C++ implementacija pravi tačno 10 skala; ostavi 10 ovde. - cdef int K = 10 - - # Alokacije kao C-kontigvne matrice - cdef np.ndarray[np.float64_t, ndim=1, mode="c"] scale_np = np.empty((K,), dtype=np.float64, order="C") - cdef np.ndarray[np.float64_t, ndim=2, mode="c"] coverage_np = np.empty((m, K), dtype=np.float64, order="C") - cdef np.ndarray[np.float64_t, ndim=1, mode="c"] maxdev_np = np.empty((m,), dtype=np.float64, order="C") - cdef np.ndarray[np.float64_t, ndim=1, mode="c"] avgdev_np = np.empty((m,), dtype=np.float64, order="C") - - # Typed memoryview-ovi nad C-kontigvnim baferima - cdef double[::1] scale_mv = scale_np - cdef double[:,::1] coverage_mv = coverage_np - cdef double[::1] maxdev_mv = maxdev_np - cdef double[::1] avgdev_mv = avgdev_np - - # Poziv C++ wrappera – prosleđuju se sirove adrese bafera - try: - self.polytope_cpp.get_sb_scaling_ratio( - tol, - min_ratio, - &scale_mv[0], - &coverage_mv[0, 0], # (m,K) row-major: i*K + k - &maxdev_mv[0], - &avgdev_mv[0] - ) - except RuntimeError as e: - # C++ baca kada nema SB/BSB uzoraka; prebaci poruku dalje. - raise RuntimeError(str(e)) - - # Vrati NumPy objekte (već su odgovarajućih dimenzija) - return scale_np, coverage_np, maxdev_np, avgdev_np - + cdef np.ndarray[np.float64_t, ndim=2] Snp = np.ascontiguousarray(S, dtype=np.float64) + cdef int d = Snp.shape[0] + cdef int N = Snp.shape[1] + cdef SBDiagnostics diag = HPolytopeCPP.sb_diagnostics(d, N, &Snp[0,0]) + return {"minESS": diag.minESS, "maxPSRF": diag.maxPSRF, "N": diag.N} # The rounding() function; as in compute_volume, more than one method is available for this step From 6d49cd9467a8d28131ddcfc4312462e3db026b95 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Iva=20Jankovi=C4=87?= Date: Mon, 10 Nov 2025 03:35:14 +0000 Subject: [PATCH 13/14] Run example --- tests/run_sb_exp.py | 205 +++++++++++++++++++++++++++++++++++++------- 1 file changed, 173 insertions(+), 32 deletions(-) diff --git a/tests/run_sb_exp.py b/tests/run_sb_exp.py index 7463c2a2..bb129619 100644 --- a/tests/run_sb_exp.py +++ b/tests/run_sb_exp.py @@ -1,50 +1,191 @@ -#!/usr/bin/python3 -import os, sys +# VolEsti (volume computation and sampling library) + +# Copyright (c) 2012-2025 Vissarion Fisikopoulos +# Copyright (c) 2018-2025 Apostolos Chalkis +# Copyright (c) 2025-2025 Iva Janković + +# Contributed and/or modified by Iva Janković, as part of Google Summer of Code 2025 program. + +import os, sys, pickle, argparse import numpy as np -from dingo import PolytopeSampler -from time import process_time -import pickle +from time import perf_counter +from volestipy import HPolytope + +def _binary_search_k(P, S_all, L0, S_last, target_ess): + """Find smallest k_rel in [1, S_last.shape[1]] where minESS >= target_ess.""" + lo, hi = 1, S_last.shape[1] + best_k_rel, best_diag = hi, None + + def pref_diag(k_rel): + return P.boundary_diag(S_all[:, : L0 + k_rel]) + while lo <= hi: + mid = (lo + hi) // 2 + dmid = pref_diag(mid) + if float(dmid.get("minESS", -1.0)) >= target_ess: + best_k_rel, best_diag = mid, dmid + hi = mid - 1 + else: + lo = mid + 1 -def sample_on_polyround_processed_polytope(p): + if best_diag is None: + best_diag = P.boundary_diag(S_all[:, : L0 + best_k_rel]) + return best_k_rel, best_diag + +def sample_on_polyround_processed_polytope(p, target_ess, n_chunk, max_total, + sampler, walk_len, burn_in_init): name = os.path.basename(p) with open(p, "rb") as f: obj = pickle.load(f) polytope = obj[0] - polyround_A = polytope.A.to_numpy() - polyround_b = polytope.b.to_numpy() + A = np.ascontiguousarray(polytope.A.to_numpy(), dtype=np.float64) + b = np.ascontiguousarray(polytope.b.to_numpy(), dtype=np.float64) + assert A.flags['C_CONTIGUOUS'] and b.flags['C_CONTIGUOUS'] - print(f"Dimensions = {polyround_A.shape[1]}, Ograničenja = {polyround_A.shape[0]}") - start = process_time() + d, m = A.shape[1], A.shape[0] + m, n = A.shape + d_meta = getattr(polytope, "dimension", None) or getattr(polytope, "d", None) + print(f"m={m}, n={n}, d(from meta)={d_meta}") - steady_states, diag = PolytopeSampler.sample_from_polytope_sb_once( - polyround_A, - polyround_b, - n=15000, - burn_in=1000, - sampler="sb", - walk_len=48, - ) + # nreflections = ceil(0.25 * d) for BSB, else 0 + if sampler.lower() == "bsb": + nreflections = int(np.ceil(0.25 * d)) + else: + nreflections = 0 - end = process_time() + print(f"Dimensions = {d}, Constraints = {m}") - print( - f"[{name}] minESS={diag['minESS']:.3f} maxPSRF={diag['maxPSRF']:.3f} " - f"N={diag['N']} phases={diag['phases']} seconds={diag['seconds']}" - ) - print(f"Total sampling time: {end - start:.2f} s") + P = HPolytope(A, b) + + S_parts, dur = [], [] + total = 0 + burn = int(burn_in_init) - out_path = f"dingo_polyround_no_multiphase_{name}.pckl" - with open(out_path, "wb") as f_out: - pickle.dump({"samples": steady_states, "diagnostics": diag}, f_out) + while True: + t0 = perf_counter() + S_chunk = P.boundary_sample( + sampler="sb", + number_of_points=int(n_chunk), + number_of_points_to_burn=int(burn), + walk_len=int(walk_len), + nreflections=int(nreflections), + ) + t1 = perf_counter() - print(f"Saved results to {out_path}") + dur.append(t1 - t0) + S_parts.append(S_chunk) + total += S_chunk.shape[1] + burn = 0 # continue without burn-in + S_all = np.concatenate(S_parts, axis=1) + diag_all = P.boundary_diag(S_all) + minESS_all = float(diag_all.get("minESS", np.nan)) + + if minESS_all >= target_ess: + # find minimal prefix in last chunk + L0 = total - S_chunk.shape[1] + k_rel, diag_at_k = _binary_search_k(P, S_all, L0, S_chunk, target_ess) + k_abs = L0 + k_rel + + # sampling time = sum of past chunks + part of last chunk + frac = k_rel / S_chunk.shape[1] + sampling_seconds = sum(dur[:-1]) + frac * dur[-1] + chunks_used = len(dur) + + print( + f"[{name}] minESS={diag_at_k['minESS']:.3f} " + f"maxPSRF={diag_at_k.get('maxPSRF', np.nan):.3f} " + f"N={int(diag_at_k.get('N', k_abs))} " + f"chunks={chunks_used} sampling_seconds={sampling_seconds:.3f}" + ) + print(f"Total sampling time: {sampling_seconds:.2f} s") + + out_path = f"dingo_polyround_no_multiphase_{name}" + with open(out_path, "wb") as f_out: + pickle.dump( + { + "samples": S_all[:, :k_abs], + "diagnostics": { + **diag_at_k, + "N_at_threshold": int(k_abs), + "chunks": chunks_used, + "sampling_seconds": sampling_seconds, + }, + "params": { + "sampler": sampler, + "walk_len": int(walk_len), + "nreflections": int(nreflections), + "target_ess": target_ess, + "n_chunk": int(n_chunk), + "max_total": int(max_total), + "burn_in_init": int(burn_in_init), + }, + }, + f_out, + ) + print(f"Saved results to {out_path}") + return + + if total >= max_total: + sampling_seconds = sum(dur) + chunks_used = len(dur) + + print( + f"[{name}] minESS={diag_all.get('minESS', np.nan):.3f} " + f"maxPSRF={diag_all.get('maxPSRF', np.nan):.3f} " + f"N={int(diag_all.get('N', S_all.shape[1]))} " + f"chunks={chunks_used} sampling_seconds={sampling_seconds:.3f}" + ) + print(f"Total sampling time: {sampling_seconds:.2f} s (target {target_ess} not reached)") + + out_path = f"dingo_polyround_no_multiphase_{name}" + with open(out_path, "wb") as f_out: + pickle.dump( + { + "samples": S_all, + "diagnostics": { + **diag_all, + "N_total_generated": int(S_all.shape[1]), + "chunks": chunks_used, + "sampling_seconds": sampling_seconds, + }, + "params": { + "sampler": sampler, + "walk_len": int(walk_len), + "nreflections": int(nreflections), + "target_ess": target_ess, + "n_chunk": int(n_chunk), + "max_total": int(max_total), + "burn_in_init": int(burn_in_init), + }, + }, + f_out, + ) + print(f"Saved results to {out_path}") + return + +def _make_parser(): + ap = argparse.ArgumentParser() + ap.add_argument("polytope_pickle", type=str) + ap.add_argument("--target-ess", type=float, default=1000.0) + ap.add_argument("--n-chunk", type=int, default=5000) + ap.add_argument("--max-total", type=int, default=200000) + ap.add_argument("--sampler", type=str, default="bsb", choices=["sb", "bsb"]) + ap.add_argument("--walk-len", type=int, default=1) + ap.add_argument("--burn-in-init", type=int, default=0) + return ap if __name__ == "__main__": - if len(sys.argv) < 2: - sys.exit("Usage: python run_bwr_exp.py ") - file_name = sys.argv[1] - sample_on_polyround_processed_polytope(file_name) + parser = _make_parser() + args = parser.parse_args() + sample_on_polyround_processed_polytope( + p=args.polytope_pickle, + target_ess=args.target_ess, + n_chunk=args.n_chunk, + max_total=args.max_total, + sampler=args.sampler, + walk_len=args.walk_len, + burn_in_init=args.burn_in_init, + ) From c0ca9ef3607ffe971d23a6e056bd00fe55f588ec Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Iva=20Jankovi=C4=87?= Date: Tue, 25 Nov 2025 01:36:02 +0000 Subject: [PATCH 14/14] (Billiard) Shake and Bake - CYTHON / bindings final --- dingo/bindings/bindings.cpp | 144 +++++++++++++++++++++++++++++++++ dingo/bindings/bindings.h | 21 +++++ dingo/volestipy.pyx | 155 ++++++++++++++++++++++++++++++++++-- 3 files changed, 313 insertions(+), 7 deletions(-) diff --git a/dingo/bindings/bindings.cpp b/dingo/bindings/bindings.cpp index af6bf810..7d759861 100644 --- a/dingo/bindings/bindings.cpp +++ b/dingo/bindings/bindings.cpp @@ -507,6 +507,44 @@ void HPolytopeCPP::get_sb_diagnostics(double* out3) const { } +void HPolytopeCPP::boundary_scaling_ratio(int d,int N,const double* samples,double tol,double min_ratio,double* scale_out,double* coverage_out,double* max_dev_out,double* avg_dev_out) const +{ + MT S(d, N); + for (int j = 0; j < N; ++j) + { + for (int i = 0; i < d; ++i) + { + S(i, j) = samples[i + j * d]; + } + } + + auto result = scaling_ratio_boundary_test(HP, S, tol, min_ratio); + const VT& scale = std::get<0>(result); + const MT& coverage = std::get<1>(result); + const VT& max_dev = std::get<2>(result); + const VT& avg_dev = std::get<3>(result); + + const int K = static_cast(scale.size()); + const int m = static_cast(coverage.rows()); + + for (int k = 0; k < K; ++k) + { + scale_out[k] = static_cast(scale[k]); + } + + for (int f = 0; f < m; ++f) + { + max_dev_out[f] = static_cast(max_dev[f]); + avg_dev_out[f] = static_cast(avg_dev[f]); + + for (int k = 0; k < K; ++k) + { + coverage_out[f * K + k] = static_cast(coverage(f, k)); + } + } +} + + ////////// Start of "rounding()" ////////// void HPolytopeCPP::apply_rounding(int rounding_method, double* new_A, double* new_b, @@ -591,3 +629,109 @@ void HPolytopeCPP::apply_rounding(int rounding_method, double* new_A, double* ne } ////////// End of "rounding()" ////////// + +////////// Known H-polytope generators wrappers ////////// + +void generate_cube_H(int dim, double scale, double* A_out, double* b_out) +{ + Hpolytope P = generate_cube(static_cast(dim), + false, // Vpoly = false + scale); + const MT& A = P.get_mat(); + const VT& b = P.get_vec(); + const int m = static_cast(A.rows()); + const int n = static_cast(A.cols()); + + for (int i = 0; i < m; ++i) { + b_out[i] = static_cast(b[i]); + for (int j = 0; j < n; ++j) { + A_out[i * n + j] = static_cast(A(i, j)); + } + } +} + +void generate_cross_H(int dim, double* A_out, double* b_out) +{ + Hpolytope P = generate_cross(static_cast(dim),false); + const MT& A = P.get_mat(); + const VT& b = P.get_vec(); + + const int m = static_cast(A.rows()); + const int n = static_cast(A.cols()); + + for (int i = 0; i < m; ++i) { + b_out[i] = static_cast(b[i]); + for (int j = 0; j < n; ++j) { + A_out[i * n + j] = static_cast(A(i, j)); + } + } +} + +void generate_simplex_H(int dim, double* A_out, double* b_out) +{ + Hpolytope P = generate_simplex(static_cast(dim),false); + const MT& A = P.get_mat(); + const VT& b = P.get_vec(); + + const int m = static_cast(A.rows()); + const int n = static_cast(A.cols()); + + for (int i = 0; i < m; ++i) { + b_out[i] = static_cast(b[i]); + for (int j = 0; j < n; ++j) { + A_out[i * n + j] = static_cast(A(i, j)); + } + } +} + +void generate_prod_simplex_H(int dim, double* A_out, double* b_out) +{ + Hpolytope P = generate_prod_simplex(static_cast(dim),false); + const MT& A = P.get_mat(); + const VT& b = P.get_vec(); + + const int m = static_cast(A.rows()); + const int n = static_cast(A.cols()); + + for (int i = 0; i < m; ++i) { + b_out[i] = static_cast(b[i]); + for (int j = 0; j < n; ++j) { + A_out[i * n + j] = static_cast(A(i, j)); + } + } +} + +void generate_skinny_cube_H(int dim, double* A_out, double* b_out) +{ + Hpolytope P = generate_skinny_cube(static_cast(dim),false); + const MT& A = P.get_mat(); + const VT& b = P.get_vec(); + + const int m = static_cast(A.rows()); + const int n = static_cast(A.cols()); + + for (int i = 0; i < m; ++i) { + b_out[i] = static_cast(b[i]); + for (int j = 0; j < n; ++j) { + A_out[i * n + j] = static_cast(A(i, j)); + } + } +} + +void generate_birkhoff_H(int n, double* A_out, double* b_out) +{ + Hpolytope P = generate_birkhoff(static_cast(n)); + const MT& A = P.get_mat(); + const VT& b = P.get_vec(); + + const int m = static_cast(A.rows()); + const int d = static_cast(A.cols()); + + for (int i = 0; i < m; ++i) { + b_out[i] = static_cast(b[i]); + for (int j = 0; j < d; ++j) { + A_out[i * d + j] = static_cast(A(i, j)); + } + } +} +////////// Ending of known H-polytope generators wrappers ////////// \ No newline at end of file diff --git a/dingo/bindings/bindings.h b/dingo/bindings/bindings.h index 0bdb058a..6eb3c92f 100644 --- a/dingo/bindings/bindings.h +++ b/dingo/bindings/bindings.h @@ -31,6 +31,7 @@ #include "sampling/parallel_mmcs.hpp" #include "diagnostics/univariate_psrf.hpp" #include "diagnostics/effective_sample_size.hpp" +#include "diagnostics/scaling_ratio.hpp" //from generate_samples, some extra headers not already included #include @@ -45,6 +46,9 @@ #include "preprocess/inscribed_ellipsoid_rounding.hpp" #include "preprocess/feasible_point.hpp" +// for creating H polytopes +#include "generators/known_polytope_generators.h" + typedef double NT; typedef Cartesian Kernel; typedef typename Kernel::Point Point; @@ -191,11 +195,28 @@ class HPolytopeCPP{ // Order: [minESS, maxPSRF, N]. void get_sb_diagnostics(double* out3) const; + // Compute boundary scaling-ratio diagnostics for a given sample buffer. + // - scale_out: length K (currently 10) scaling factors + // - coverage_out: m x K coverage matrix, stored row-major (facet-major) + // - max_dev_out: length m, maximum deviation per facet (in %) + // - avg_dev_out: length m, average deviation per facet (in %) + void boundary_scaling_ratio(int d,int N,const double* samples,double tol,double min_ratio,double* scale_out,double* coverage_out,double* max_dev_out,double* avg_dev_out) const; + private: SBDiagnostics sb_diag_; MT sb_samples_; }; +// Known H-polytopes generators + +void generate_cube_H(int dim, double scale, double* A_out, double* b_out); +void generate_cross_H(int dim, double* A_out, double* b_out); +void generate_simplex_H(int dim, double* A_out, double* b_out); +void generate_prod_simplex_H(int dim, double* A_out, double* b_out); +void generate_skinny_cube_H(int dim, double* A_out, double* b_out); +void generate_birkhoff_H(int n, double* A_out, double* b_out); + + #endif diff --git a/dingo/volestipy.pyx b/dingo/volestipy.pyx index e77182be..40112251 100644 --- a/dingo/volestipy.pyx +++ b/dingo/volestipy.pyx @@ -44,6 +44,13 @@ def get_time_seed(): # Get classes from the bindings.h file cdef extern from "bindings.h": + void generate_cube_H(int dim, double scale, double* A_out, double* b_out) + void generate_cross_H(int dim, double* A_out, double* b_out) + void generate_simplex_H(int dim, double* A_out, double* b_out) + void generate_prod_simplex_H(int dim, double* A_out, double* b_out) + void generate_skinny_cube_H(int dim, double* A_out, double* b_out) + void generate_birkhoff_H(int n, double* A_out, double* b_out) + cdef struct SBDiagnostics: double minESS double maxPSRF @@ -79,16 +86,12 @@ 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); - int apply_boundary_sampling(int walk_len, - int number_of_points, - int number_of_points_to_burn, - const char* sampler, - int nreflections, - double* samples) + int apply_boundary_sampling(int walk_len,int number_of_points,int number_of_points_to_burn,const char* sampler,int nreflections,double* samples) void set_sb_state_from_buffer(int d, int N, const double* samples, SBDiagnostics diag) void get_sb_samples(double* out) void get_sb_diagnostics(double* out5) + void boundary_scaling_ratio(int d,int N,const double* samples,double tol,double min_ratio,double* scale_out,double* coverage_out,double* max_dev_out,double* avg_dev_out) const @staticmethod SBDiagnostics sb_diagnostics(int d, int N, const double* samples) @@ -218,7 +221,36 @@ cdef class HPolytope: cdef int N = Snp.shape[1] cdef SBDiagnostics diag = HPolytopeCPP.sb_diagnostics(d, N, &Snp[0,0]) return {"minESS": diag.minESS, "maxPSRF": diag.maxPSRF, "N": diag.N} - + + def boundary_scaling_ratio(self, S, double tol=1e-10, double min_ratio=0.01): + """ + Compute boundary scaling-ratio diagnostics from a sample matrix S. + Returns + scale : np.ndarray, shape (K,) + coverage : np.ndarray, shape (m, K) + max_dev : np.ndarray, shape (m,) + avg_dev : np.ndarray, shape (m,) + """ + cdef np.ndarray[np.float64_t, ndim=2, mode="fortran"] Snp = \ + np.array(S, dtype=np.float64, order="F") + + cdef int d = Snp.shape[0] + cdef int N = Snp.shape[1] + cdef int m = self._A.shape[0] + cdef int K = 10 # must match C++ scale(10) + + cdef np.ndarray[np.float64_t, ndim=1] scale = \ + np.zeros(K, dtype=np.float64) + cdef np.ndarray[np.float64_t, ndim=1] max_dev = \ + np.zeros(m, dtype=np.float64) + cdef np.ndarray[np.float64_t, ndim=1] avg_dev = \ + np.zeros(m, dtype=np.float64) + cdef np.ndarray[np.float64_t, ndim=2] coverage = \ + np.zeros((m, K), dtype=np.float64, order="C") + + self.polytope_cpp.boundary_scaling_ratio(d,N, &Snp[0, 0],tol,min_ratio,&scale[0],&coverage[0, 0],&max_dev[0],&avg_dev[0]) + + return (np.asarray(scale),np.asarray(coverage),np.asarray(max_dev),np.asarray(avg_dev),) # 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): @@ -298,3 +330,112 @@ cdef class HPolytope: def dimension(self): return self._A.shape[1] + +def create_cube(int dim, double scale=1.0): + """ + Create an H-polytope for a hypercube in R^dim with side length 2*scale + centered at the origin: [-scale, scale]^dim. + + """ + cdef int m = 2 * dim + cdef int n = dim + + cdef np.ndarray[np.float64_t, ndim=2] A = \ + np.zeros((m, n), dtype=np.float64) + cdef np.ndarray[np.float64_t, ndim=1] b = \ + np.zeros(m, dtype=np.float64) + + generate_cube_H(dim, scale, &A[0, 0], &b[0]) + + return np.asarray(A), np.asarray(b) + + +def create_crosspolytope(int dim): + """ + Create an H-polytope for the crosspolytope (l1-ball) in R^dim. + + """ + cdef int m = 1 << dim # 2^(dim) + cdef int n = dim + + cdef np.ndarray[np.float64_t, ndim=2] A = \ + np.zeros((m, n), dtype=np.float64) + cdef np.ndarray[np.float64_t, ndim=1] b = \ + np.zeros(m, dtype=np.float64) + + generate_cross_H(dim, &A[0, 0], &b[0]) + + return np.asarray(A), np.asarray(b) + + +def create_simplex(int dim): + """ + Create an H-polytope for a standard simplex in R^dim. + + """ + cdef int m = dim + 1 + cdef int n = dim + + cdef np.ndarray[np.float64_t, ndim=2] A = \ + np.zeros((m, n), dtype=np.float64) + cdef np.ndarray[np.float64_t, ndim=1] b = \ + np.zeros(m, dtype=np.float64) + + generate_simplex_H(dim, &A[0, 0], &b[0]) + + return np.asarray(A), np.asarray(b) + + +def create_prod_simplex(int dim): + """ + Create an H-polytope for the product of two dim-dimensional simplices. + + """ + cdef int m = 2 * dim + 2 + cdef int n = 2 * dim + + cdef np.ndarray[np.float64_t, ndim=2] A = \ + np.zeros((m, n), dtype=np.float64) + cdef np.ndarray[np.float64_t, ndim=1] b = \ + np.zeros(m, dtype=np.float64) + + generate_prod_simplex_H(dim, &A[0, 0], &b[0]) + + return np.asarray(A), np.asarray(b) + + +def create_skinny_cube(int dim): + """ + Create an H-polytope for a 'skinny cube' in R^dim, + where the first coordinate is scaled differently. + + """ + cdef int m = 2 * dim + cdef int n = dim + + cdef np.ndarray[np.float64_t, ndim=2] A = \ + np.zeros((m, n), dtype=np.float64) + cdef np.ndarray[np.float64_t, ndim=1] b = \ + np.zeros(m, dtype=np.float64) + + generate_skinny_cube_H(dim, &A[0, 0], &b[0]) + + return np.asarray(A), np.asarray(b) + + +def create_birkhoff(int n): + """ + Creates an H-polytope for the Birkhoff polytope of size n. + + """ + cdef int m = n * n + cdef int d = n * n - 2 * n + 1 + + cdef np.ndarray[np.float64_t, ndim=2] A = \ + np.zeros((m, d), dtype=np.float64) + cdef np.ndarray[np.float64_t, ndim=1] b = \ + np.zeros(m, dtype=np.float64) + + generate_birkhoff_H(n, &A[0, 0], &b[0]) + + return np.asarray(A), np.asarray(b)