Skip to content

Fixing inconsistent results when using seed in parmest examples #3621

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
blnicho merged 21 commits into from
Jul 1, 2025
Merged

Fixing inconsistent results when using seed in parmest examples #3621

Show file tree
Hide file tree
Changes from 18 commits
Commits
Show all changes
21 commits
Select commit Hold shift + click to select a range
e9308d2
Added locations of fixes needed.
sscini May 21, 2025
34849d7
Merge branch 'Pyomo:main' into fix-SSE-w-all-theta_est
sscini May 28, 2025
effdc8c
Create mwe_bootstrap_error_SSE.py
sscini May 28, 2025
82b4a7e
Moved temp error file to correct folder
sscini May 28, 2025
3c97998
Update mwe_bootstrap_error_SSE.py
sscini May 29, 2025
72e9974
Commented out plotting code for output simplicity.
sscini May 29, 2025
b1e1181
Renamed example files, and added potential solution to parmest.py
sscini May 29, 2025
d268613
Update parmest.py
sscini May 29, 2025
ab2e67b
Update for discussion in dev meeting
sscini Jun 2, 2025
e75f03b
Merge branch 'Pyomo:main' into fix-SSE-w-all-theta_est
sscini Jun 6, 2025
ce310ef
Added seed to 3 examples, removed temp.py files
sscini Jun 6, 2025
279036f
Removed extra spaces
sscini Jun 6, 2025
f0d74d5
Seed is now available option for all random generators I could find
sscini Jun 6, 2025
1fbe06c
Confirmed tests run with changes
sscini Jun 6, 2025
c19d59f
Ran black on parmest folder
sscini Jun 23, 2025
0da9794
Merge branch 'Pyomo:main' into Example-seed-issue
sscini Jun 23, 2025
a1c1073
Update test_parmest.py
sscini Jun 24, 2025
f760b8f
Merge branch 'Pyomo:main' into Example-seed-issue
sscini Jun 24, 2025
ab30261
Update test_parmest.py
sscini Jun 26, 2025
466e3b4
Ran black, all tests pass, final lookover
sscini Jun 26, 2025
d7db8d3
Merge branch 'main' into Example-seed-issue
blnicho Jun 30, 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
2 changes: 1 addition & 1 deletion pyomo/contrib/parmest/examples/reactor_design/bootstrap_example.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,7 @@ def main():
obj, theta = pest.theta_est()

# Parameter estimation with bootstrap resampling
bootstrap_theta = pest.theta_est_bootstrap(50)
bootstrap_theta = pest.theta_est_bootstrap(50, seed=524)

# Plot results
parmest.graphics.pairwise_plot(bootstrap_theta, title="Bootstrap theta")
Expand Down
4 changes: 2 additions & 2 deletions pyomo/contrib/parmest/examples/reactor_design/confidence_region_example.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,11 +39,11 @@ def main():
obj, theta = pest.theta_est()

# Bootstrapping
bootstrap_theta = pest.theta_est_bootstrap(10)
bootstrap_theta = pest.theta_est_bootstrap(10, seed=524)
print(bootstrap_theta)

# Confidence region test
CR = pest.confidence_region_test(bootstrap_theta, "MVN", [0.5, 0.75, 1.0])
CR = pest.confidence_region_test(bootstrap_theta, "MVN", [0.5, 0.75, 1.0], seed=524)
print(CR)


Expand Down
9 changes: 7 additions & 2 deletions pyomo/contrib/parmest/examples/reactor_design/leaveNout_example.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,9 @@ def main():
file_name = abspath(join(file_dirname, "reactor_data.csv"))
data = pd.read_csv(file_name)

seed = 524 # Set seed for reproducibility
np.random.seed(seed) # Set seed for reproducibility

# Create more data for the example
N = 50
df_std = data.std().to_frame().transpose()
Expand All @@ -50,10 +53,10 @@ def main():
### Parameter estimation with 'leave-N-out'
# Example use case: For each combination of data where one data point is left
# out, estimate theta
lNo_theta = pest.theta_est_leaveNout(1)
lNo_theta = pest.theta_est_leaveNout(1, seed=524)
print(lNo_theta.head())

parmest.graphics.pairwise_plot(lNo_theta, theta)
parmest.graphics.pairwise_plot(lNo_theta, theta, seed=524)

### Leave one out/boostrap analysis
# Example use case: leave 25 data points out, run 20 bootstrap samples with the
Expand Down Expand Up @@ -81,6 +84,8 @@ def main():
alpha,
["MVN"],
title="Alpha: " + str(alpha) + ", " + str(theta_est_N.loc[0, alpha]),
seed=524 + i, # setting the seed for testing repeatability,
# for typical use cases, this should not be set
)

# Extract the percent of points that are within the alpha region
Expand Down
22 changes: 16 additions & 6 deletions pyomo/contrib/parmest/graphics.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -218,6 +218,7 @@ def pairwise_plot(
add_obj_contour=True,
add_legend=True,
filename=None,
seed=None,
):
"""
Plot pairwise relationship for theta values, and optionally alpha-level
Expand Down Expand Up @@ -283,6 +284,9 @@ def pairwise_plot(
assert isinstance(add_obj_contour, bool)
assert isinstance(filename, (type(None), str))

if seed is not None:
np.random.seed(seed)

# If theta_values is a tuple containing (mean, cov, n), create a DataFrame of values
if isinstance(theta_values, tuple):
assert len(theta_values) == 3
Expand All @@ -292,7 +296,7 @@ def pairwise_plot(
if isinstance(mean, dict):
mean = pd.Series(mean)
theta_names = mean.index
mvn_dist = stats.multivariate_normal(mean, cov)
mvn_dist = stats.multivariate_normal(mean, cov, seed=seed)
theta_values = pd.DataFrame(
mvn_dist.rvs(n, random_state=1), columns=theta_names
)
Expand Down Expand Up @@ -402,7 +406,7 @@ def pairwise_plot(
)

elif dist == "MVN":
mvn_dist = fit_mvn_dist(thetas)
mvn_dist = fit_mvn_dist(thetas, seed=seed)
Z = mvn_dist.pdf(thetas)
score = stats.scoreatpercentile(Z, (1 - alpha) * 100)
g.map_offdiag(
Expand All @@ -419,7 +423,7 @@ def pairwise_plot(
)

elif dist == "KDE":
kde_dist = fit_kde_dist(thetas)
kde_dist = fit_kde_dist(thetas, seed=seed)
Z = kde_dist.pdf(thetas.transpose())
score = stats.scoreatpercentile(Z, (1 - alpha) * 100)
g.map_offdiag(
Expand Down Expand Up @@ -512,7 +516,7 @@ def fit_rect_dist(theta_values, alpha):
return lower_bound, upper_bound


def fit_mvn_dist(theta_values):
def fit_mvn_dist(theta_values, seed=None):
"""
Fit a multivariate normal distribution to theta values

Expand All @@ -527,13 +531,17 @@ def fit_mvn_dist(theta_values):
"""
assert isinstance(theta_values, pd.DataFrame)

if seed is not None:
np.random.seed(seed)

dist = stats.multivariate_normal(
theta_values.mean(), theta_values.cov(), allow_singular=True
theta_values.mean(), theta_values.cov(), allow_singular=True, seed=seed
)

return dist


def fit_kde_dist(theta_values):
def fit_kde_dist(theta_values, seed=None):
"""
Fit a Gaussian kernel-density distribution to theta values

Expand All @@ -547,6 +555,8 @@ def fit_kde_dist(theta_values):
scipy.stats.gaussian_kde distribution
"""
assert isinstance(theta_values, pd.DataFrame)
if seed is not None:
np.random.seed(seed)

dist = stats.gaussian_kde(theta_values.transpose().values)

Expand Down
21 changes: 14 additions & 7 deletions pyomo/contrib/parmest/parmest.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -407,7 +407,6 @@ def _create_parmest_model(self, experiment_number):

# Add objective function (optional)
if self.obj_function:

# Check for component naming conflicts
reserved_names = [
'Total_Cost_Objective',
Expand Down Expand Up @@ -818,7 +817,7 @@ def _Q_at_theta(self, thetavals, initialize_parmest_model=False):

return retval, thetavals, WorstStatus

def _get_sample_list(self, samplesize, num_samples, replacement=True):
def _get_sample_list(self, samplesize, num_samples, replacement=True, seed=None):
samplelist = list()

scenario_numbers = list(range(len(self.exp_list)))
Expand All @@ -835,6 +834,8 @@ def _get_sample_list(self, samplesize, num_samples, replacement=True):
while (unique_samples <= len(self._return_theta_names())) and (
not duplicate
):
# if seed is not None:
# np.random.seed(seed) # set seed for reproducibility
sample = np.random.choice(
scenario_numbers, samplesize, replace=replacement
)
Expand Down Expand Up @@ -1036,7 +1037,9 @@ def theta_est_leaveNout(
if seed is not None:
np.random.seed(seed)

global_list = self._get_sample_list(samplesize, lNo_samples, replacement=False)
global_list = self._get_sample_list(
samplesize, lNo_samples, replacement=False, seed=seed
)

task_mgr = utils.ParallelTaskManager(len(global_list))
local_list = task_mgr.global_to_local_data(global_list)
Expand Down Expand Up @@ -1123,13 +1126,14 @@ def leaveNout_bootstrap_test(

obj, theta = self.theta_est()

bootstrap_theta = self.theta_est_bootstrap(bootstrap_samples)
bootstrap_theta = self.theta_est_bootstrap(bootstrap_samples, seed=seed)

training, test = self.confidence_region_test(
bootstrap_theta,
distribution=distribution,
alphas=alphas,
test_theta_values=theta,
seed=seed,
)

results.append((sample, test, training))
Expand Down Expand Up @@ -1287,7 +1291,7 @@ def likelihood_ratio_test(
return LR

def confidence_region_test(
self, theta_values, distribution, alphas, test_theta_values=None
self, theta_values, distribution, alphas, test_theta_values=None, seed=None
):
"""
Confidence region test to determine if theta values are within a
Expand Down Expand Up @@ -1341,6 +1345,9 @@ def confidence_region_test(
if test_theta_values is not None:
test_result = test_theta_values.copy()

if seed is not None:
np.random.seed(seed)

for a in alphas:
if distribution == 'Rect':
lb, ub = graphics.fit_rect_dist(theta_values, a)
Expand All @@ -1355,7 +1362,7 @@ def confidence_region_test(
).all(axis=1)

elif distribution == 'MVN':
dist = graphics.fit_mvn_dist(theta_values)
dist = graphics.fit_mvn_dist(theta_values, seed=seed)
Z = dist.pdf(theta_values)
score = scipy.stats.scoreatpercentile(Z, (1 - a) * 100)
training_results[a] = Z >= score
Expand All @@ -1366,7 +1373,7 @@ def confidence_region_test(
test_result[a] = Z >= score

elif distribution == 'KDE':
dist = graphics.fit_kde_dist(theta_values)
dist = graphics.fit_kde_dist(theta_values, seed=seed)
Z = dist.pdf(theta_values.transpose())
score = scipy.stats.scoreatpercentile(Z, (1 - a) * 100)
training_results[a] = Z >= score
Expand Down
24 changes: 18 additions & 6 deletions pyomo/contrib/parmest/tests/test_parmest.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -33,6 +33,10 @@
pynumero_ASL_available = AmplInterface.available()
testdir = this_file_dir()

# Set the global seed for reproducibility
_RANDOM_SEED_FOR_TESTING = 524
np.random.seed(_RANDOM_SEED_FOR_TESTING) # Set seed for reproducibility

@mrmundt mrmundt Jun 24, 2025

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You'll probably need to put a try ... except around this. You are now making numpy a hard requirement.

@adowling2 adowling2 Jun 24, 2025

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is a good suggestion.

@sscini Do the bootstrap, leave_n_out, ... functionality require numpy or scipy?

@sscini sscini Jun 24, 2025

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@adowling2 Agreed, I will add a try, except. Bootstrap and leave_N_out rely on numpy, graphics.pairwise_plot, and other graphics functions use scipy.

@sscini sscini Jun 24, 2025

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@mrmundt Would the "except" in this case be an import error for numpy? The other option would be adding the seed within each of the test classes, but if the try except works that's a lot less additions.

@mrmundt mrmundt Jun 24, 2025 •
edited
Loading

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think, considering that the random seed is necessary for these tests to run, that it would be best to do something like:

if numpy_available:
    np.random.seed(...)
else:
    raise unittest.SkipTest('Numpy is not available')

You can import numpy_available in line 25:

from pyomo.common.dependencies import numpy as np, numpy_available, pandas as pd, scipy, matplotlib

Caveat: This will make ALL parmest tests unavailable, so don't do that if you still want some tests to run.

@sscini sscini Jun 25, 2025

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@adowling2 Thoughts on this? Numpy is present in a lot of the tests here, so I think having the numpy_available would not be horrible to add.

@adowling2 adowling2 Jun 25, 2025

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It would probably be best practice to explicitly set the random seed in only the tests that need it. This gets around the caveat that @mrmundt points out. It also makes it super clear to future developers which tests/features require numpy (or scipy).

@sscini sscini Jun 25, 2025

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@adowling2 @mrmundt Thank you both for suggestions. So final action item I interpreted will be leaving the number (line 37) at the top of the file, and then adding the np.random.seed(...) to the tests that require it at the class level.



@unittest.skipIf(
not parmest.parmest_available,
Expand Down Expand Up @@ -93,7 +97,9 @@ def test_bootstrap(self):
objval, thetavals = self.pest.theta_est()

num_bootstraps = 10
theta_est = self.pest.theta_est_bootstrap(num_bootstraps, return_samples=True)
theta_est = self.pest.theta_est_bootstrap(
num_bootstraps, return_samples=True, seed=_RANDOM_SEED_FOR_TESTING
)

num_samples = theta_est["samples"].apply(len)
self.assertEqual(len(theta_est.index), 10)
Expand All @@ -109,9 +115,15 @@ def test_bootstrap(self):
self.assertEqual(CR[0.75].sum(), 7)
self.assertEqual(CR[1.0].sum(), 10) # all true

graphics.pairwise_plot(theta_est)
graphics.pairwise_plot(theta_est, thetavals)
graphics.pairwise_plot(theta_est, thetavals, 0.8, ["MVN", "KDE", "Rect"])
graphics.pairwise_plot(theta_est, seed=_RANDOM_SEED_FOR_TESTING)
graphics.pairwise_plot(theta_est, thetavals, seed=_RANDOM_SEED_FOR_TESTING)
graphics.pairwise_plot(
theta_est,
thetavals,
0.8,
["MVN", "KDE", "Rect"],
seed=_RANDOM_SEED_FOR_TESTING,
)

@unittest.skipIf(
not graphics.imports_available, "parmest.graphics imports are unavailable"
Expand Down Expand Up @@ -140,7 +152,7 @@ def test_leaveNout(self):
self.assertTrue(lNo_theta.shape == (6, 2))

results = self.pest.leaveNout_bootstrap_test(
1, None, 3, "Rect", [0.5, 1.0], seed=5436
1, None, 3, "Rect", [0.5, 1.0], seed=_RANDOM_SEED_FOR_TESTING
)
self.assertEqual(len(results), 6) # 6 lNo samples
i = 1
Expand Down Expand Up @@ -1140,7 +1152,7 @@ def test_leaveNout(self):
self.assertTrue(lNo_theta.shape == (6, 2))

results = self.pest.leaveNout_bootstrap_test(
1, None, 3, "Rect", [0.5, 1.0], seed=5436
1, None, 3, "Rect", [0.5, 1.0], seed=_RANDOM_SEED_FOR_TESTING
)
self.assertTrue(len(results) == 6) # 6 lNo samples
i = 1
Expand Down
Loading