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

pyomo.doe adding more verbose output for sensitivity analysis #3525

Draft
wants to merge 9 commits into
base: main
Choose a base branch
from
29 changes: 21 additions & 8 deletions pyomo/contrib/doe/doe.py
Original file line number Diff line number Diff line change
Expand Up @@ -1389,7 +1389,7 @@ def update_FIM_prior(self, model=None, FIM=None):

self.logger.info("FIM prior has been updated.")

# ToDo: Add an update function for the parameter values? --> closed loop parameter estimation?
# TODO: Add an update function for the parameter values? --> closed loop parameter estimation?
# Or leave this to the user?????
def update_unknown_parameter_values(self, model=None, param_vals=None):
raise NotImplementedError(
Expand Down Expand Up @@ -1448,15 +1448,19 @@ def compute_FIM_full_factorial(
"Design ranges keys must be a subset of experimental design names."
)

# ToDo: Add more objective types? i.e., modified-E; G-opt; V-opt; etc?
# ToDo: Also, make this a result object, or more user friendly.
# TODO: Add more objective types? i.e., modified-E; G-opt; V-opt; etc?
# TODO: Also, make this a result object, or more user friendly.
fim_factorial_results = {k.name: [] for k, v in model.experiment_inputs.items()}
fim_factorial_results.update(
{
"log10 D-opt": [],
"log10 A-opt": [],
"log10 E-opt": [],
"log10 ME-opt": [],
"eigval_min": [],
"eigval_max": [],
"det_FIM": [],
"trace_FIM": [],
"solve_time": [],
}
)
Expand Down Expand Up @@ -1521,17 +1525,22 @@ def compute_FIM_full_factorial(
# Compute and record metrics on FIM
D_opt = np.log10(np.linalg.det(FIM))
A_opt = np.log10(np.trace(FIM))
E_vals, E_vecs = np.linalg.eig(FIM) # Grab eigenvalues
E_ind = np.argmin(E_vals.real) # Grab index of minima to check imaginary
E_vals =np.linalg.eigvalsh(FIM) # np.linalg.eigvalsh(FIM) # Grab eigenvalues

#Shuvo: ????????????????????????
# Probably not required?
E_ind = 0 # np.argmin(E_vals.real) # Grab index of minima to check imaginary
# Warn the user if there is a ``large`` imaginary component (should not be)
if abs(E_vals.imag[E_ind]) > 1e-8:
self.logger.warning(
"Eigenvalue has imaginary component greater than 1e-6, contact developers if this issue persists."
"Eigenvalue has imaginary component greater than 1e-8, contact developers if this issue persists."
Copy link
Member

Choose a reason for hiding this comment

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

nit: this is why "embedded constants are evil". An improved solution would be to declare a module-level constant (maybe IMAG_THRESHOLD) and then use that both in the test and in the message.

)
#????????????????????????

# If the real value is less than or equal to zero, set the E_opt value to nan
if E_vals.real[E_ind] <= 0:
E_opt = np.nan
E_opt = np.nan # Shuvo: Also, instead of using nan, can we not use `E_opt = np.log10(abs(E_vals[0]))`
# Shuvo: if the value is very small, e.g. 1e-10?
else:
E_opt = np.log10(E_vals.real[E_ind])

Expand All @@ -1545,8 +1554,12 @@ def compute_FIM_full_factorial(
fim_factorial_results["log10 A-opt"].append(A_opt)
fim_factorial_results["log10 E-opt"].append(E_opt)
fim_factorial_results["log10 ME-opt"].append(ME_opt)
fim_factorial_results["eigval_min"].append(E_vals[0])
fim_factorial_results["eigval_max"].append(E_vals[-1])
fim_factorial_results["det_FIM"].append(np.linalg.det(FIM))
fim_factorial_results["trace_FIM"].append(np.trace(FIM))
fim_factorial_results["solve_time"].append(time_set[-1])

self.fim_factorial_results = fim_factorial_results

return self.fim_factorial_results
Expand Down