Implementing multistart version of theta_est using multiple sampling methods #3575
Conversation
|
@djlaky @adowling2 Please provide early feedback. |
|
Dynamic saving using flush, add. |
pyomo/contrib/parmest/parmest.py
Outdated
| # # If only one restart, return an empty list | ||
| # return [] | ||
|
|
||
| # return {theta_names[i]: initial_theta[i] for i in range(len(theta_names))} |
There was a problem hiding this comment.
We discussed adding a "dataframe" sampling method that uses multistart points defined by the user. This is helpful if we want to try the same set of multistart points for multiple experiments.
| "Multistart is not supported in the deprecated parmest interface") | ||
| ) | ||
|
|
||
| assert isinstance(n_restarts, int) |
There was a problem hiding this comment.
Please look at other Pyomo code fgor exampels of throwing exceptions
There was a problem hiding this comment.
Agree with @adowling2 here, you need to throw an exception so you can test the exception is caught.
pyomo/contrib/parmest/parmest.py
Outdated
| ) | ||
|
|
||
|
|
||
| results = [] |
There was a problem hiding this comment.
It might make more sense to create a dataframe and then add rows as you go. Or you could preallocate the dataframe size because you know how many restarts.
There was a problem hiding this comment.
You could even have your generate_samples function generate this empty dataframe.
|
Extend existing tests for parmest to include multistart, add. |
|
Models provided need to include bounds, add exception |
| upper_bound = np.array([parmest_model.find_component(name).ub for name in theta_names]) | ||
| # Check if the lower and upper bounds are defined | ||
| if np.any(np.isnan(lower_bound)) or np.any(np.isnan(upper_bound)): | ||
| raise ValueError( |
There was a problem hiding this comment.
You probably already know this, but you will need to check all the errors are raised when expected.
| ) | ||
|
|
||
| if self.method == "random": | ||
| np.random.seed(seed) |
There was a problem hiding this comment.
Do you want to skip setting the random seed if seed=None (default)?
There was a problem hiding this comment.
The default is none for all the functions I use that set seed, so if it receives seed = None, it would work as expected. Would skipping it still be best practice?
pyomo/contrib/parmest/parmest.py
Outdated
| elif self.method == "latin_hypercube": | ||
| # Generate theta values using Latin hypercube sampling | ||
| sampler = scipy.stats.qmc.LatinHypercube(d=len(theta_names), seed=seed) | ||
| samples = sampler.random(n=self.n_restarts+1)[1:] # Skip the first sample |
There was a problem hiding this comment.
Why are you skipping the first sample? Please explain in the comments.
There was a problem hiding this comment.
I will add a comment in code to explain as well. The first sample generated using qmc.sobol is always the origin (zero vector). I thought logic applied to all qmc methods, but no only sobol. So to get nonzero points, you need to skip first sample
| # Generate theta values using Latin hypercube sampling | ||
| sampler = scipy.stats.qmc.LatinHypercube(d=len(theta_names), seed=seed) | ||
| samples = sampler.random(n=self.n_restarts+1)[1:] # Skip the first sample | ||
| theta_vals_multistart = np.array([lower_bound + (upper_bound - lower_bound) * theta for theta in samples]) |
There was a problem hiding this comment.
Add a comment to explain what you are doing here. I suspect that this function returns an LHS refined on [0, 1] and you are rescaling. But it is best to be explicit.
| @@ -921,6 +1020,116 @@ def theta_est( | |||
| cov_n=cov_n, | |||
| ) | |||
|
|
|||
| def theta_est_multistart( | |||
| self, | |||
| buffer=10, | |||
There was a problem hiding this comment.
Need to explain the buffer in the doc string.
pyomo/contrib/parmest/parmest.py
Outdated
| "Multistart is not supported in the deprecated parmest interface" | ||
| ) | ||
|
|
||
| assert isinstance(self.n_restarts, int) |
There was a problem hiding this comment.
Replace all of these with more descriptive error messages. Remember that we need tests for each error message.
| # This method will be used to generate the initial theta values for multistart | ||
| # optimization. It will take the theta names and the initial theta values | ||
| # and return a dictionary of theta names and their corresponding values. | ||
| def _generate_initial_theta(self, parmest_model, seed=None, n_restarts=None, multistart_sampling_method=None, user_provided=None): |
There was a problem hiding this comment.
All methods (even private) should have a description of each input, what type is expected, and what it's for in a brief yet coherent sentence.
There was a problem hiding this comment.
Also, it's probably a good idea to make sure each of the inputs are valid, in some sense, when called. For instance, if parmest_model is not a model in a specific form, there may be an error. We didn't require a model input for DoE, instead made it None and handled the case where this was sent.
| # optimization. It will take the theta names and the initial theta values | ||
| # and return a dictionary of theta names and their corresponding values. | ||
| def _generate_initial_theta(self, parmest_model, seed=None, n_restarts=None, multistart_sampling_method=None, user_provided=None): | ||
| if n_restarts == 1: |
There was a problem hiding this comment.
I like just sending a warning, and not returning. For example, n_restarts might be 1 by default. You should check if n_restarts is an int as well. Then, if n_restarts is 1, you should send a warning that the tool is intended for this number to be greater than one and solve as normal.
|
|
||
| # Get the theta names and initial theta values | ||
| theta_names = self._return_theta_names() | ||
| initial_theta = [parmest_model.find_component(name)() for name in theta_names] |
There was a problem hiding this comment.
Is it better to use the suffix for this? The suffix value shouldn't change, but the theta value may if the model has been solved for some reason. I don't know if this is a potential issue but I think that grabbing these values from the suffixes would be more dummy-proof.
| return print("No multistart optimization needed. Please use normal theta_est()") | ||
|
|
||
| # Get the theta names and initial theta values | ||
| theta_names = self._return_theta_names() |
There was a problem hiding this comment.
Is this strings? We went through the effort of trying to avoid using strings to do tasks because they are fragile. Is it possible to stick with suffixes and use the actual component/component ID?
| if multistart_sampling_method == "random": | ||
| np.random.seed(seed) | ||
| # Generate random theta values | ||
| theta_vals_multistart = np.random.uniform(lower_bound, upper_bound, size=len(theta_names)) |
There was a problem hiding this comment.
Is there any merit in allowing the user to pull from a non-uniform distribution? Unsure what extra machinery is required but maybe for a future PR/enhancement.
| "The length of theta_names and initial_theta must be the same." | ||
| ) | ||
|
|
||
| if multistart_sampling_method == "random": |
There was a problem hiding this comment.
I don't necessarily like that the first method is called "random". Aren't most of the methods inherently based on random sampling? It should probably be more explicit, like uniform sampling, or random uniform, or something like this.
| samples = sampler.random(n=n_restarts+1)[1:] | ||
| theta_vals_multistart = np.array([lower_bound + (upper_bound - lower_bound) * theta for theta in samples]) | ||
|
|
||
| elif multistart_sampling_method == "user_provided": |
There was a problem hiding this comment.
I think you probably should have a comment to start each different method that describes what you're doing. It is not entirely obvious what's going on when you haven't yet read the code.
|
|
||
| if isinstance(user_provided, np.ndarray): | ||
| # Check if the user provided numpy array has the same number of rows as the number of restarts | ||
| if user_provided.shape[0] != n_restarts: |
There was a problem hiding this comment.
For instance (from last comment), it is unclear that a 2-D numpy array is provided where the first dimension is the number of restarts and the second dimension in the number of theta names, being a N x M matrix representing the N instances of M parameters to be "multi-started".
There was a problem hiding this comment.
This is especially problematic with N = M, as the users would get results back and not know that they have made a mistake.
|
|
||
| # Add the output info values to the dataframe, starting values as nan | ||
| for i in range(len(theta_names)): | ||
| df_multistart[f'converged_{theta_names[i]}'] = np.nan |
There was a problem hiding this comment.
Are all string characters legal for pandas names? I would think so but this line seems (maybe) dangerous? For instance, what if we start getting into block structuring with multi-index parameters? We should make sure there is a test to ensure the system is robust. I believe I have this on the back burner to make one for Pyomo.DoE as well.
| theta_sampling_method: string, optional | ||
| Method used to sample theta values. Options are "random", "latin_hypercube", or "sobol". | ||
| Default is "random". | ||
| solver: string, optional |
There was a problem hiding this comment.
Is there a way to pass custom solver options? I'm not super familiar with this but if this is called separate of the standard parmest call, you need to make sure the user can specify their own solver options.
|
|
||
| # Store the results in a list or DataFrame | ||
| # depending on the number of restarts | ||
| results_df.iloc[i, len(theta_names):len(theta_names) + len(theta_names)] = converged_theta_vals[i, :] |
There was a problem hiding this comment.
It would be good to have the general format of the pandas frame described in text so lines like this can be checked by reviewers. Currently, it is very hard to check to make sure this is correct without running code ourselves.
Fixes # .
Summary/Motivation:
Currently, the optimization is only done from a single initial value. This implementation adds the ability to specify multiple initial values using selected sampling techniques: from a random uniform distribution, using Latin Hypercube Sampling, or using Sobol Quasi-Monte Carlo sampling.
Changes proposed in this PR:
TODO before converting from draft:
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: