-
Notifications
You must be signed in to change notification settings - Fork 563
Pyomo.DoE: Correct A optimality #3803
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.
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
base: main
Are you sure you want to change the base?
Conversation
…ample did not solve to the desired point using A-optimality. However, when initialized to small value `1e-5`, it solved to the optimal point.
…scheme, the MM model solved without initialization from outside.
…`doe/tests`` files to reflect the changes.
…n `doe/test_utils.py`
…erstanding in ``doe/test_utils.py``
…and `_Cholesky_option = False`. Add flag for functions are the only used in `trace` calculation.
|
@adowling2 @blnicho This PR is ready for initial review. |
Replace `model` with `m` in `cholesky_LLinv_imp()` Co-authored-by: Bethany Nicholson <blnicho@users.noreply.github.com>
blnicho
left a comment
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@smondal13 thanks for addressing my initial review comments. I think these changes look reasonable so far but would like to see tests for the corrected A-optimality before I'll approve the PR.
… and add `doe_example..py` based on the rooney_biegler.py
…gler/doe_example.py`
Codecov Report❌ Patch coverage is Additional details and impacted files@@ Coverage Diff @@
## main #3803 +/- ##
==========================================
+ Coverage 89.41% 89.46% +0.04%
==========================================
Files 909 907 -2
Lines 105579 105621 +42
==========================================
+ Hits 94408 94493 +85
+ Misses 11171 11128 -43
Flags with carried forward coverage won't be shown. Click here to find out more. ☔ View full report in Codecov by Sentry. 🚀 New features to boost your workflow:
|
…oe/examples/rooney_biegler_doe_example.py`
adowling2
left a comment
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Started to leave review comments. I will add more comments soon.
pyomo/contrib/doe/doe.py
Outdated
|
|
||
|
|
||
| class ObjectiveLib(Enum): | ||
| determinant = "determinant" |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Let's add the following comments after each line
# det(FIM), D-optimality
# trace(inv(FIM)), A-optimalty
# trace(FIM), "psuedo"-A-optimality
# min(eig(FIM)), E-optimality
# cond(FIM), ME-optimality
# constant zero, useful for initialization and debugging
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Added comments
adowling2
left a comment
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
From a computational perspective, pseudo-A-optimality, i.e., trace(FIM), is the easiest objective to solve. I suggest we do a majority of the testing of permutations of DoE options with pseudo-A-optimality. I suggest that we then do a fewer number of tests with D-optimality (both Cholseky and GreyBox) and A-optimality (both Cholseky and GreyBox).
| rooney_biegler_doe = DesignOfExperiments( | ||
| experiment=rooney_biegler_experiment, | ||
| objective_option="determinant", | ||
| tee=tee, | ||
| prior_FIM=FIM, | ||
| improve_cholesky_roundoff_error=True, | ||
| ) | ||
|
|
||
| if optimize_experiment_D: | ||
| # D-Optimality | ||
| rooney_biegler_doe_D = DesignOfExperiments( | ||
| experiment=rooney_biegler_experiment, | ||
| objective_option="determinant", | ||
| tee=tee, | ||
| prior_FIM=FIM, | ||
| improve_cholesky_roundoff_error=True, | ||
| ) | ||
| rooney_biegler_doe_D.run_doe() |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Why do you need both rooney_biegler_doe and rooney_biegler_doe_D? Seems like they are identical instances of the DesignOfExperiments class. Why not just call run_doe directly on rooney_biegler_doe when optimize_experiment_D is True?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
That was repeated. It's corrected now.
| opt_D = results_container["optimization"]["D"] | ||
| opt_A = results_container["optimization"]["A"] |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This will fail if either of optimize_experiment_A or optimize_experiment_D are set to False. You should either guard this code against different options sent to the run_rooney_biegler_doe function or remove the options and have this example always run DoE D/A optimization.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I have removed the plotting function. It is to visually confirm the result. It is not required.
| A_opt_value = run_rooney_biegler_doe(optimize_experiment_A=True)[ | ||
| "optimization" | ||
| ]["A"]["value"] | ||
| A_opt_design_value = run_rooney_biegler_doe(optimize_experiment_A=True)[ | ||
| "optimization" | ||
| ]["A"]["design"][0] |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Calling run_rooney_biegler_doe twice here seems inefficient. Why not save the return value in a local variable? Then you can extract any of the values in the saved results
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yeah, it was inefficient. I have saved the return as a variable and reused it.
| @unittest.skipUnless(pandas_available, "test requires pandas") | ||
| @unittest.skipUnless(matplotlib_available, "test requires matplotlib") |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The run_rooney_biegler_doe function has an option to plot or not. So why not have this test change the value of plot_optimal_design depending on the value of matplotlib_available? Then you can still test most of the example even when matplotlib isn't available.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yes, I have done that.
| # Switch backend to 'Agg' to prevent GUI windows from popping up | ||
| plt.switch_backend('Agg') |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Is this still needed? Didn't you set the backend to 'Agg` up on L28?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Not needed, I have removed that
| fd_method = "central" | ||
| obj_used = "trace" | ||
|
|
||
| experiment = run_rooney_biegler_doe()["experiment"] |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Could you rework this to use FullReactorExperiment like the other tests in this file?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Since we talked about using inexpensive tests during the IDAES summit in Pittsburgh, I have added simpler rooney_biegler example here for the test instead of FullReactorExperiment, which is more complex and sometimes fails to converge. Eventually, we want to make all tests inexpensive (or comparatively inexpensive) for Pyomo.DoE. However, it does not mean all tests that use FullReactorExperiment are expensive. I can add a test using FullReactorExperiment. Let me know whether I should add FullReactorExperiment.
| params = list(model.parameter_names) | ||
|
|
||
| # Check cov_trace initialization | ||
| cov_trace_from_fim_inv = sum(pyo.value(model.fim_inv[j, j]) for j in params) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The expected value in a test should be a hard-coded number not a repeat of the implementation.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
A hard-coded value is added.
| val = sum( | ||
| pyo.value(model.L[c, params[k]]) | ||
| * pyo.value(model.L_inv[params[k], d]) | ||
| for k in range(len(params)) | ||
| ) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
| val = sum( | |
| pyo.value(model.L[c, params[k]]) | |
| * pyo.value(model.L_inv[params[k], d]) | |
| for k in range(len(params)) | |
| ) | |
| val = pyo.value(sum(model.L[c, params[k]] * model.L_inv[params[k], d] | |
| for k in range(len(params)) | |
| )) |
I'm pretty sure you can simplify this slightly by just calling pyo.value once. (Note: if you accept this change through the GitHub GUI the line indentation is off)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Added manually.
… and remove redundant/bad code.
|
@smondal13 - Gurobi just released a new version that supports Python 3.14, and we are seeing failures for our gurobi tests because of it. I cancelled your jobs because I have a fix already, just need to get it in. TL;DR - if you see 3.14 failures in gurobi, not your fault, we're on it. |
… change-a-optimality
Fixes
Summary/Motivation:
In the old definition of A-optimality (we will call it pseudo A-optimality from now on), DoE computed the trace of the Fisher Information Matrix (FIM), as few studies used that definition. However, that is not correct, and new literature mentions the trace of the inverse of the FIM as A-optimality. This PR has adopted this new definition of A-optimality in DoE.
Changes proposed in this PR:
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: