Skip to content

Conversation

@sscini
Copy link
Contributor

@sscini sscini commented Nov 20, 2025

Fixes # .

Summary/Motivation:

_Q_opt, the main model building and solving function within parmest, is still heavily dependent on MSISPPY, and the code is becoming outdated and difficult to interpret. The goal of this PR is to redesign _Q_opt to work without this dependence, retaining all the current functionality

Changes proposed in this PR:

-_Q_opt redesign, using a block structure, with minimal changes needed to functions that utilize it.

Legal Acknowledgement

By contributing to this software project, I have read the contribution guide and agree to the following terms and conditions for my contribution:

  1. I agree my contributions are submitted under the BSD license.
  2. I represent I am authorized to make the contributions and grant the license. If my employer has rights to intellectual property that includes these contributions, I represent that I have received permission to make contributions and grant the required license on behalf of that employer.

@sscini
Copy link
Contributor Author

sscini commented Nov 20, 2025

@djlaky @adowling2 Here are my initial thoughts on the redesign. Please provide feedback on code layout to this point.

@sscini
Copy link
Contributor Author

sscini commented Nov 20, 2025

Dan and I had a good conversation today about the functional things needed in the redesign, and how to go about adding components to the overall block and sub-models. Working to implement changes, closely related to what is done currently in doe's _generate_scenario_blocks.

Copy link
Member

@adowling2 adowling2 left a comment

Choose a reason for hiding this comment

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

Nice progress

@sscini
Copy link
Contributor Author

sscini commented Nov 25, 2025

@adowling2 @djlaky Added case for bootstrap, now works with theta_est, theta_est_bootstrap, and theta_est_leaveNout if I replace _Q_opt with _Q_opt_blocks. Please review.

I am aware of the duplicated code, but when I attempted to unify them using n_scenarios = len(self.exp_list) or len(bootlist), I am getting an error for bootstrap I am still working out.

Would we want to fully transition to using theta_est for obj and theta, and cov_est for covariance and add an error/warning? Or should I make it work to calculate cov if calc_cov=True?

Thanks!

@codecov
Copy link

codecov bot commented Nov 30, 2025

Codecov Report

❌ Patch coverage is 4.76190% with 100 lines in your changes missing coverage. Please review.
✅ Project coverage is 89.31%. Comparing base (13facca) to head (d91ce3f).
⚠️ Report is 21 commits behind head on main.

Files with missing lines Patch % Lines
pyomo/contrib/parmest/parmest.py 4.76% 100 Missing ⚠️
Additional details and impacted files
@@            Coverage Diff             @@
##             main    #3789      +/-   ##
==========================================
- Coverage   89.40%   89.31%   -0.10%     
==========================================
  Files         909      909              
  Lines      105541   105646     +105     
==========================================
- Hits        94364    94360       -4     
- Misses      11177    11286     +109     
Flag Coverage Δ
builders 29.07% <4.76%> (-0.04%) ⬇️
default 85.94% <4.76%> (?)
expensive 35.73% <4.76%> (?)
linux 86.63% <4.76%> (-2.55%) ⬇️
linux_other 86.63% <4.76%> (-0.10%) ⬇️
osx 82.78% <4.76%> (-0.09%) ⬇️
win 84.88% <4.76%> (-0.09%) ⬇️
win_other 84.88% <4.76%> (-0.09%) ⬇️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:
  • ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

@sscini
Copy link
Contributor Author

sscini commented Dec 9, 2025

@adowling2 @djlaky Should I tag larger team on this? Please review when available. Thanks!

@sscini sscini requested a review from adowling2 January 19, 2026 18:22
Copy link
Member

@adowling2 adowling2 left a comment

Choose a reason for hiding this comment

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

@sscini Here is my first pass at the review. A handful of places would benefit from more verbose comments to explain the logic or assumptions.

I wonder if there is a better way besides getattr to access variable values. I am thinking how we access values in Pyomo.DoE for variables that are in a suffix.

@slilonfe5, please take a look. I think this would benefit from a conversation with you to integrate the different methods to estimate the parameter covariance.

model.x2 = Param(initialize=float(data['x2']))

# Rate constants
# @Reviewers: Can we switch this to explicitly defining which parameters are to be
Copy link
Member

Choose a reason for hiding this comment

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

Yes, I think this example could/should be reworked to fully embrace the Experiment class.


# Parameter estimation with covariance
obj, theta, cov = pest.theta_est(calc_cov=True, cov_n=17)
obj, theta, cov = pest.theta_est(calc_cov=True, cov_n=19)
Copy link
Member

Choose a reason for hiding this comment

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

Just to satisfy my curiosity, is 17 incorrect and 19 correct? Is this a mistake in the example that is being fixed?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

From the data in the csv file, it should be 19. There are 19 experiments and 4 experiment outputs. The experiment counter (_count_total_experiments) is currently returning 76, which is why I have mentioned making it more robust and commented out the related assertion at this time.

# Parameter estimation with covariance
obj, theta, cov = pest.theta_est(calc_cov=True, cov_n=17)
obj, theta, cov = pest.theta_est(calc_cov=True, cov_n=19)
print(obj)
Copy link
Member

Choose a reason for hiding this comment

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

While we are editing this, let's print out something such as "print("Least squares objective = ",obj)", etc. to help a new user.

) # min^-1
model.k3 = pyo.Param(
initialize=1.0 / 6000.0, within=pyo.PositiveReals, mutable=True
# Rate constants, make unknown parameters variables
Copy link
Member

Choose a reason for hiding this comment

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

Should this example get refactored to use the Experiment class? Is this something @smondal13 is thinking about for a separate PR?

#### Redesign with Experiment class Dec 2023

# Options for using mpi-sppy or local EF only used in the deprecatedEstimator class
# TODO: move use_mpisppy to a Pyomo configuration option
Copy link
Member

Choose a reason for hiding this comment

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

Does this need to get cleaned up? Is this for the deprecated interface only?

Copy link
Contributor Author

@sscini sscini Jan 20, 2026

Choose a reason for hiding this comment

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

Will be cleaned up. Yes, currently the use_mpisppy option is only included in deprecated interface. This was added before I started this PR.

for exp_i in exp_blocks:
vals = {}
for var in return_values:
exp_i_var = exp_i.find_component(str(var))
Copy link
Member

Choose a reason for hiding this comment

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

What is this doing? More comments would help? Is it assuming specific things about the structure of the parmest model? At a minimum, I would recommend documenting any assumptions in comments at least.

var_values = pd.DataFrame(var_values)

# Calculate covariance if requested using cov_est()
if calc_cov is not NOTSET and calc_cov:
Copy link
Member

Choose a reason for hiding this comment

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

When would calc_cov be assigned NOTSET?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

[Addresses next few comments]. Need to confirm with @blnicho and team which tests need to be kept and enforced with these changes. If we call cov_est or _cov_at_theta here, then we are restraining the user to only using variables for unknown parameters, which makes multiple of the existing tests and examples fail. The inverse reduced hessian method was what @slilonfe5 implemented to calculate the covariance if a user still wishes to use theta_est for covariance. Previous conversations indicated these features should be preserved from my understanding.
Will also add more comments

Copy link
Member

Choose a reason for hiding this comment

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

If the user doesn't pass the calc_cov and cov_n arguments in theta_est, then calc_cov is NOTSET and cov_n is NOTSET.

return obj_val, theta_vals, cov
elif calc_cov is NOTSET or not calc_cov:
return obj_val, theta_vals
measurement_var = self.obj_value / (
Copy link
Member

Choose a reason for hiding this comment

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

How does this work with the changes @slilonfe5 introduced to define different objectives?

Copy link
Member

Choose a reason for hiding this comment

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

@sscini, please see my comments above for answers

return obj_val, theta_vals, var_values, cov
elif calc_cov is NOTSET or not calc_cov:
return obj_val, theta_vals, var_values
solve_result, inv_red_hes = (
Copy link
Member

Choose a reason for hiding this comment

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

Should this call the _cov_at_theta function?

Copy link
Member

Choose a reason for hiding this comment

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

I agree @adowling2. @sscini, to get the covariance matrix when the user employs the old method of theta_est, instead of duplicating code, you should get the covariance matrix from cov_est which uses _cov_at_theta function (note _cov_at_theta is not related to and does perform similar tasks as _Q_at_theta, you can view or rename it as _cov_opt). E.g., cov = cov_est(method="reduced_hessian"). For covariance calculation using the legacy theta_est, use the "reduced_hessian" method.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

@slilonfe5 I had this implemented. Main concern here was this caused model variants tests to fail (e.g. Params, indexedParams) as we discussed. I will look at extracting the values from the parmest model and review the change with you to fix this.

Copy link
Member

Choose a reason for hiding this comment

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

@sscini I'm curious as to why they are failing, given that previously, when we've been using the "reduced_hessian" method to compute covariance using theta_est, those tests did not fail. If it's methods like "finite_difference", I will understand that the parameters need to be defined with Var.

Copy link
Member

Choose a reason for hiding this comment

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

If I run the tests on the current version of parmest that I have, those tests do not fail, and they also use the "reduced_hessian" method. Maybe something else is causing them to fail.

1
) # initialization performed using just 1 set of theta values

# print("DEBUG objective_at_theta_blocks")
Copy link
Member

Choose a reason for hiding this comment

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

Remove these commented out print statements before we merge. They are okay to keep while debugging / working on other feedback.

@slilonfe5
Copy link
Member

@sscini Here is my first pass at the review. A handful of places would benefit from more verbose comments to explain the logic or assumptions.

I wonder if there is a better way besides getattr to access variable values. I am thinking how we access values in Pyomo.DoE for variables that are in a suffix.

@slilonfe5, please take a look. I think this would benefit from a conversation with you to integrate the different methods to estimate the parameter covariance.

Hello @adowling2, I agree. Stephen and I have been talking a lot about integrating the covariance estimation, as well as help with understanding the parmest test file. I will review and leave comments and suggestions.

@sscini
Copy link
Contributor Author

sscini commented Jan 20, 2026

Notes to discuss in 01/20/26 meeting:

Does the use of parameters and indexed parameters still need to be preserved?
The _count_total_experiments function: Suggestions for path forward.

@slilonfe5
Copy link
Member

Notes to discuss in 01/20/26 meeting:

Does the use of parameters and indexed parameters still need to be preserved? The _count_total_experiments function: Suggestions for path forward.

Hello @sscini. I think you want your code to be robust to those old parmest examples, we can discuss later. I will leave suggestions on the _count_total_experiments function.

model.exp_scenarios[i].transfer_attributes_from(parmest_model)
# model.exp_scenarios[i].pprint()

# Transfer all the unknown parameters to the parent model
Copy link
Member

Choose a reason for hiding this comment

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

Are the unknown_parameters not already defined in the Experiment class? I don't understand what you're trying to do here

for index in var:
val = ThetaVals[name][index]
var[index].set_value(val)
if fix_theta:
Copy link
Member

Choose a reason for hiding this comment

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

I don't understand why we are fixing theta? Isn't theta what we are solving for in the minimization problem?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I combined _Q_opt and _Q_at_theta. _Q_opt estimates, _Q_at_theta returns objective with fixed theta values (objective_at_theta). Previously, _Q_at_theta built each scenario sequentially to avoid building extended form model, and divided by number of experiments. Now it just fixes values in extended form. It is a lot less costly, and actually made objective_at_theta a lot faster. This also lets me use it, when unfixed, for multistart estimation (next PR). I will definitely add more comments about this.

obj_val = pyo.value(ef.EF_Obj)
self.obj_value = obj_val
self.estimated_theta = theta_vals
model.Obj = pyo.Objective(rule=total_obj, sense=pyo.minimize)
Copy link
Member

Choose a reason for hiding this comment

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

If theta is fixed, I don't get how we are finding the optimal theta here.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Fixed theta, returns objective, previously used for model initialization. Will add more comments for logic.

cov, index=theta_vals.keys(), columns=theta_vals.keys()
)
# Solve model
solve_result = sol.solve(model, tee=self.tee)
Copy link
Member

Choose a reason for hiding this comment

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

You need to check the solve_result status here. You need the check to be global

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Since I am combining, the global check is below following if not fix_theta:, where the rest of the normal _Q_opt workflow takes place. The if fix_theta is meant to separate _Q_at_theta logic, which reports the worst solver status back to the user, as long as it is not infeasible.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

Projects

Development

Successfully merging this pull request may close these issues.

6 participants