Skip to content
Draft
Changes from 9 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
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()
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