Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
52 changes: 29 additions & 23 deletions pyomo/contrib/parmest/parmest.py
Original file line number Diff line number Diff line change
Expand Up @@ -1074,16 +1074,13 @@ def __init__(
assert isinstance(experiment_list, list)
self.exp_list = experiment_list

# get the number of experiments
self.number_exp = _count_total_experiments(self.exp_list)

# check if the experiment has a ``get_labeled_model`` function
model = _get_labeled_model(self.exp_list[0])

# check if the model has all the required suffixes
_check_model_labels(model)

# populate keyword argument options
# check if the objective function is supplied correctly
if isinstance(obj_function, str):
try:
self.obj_function = ObjectiveType(obj_function)
Expand Down Expand Up @@ -1113,6 +1110,18 @@ def __init__(
f"{[e.value for e in RegularizationType]}."
)

# get the number of data points
# this is used to compute the covariance matrix for the case
# when the measurement errors are not supplied by the user
get_measurement_error = getattr(model, "measurement_error", None)

all_unknown_errors = get_measurement_error is None or all(
model.measurement_error[y_hat] is None for y_hat in model.experiment_outputs
)

if self.obj_function == ObjectiveType.SSE and all_unknown_errors:
self.number_exp = _count_total_experiments(self.exp_list)

self.tee = tee
self.diagnostic_mode = diagnostic_mode
self.solver_options = solver_options
Expand Down Expand Up @@ -1691,7 +1700,11 @@ def _cov_at_theta(self, method, solver, step):

# fix the value of the unknown parameters to the estimated values
for param in model.unknown_parameters:
param.fix(self.estimated_theta[param.name])
if param.is_indexed():
for idx in param:
param[idx].fix(self.estimated_theta[param[idx].name])
else:
param.fix(self.estimated_theta[param.name])

# re-solve the model with the estimated parameters
results = pyo.SolverFactory(solver).solve(model, tee=self.tee)
Expand Down Expand Up @@ -1720,12 +1733,6 @@ def _cov_at_theta(self, method, solver, step):
f"The sum of squared errors at the estimated parameter(s) is: {sse}"
)

# Number of data points considered
n = self.number_exp

# Extract the number of fitted parameters
l = len(self.estimated_theta)

"""Calculate covariance assuming experimental observation errors are
independent and follow a Gaussian distribution with constant variance.

Expand Down Expand Up @@ -1763,6 +1770,17 @@ def _cov_at_theta(self, method, solver, step):

# check if the user supplied the values of the measurement errors
if all(item is None for item in meas_error):
# Number of data points considered
n = self.number_exp

# Extract the number of fitted parameters
l = len(self.estimated_theta)

assert n > l, (
"The number of datapoints must be greater than the "
"number of parameters to estimate."
)

if cov_method == CovarianceMethod.reduced_hessian:
# in the "reduced_hessian" method, use the objective value
# to calculate the measurement error variance because this
Expand Down Expand Up @@ -2027,18 +2045,6 @@ def cov_est(self, method="finite_difference", solver="ipopt", step=1e-3):
if not isinstance(step, float):
raise TypeError("Expected a float for the step, e.g., 1e-2")

# number of unknown parameters
num_unknowns = max(
[
len(_expanded_unknown_parameter_info(experiment.get_labeled_model())[0])
for experiment in self.exp_list
]
)
assert self.number_exp > num_unknowns, (
"The number of datapoints must be greater than the "
"number of parameters to estimate."
)

return self._cov_at_theta(method=method, solver=solver, step=step)

def theta_est_bootstrap(
Expand Down
3 changes: 3 additions & 0 deletions pyomo/contrib/parmest/tests/test_parmest.py
Original file line number Diff line number Diff line change
Expand Up @@ -2272,9 +2272,12 @@ def test_indexed_unknown_parameter_names_are_expanded_consistently(self):

self.assertEqual(pest._return_theta_names(), ["theta[a]", "theta[b]"])

@unittest.skipUnless(ipopt_available, "Test requires ipopt")
def test_cov_est_counts_expanded_indexed_unknown_parameters(self):
pest = _build_indexed_theta_estimator([(1.0, 2.0), (2.0, 4.0)])

obj, theta = pest.theta_est()

with self.assertRaisesRegex(
AssertionError,
"The number of datapoints must be greater than the number of parameters to estimate.",
Expand Down