Skip to content

Remove mmcs option for generate_steady_states_no_multiphase and fix a warning #122

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
vissarion merged 2 commits into from
Feb 1, 2026
Merged

Remove mmcs option for generate_steady_states_no_multiphase and fix a warning #122

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
8fcc8d0
Remove mmcs option from apply_sampling and generate_steady_states_no_…
vissarion Dec 9, 2025
04e8376
Remove deprecated warning
vissarion Dec 9, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
23 changes: 13 additions & 10 deletions dingo/PolytopeSampler.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 == []
len(self._A) == 0
or len(self._b) == 0
or len(self._N) == 0
or len(self._N_shift) == 0
or len(self._T) == 0
or len(self._T_shift) == 0
):

(
Expand Down Expand Up @@ -161,7 +161,8 @@ def generate_steady_states(
return steady_states

def generate_steady_states_no_multiphase(
self, method = 'billiard_walk', n=1000, burn_in=0, thinning=1, variance=1.0, bias_vector=None, ess=1000
self, method = 'billiard_walk', n=1000, burn_in=0, thinning=1,
variance=1.0, bias_vector=None
):
"""A member function to sample steady states.

Expand All @@ -171,6 +172,8 @@ def generate_steady_states_no_multiphase(
burn_in -- the number of points to burn before sampling
thinning -- the walk length of the chain
"""
if method == "mmcs":
raise ValueError("method cannot be 'mmcs'. Use generate_steady_states() instead.")

self.get_polytope()

Expand All @@ -181,7 +184,7 @@ def generate_steady_states_no_multiphase(
else:
bias_vector = bias_vector.astype('float64')

samples = P.generate_samples(method.encode('utf-8'), n, burn_in, thinning, variance, bias_vector, self._parameters["solver"], ess)
samples = P.generate_samples(method.encode('utf-8'), n, burn_in, thinning, variance, bias_vector, self._parameters["solver"])
samples_T = samples.T

steady_states = map_samples_to_steady_states(
Expand Down Expand Up @@ -216,7 +219,7 @@ def sample_from_polytope(

@staticmethod
def sample_from_polytope_no_multiphase(
A, b, method = 'billiard_walk', n=1000, burn_in=0, thinning=1, variance=1.0, bias_vector=None, solver=None, ess=1000
A, b, method = 'billiard_walk', n=1000, burn_in=0, thinning=1, variance=1.0, bias_vector=None, solver=None
):
"""A static function to sample from a full dimensional polytope with an MCMC method.

Expand All @@ -235,7 +238,7 @@ def sample_from_polytope_no_multiphase(

P = HPolytope(A, b)

samples = P.generate_samples(method.encode('utf-8'), n, burn_in, thinning, variance, bias_vector, solver, ess)
samples = P.generate_samples(method.encode('utf-8'), n, burn_in, thinning, variance, bias_vector, solver)

samples_T = samples.T
return samples_T
Expand Down
24 changes: 7 additions & 17 deletions dingo/bindings/bindings.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -91,8 +91,7 @@ double HPolytopeCPP::apply_sampling(int walk_len,
double radius,
double* samples,
double variance_value,
double* bias_vector_,
int ess){
double* bias_vector_){

RNGType rng(HP.dimension());
HP.normalize();
Expand Down Expand Up @@ -134,13 +133,6 @@ double HPolytopeCPP::apply_sampling(int walk_len,
} else if (strcmp(method, "vaidya_walk") == 0) { // vaidya walk
uniform_sampling<VaidyaWalk>(rand_points, HP, rng, walk_len, number_of_points,
starting_point, number_of_points_to_burn);
} else if (strcmp(method, "mmcs") == 0) { // vaidya walk
MT S;
int total_ess;
//TODO: avoid passing polytopes as non-const references
const Hpolytope HP_const = HP;
mmcs(HP_const, ess, S, total_ess, walk_len, rng);
samples = S.data();
} else if (strcmp(method, "gaussian_hmc_walk") == 0) { // Gaussian sampling with exact HMC walk
NT a = NT(1)/(NT(2)*variance);
gaussian_sampling<GaussianHamiltonianMonteCarloExactWalk>(rand_points, HP, rng, walk_len, number_of_points, a,
Expand Down Expand Up @@ -170,14 +162,12 @@ double HPolytopeCPP::apply_sampling(int walk_len,
throw std::runtime_error("This function must not be called.");
}

if (strcmp(method, "mmcs") != 0) {
// The following block of code allows us to copy the sampled points
auto n_si=0;
for (auto it_s = rand_points.cbegin(); it_s != rand_points.cend(); it_s++){
for (auto i = 0; i != it_s->dimension(); i++){
samples[n_si++] = (*it_s)[i];
}
}
// The following block of code allows us to copy the sampled points
auto n_si=0;
for (auto it_s = rand_points.cbegin(); it_s != rand_points.cend(); it_s++){
for (auto i = 0; i != it_s->dimension(); i++){
samples[n_si++] = (*it_s)[i];
}
}
return 0.0;
}
Expand Down
2 changes: 1 addition & 1 deletion dingo/bindings/bindings.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -141,7 +141,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);

void mmcs_initialize(int d, int ess, bool psrf_check, bool parallelism, int num_threads);

Expand Down
6 changes: 3 additions & 3 deletions dingo/volestipy.pyx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -55,7 +55,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)

# 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);
Expand Down Expand Up @@ -119,7 +119,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):

n_variables = self._A.shape[1]
cdef double[:,::1] samples = np.zeros((number_of_points, n_variables), dtype = np.float64, order = "C")
Expand All @@ -133,7 +133,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])
return np.asarray(samples)


Expand Down
42 changes: 16 additions & 26 deletions tests/sampling_no_multiphase.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
# dingo is part of GeomScale project

# Copyright (c) 2022 Apostolos Chalkis
# Copyright (c) 2022 Vissarion Fisikopoulos
# Copyright (c) 2022-2025 Vissarion Fisikopoulos
# Copyright (c) 2022 Haris Zafeiropoulos
# Copyright (c) 2024 Ke Shi

Expand All @@ -11,41 +11,33 @@
import unittest
import os
import sys
import numpy as np
from dingo import MetabolicNetwork, PolytopeSampler
from dingo.pyoptinterface_based_impl import set_default_solver

def test_shape_and_non_zero(steady_states, shape0, shape1, testing_class):
testing_class.assertTrue( steady_states.shape[0] == shape0 )
testing_class.assertTrue( steady_states.shape[1] == shape1 )
testing_class.assertTrue( np.any(np.abs(steady_states) > 1e-12) )

def sampling(model, testing_class):
sampler = PolytopeSampler(model)

#Gaussian hmc sampling
steady_states = sampler.generate_steady_states_no_multiphase(method = 'mmcs', ess=1000)

testing_class.assertTrue( steady_states.shape[0] == 95 )
testing_class.assertTrue( steady_states.shape[1] == 1000 )

#Gaussian hmc sampling
steady_states = sampler.generate_steady_states_no_multiphase(method = 'gaussian_hmc_walk', n=500)

testing_class.assertTrue( steady_states.shape[0] == 95 )
testing_class.assertTrue( steady_states.shape[1] == 500 )
steady_states = sampler.generate_steady_states_no_multiphase(method = 'gaussian_hmc_walk', n=1000)
test_shape_and_non_zero(steady_states, 95, 1000, testing_class)

#exponential hmc sampling
steady_states = sampler.generate_steady_states_no_multiphase(method = 'exponential_hmc_walk', n=500, variance=50)

testing_class.assertTrue( steady_states.shape[0] == 95 )
testing_class.assertTrue( steady_states.shape[1] == 500 )

steady_states = sampler.generate_steady_states_no_multiphase(method = 'exponential_hmc_walk', n=1000, variance=50)
test_shape_and_non_zero(steady_states, 95, 1000, testing_class)

#hmc sampling with Gaussian distribution
steady_states = sampler.generate_steady_states_no_multiphase(method = 'hmc_leapfrog_gaussian', n=500)

testing_class.assertTrue( steady_states.shape[0] == 95 )
testing_class.assertTrue( steady_states.shape[1] == 500 )
steady_states = sampler.generate_steady_states_no_multiphase(method = 'hmc_leapfrog_gaussian', n=1000)
test_shape_and_non_zero(steady_states, 95, 1000, testing_class)

#hmc sampling with exponential distribution
steady_states = sampler.generate_steady_states_no_multiphase(method = 'hmc_leapfrog_exponential', n=500, variance=50)

testing_class.assertTrue( steady_states.shape[0] == 95 )
testing_class.assertTrue( steady_states.shape[1] == 500 )
steady_states = sampler.generate_steady_states_no_multiphase(method = 'hmc_leapfrog_exponential', n=1000, variance=50)
test_shape_and_non_zero(steady_states, 95, 1000, testing_class)

#steady_states[12].mean() seems to have a lot of discrepancy between experiments, so we won't check the mean for now
#self.assertTrue( abs( steady_states[12].mean() - 2.504 ) < 1e-03 )
Expand All @@ -70,8 +62,6 @@ def test_sample_sbml(self):
model = MetabolicNetwork.from_sbml( input_file_sbml )
sampling(model, self)



if __name__ == "__main__":
if len(sys.argv) > 1:
set_default_solver(sys.argv[1])
Expand Down