Skip to content
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

Adding regularization feature in parmest as an option to the default SSE objective #3550

Draft
wants to merge 10 commits into
base: main
Choose a base branch
from
32 changes: 31 additions & 1 deletion pyomo/contrib/parmest/parmest.py
Original file line number Diff line number Diff line change
Expand Up @@ -234,6 +234,21 @@ def SSE(model):
expr = sum((y - y_hat) ** 2 for y, y_hat in model.experiment_outputs.items())
return expr

def regularize_term(model, prior_FIM, theta_ref):
"""
Regularization term for the objective function, which is used to penalize deviation from a
reference theta
(theta - theta_ref).transpose() * prior_FIM * (theta - theta_ref)

theta_ref: Reference parameter value, element of matrix
prior_FIM: Fisher Information Matrix from prior experimental design, matrix
theta: Parameter value, matrix

Added to SSE objective function
"""
expr = ((theta - theta_ref).transpose() * prior_FIM * (theta - theta_ref) for theta in model.unknown_parameters.items())
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 run? My intuition is that you need to write out the matrix multiplication and cannot use matrix multiplication.

return expr


class Estimator(object):
"""
Expand Down Expand Up @@ -270,6 +285,8 @@ def __init__(
self,
experiment_list,
obj_function=None,
prior_FIM=None,
theta_ref=None,
tee=False,
diagnostic_mode=False,
solver_options=None,
Expand All @@ -296,6 +313,8 @@ def __init__(

# populate keyword argument options
self.obj_function = obj_function
self.prior_FIM = prior_FIM
self.theta_ref = theta_ref
self.tee = tee
self.diagnostic_mode = diagnostic_mode
self.solver_options = solver_options
Expand Down Expand Up @@ -427,7 +446,18 @@ def _create_parmest_model(self, experiment_number):
# TODO, this needs to be turned into an enum class of options that still support
# custom functions
if self.obj_function == 'SSE':
second_stage_rule = SSE

if self.prior_FIM and self.theta_ref is not None:
# Regularize the objective function
second_stage_rule = SSE + regularize_term(prior_FIM = self.prior_FIM, theta_ref = self.theta_ref)
elif self.prior_FIM:
theta_ref = model.unknown_parameters.values()
Copy link
Member

Choose a reason for hiding this comment

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

Why is this line needed?

second_stage_rule = SSE + regularize_term(prior_FIM = self.prior_FIM, theta_ref = self.theta_ref)

else:
# Sum of squared errors
second_stage_rule = SSE

else:
# A custom function uses model.experiment_outputs as data
second_stage_rule = self.obj_function
Expand Down