From ad3c416df327e970048f8f770cf18daa7e3b4d38 Mon Sep 17 00:00:00 2001 From: Jeffrey Larson Date: Thu, 1 Feb 2024 15:23:09 -0600 Subject: [PATCH 01/24] A starting branch having modular model building/updating in pounders --- pounders/py/formquad.py | 89 ++++++++++++++++++++++++++++++++++++++ pounders/py/pounders.py | 94 ++++++++--------------------------------- 2 files changed, 107 insertions(+), 76 deletions(-) diff --git a/pounders/py/formquad.py b/pounders/py/formquad.py index 2299fc57..e2ffbd41 100644 --- a/pounders/py/formquad.py +++ b/pounders/py/formquad.py @@ -2,6 +2,8 @@ import scipy.linalg from .phi2eval import phi2eval +from .bmpts import bmpts + # from .flipFirstRow import flipFirstRow # from .flipSignQ import flipSignQ @@ -167,3 +169,90 @@ def formquad(X, F, delta, xk_in, np_max, Pars, vf): G = G / delta return [Mdir, mp, valid, G, H, Mind] + +def build_formquad_models(Cres, Hres, Model, delta, n, Low, Upp, X, F, hF, nf, xk_in, nf_max, Ffun, hfun, printf): + Res = np.zeros(np.shape(F)) + + # 1a. Compute the interpolation set. + D = X[: nf + 1] - X[xk_in] + Res[: nf + 1, :] = (F[: nf + 1, :] - Cres) - np.diagonal(0.5 * D @ (np.tensordot(D, Hres, axes=1))).T + [Mdir, mp, valid, Gres, Hresdel, Mind] = formquad(X[0 : nf + 1, :], Res[0 : nf + 1, :], delta, xk_in, Model["np_max"], Model["Par"], 0) + if mp < n: + [Mdir, mp] = bmpts(X[xk_in], Mdir[0 : n - mp, :], Low, Upp, delta, Model["Par"][2]) + for i in range(int(min(n - mp, nf_max - (nf + 1)))): + nf += 1 + X[nf] = np.minimum(Upp, np.maximum(Low, X[xk_in] + Mdir[i, :])) + F[nf] = Ffun(X[nf]) + if np.any(np.isnan(F[nf])): + raise NanValueError("NaN encountered") + hF[nf] = hfun(F[nf]) + if printf: + print("%4i Geometry point %11.5e\n" % (nf, hF[nf])) + D = Mdir[i, :] + Res[nf, :] = (F[nf, :] - Cres) - 0.5 * D @ np.tensordot(D.T, Hres, 1) + if nf + 1 >= nf_max: + raise MaxEvalError("All done") + [_, mp, valid, Gres, Hresdel, Mind] = formquad(X[0 : nf + 1, :], Res[0 : nf + 1, :], delta, xk_in, Model["np_max"], Model["Par"], False) + if mp < n: + raise ModelBuildingError("yeah") + + # 1b. Update the quadratic model + Cres = F[xk_in] + Hres = Hres + Hresdel + + return Cres, Gres, Hres, nf, X, F, hF, Mind, valid + + +def formquad_model_improvement(nf, nf_max, valid, rho, eta_1, X, F, hF, delta, xk_in, Model, Cres, Gres, Hres, combinemodels, n, Low, Upp, Ffun, hfun, printf): + Res = np.zeros(np.shape(F)) + + # Need to check because model may be valid after Xsp evaluation + [Mdir, mp, valid, _, _, _] = formquad(X[: nf + 1, :], F[: nf + 1, :], delta, xk_in, Model["np_max"], Model["Par"], 1) + if not valid: # ! One strategy for choosing model-improving point: + # Update model (exists because delta & xk_in unchanged) + D = X[: nf + 1] - X[xk_in] + Res[: nf + 1, :] = (F[: nf + 1, :] - Cres) - np.diagonal(0.5 * D @ (np.tensordot(D, Hres, axes=1))).T + [_, _, valid, Gres, Hresdel, Mind] = formquad(X[: nf + 1, :], Res[: nf + 1, :], delta, xk_in, Model["np_max"], Model["Par"], False) + Hres = Hres + Hresdel + # Update for modelimp; Cres unchanged b/c xk_in unchanged + G, H = combinemodels(Cres, Gres, Hres) + # Evaluate model-improving points to pick best one + # May eventually want to normalize Mdir first for infty norm + # Plus directions + [Mdir1, mp1] = bmpts(X[xk_in], Mdir[0 : n - mp, :], Low, Upp, delta, Model["Par"][2]) + for i in range(n - mp1): + D = Mdir1[i, :] + Res[i, 0] = D @ (G + 0.5 * H @ D.T) + b = np.argmin(Res[: n - mp1, 0:1]) + a1 = np.min(Res[: n - mp1, 0:1]) + Xsp = Mdir1[b, :] + # Minus directions + [Mdir1, mp2] = bmpts(X[xk_in], -Mdir[0 : n - mp, :], Low, Upp, delta, Model["Par"][2]) + for i in range(n - mp2): + D = Mdir1[i, :] + Res[i, 0] = D @ (G + 0.5 * H @ D.T) + b = np.argmin(Res[: n - mp2, 0:1]) + a2 = np.min(Res[: n - mp2, 0:1]) + if a2 < a1: + Xsp = Mdir1[b, :] + nf += 1 + X[nf] = np.minimum(Upp, np.maximum(Low, X[xk_in] + Xsp)) # Temp safeguard + F[nf] = Ffun(X[nf]) + if np.any(np.isnan(F[nf])): + X, F, hF, flag = prepare_outputs_before_return(X, F, hF, nf, -3) + return X, F, hF, flag, xk_in + hF[nf] = hfun(F[nf]) + if printf: + print("%4i Model point %11.5e\n" % (nf, hF[nf])) + if hF[nf] < hF[xk_in]: # ! Eventually check stuff decrease here + if printf: + print("**improvement from model point****") + # Update model to reflect new base point + D = X[nf] - X[xk_in] + xk_in = nf # Change current center + Cres = F[xk_in] + # Don't actually use + # for j in range(m): + # Gres[:, j] = Gres[:, j] + Hres[:, :, j] @ D.T + + return Cres, Gres, Hres, nf, X, F, hF, valid, xk_in diff --git a/pounders/py/pounders.py b/pounders/py/pounders.py index 0feba490..b0f96650 100644 --- a/pounders/py/pounders.py +++ b/pounders/py/pounders.py @@ -5,9 +5,14 @@ from .bmpts import bmpts from .bqmin import bqmin from .checkinputss import checkinputss -from .formquad import formquad +from .formquad import formquad, build_formquad_models, formquad_model_improvement from .prepare_outputs_before_return import prepare_outputs_before_return +class NanValueError(Exception): + pass + +class ModelBuildingError(Exception): + pass def _default_model_par_values(n): par = np.zeros(4) @@ -185,39 +190,22 @@ def pounders(Ffun, X_0, n, nf_max, g_tol, delta_0, m, Low, Upp, Prior=None, Opti nf_max = nf_max + nfs for i in range(nf + 1): hF[i] = hfun(F[i]) - Res = np.zeros(np.shape(F)) Cres = F[xk_in] Hres = np.zeros((n, n, m)) ng = np.nan # Needed for early termination, e.g., if a model is never built while nf + 1 < nf_max: - # 1a. Compute the interpolation set. - D = X[: nf + 1] - X[xk_in] - Res[: nf + 1, :] = (F[: nf + 1, :] - Cres) - np.diagonal(0.5 * D @ (np.tensordot(D, Hres, axes=1))).T - [Mdir, mp, valid, Gres, Hresdel, Mind] = formquad(X[0 : nf + 1, :], Res[0 : nf + 1, :], delta, xk_in, Model["np_max"], Model["Par"], 0) - if mp < n: - [Mdir, mp] = bmpts(X[xk_in], Mdir[0 : n - mp, :], Low, Upp, delta, Model["Par"][2]) - for i in range(int(min(n - mp, nf_max - (nf + 1)))): - nf += 1 - X[nf] = np.minimum(Upp, np.maximum(Low, X[xk_in] + Mdir[i, :])) - F[nf] = Ffun(X[nf]) - if np.any(np.isnan(F[nf])): - X, F, hF, flag = prepare_outputs_before_return(X, F, hF, nf, -3) - return X, F, hF, flag, xk_in - hF[nf] = hfun(F[nf]) - if printf: - print("%4i Geometry point %11.5e\n" % (nf, hF[nf])) - D = Mdir[i, :] - Res[nf, :] = (F[nf, :] - Cres) - 0.5 * D @ np.tensordot(D.T, Hres, 1) - if nf + 1 >= nf_max: - break - [_, mp, valid, Gres, Hresdel, Mind] = formquad(X[0 : nf + 1, :], Res[0 : nf + 1, :], delta, xk_in, Model["np_max"], Model["Par"], False) - if mp < n: + try: + [Cres, Gres, Hres, nf, X, F, hF, Mind, valid] = build_formquad_models(Cres, Hres, Model, delta, n, Low, Upp, X, F, hF, nf, xk_in, nf_max, Ffun, hfun, printf) + except NanValueError: + + X, F, hF, flag = prepare_outputs_before_return(X, F, hF, nf, -3) + return X, F, hF, flag, xk_in + + except ModelBuildingError: + X, F, hF, flag = prepare_outputs_before_return(X, F, hF, nf, -5) return X, F, hF, flag, xk_in - # 1b. Update the quadratic model - Cres = F[xk_in] - Hres = Hres + Hresdel c = hF[xk_in] G, H = combinemodels(Cres, Gres, Hres) ind_Lownotbinding = (X[xk_in] > Low) * (G.T > 0) @@ -260,7 +248,7 @@ def pounders(Ffun, X_0, n, nf_max, g_tol, delta_0, m, Low, Upp, Prior=None, Opti if nf + 1 >= nf_max: break # Recalculate gradient based on a MFN model - [_, _, valid, Gres, Hres, Mind] = formquad(X[: nf + 1, :], F[: nf + 1, :], delta, xk_in, Model["np_max"], Model["Par"], 0) + [_, _, valid, Gres, Hres, _] = formquad(X[: nf + 1, :], F[: nf + 1, :], delta, xk_in, Model["np_max"], Model["Par"], 0) G, H = combinemodels(Cres, Gres, Hres) ind_Lownotbinding = (X[xk_in] > Low) * (G.T > 0) ind_Uppnotbinding = (X[xk_in] < Upp) * (G.T < 0) @@ -335,54 +323,8 @@ def pounders(Ffun, X_0, n, nf_max, g_tol, delta_0, m, Low, Upp, Prior=None, Opti print("Warning: skipping sp soln!-----------") # 5. Evaluate a model-improving point if necessary if not valid and (nf + 1 < nf_max) and (rho < eta_1): # Implies xk_in, delta unchanged - # Need to check because model may be valid after Xsp evaluation - [Mdir, mp, valid, _, _, _] = formquad(X[: nf + 1, :], F[: nf + 1, :], delta, xk_in, Model["np_max"], Model["Par"], 1) - if not valid: # ! One strategy for choosing model-improving point: - # Update model (exists because delta & xk_in unchanged) - D = X[: nf + 1] - X[xk_in] - Res[: nf + 1, :] = (F[: nf + 1, :] - Cres) - np.diagonal(0.5 * D @ (np.tensordot(D, Hres, axes=1))).T - [_, _, valid, Gres, Hresdel, Mind] = formquad(X[: nf + 1, :], Res[: nf + 1, :], delta, xk_in, Model["np_max"], Model["Par"], False) - Hres = Hres + Hresdel - # Update for modelimp; Cres unchanged b/c xk_in unchanged - G, H = combinemodels(Cres, Gres, Hres) - # Evaluate model-improving points to pick best one - # May eventually want to normalize Mdir first for infty norm - # Plus directions - [Mdir1, mp1] = bmpts(X[xk_in], Mdir[0 : n - mp, :], Low, Upp, delta, Model["Par"][2]) - for i in range(n - mp1): - D = Mdir1[i, :] - Res[i, 0] = D @ (G + 0.5 * H @ D.T) - b = np.argmin(Res[: n - mp1, 0:1]) - a1 = np.min(Res[: n - mp1, 0:1]) - Xsp = Mdir1[b, :] - # Minus directions - [Mdir1, mp2] = bmpts(X[xk_in], -Mdir[0 : n - mp, :], Low, Upp, delta, Model["Par"][2]) - for i in range(n - mp2): - D = Mdir1[i, :] - Res[i, 0] = D @ (G + 0.5 * H @ D.T) - b = np.argmin(Res[: n - mp2, 0:1]) - a2 = np.min(Res[: n - mp2, 0:1]) - if a2 < a1: - Xsp = Mdir1[b, :] - nf += 1 - X[nf] = np.minimum(Upp, np.maximum(Low, X[xk_in] + Xsp)) # Temp safeguard - F[nf] = Ffun(X[nf]) - if np.any(np.isnan(F[nf])): - X, F, hF, flag = prepare_outputs_before_return(X, F, hF, nf, -3) - return X, F, hF, flag, xk_in - hF[nf] = hfun(F[nf]) - if printf: - print("%4i Model point %11.5e\n" % (nf, hF[nf])) - if hF[nf] < hF[xk_in]: # ! Eventually check stuff decrease here - if printf: - print("**improvement from model point****") - # Update model to reflect new base point - D = X[nf] - X[xk_in] - xk_in = nf # Change current center - Cres = F[xk_in] - # Don't actually use - for j in range(m): - Gres[:, j] = Gres[:, j] + Hres[:, :, j] @ D.T + [Cres, Gres, Hres, nf, X, F, hF, valid, xk_in] = formquad_model_improvement(nf, nf_max, valid, rho, eta_1, X, F, hF, delta, xk_in, Model, Cres, Gres, Hres, combinemodels, n, Low, Upp, Ffun, hfun, printf) + if printf: print("Number of function evals exceeded") flag = ng From f7a16480a39f7895950748f82ab1b0cd5ed7fe98 Mon Sep 17 00:00:00 2001 From: Jeffrey Larson Date: Thu, 1 Feb 2024 15:30:49 -0600 Subject: [PATCH 02/24] Some missing arguments --- pounders/py/formquad.py | 9 +++++++-- pounders/py/pounders.py | 4 ++-- 2 files changed, 9 insertions(+), 4 deletions(-) diff --git a/pounders/py/formquad.py b/pounders/py/formquad.py index e2ffbd41..d82929bc 100644 --- a/pounders/py/formquad.py +++ b/pounders/py/formquad.py @@ -8,6 +8,11 @@ # from .flipFirstRow import flipFirstRow # from .flipSignQ import flipSignQ +class NanValueError(Exception): + pass + +class ModelBuildingError(Exception): + pass def formquad(X, F, delta, xk_in, np_max, Pars, vf): """ @@ -200,7 +205,7 @@ def build_formquad_models(Cres, Hres, Model, delta, n, Low, Upp, X, F, hF, nf, x Cres = F[xk_in] Hres = Hres + Hresdel - return Cres, Gres, Hres, nf, X, F, hF, Mind, valid + return Cres, Gres, Hres, nf, X, F, hF, Mind, valid, mp def formquad_model_improvement(nf, nf_max, valid, rho, eta_1, X, F, hF, delta, xk_in, Model, Cres, Gres, Hres, combinemodels, n, Low, Upp, Ffun, hfun, printf): @@ -255,4 +260,4 @@ def formquad_model_improvement(nf, nf_max, valid, rho, eta_1, X, F, hF, delta, x # for j in range(m): # Gres[:, j] = Gres[:, j] + Hres[:, :, j] @ D.T - return Cres, Gres, Hres, nf, X, F, hF, valid, xk_in + return Cres, Gres, Hres, nf, X, F, hF, valid, xk_in, mp diff --git a/pounders/py/pounders.py b/pounders/py/pounders.py index b0f96650..b88674c1 100644 --- a/pounders/py/pounders.py +++ b/pounders/py/pounders.py @@ -195,7 +195,7 @@ def pounders(Ffun, X_0, n, nf_max, g_tol, delta_0, m, Low, Upp, Prior=None, Opti ng = np.nan # Needed for early termination, e.g., if a model is never built while nf + 1 < nf_max: try: - [Cres, Gres, Hres, nf, X, F, hF, Mind, valid] = build_formquad_models(Cres, Hres, Model, delta, n, Low, Upp, X, F, hF, nf, xk_in, nf_max, Ffun, hfun, printf) + [Cres, Gres, Hres, nf, X, F, hF, Mind, valid, mp] = build_formquad_models(Cres, Hres, Model, delta, n, Low, Upp, X, F, hF, nf, xk_in, nf_max, Ffun, hfun, printf) except NanValueError: X, F, hF, flag = prepare_outputs_before_return(X, F, hF, nf, -3) @@ -323,7 +323,7 @@ def pounders(Ffun, X_0, n, nf_max, g_tol, delta_0, m, Low, Upp, Prior=None, Opti print("Warning: skipping sp soln!-----------") # 5. Evaluate a model-improving point if necessary if not valid and (nf + 1 < nf_max) and (rho < eta_1): # Implies xk_in, delta unchanged - [Cres, Gres, Hres, nf, X, F, hF, valid, xk_in] = formquad_model_improvement(nf, nf_max, valid, rho, eta_1, X, F, hF, delta, xk_in, Model, Cres, Gres, Hres, combinemodels, n, Low, Upp, Ffun, hfun, printf) + [Cres, Gres, Hres, nf, X, F, hF, valid, xk_in, mp] = formquad_model_improvement(nf, nf_max, valid, rho, eta_1, X, F, hF, delta, xk_in, Model, Cres, Gres, Hres, combinemodels, n, Low, Upp, Ffun, hfun, printf) if printf: print("Number of function evals exceeded") From 52ce9307e47f9eb1bc509f789cd01c9d7d79d35c Mon Sep 17 00:00:00 2001 From: Jeffrey Larson Date: Thu, 1 Feb 2024 15:47:43 -0600 Subject: [PATCH 03/24] Some missing arguments --- pounders/py/pounders.py | 15 +++------------ 1 file changed, 3 insertions(+), 12 deletions(-) diff --git a/pounders/py/pounders.py b/pounders/py/pounders.py index b88674c1..f402554c 100644 --- a/pounders/py/pounders.py +++ b/pounders/py/pounders.py @@ -5,15 +5,9 @@ from .bmpts import bmpts from .bqmin import bqmin from .checkinputss import checkinputss -from .formquad import formquad, build_formquad_models, formquad_model_improvement +from .formquad import formquad, build_formquad_models, formquad_model_improvement, NanValueError, ModelBuildingError from .prepare_outputs_before_return import prepare_outputs_before_return -class NanValueError(Exception): - pass - -class ModelBuildingError(Exception): - pass - def _default_model_par_values(n): par = np.zeros(4) par[0] = np.sqrt(n) @@ -197,14 +191,11 @@ def pounders(Ffun, X_0, n, nf_max, g_tol, delta_0, m, Low, Upp, Prior=None, Opti try: [Cres, Gres, Hres, nf, X, F, hF, Mind, valid, mp] = build_formquad_models(Cres, Hres, Model, delta, n, Low, Upp, X, F, hF, nf, xk_in, nf_max, Ffun, hfun, printf) except NanValueError: - X, F, hF, flag = prepare_outputs_before_return(X, F, hF, nf, -3) return X, F, hF, flag, xk_in - except ModelBuildingError: - - X, F, hF, flag = prepare_outputs_before_return(X, F, hF, nf, -5) - return X, F, hF, flag, xk_in + X, F, hF, flag = prepare_outputs_before_return(X, F, hF, nf, -5) + return X, F, hF, flag, xk_in c = hF[xk_in] G, H = combinemodels(Cres, Gres, Hres) From cecf0f2217852317478cf76a3ee9675ca48cb6fb Mon Sep 17 00:00:00 2001 From: Jeffrey Larson Date: Thu, 1 Feb 2024 15:48:00 -0600 Subject: [PATCH 04/24] formatting --- pounders/py/formquad.py | 6 +++++- pounders/py/pounders.py | 5 ++++- 2 files changed, 9 insertions(+), 2 deletions(-) diff --git a/pounders/py/formquad.py b/pounders/py/formquad.py index d82929bc..f7d64956 100644 --- a/pounders/py/formquad.py +++ b/pounders/py/formquad.py @@ -8,12 +8,15 @@ # from .flipFirstRow import flipFirstRow # from .flipSignQ import flipSignQ + class NanValueError(Exception): pass + class ModelBuildingError(Exception): pass + def formquad(X, F, delta, xk_in, np_max, Pars, vf): """ formquad(X, F, delta, xk_in, np_max, Pars, vf) -> [Mdir, np, valid, G, H, Mind] @@ -175,6 +178,7 @@ def formquad(X, F, delta, xk_in, np_max, Pars, vf): return [Mdir, mp, valid, G, H, Mind] + def build_formquad_models(Cres, Hres, Model, delta, n, Low, Upp, X, F, hF, nf, xk_in, nf_max, Ffun, hfun, printf): Res = np.zeros(np.shape(F)) @@ -260,4 +264,4 @@ def formquad_model_improvement(nf, nf_max, valid, rho, eta_1, X, F, hF, delta, x # for j in range(m): # Gres[:, j] = Gres[:, j] + Hres[:, :, j] @ D.T - return Cres, Gres, Hres, nf, X, F, hF, valid, xk_in, mp + return Cres, Gres, Hres, nf, X, F, hF, valid, xk_in, mp diff --git a/pounders/py/pounders.py b/pounders/py/pounders.py index f402554c..4a834d06 100644 --- a/pounders/py/pounders.py +++ b/pounders/py/pounders.py @@ -8,6 +8,7 @@ from .formquad import formquad, build_formquad_models, formquad_model_improvement, NanValueError, ModelBuildingError from .prepare_outputs_before_return import prepare_outputs_before_return + def _default_model_par_values(n): par = np.zeros(4) par[0] = np.sqrt(n) @@ -314,7 +315,9 @@ def pounders(Ffun, X_0, n, nf_max, g_tol, delta_0, m, Low, Upp, Prior=None, Opti print("Warning: skipping sp soln!-----------") # 5. Evaluate a model-improving point if necessary if not valid and (nf + 1 < nf_max) and (rho < eta_1): # Implies xk_in, delta unchanged - [Cres, Gres, Hres, nf, X, F, hF, valid, xk_in, mp] = formquad_model_improvement(nf, nf_max, valid, rho, eta_1, X, F, hF, delta, xk_in, Model, Cres, Gres, Hres, combinemodels, n, Low, Upp, Ffun, hfun, printf) + [Cres, Gres, Hres, nf, X, F, hF, valid, xk_in, mp] = formquad_model_improvement( + nf, nf_max, valid, rho, eta_1, X, F, hF, delta, xk_in, Model, Cres, Gres, Hres, combinemodels, n, Low, Upp, Ffun, hfun, printf + ) if printf: print("Number of function evals exceeded") From 59c350f1be2e8abf8c8cf83320a090db7c1cfaec Mon Sep 17 00:00:00 2001 From: Jeffrey Larson Date: Thu, 1 Feb 2024 15:50:57 -0600 Subject: [PATCH 05/24] Another exception --- pounders/py/formquad.py | 7 +++++-- pounders/py/pounders.py | 12 +++++++++--- 2 files changed, 14 insertions(+), 5 deletions(-) diff --git a/pounders/py/formquad.py b/pounders/py/formquad.py index f7d64956..bda66ee5 100644 --- a/pounders/py/formquad.py +++ b/pounders/py/formquad.py @@ -13,6 +13,10 @@ class NanValueError(Exception): pass +class MaxEvalError(Exception): + pass + + class ModelBuildingError(Exception): pass @@ -248,8 +252,7 @@ def formquad_model_improvement(nf, nf_max, valid, rho, eta_1, X, F, hF, delta, x X[nf] = np.minimum(Upp, np.maximum(Low, X[xk_in] + Xsp)) # Temp safeguard F[nf] = Ffun(X[nf]) if np.any(np.isnan(F[nf])): - X, F, hF, flag = prepare_outputs_before_return(X, F, hF, nf, -3) - return X, F, hF, flag, xk_in + raise NanValueError("NaN encountered") hF[nf] = hfun(F[nf]) if printf: print("%4i Model point %11.5e\n" % (nf, hF[nf])) diff --git a/pounders/py/pounders.py b/pounders/py/pounders.py index 4a834d06..488fa3ec 100644 --- a/pounders/py/pounders.py +++ b/pounders/py/pounders.py @@ -197,6 +197,8 @@ def pounders(Ffun, X_0, n, nf_max, g_tol, delta_0, m, Low, Upp, Prior=None, Opti except ModelBuildingError: X, F, hF, flag = prepare_outputs_before_return(X, F, hF, nf, -5) return X, F, hF, flag, xk_in + except MaxEvalError: + break c = hF[xk_in] G, H = combinemodels(Cres, Gres, Hres) @@ -315,9 +317,13 @@ def pounders(Ffun, X_0, n, nf_max, g_tol, delta_0, m, Low, Upp, Prior=None, Opti print("Warning: skipping sp soln!-----------") # 5. Evaluate a model-improving point if necessary if not valid and (nf + 1 < nf_max) and (rho < eta_1): # Implies xk_in, delta unchanged - [Cres, Gres, Hres, nf, X, F, hF, valid, xk_in, mp] = formquad_model_improvement( - nf, nf_max, valid, rho, eta_1, X, F, hF, delta, xk_in, Model, Cres, Gres, Hres, combinemodels, n, Low, Upp, Ffun, hfun, printf - ) + try: + [Cres, Gres, Hres, nf, X, F, hF, valid, xk_in, mp] = formquad_model_improvement( + nf, nf_max, valid, rho, eta_1, X, F, hF, delta, xk_in, Model, Cres, Gres, Hres, combinemodels, n, Low, Upp, Ffun, hfun, printf + ) + except NanValueError: + X, F, hF, flag = prepare_outputs_before_return(X, F, hF, nf, -3) + return X, F, hF, flag, xk_in if printf: print("Number of function evals exceeded") From 6c6503fe534d07cd43dc99b42a3be81bbe2a1f22 Mon Sep 17 00:00:00 2001 From: Jeffrey Larson Date: Thu, 1 Feb 2024 16:09:19 -0600 Subject: [PATCH 06/24] import --- pounders/py/formquad.py | 3 +-- pounders/py/pounders.py | 2 +- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/pounders/py/formquad.py b/pounders/py/formquad.py index bda66ee5..46172ead 100644 --- a/pounders/py/formquad.py +++ b/pounders/py/formquad.py @@ -1,9 +1,8 @@ import numpy as np import scipy.linalg -from .phi2eval import phi2eval from .bmpts import bmpts - +from .phi2eval import phi2eval # from .flipFirstRow import flipFirstRow # from .flipSignQ import flipSignQ diff --git a/pounders/py/pounders.py b/pounders/py/pounders.py index 488fa3ec..4dbd09a5 100644 --- a/pounders/py/pounders.py +++ b/pounders/py/pounders.py @@ -5,7 +5,7 @@ from .bmpts import bmpts from .bqmin import bqmin from .checkinputss import checkinputss -from .formquad import formquad, build_formquad_models, formquad_model_improvement, NanValueError, ModelBuildingError +from .formquad import formquad, build_formquad_models, formquad_model_improvement, NanValueError, ModelBuildingError, MaxEvalError from .prepare_outputs_before_return import prepare_outputs_before_return From 88f79119961d2aee00ab597fae10000510c0eabc Mon Sep 17 00:00:00 2001 From: Jeffrey Larson Date: Mon, 5 Feb 2024 15:55:06 -0600 Subject: [PATCH 07/24] Moving more model improvement to pounders --- pounders/py/formquad.py | 133 +++++++++++----------------------------- pounders/py/pounders.py | 78 +++++++++++++++++------ 2 files changed, 96 insertions(+), 115 deletions(-) diff --git a/pounders/py/formquad.py b/pounders/py/formquad.py index 46172ead..3efe727f 100644 --- a/pounders/py/formquad.py +++ b/pounders/py/formquad.py @@ -8,19 +8,8 @@ # from .flipSignQ import flipSignQ -class NanValueError(Exception): - pass - -class MaxEvalError(Exception): - pass - - -class ModelBuildingError(Exception): - pass - - -def formquad(X, F, delta, xk_in, np_max, Pars, vf): +def formquad(X, F, delta, xk_in, np_max, Pars, vf, H_flag=False, Old_H=[]): """ formquad(X, F, delta, xk_in, np_max, Pars, vf) -> [Mdir, np, valid, G, H, Mind] Computes the parameters for m quadratics @@ -35,13 +24,17 @@ def formquad(X, F, delta, xk_in, np_max, Pars, vf): X [dbl] [nf-by-n] Locations of evaluated points F [dbl] [nf-by-m] Function values of evaluated points delta [dbl] Positive trust region radius - xk_in [int] Index in (X and F) of the current center - np_max [int] Max # interpolation points (>=n+1) (.5*(n+1)*(n+2)) + xk_in [int] Index in (X and F) of the current center + np_max [int] Max # interpolation points (>=n+1) (.5*(n+1)*(n+2)) Pars[0] [dbl] delta multiplier for checking validity Pars[1] [dbl] delta multiplier for all interpolation points Pars[2] [dbl] Pivot threshold for validity Pars[3] [dbl] Pivot threshold for additional points (.001) vf [log] Flag indicating you just want to check model validity + H_flag [int] Flag indicating type of Hessians to return + 0: Minimum change in Frobenius norm Hessians + 1: Minimum-norm Hessians + --OUTPUTS---------------------------------------------------------------- Mdir [dbl] [(n-np+1)-by-n] Unit directions to improve model np [int] Number of interpolation points (=length(Mind)) @@ -179,91 +172,37 @@ def formquad(X, F, delta, xk_in, np_max, Pars, vf): H = H / (delta**2) G = G / delta - return [Mdir, mp, valid, G, H, Mind] - - -def build_formquad_models(Cres, Hres, Model, delta, n, Low, Upp, X, F, hF, nf, xk_in, nf_max, Ffun, hfun, printf): - Res = np.zeros(np.shape(F)) - - # 1a. Compute the interpolation set. - D = X[: nf + 1] - X[xk_in] - Res[: nf + 1, :] = (F[: nf + 1, :] - Cres) - np.diagonal(0.5 * D @ (np.tensordot(D, Hres, axes=1))).T - [Mdir, mp, valid, Gres, Hresdel, Mind] = formquad(X[0 : nf + 1, :], Res[0 : nf + 1, :], delta, xk_in, Model["np_max"], Model["Par"], 0) - if mp < n: - [Mdir, mp] = bmpts(X[xk_in], Mdir[0 : n - mp, :], Low, Upp, delta, Model["Par"][2]) - for i in range(int(min(n - mp, nf_max - (nf + 1)))): - nf += 1 - X[nf] = np.minimum(Upp, np.maximum(Low, X[xk_in] + Mdir[i, :])) - F[nf] = Ffun(X[nf]) - if np.any(np.isnan(F[nf])): - raise NanValueError("NaN encountered") - hF[nf] = hfun(F[nf]) - if printf: - print("%4i Geometry point %11.5e\n" % (nf, hF[nf])) - D = Mdir[i, :] - Res[nf, :] = (F[nf, :] - Cres) - 0.5 * D @ np.tensordot(D.T, Hres, 1) - if nf + 1 >= nf_max: - raise MaxEvalError("All done") - [_, mp, valid, Gres, Hresdel, Mind] = formquad(X[0 : nf + 1, :], Res[0 : nf + 1, :], delta, xk_in, Model["np_max"], Model["Par"], False) - if mp < n: - raise ModelBuildingError("yeah") - - # 1b. Update the quadratic model - Cres = F[xk_in] - Hres = Hres + Hresdel - - return Cres, Gres, Hres, nf, X, F, hF, Mind, valid, mp + if H_flag: + H = Old_H + H + return [Mdir, mp, valid, G, H, Mind] -def formquad_model_improvement(nf, nf_max, valid, rho, eta_1, X, F, hF, delta, xk_in, Model, Cres, Gres, Hres, combinemodels, n, Low, Upp, Ffun, hfun, printf): - Res = np.zeros(np.shape(F)) - # Need to check because model may be valid after Xsp evaluation - [Mdir, mp, valid, _, _, _] = formquad(X[: nf + 1, :], F[: nf + 1, :], delta, xk_in, Model["np_max"], Model["Par"], 1) - if not valid: # ! One strategy for choosing model-improving point: - # Update model (exists because delta & xk_in unchanged) - D = X[: nf + 1] - X[xk_in] - Res[: nf + 1, :] = (F[: nf + 1, :] - Cres) - np.diagonal(0.5 * D @ (np.tensordot(D, Hres, axes=1))).T - [_, _, valid, Gres, Hresdel, Mind] = formquad(X[: nf + 1, :], Res[: nf + 1, :], delta, xk_in, Model["np_max"], Model["Par"], False) - Hres = Hres + Hresdel - # Update for modelimp; Cres unchanged b/c xk_in unchanged - G, H = combinemodels(Cres, Gres, Hres) - # Evaluate model-improving points to pick best one - # May eventually want to normalize Mdir first for infty norm - # Plus directions - [Mdir1, mp1] = bmpts(X[xk_in], Mdir[0 : n - mp, :], Low, Upp, delta, Model["Par"][2]) - for i in range(n - mp1): - D = Mdir1[i, :] - Res[i, 0] = D @ (G + 0.5 * H @ D.T) - b = np.argmin(Res[: n - mp1, 0:1]) - a1 = np.min(Res[: n - mp1, 0:1]) +def formquad_model_improvement(x_k, Cres, Gres, Hres, Mdir, mp, Low, Upp, delta, Model, combinemodels): + n = len(x_k) + + # Update for modelimp; Cres unchanged b/c xk_in unchanged + G, H = combinemodels(Cres, Gres, Hres) + # Evaluate model-improving points to pick best one + # May eventually want to normalize Mdir first for infty norm + # Plus directions + [Mdir1, mp1] = bmpts(x_k, Mdir[0 : n - mp, :], Low, Upp, delta, Model["Par"][2]) + Res = np.zeros((n - mp1, 1)) + for i in range(n - mp1): + D = Mdir1[i, :] + Res[i,0] = D @ (G + 0.5 * H @ D.T) + b = np.argmin(Res[: n - mp1, 0:1]) + a1 = np.min(Res[: n - mp1, 0:1]) + Xsp = Mdir1[b, :] + # Minus directions + [Mdir1, mp2] = bmpts(x_k, -Mdir[0 : n - mp, :], Low, Upp, delta, Model["Par"][2]) + Res = np.zeros((n - mp2, 1)) + for i in range(n - mp2): + D = Mdir1[i, :] + Res[i,0] = D @ (G + 0.5 * H @ D.T) + b = np.argmin(Res[: n - mp2, 0:1]) + a2 = np.min(Res[: n - mp2, 0:1]) + if a2 < a1: Xsp = Mdir1[b, :] - # Minus directions - [Mdir1, mp2] = bmpts(X[xk_in], -Mdir[0 : n - mp, :], Low, Upp, delta, Model["Par"][2]) - for i in range(n - mp2): - D = Mdir1[i, :] - Res[i, 0] = D @ (G + 0.5 * H @ D.T) - b = np.argmin(Res[: n - mp2, 0:1]) - a2 = np.min(Res[: n - mp2, 0:1]) - if a2 < a1: - Xsp = Mdir1[b, :] - nf += 1 - X[nf] = np.minimum(Upp, np.maximum(Low, X[xk_in] + Xsp)) # Temp safeguard - F[nf] = Ffun(X[nf]) - if np.any(np.isnan(F[nf])): - raise NanValueError("NaN encountered") - hF[nf] = hfun(F[nf]) - if printf: - print("%4i Model point %11.5e\n" % (nf, hF[nf])) - if hF[nf] < hF[xk_in]: # ! Eventually check stuff decrease here - if printf: - print("**improvement from model point****") - # Update model to reflect new base point - D = X[nf] - X[xk_in] - xk_in = nf # Change current center - Cres = F[xk_in] - # Don't actually use - # for j in range(m): - # Gres[:, j] = Gres[:, j] + Hres[:, :, j] @ D.T - return Cres, Gres, Hres, nf, X, F, hF, valid, xk_in, mp + return Xsp diff --git a/pounders/py/pounders.py b/pounders/py/pounders.py index 4dbd09a5..8cbd8742 100644 --- a/pounders/py/pounders.py +++ b/pounders/py/pounders.py @@ -5,7 +5,7 @@ from .bmpts import bmpts from .bqmin import bqmin from .checkinputss import checkinputss -from .formquad import formquad, build_formquad_models, formquad_model_improvement, NanValueError, ModelBuildingError, MaxEvalError +from .formquad import formquad, formquad_model_improvement from .prepare_outputs_before_return import prepare_outputs_before_return @@ -109,9 +109,12 @@ def pounders(Ffun, X_0, n, nf_max, g_tol, delta_0, m, Low, Upp, Prior=None, Opti if Model is None: Model = {} + Model["H_flag"] = True Model["Par"] = _default_model_par_values(n) Model["np_max"] = _default_model_np_max(n) else: + if "H_flag" not in Model: + Model["H_flag"] = True if "Par" not in Model: Model["Par"] = _default_model_par_values(n) if "np_max" not in Model: @@ -185,20 +188,36 @@ def pounders(Ffun, X_0, n, nf_max, g_tol, delta_0, m, Low, Upp, Prior=None, Opti nf_max = nf_max + nfs for i in range(nf + 1): hF[i] = hfun(F[i]) + Res = np.zeros(np.shape(F)) Cres = F[xk_in] Hres = np.zeros((n, n, m)) ng = np.nan # Needed for early termination, e.g., if a model is never built while nf + 1 < nf_max: - try: - [Cres, Gres, Hres, nf, X, F, hF, Mind, valid, mp] = build_formquad_models(Cres, Hres, Model, delta, n, Low, Upp, X, F, hF, nf, xk_in, nf_max, Ffun, hfun, printf) - except NanValueError: - X, F, hF, flag = prepare_outputs_before_return(X, F, hF, nf, -3) - return X, F, hF, flag, xk_in - except ModelBuildingError: - X, F, hF, flag = prepare_outputs_before_return(X, F, hF, nf, -5) - return X, F, hF, flag, xk_in - except MaxEvalError: - break + # 1a. Compute displacements + D = X[: nf + 1] - X[xk_in] + Res[: nf + 1, :] = (F[: nf + 1, :] - Cres) - np.diagonal(0.5 * D @ (np.tensordot(D, Hres, axes=1))).T + + # Request directions for points to be evaluated + [Mdir, mp, valid, Gres, Hresdel, Mind] = formquad(X[0 : nf + 1, :], Res[0 : nf + 1, :], delta, xk_in, Model["np_max"], Model["Par"], False, Model["H_flag"], Hres) + + if mp < n: + [Mdir, mp] = bmpts(X[xk_in], Mdir[0 : n - mp, :], Low, Upp, delta, Model["Par"][2]) + for i in range(int(min(n - mp, nf_max - (nf + 1)))): + nf += 1 + X[nf] = np.minimum(Upp, np.maximum(Low, X[xk_in] + Mdir[i, :])) + F[nf] = Ffun(X[nf]) + if np.any(np.isnan(F[nf])): + X, F, hF, flag = prepare_outputs_before_return(X, F, hF, nf, -3) + return X, F, hF, flag, xk_in + hF[nf] = hfun(F[nf]) + if printf: + print("%4i Geometry point %11.5e\n" % (nf, hF[nf])) + D = Mdir[i, :] + Res[nf, :] = (F[nf, :] - Cres) - 0.5 * D @ np.tensordot(D.T, Hres, 1) + if nf + 1 >= nf_max: + break + + [_, mp, valid, Gres, Hresdel, Mind] = formquad(X[0 : nf + 1, :], Res[0 : nf + 1, :], delta, xk_in, Model["np_max"], Model["Par"], False, Model["H_flag"], Hres) c = hF[xk_in] G, H = combinemodels(Cres, Gres, Hres) @@ -317,13 +336,36 @@ def pounders(Ffun, X_0, n, nf_max, g_tol, delta_0, m, Low, Upp, Prior=None, Opti print("Warning: skipping sp soln!-----------") # 5. Evaluate a model-improving point if necessary if not valid and (nf + 1 < nf_max) and (rho < eta_1): # Implies xk_in, delta unchanged - try: - [Cres, Gres, Hres, nf, X, F, hF, valid, xk_in, mp] = formquad_model_improvement( - nf, nf_max, valid, rho, eta_1, X, F, hF, delta, xk_in, Model, Cres, Gres, Hres, combinemodels, n, Low, Upp, Ffun, hfun, printf - ) - except NanValueError: - X, F, hF, flag = prepare_outputs_before_return(X, F, hF, nf, -3) - return X, F, hF, flag, xk_in + # Need to check because model may be valid after Xsp evaluation + [Mdir, mp, valid, _, _, _] = formquad(X[: nf + 1, :], F[: nf + 1, :], delta, xk_in, Model["np_max"], Model["Par"], True) + + if not valid: # ! One strategy for choosing model-improving point: + # Update model (exists because delta & xk_in unchanged) + D = X[: nf + 1] - X[xk_in] + Res[: nf + 1, :] = (F[: nf + 1, :] - Cres) - np.diagonal(0.5 * D @ (np.tensordot(D, Hres, axes=1))).T + + [_, _, valid, Gres, Hresdel, Mind] = formquad(X[: nf + 1, :], Res[: nf + 1, :], delta, xk_in, Model["np_max"], Model["Par"], False, Model["H_flag"], Hres) + + Xsp = formquad_model_improvement(X[xk_in], Cres, Gres, Hres, Mdir, mp, Low, Upp, delta, Model, combinemodels) + + nf += 1 + X[nf] = np.minimum(Upp, np.maximum(Low, X[xk_in] + Xsp)) # Temp safeguard + F[nf] = Ffun(X[nf]) + if np.any(np.isnan(F[nf])): + raise NanValueError("NaN encountered") + hF[nf] = hfun(F[nf]) + if printf: + print("%4i Model point %11.5e\n" % (nf, hF[nf])) + if hF[nf] < hF[xk_in]: # ! Eventually check stuff decrease here + if printf: + print("**improvement from model point****") + # Update model to reflect new base point + D = X[nf] - X[xk_in] + xk_in = nf # Change current center + Cres = F[xk_in] + # Don't actually use + # for j in range(m): + # Gres[:, j] = Gres[:, j] + Hres[:, :, j] @ D.T if printf: print("Number of function evals exceeded") From 46b933fd5585fff66d21c3e129c3469cd411151c Mon Sep 17 00:00:00 2001 From: Jeffrey Larson Date: Mon, 5 Feb 2024 15:55:20 -0600 Subject: [PATCH 08/24] format --- pounders/py/formquad.py | 5 ++--- pounders/py/pounders.py | 2 +- 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/pounders/py/formquad.py b/pounders/py/formquad.py index 3efe727f..fa93a022 100644 --- a/pounders/py/formquad.py +++ b/pounders/py/formquad.py @@ -8,7 +8,6 @@ # from .flipSignQ import flipSignQ - def formquad(X, F, delta, xk_in, np_max, Pars, vf, H_flag=False, Old_H=[]): """ formquad(X, F, delta, xk_in, np_max, Pars, vf) -> [Mdir, np, valid, G, H, Mind] @@ -190,7 +189,7 @@ def formquad_model_improvement(x_k, Cres, Gres, Hres, Mdir, mp, Low, Upp, delta, Res = np.zeros((n - mp1, 1)) for i in range(n - mp1): D = Mdir1[i, :] - Res[i,0] = D @ (G + 0.5 * H @ D.T) + Res[i, 0] = D @ (G + 0.5 * H @ D.T) b = np.argmin(Res[: n - mp1, 0:1]) a1 = np.min(Res[: n - mp1, 0:1]) Xsp = Mdir1[b, :] @@ -199,7 +198,7 @@ def formquad_model_improvement(x_k, Cres, Gres, Hres, Mdir, mp, Low, Upp, delta, Res = np.zeros((n - mp2, 1)) for i in range(n - mp2): D = Mdir1[i, :] - Res[i,0] = D @ (G + 0.5 * H @ D.T) + Res[i, 0] = D @ (G + 0.5 * H @ D.T) b = np.argmin(Res[: n - mp2, 0:1]) a2 = np.min(Res[: n - mp2, 0:1]) if a2 < a1: diff --git a/pounders/py/pounders.py b/pounders/py/pounders.py index 8cbd8742..0db1d078 100644 --- a/pounders/py/pounders.py +++ b/pounders/py/pounders.py @@ -113,7 +113,7 @@ def pounders(Ffun, X_0, n, nf_max, g_tol, delta_0, m, Low, Upp, Prior=None, Opti Model["Par"] = _default_model_par_values(n) Model["np_max"] = _default_model_np_max(n) else: - if "H_flag" not in Model: + if "H_flag" not in Model: Model["H_flag"] = True if "Par" not in Model: Model["Par"] = _default_model_par_values(n) From 76a05c54f305b2a34c6b586f81b01402adff05ee Mon Sep 17 00:00:00 2001 From: Jeffrey Larson Date: Mon, 5 Feb 2024 15:56:21 -0600 Subject: [PATCH 09/24] Reverting nan catching --- pounders/py/pounders.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pounders/py/pounders.py b/pounders/py/pounders.py index 0db1d078..6cb9cbf8 100644 --- a/pounders/py/pounders.py +++ b/pounders/py/pounders.py @@ -352,7 +352,8 @@ def pounders(Ffun, X_0, n, nf_max, g_tol, delta_0, m, Low, Upp, Prior=None, Opti X[nf] = np.minimum(Upp, np.maximum(Low, X[xk_in] + Xsp)) # Temp safeguard F[nf] = Ffun(X[nf]) if np.any(np.isnan(F[nf])): - raise NanValueError("NaN encountered") + X, F, hF, flag = prepare_outputs_before_return(X, F, hF, nf, -3) + return X, F, flag, xk_in hF[nf] = hfun(F[nf]) if printf: print("%4i Model point %11.5e\n" % (nf, hF[nf])) From c9354312a178c32f0190a1a22ec5bb7d9c54f1de Mon Sep 17 00:00:00 2001 From: Jeffrey Larson Date: Mon, 5 Feb 2024 15:59:27 -0600 Subject: [PATCH 10/24] Reverting nan catching --- pounders/py/pounders.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/pounders/py/pounders.py b/pounders/py/pounders.py index 6cb9cbf8..1dc21c82 100644 --- a/pounders/py/pounders.py +++ b/pounders/py/pounders.py @@ -218,6 +218,9 @@ def pounders(Ffun, X_0, n, nf_max, g_tol, delta_0, m, Low, Upp, Prior=None, Opti break [_, mp, valid, Gres, Hresdel, Mind] = formquad(X[0 : nf + 1, :], Res[0 : nf + 1, :], delta, xk_in, Model["np_max"], Model["Par"], False, Model["H_flag"], Hres) + if mp < n: + X, F, hF, flag = prepare_outputs_before_return(X, F, hF, nf, -5) + return X, F, hF, flag, xk_in c = hF[xk_in] G, H = combinemodels(Cres, Gres, Hres) @@ -254,7 +257,7 @@ def pounders(Ffun, X_0, n, nf_max, g_tol, delta_0, m, Low, Upp, Prior=None, Opti F[nf] = Ffun(X[nf]) if np.any(np.isnan(F[nf])): X, F, hF, flag = prepare_outputs_before_return(X, F, hF, nf, -3) - return X, F, flag, xk_in + return X, F, hF, flag, xk_in hF[nf] = hfun(F[nf]) if printf: print("%4i Critical point %11.5e\n" % (nf, hF[nf])) @@ -353,7 +356,7 @@ def pounders(Ffun, X_0, n, nf_max, g_tol, delta_0, m, Low, Upp, Prior=None, Opti F[nf] = Ffun(X[nf]) if np.any(np.isnan(F[nf])): X, F, hF, flag = prepare_outputs_before_return(X, F, hF, nf, -3) - return X, F, flag, xk_in + return X, F, hF, flag, xk_in hF[nf] = hfun(F[nf]) if printf: print("%4i Model point %11.5e\n" % (nf, hF[nf])) From e26bcdf3afa8c3104ec736eb62b372d6847902c1 Mon Sep 17 00:00:00 2001 From: Jeffrey Larson Date: Mon, 5 Feb 2024 16:02:10 -0600 Subject: [PATCH 11/24] Decreasing diffs --- pounders/py/pounders.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pounders/py/pounders.py b/pounders/py/pounders.py index 1dc21c82..6eb473b2 100644 --- a/pounders/py/pounders.py +++ b/pounders/py/pounders.py @@ -368,8 +368,8 @@ def pounders(Ffun, X_0, n, nf_max, g_tol, delta_0, m, Low, Upp, Prior=None, Opti xk_in = nf # Change current center Cres = F[xk_in] # Don't actually use - # for j in range(m): - # Gres[:, j] = Gres[:, j] + Hres[:, :, j] @ D.T + for j in range(m): + Gres[:, j] = Gres[:, j] + Hres[:, :, j] @ D.T if printf: print("Number of function evals exceeded") From ed19f7d28cddf7dbcacbfa2a59946f6708169412 Mon Sep 17 00:00:00 2001 From: Jeffrey Larson Date: Fri, 16 Feb 2024 08:23:57 -0600 Subject: [PATCH 12/24] Hres retrns --- pounders/py/pounders.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pounders/py/pounders.py b/pounders/py/pounders.py index 6eb473b2..21a5d460 100644 --- a/pounders/py/pounders.py +++ b/pounders/py/pounders.py @@ -198,7 +198,7 @@ def pounders(Ffun, X_0, n, nf_max, g_tol, delta_0, m, Low, Upp, Prior=None, Opti Res[: nf + 1, :] = (F[: nf + 1, :] - Cres) - np.diagonal(0.5 * D @ (np.tensordot(D, Hres, axes=1))).T # Request directions for points to be evaluated - [Mdir, mp, valid, Gres, Hresdel, Mind] = formquad(X[0 : nf + 1, :], Res[0 : nf + 1, :], delta, xk_in, Model["np_max"], Model["Par"], False, Model["H_flag"], Hres) + [Mdir, mp, valid, Gres, Hres, Mind] = formquad(X[0 : nf + 1, :], Res[0 : nf + 1, :], delta, xk_in, Model["np_max"], Model["Par"], False, Model["H_flag"], Hres) if mp < n: [Mdir, mp] = bmpts(X[xk_in], Mdir[0 : n - mp, :], Low, Upp, delta, Model["Par"][2]) @@ -217,7 +217,7 @@ def pounders(Ffun, X_0, n, nf_max, g_tol, delta_0, m, Low, Upp, Prior=None, Opti if nf + 1 >= nf_max: break - [_, mp, valid, Gres, Hresdel, Mind] = formquad(X[0 : nf + 1, :], Res[0 : nf + 1, :], delta, xk_in, Model["np_max"], Model["Par"], False, Model["H_flag"], Hres) + [_, mp, valid, Gres, Hres, Mind] = formquad(X[0 : nf + 1, :], Res[0 : nf + 1, :], delta, xk_in, Model["np_max"], Model["Par"], False, Model["H_flag"], Hres) if mp < n: X, F, hF, flag = prepare_outputs_before_return(X, F, hF, nf, -5) return X, F, hF, flag, xk_in @@ -347,7 +347,7 @@ def pounders(Ffun, X_0, n, nf_max, g_tol, delta_0, m, Low, Upp, Prior=None, Opti D = X[: nf + 1] - X[xk_in] Res[: nf + 1, :] = (F[: nf + 1, :] - Cres) - np.diagonal(0.5 * D @ (np.tensordot(D, Hres, axes=1))).T - [_, _, valid, Gres, Hresdel, Mind] = formquad(X[: nf + 1, :], Res[: nf + 1, :], delta, xk_in, Model["np_max"], Model["Par"], False, Model["H_flag"], Hres) + [_, _, valid, Gres, Hres, Mind] = formquad(X[: nf + 1, :], Res[: nf + 1, :], delta, xk_in, Model["np_max"], Model["Par"], False, Model["H_flag"], Hres) Xsp = formquad_model_improvement(X[xk_in], Cres, Gres, Hres, Mdir, mp, Low, Upp, delta, Model, combinemodels) From 0dc30d2a3c56210c4deeaa73aa439024d756433b Mon Sep 17 00:00:00 2001 From: Jeffrey Larson Date: Mon, 18 Mar 2024 15:11:38 -0500 Subject: [PATCH 13/24] not setting H to be empty --- pounders/py/formquad.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pounders/py/formquad.py b/pounders/py/formquad.py index fa93a022..b0a06cc9 100644 --- a/pounders/py/formquad.py +++ b/pounders/py/formquad.py @@ -98,7 +98,7 @@ def formquad(X, F, delta, xk_in, np_max, Pars, vf, H_flag=False, Old_H=[]): elif aff == 1 and mp < n: # Need to evaluate more points, then recall Mdir = Q[:, mp:n].T # Output Geometry directions G = np.empty(shape=(0, 0)) - H = np.empty(shape=(0, 0)) + # H = np.empty(shape=(0, 0)) return [Mdir, mp, valid, G, H, Mind] elif aff == 0: # Output model-improving directions Mdir = Q[:, mp:n].T # Will be empty if mp=n From 72241574448efc04f20e45b50deac25d98bb82b2 Mon Sep 17 00:00:00 2001 From: Jeffrey Larson Date: Thu, 18 Apr 2024 13:27:52 -0500 Subject: [PATCH 14/24] Whitespace --- pounders/py/pounders.py | 5 ----- 1 file changed, 5 deletions(-) diff --git a/pounders/py/pounders.py b/pounders/py/pounders.py index 21a5d460..b17e0759 100644 --- a/pounders/py/pounders.py +++ b/pounders/py/pounders.py @@ -196,8 +196,6 @@ def pounders(Ffun, X_0, n, nf_max, g_tol, delta_0, m, Low, Upp, Prior=None, Opti # 1a. Compute displacements D = X[: nf + 1] - X[xk_in] Res[: nf + 1, :] = (F[: nf + 1, :] - Cres) - np.diagonal(0.5 * D @ (np.tensordot(D, Hres, axes=1))).T - - # Request directions for points to be evaluated [Mdir, mp, valid, Gres, Hres, Mind] = formquad(X[0 : nf + 1, :], Res[0 : nf + 1, :], delta, xk_in, Model["np_max"], Model["Par"], False, Model["H_flag"], Hres) if mp < n: @@ -216,7 +214,6 @@ def pounders(Ffun, X_0, n, nf_max, g_tol, delta_0, m, Low, Upp, Prior=None, Opti Res[nf, :] = (F[nf, :] - Cres) - 0.5 * D @ np.tensordot(D.T, Hres, 1) if nf + 1 >= nf_max: break - [_, mp, valid, Gres, Hres, Mind] = formquad(X[0 : nf + 1, :], Res[0 : nf + 1, :], delta, xk_in, Model["np_max"], Model["Par"], False, Model["H_flag"], Hres) if mp < n: X, F, hF, flag = prepare_outputs_before_return(X, F, hF, nf, -5) @@ -346,7 +343,6 @@ def pounders(Ffun, X_0, n, nf_max, g_tol, delta_0, m, Low, Upp, Prior=None, Opti # Update model (exists because delta & xk_in unchanged) D = X[: nf + 1] - X[xk_in] Res[: nf + 1, :] = (F[: nf + 1, :] - Cres) - np.diagonal(0.5 * D @ (np.tensordot(D, Hres, axes=1))).T - [_, _, valid, Gres, Hres, Mind] = formquad(X[: nf + 1, :], Res[: nf + 1, :], delta, xk_in, Model["np_max"], Model["Par"], False, Model["H_flag"], Hres) Xsp = formquad_model_improvement(X[xk_in], Cres, Gres, Hres, Mdir, mp, Low, Upp, delta, Model, combinemodels) @@ -370,7 +366,6 @@ def pounders(Ffun, X_0, n, nf_max, g_tol, delta_0, m, Low, Upp, Prior=None, Opti # Don't actually use for j in range(m): Gres[:, j] = Gres[:, j] + Hres[:, :, j] @ D.T - if printf: print("Number of function evals exceeded") flag = ng From aa63bda9b04309ff631ef1f1f683e298488b8a61 Mon Sep 17 00:00:00 2001 From: Jeffrey Larson Date: Thu, 18 Apr 2024 13:43:41 -0500 Subject: [PATCH 15/24] Doc fix --- pounders/py/formquad.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pounders/py/formquad.py b/pounders/py/formquad.py index b0a06cc9..af5127dd 100644 --- a/pounders/py/formquad.py +++ b/pounders/py/formquad.py @@ -31,8 +31,8 @@ def formquad(X, F, delta, xk_in, np_max, Pars, vf, H_flag=False, Old_H=[]): Pars[3] [dbl] Pivot threshold for additional points (.001) vf [log] Flag indicating you just want to check model validity H_flag [int] Flag indicating type of Hessians to return - 0: Minimum change in Frobenius norm Hessians - 1: Minimum-norm Hessians + 0: Minimum-norm Hessians + 1: Minimum change in Frobenius norm Hessians --OUTPUTS---------------------------------------------------------------- Mdir [dbl] [(n-np+1)-by-n] Unit directions to improve model From e4fb842c2706a78453d04a31950f385310602aca Mon Sep 17 00:00:00 2001 From: Jeffrey Larson Date: Thu, 18 Apr 2024 13:50:14 -0500 Subject: [PATCH 16/24] Easier to read --- pounders/py/pounders.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pounders/py/pounders.py b/pounders/py/pounders.py index b17e0759..24fbdb64 100644 --- a/pounders/py/pounders.py +++ b/pounders/py/pounders.py @@ -242,10 +242,10 @@ def pounders(Ffun, X_0, n, nf_max, g_tol, delta_0, m, Low, Upp, Prior=None, Opti jerr[i, j] = (Cres[j] - F[Mind[i], j]) + D @ (Gres[:, j] + 0.5 * Hres[:, :, j] @ D) print(jerr) # input("Enter a key and press Enter to continue\n") - Don't uncomment when using Pytest with test_pounders.py - # 2. Critically test invoked if the projected model gradient is small + # 2. Criticality test invoked if the projected model gradient is small if ng < g_tol: delta = max(g_tol, np.max(np.abs(X[xk_in])) * eps) - [Mdir, _, valid, _, _, _] = formquad(X[: nf + 1, :], F[: nf + 1, :], delta, xk_in, Model["np_max"], Model["Par"], 1) + [Mdir, _, valid, _, _, _] = formquad(X[: nf + 1, :], F[: nf + 1, :], delta, xk_in, Model["np_max"], Model["Par"], True) # Not passing in H_flag, so we are starting from scratch to make sure our geometry is good (i.e., not concerned with past modesl) so we don't pass in old_H, even if that Model["H_flag"] is True. if not valid: [Mdir, mp] = bmpts(X[xk_in], Mdir, Low, Upp, delta, Model["Par"][2]) for i in range(min(n - mp, nf_max - (nf + 1))): @@ -261,7 +261,7 @@ def pounders(Ffun, X_0, n, nf_max, g_tol, delta_0, m, Low, Upp, Prior=None, Opti if nf + 1 >= nf_max: break # Recalculate gradient based on a MFN model - [_, _, valid, Gres, Hres, _] = formquad(X[: nf + 1, :], F[: nf + 1, :], delta, xk_in, Model["np_max"], Model["Par"], 0) + [_, _, valid, Gres, Hres, _] = formquad(X[: nf + 1, :], F[: nf + 1, :], delta, xk_in, Model["np_max"], Model["Par"], False) G, H = combinemodels(Cres, Gres, Hres) ind_Lownotbinding = (X[xk_in] > Low) * (G.T > 0) ind_Uppnotbinding = (X[xk_in] < Upp) * (G.T < 0) From da05502742dff519701319ec643558e6f673d0bd Mon Sep 17 00:00:00 2001 From: Jeffrey Larson Date: Thu, 18 Apr 2024 13:52:06 -0500 Subject: [PATCH 17/24] comment --- pounders/py/pounders.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/pounders/py/pounders.py b/pounders/py/pounders.py index 24fbdb64..37a602c3 100644 --- a/pounders/py/pounders.py +++ b/pounders/py/pounders.py @@ -245,7 +245,11 @@ def pounders(Ffun, X_0, n, nf_max, g_tol, delta_0, m, Low, Upp, Prior=None, Opti # 2. Criticality test invoked if the projected model gradient is small if ng < g_tol: delta = max(g_tol, np.max(np.abs(X[xk_in])) * eps) - [Mdir, _, valid, _, _, _] = formquad(X[: nf + 1, :], F[: nf + 1, :], delta, xk_in, Model["np_max"], Model["Par"], True) # Not passing in H_flag, so we are starting from scratch to make sure our geometry is good (i.e., not concerned with past modesl) so we don't pass in old_H, even if that Model["H_flag"] is True. + # Not passing in H_flag, because we are starting from scratch and + # we want to make sure our geometry is good (i.e., not concerned + # with past models) so we don't pass in old_H, even if + # Model["H_flag"] is True. + [Mdir, _, valid, _, _, _] = formquad(X[: nf + 1, :], F[: nf + 1, :], delta, xk_in, Model["np_max"], Model["Par"], True) if not valid: [Mdir, mp] = bmpts(X[xk_in], Mdir, Low, Upp, delta, Model["Par"][2]) for i in range(min(n - mp, nf_max - (nf + 1))): From bd2798025537499ee6c96f1492411e953bc2b2a5 Mon Sep 17 00:00:00 2001 From: Jeffrey Larson Date: Mon, 22 Apr 2024 10:50:48 -0500 Subject: [PATCH 18/24] black --- pounders/py/pounders.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pounders/py/pounders.py b/pounders/py/pounders.py index 37a602c3..ad2be008 100644 --- a/pounders/py/pounders.py +++ b/pounders/py/pounders.py @@ -249,7 +249,7 @@ def pounders(Ffun, X_0, n, nf_max, g_tol, delta_0, m, Low, Upp, Prior=None, Opti # we want to make sure our geometry is good (i.e., not concerned # with past models) so we don't pass in old_H, even if # Model["H_flag"] is True. - [Mdir, _, valid, _, _, _] = formquad(X[: nf + 1, :], F[: nf + 1, :], delta, xk_in, Model["np_max"], Model["Par"], True) + [Mdir, _, valid, _, _, _] = formquad(X[: nf + 1, :], F[: nf + 1, :], delta, xk_in, Model["np_max"], Model["Par"], True) if not valid: [Mdir, mp] = bmpts(X[xk_in], Mdir, Low, Upp, delta, Model["Par"][2]) for i in range(min(n - mp, nf_max - (nf + 1))): From 7cc3a9bb0f7235f2b3e33f8eb9cf25f07f1935cb Mon Sep 17 00:00:00 2001 From: Matt Menickelly Date: Wed, 24 Apr 2024 15:29:27 -0500 Subject: [PATCH 19/24] model_builder is now part of the Model struct and generalizes formquad --- pounders/m/pounders.m | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/pounders/m/pounders.m b/pounders/m/pounders.m index 8f696030..2ab54596 100644 --- a/pounders/m/pounders.m +++ b/pounders/m/pounders.m @@ -66,6 +66,10 @@ Options.printf = 0; % Don't print by default end +if ~isfield(Model, 'model_builder') + Model.model_builder = @formquad; +end + if ~isfield(Model, 'np_max') Model.np_max = 2 * n + 1; end @@ -171,7 +175,7 @@ Res(i, :) = F(i, :) - Cres - .5 * D * reshape(D * reshape(Hres, n, m * n), n, m); end [Mdir, np, valid, Gres, Hresdel, Mind] = ... - formquad(X(1:nf, :), Res(1:nf, :), delta, xk_in, np_max, Model.Par, 0); + Model.model_builder(X(1:nf, :), Res(1:nf, :), delta, xk_in, np_max, Model.Par, 0); if np < n % Must obtain and evaluate bounded geometry points [Mdir, np] = bmpts(X(xk_in, :), Mdir(1:n - np, :), Low, Upp, delta, Model.Par(3)); for i = 1:min(n - np, nf_max - nf) @@ -193,7 +197,7 @@ break end [~, np, valid, Gres, Hresdel, Mind] = ... - formquad(X(1:nf, :), Res(1:nf, :), delta, xk_in, np_max, Model.Par, 0); + Model.model_builder(X(1:nf, :), Res(1:nf, :), delta, xk_in, np_max, Model.Par, 0); if np < n [X, F, hF, flag] = prepare_outputs_before_return(X, F, hF, nf, -5); return @@ -233,7 +237,7 @@ % Check to see if the model is valid within a region of size g_tol delta = max(g_tol, max(abs(X(xk_in, :))) * eps); % Safety for tiny g_tol values [Mdir, ~, valid] = ... - formquad(X(1:nf, :), F(1:nf, :), delta, xk_in, np_max, Model.Par, 1); + Model.model_builder(X(1:nf, :), F(1:nf, :), delta, xk_in, np_max, Model.Par, 1); if ~valid % Make model valid in this small region [Mdir, np] = bmpts(X(xk_in, :), Mdir, Low, Upp, delta, Model.Par(3)); for i = 1:min(n - np, nf_max - nf) @@ -254,7 +258,7 @@ end % Recalculate gradient based on a MFN model [~, ~, valid, Gres, Hres, Mind] = ... - formquad(X(1:nf, :), F(1:nf, :), delta, xk_in, np_max, Model.Par, 0); + Model.model_builder(X(1:nf, :), F(1:nf, :), delta, xk_in, np_max, Model.Par, 0); [G, H] = combinemodels(Cres, Gres, Hres); ind_Lnotbinding = and(X(xk_in, :) > Low, G' > 0); ind_Unotbinding = and(X(xk_in, :) < Upp, G' < 0); @@ -357,7 +361,7 @@ if ~(valid) && (nf < nf_max) && (rho < eta_1) % Implies xk_in & delta unchanged % Need to check because model may be valid after Xsp evaluation [Mdir, np, valid] = ... - formquad(X(1:nf, :), F(1:nf, :), delta, xk_in, np_max, Model.Par, 1); + Model.model_builder(X(1:nf, :), F(1:nf, :), delta, xk_in, np_max, Model.Par, 1); if ~(valid) % ! One strategy for choosing model-improving point: % Update model (exists because delta & xk_in unchanged) for i = 1:nf @@ -365,7 +369,7 @@ Res(i, :) = F(i, :) - Cres - .5 * D * reshape(D * reshape(Hres, n, m * n), n, m); end [~, ~, valid, Gres, Hresdel, Mind] = ... - formquad(X(1:nf, :), Res(1:nf, :), delta, xk_in, np_max, Model.Par, 0); + Model.model_builder(X(1:nf, :), Res(1:nf, :), delta, xk_in, np_max, Model.Par, 0); Hres = Hres + Hresdel; % Update for modelimp; Cres unchanged b/c xk_in unchanged [G, H] = combinemodels(Cres, Gres, Hres); From c6c257f12ba8ea24bd16a5e125f3ab9afb690ad8 Mon Sep 17 00:00:00 2001 From: Matt Menickelly Date: Wed, 24 Apr 2024 16:16:42 -0500 Subject: [PATCH 20/24] Still not broken, comments to come in PR --- pounders/m/pounders.m | 131 ++++++++++++++++++++++++++++-------------- 1 file changed, 88 insertions(+), 43 deletions(-) diff --git a/pounders/m/pounders.m b/pounders/m/pounders.m index 2ab54596..618cc6a5 100644 --- a/pounders/m/pounders.m +++ b/pounders/m/pounders.m @@ -15,6 +15,7 @@ Prior.nfs = 0; Prior.X_init = []; Prior.F_init = []; + Prior.aux_init = {}; Prior.xk_in = 1; end @@ -70,6 +71,10 @@ Model.model_builder = @formquad; end +if ~isfield(Model, 'Ffun_nargout') + Model.Ffun_nargout = 1; +end + if ~isfield(Model, 'np_max') Model.np_max = 2 * n + 1; end @@ -98,7 +103,7 @@ checkinputss(Ffun, X_0, n, Model.np_max, nf_max, g_tol, delta, nfs, m, Prior.F_init, Prior.xk_in, Low, Upp); if flag == -1 % Problem with the input X = []; - F = []; + eval_data = {}; hF = []; return end @@ -131,36 +136,47 @@ if nfs == 0 % Need to do the first evaluation X = [X_0; zeros(nf_max - 1, n)]; % Stores the point locations - F = zeros(nf_max, m); % Stores the function values + eval_data.F = zeros(nf_max, m); % Stores the function values hF = zeros(nf_max, 1); % Stores the sum of squares of evaluated points nf = 1; - F0 = Ffun(X(nf, :)); + switch Model.Ffun_nargout + case 1 + F0 = Ffun(X(nf, :)); + case 2 + [F0, aux0] = Ffun(X(nf, :)); + end if length(F0) ~= m disp(' Error: F0 does not contain the right number of residuals'); flag = -1; return end - F(nf, :) = F0; - if any(isnan(F(nf, :))) - [X, F, hF, flag] = prepare_outputs_before_return(X, F, hF, nf, -3); + eval_data.F(nf, :) = F0; + if Model.Ffun_nargout == 2 + eval_data.aux{nf} = aux0; + end + if any(isnan(eval_data.F(nf, :))) + [X, F, hF, flag] = prepare_outputs_before_return(X, eval_data.F, hF, nf, -3); return end if printf - fprintf('%4i Initial point %11.5e\n', nf, hfun(F(nf, :))); + fprintf('%4i Initial point %11.5e\n', nf, hfun(eval_data.F(nf, :))); end else % Have other function values around - X = [X_0(1:nfs, :); zeros(nf_max, n)]; % Stores the point locations - F = [F0(1:nfs, :); zeros(nf_max, m)]; % Stores the function values + X = [Prior.X_init(1:nfs, :); zeros(nf_max, n)]; % Stores the point locations + eval_data.F = [Prior.F_init(1:nfs, :); zeros(nf_max, m)]; % Stores the function values + if Model.Ffun_nargout == 2 + eval_data.aux{1:nfs} = Prior.aux_init{1:nfs}; + end hF = zeros(nf_max + nfs, 1); % Stores the sum of squares of evaluated points nf = nfs; nf_max = nf_max + nfs; end for i = 1:nf - hF(i) = hfun(F(i, :)); + hF(i) = hfun(eval_data.F(i, :)); end -Res = zeros(size(F)); % Stores the residuals for model updates -Cres = F(xk_in, :); +Res = zeros(size(eval_data.F)); % Stores the residuals for model updates +Cres = eval_data.F(xk_in, :); Hres = zeros(n, n, m); ng = NaN; % Needed for early termination, e.g., if a model is never built % H = zeros(n); G = zeros(n,1); c = hF(xk_in); @@ -172,7 +188,7 @@ % 1a. Compute the interpolation set. for i = 1:nf D = X(i, :) - X(xk_in, :); - Res(i, :) = F(i, :) - Cres - .5 * D * reshape(D * reshape(Hres, n, m * n), n, m); + Res(i, :) = eval_data.F(i, :) - Cres - .5 * D * reshape(D * reshape(Hres, n, m * n), n, m); end [Mdir, np, valid, Gres, Hresdel, Mind] = ... Model.model_builder(X(1:nf, :), Res(1:nf, :), delta, xk_in, np_max, Model.Par, 0); @@ -181,17 +197,24 @@ for i = 1:min(n - np, nf_max - nf) nf = nf + 1; X(nf, :) = min(Upp, max(Low, X(xk_in, :) + Mdir(i, :))); % Temp safeguard - F(nf, :) = Ffun(X(nf, :)); - if any(isnan(F(nf, :))) - [X, F, hF, flag] = prepare_outputs_before_return(X, F, hF, nf, -3); + switch Model.Ffun_nargout + case 1 + eval_data.F(nf, :) = Ffun(X(nf, :)); + case 2 + [F_nf, aux_nf] = Ffun(X(nf, :)); + eval_data.F(nf, :) = F_nf; + eval_data.aux{nf} = aux_nf; + end + if any(isnan(eval_data.F(nf, :))) + [X, F, hF, flag] = prepare_outputs_before_return(X, eval_data.F, hF, nf, -3); return end - hF(nf) = hfun(F(nf, :)); + hF(nf) = hfun(eval_data.F(nf, :)); if printf fprintf('%4i Geometry point %11.5e\n', nf, hF(nf)); end D = Mdir(i, :); - Res(nf, :) = F(nf, :) - Cres - .5 * D * reshape(D * reshape(Hres, n, m * n), n, m); + Res(nf, :) = eval_data.F(nf, :) - Cres - .5 * D * reshape(D * reshape(Hres, n, m * n), n, m); end if nf >= nf_max break @@ -199,13 +222,13 @@ [~, np, valid, Gres, Hresdel, Mind] = ... Model.model_builder(X(1:nf, :), Res(1:nf, :), delta, xk_in, np_max, Model.Par, 0); if np < n - [X, F, hF, flag] = prepare_outputs_before_return(X, F, hF, nf, -5); + [X, F, hF, flag] = prepare_outputs_before_return(X, eval_data.F, hF, nf, -5); return end end % 1b. Update the quadratic model - Cres = F(xk_in, :); + Cres = eval_data.F(xk_in, :); Hres = Hres + Hresdel; c = hF(xk_in); [G, H] = combinemodels(Cres, Gres, Hres); @@ -225,7 +248,7 @@ for i = 1:size(Mind, 1) D = (X(Mind(i), :) - X(xk_in, :)); for j = 1:m - jerr(i, j) = (Cres(j) - F(Mind(i), j)) + D * (Gres(:, j) + .5 * Hres(:, :, j) * D'); + jerr(i, j) = (Cres(j) - eval_data.F(Mind(i), j)) + D * (Gres(:, j) + .5 * Hres(:, :, j) * D'); end end disp(jerr); @@ -237,18 +260,25 @@ % Check to see if the model is valid within a region of size g_tol delta = max(g_tol, max(abs(X(xk_in, :))) * eps); % Safety for tiny g_tol values [Mdir, ~, valid] = ... - Model.model_builder(X(1:nf, :), F(1:nf, :), delta, xk_in, np_max, Model.Par, 1); + Model.model_builder(X(1:nf, :), eval_data.F(1:nf, :), delta, xk_in, np_max, Model.Par, 1); if ~valid % Make model valid in this small region [Mdir, np] = bmpts(X(xk_in, :), Mdir, Low, Upp, delta, Model.Par(3)); for i = 1:min(n - np, nf_max - nf) nf = nf + 1; X(nf, :) = min(Upp, max(Low, X(xk_in, :) + Mdir(i, :))); % Temp safeg. - F(nf, :) = Ffun(X(nf, :)); - if any(isnan(F(nf, :))) - [X, F, hF, flag] = prepare_outputs_before_return(X, F, hF, nf, -3); + switch Model.Ffun_nargout + case 1 + eval_data.F(nf, :) = Ffun(X(nf, :)); + case 2 + [F_nf, aux_nf] = Ffun(X(nf, :)); + eval_data.F(nf, :) = F_nf; + eval_data.aux{nf} = aux_nf; + end + if any(isnan(eval_data.F(nf, :))) + [X, F, hF, flag] = prepare_outputs_before_return(X, eval_data.F, hF, nf, -3); return end - hF(nf) = hfun(F(nf, :)); + hF(nf) = hfun(eval_data.F(nf, :)); if printf fprintf('%4i Critical point %11.5e\n', nf, hF(nf)); end @@ -258,14 +288,14 @@ end % Recalculate gradient based on a MFN model [~, ~, valid, Gres, Hres, Mind] = ... - Model.model_builder(X(1:nf, :), F(1:nf, :), delta, xk_in, np_max, Model.Par, 0); + Model.model_builder(X(1:nf, :), eval_data.F(1:nf, :), delta, xk_in, np_max, Model.Par, 0); [G, H] = combinemodels(Cres, Gres, Hres); ind_Lnotbinding = and(X(xk_in, :) > Low, G' > 0); ind_Unotbinding = and(X(xk_in, :) < Upp, G' < 0); ng = norm(G .* (ind_Lnotbinding + ind_Unotbinding)'); end if ng < g_tol % We trust the small gradient norm and return - [X, F, hF, flag] = prepare_outputs_before_return(X, F, hF, nf, 0); + [X, F, hF, flag] = prepare_outputs_before_return(X, eval_data.F, hF, nf, 0); return end end @@ -278,7 +308,7 @@ elseif spsolver == 2 % Arnold Neumaier's minq5 [Xsp, mdec, minq_err] = minqsw(0, G, H, Lows', Upps', 0, zeros(n, 1)); if minq_err < 0 - [X, F, hF, flag] = prepare_outputs_before_return(X, F, hF, nf, -4); + [X, F, hF, flag] = prepare_outputs_before_return(X, eval_data.F, hF, nf, -4); return end @@ -313,24 +343,31 @@ end if mdec == 0 && valid && all(Xsp == X(xk_in, :)) - [X, F, hF, flag] = prepare_outputs_before_return(X, F, hF, nf, -2); + [X, F, hF, flag] = prepare_outputs_before_return(X, eval_data.F, hF, nf, -2); return end nf = nf + 1; X(nf, :) = Xsp; - F(nf, :) = Ffun(X(nf, :)); - if any(isnan(F(nf, :))) - [X, F, hF, flag] = prepare_outputs_before_return(X, F, hF, nf, -3); + switch Model.Ffun_nargout + case 1 + eval_data.F(nf, :) = Ffun(X(nf, :)); + case 2 + [F_nf, aux_nf] = Ffun(X(nf, :)); + eval_data.F(nf, :) = F_nf; + eval_data.aux{nf} = aux_nf; + end + if any(isnan(eval_data.F(nf, :))) + [X, F, hF, flag] = prepare_outputs_before_return(X, eval_data.F, hF, nf, -3); return end - hF(nf) = hfun(F(nf, :)); + hF(nf) = hfun(eval_data.F(nf, :)); if mdec ~= 0 rho = (hF(nf) - hF(xk_in)) / mdec; else % Note: this conditional only occurs when model is valid if hF(nf) == hF(xk_in) - [X, F, hF, flag] = prepare_outputs_before_return(X, F, hF, nf, -2); + [X, F, hF, flag] = prepare_outputs_before_return(X, eval_data.F, hF, nf, -2); return else rho = inf * sign(hF(nf) - hF(xk_in)); @@ -340,7 +377,7 @@ % 4a. Update the center if (rho >= eta_1) || ((rho > 0) && (valid)) % Update model to reflect new center - Cres = F(xk_in, :); + Cres = eval_data.F(xk_in, :); xk_in = nf; % Change current center end @@ -361,12 +398,12 @@ if ~(valid) && (nf < nf_max) && (rho < eta_1) % Implies xk_in & delta unchanged % Need to check because model may be valid after Xsp evaluation [Mdir, np, valid] = ... - Model.model_builder(X(1:nf, :), F(1:nf, :), delta, xk_in, np_max, Model.Par, 1); + Model.model_builder(X(1:nf, :), eval_data.F(1:nf, :), delta, xk_in, np_max, Model.Par, 1); if ~(valid) % ! One strategy for choosing model-improving point: % Update model (exists because delta & xk_in unchanged) for i = 1:nf D = (X(i, :) - X(xk_in, :)); - Res(i, :) = F(i, :) - Cres - .5 * D * reshape(D * reshape(Hres, n, m * n), n, m); + Res(i, :) = eval_data.F(i, :) - Cres - .5 * D * reshape(D * reshape(Hres, n, m * n), n, m); end [~, ~, valid, Gres, Hresdel, Mind] = ... Model.model_builder(X(1:nf, :), Res(1:nf, :), delta, xk_in, np_max, Model.Par, 0); @@ -397,12 +434,19 @@ nf = nf + 1; X(nf, :) = min(Upp, max(Low, X(xk_in, :) + Xsp)); % Temp safeguard - F(nf, :) = Ffun(X(nf, :)); - if any(isnan(F(nf, :))) - [X, F, hF, flag] = prepare_outputs_before_return(X, F, hF, nf, -3); + switch Model.Ffun_nargout + case 1 + eval_data.F(nf, :) = Ffun(X(nf, :)); + case 2 + [F_nf, aux_nf] = Ffun(X(nf, :)); + eval_data.F(nf, :) = F_nf; + eval_data.aux{nf} = aux_nf; + end + if any(isnan(eval_data.F(nf, :))) + [X, F, hF, flag] = prepare_outputs_before_return(X, eval_data.F, hF, nf, -3); return end - hF(nf) = hfun(F(nf, :)); + hF(nf) = hfun(eval_data.F(nf, :)); if printf fprintf('%4i Model point %11.5e\n', nf, hF(nf)); end @@ -413,7 +457,7 @@ % Update model to reflect new base point D = (X(nf, :) - X(xk_in, :)); xk_in = nf; % Change current center - Cres = F(xk_in, :); + Cres = eval_data.F(xk_in, :); % Don't actually use: for j = 1:m Gres(:, j) = Gres(:, j) + Hres(:, :, j) * D'; @@ -426,3 +470,4 @@ disp('Number of function evals exceeded'); end flag = ng; +F = eval_data.F; From 4f67bf3732abab075842daa16960e7c38e4a6217 Mon Sep 17 00:00:00 2001 From: Jeffrey Larson Date: Mon, 29 Apr 2024 11:56:49 -0500 Subject: [PATCH 21/24] assert --- pounders/m/pounders.m | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/pounders/m/pounders.m b/pounders/m/pounders.m index 618cc6a5..23ac3e6d 100644 --- a/pounders/m/pounders.m +++ b/pounders/m/pounders.m @@ -15,7 +15,7 @@ Prior.nfs = 0; Prior.X_init = []; Prior.F_init = []; - Prior.aux_init = {}; + Prior.aux_init = {}; Prior.xk_in = 1; end @@ -73,6 +73,8 @@ if ~isfield(Model, 'Ffun_nargout') Model.Ffun_nargout = 1; +else + assert(Model.Ffun_nargout == 1 || Model.Ffun_nargout == 2, "Ffun_nargout must be 1 or 2"); end if ~isfield(Model, 'np_max') From cce223a3afb57dcc6fcd33e9c39c6ed0511a2509 Mon Sep 17 00:00:00 2001 From: Matt Menickelly Date: Wed, 22 May 2024 11:08:55 -0500 Subject: [PATCH 22/24] Created model improvement function specific to formquad --- pounders/m/formquad_model_improvement.m | 30 +++++++++++++++++++++++++ pounders/m/pounders.m | 29 +++++------------------- 2 files changed, 36 insertions(+), 23 deletions(-) create mode 100644 pounders/m/formquad_model_improvement.m diff --git a/pounders/m/formquad_model_improvement.m b/pounders/m/formquad_model_improvement.m new file mode 100644 index 00000000..ff4b3158 --- /dev/null +++ b/pounders/m/formquad_model_improvement.m @@ -0,0 +1,30 @@ +function Xsp = formquad_model_improvement(x_k, Cres, Gres, Hres, Mdir, np, Low, Upp, delta, Model, combinemodels) + + n = length(x_k); + + % Update for modelimp; Cres unchanged b/c xk_in unchanged + [G, H] = combinemodels(Cres, Gres, Hres); + + % Evaluate model-improving points to pick best one + % ! May eventually want to normalize Mdir first for infty norm + % Plus directions + [Mdir1, np1] = bmpts(x_k, Mdir(1:n - np, :), Low, Upp, delta, Model.Par(3)); + for i = 1:n - np1 + D = Mdir1(i, :); + Res(i, 1) = D * (G + .5 * H * D'); + end + [a1, b] = min(Res(1:n - np1, 1)); + Xsp = Mdir1(b, :); + % Minus directions + [Mdir1, np2] = bmpts(x_k, -Mdir(1:n - np, :), Low, Upp, delta, Model.Par(3)); + for i = 1:n - np2 + D = Mdir1(i, :); + Res(i, 1) = D * (G + .5 * H * D'); + end + [a2, b] = min(Res(1:n - np2, 1)); + if a2 < a1 + Xsp = Mdir1(b, :); + end + +end + diff --git a/pounders/m/pounders.m b/pounders/m/pounders.m index 23ac3e6d..f6a33f7a 100644 --- a/pounders/m/pounders.m +++ b/pounders/m/pounders.m @@ -71,6 +71,10 @@ Model.model_builder = @formquad; end +if ~isfield(Model, 'model_improver') + Model.model_improver = @formquad_model_improvement; +end + if ~isfield(Model, 'Ffun_nargout') Model.Ffun_nargout = 1; else @@ -410,30 +414,9 @@ [~, ~, valid, Gres, Hresdel, Mind] = ... Model.model_builder(X(1:nf, :), Res(1:nf, :), delta, xk_in, np_max, Model.Par, 0); Hres = Hres + Hresdel; - % Update for modelimp; Cres unchanged b/c xk_in unchanged - [G, H] = combinemodels(Cres, Gres, Hres); - - % Evaluate model-improving points to pick best one - % ! May eventually want to normalize Mdir first for infty norm - % Plus directions - [Mdir1, np1] = bmpts(X(xk_in, :), Mdir(1:n - np, :), Low, Upp, delta, Model.Par(3)); - for i = 1:n - np1 - D = Mdir1(i, :); - Res(i, 1) = D * (G + .5 * H * D'); - end - [a1, b] = min(Res(1:n - np1, 1)); - Xsp = Mdir1(b, :); - % Minus directions - [Mdir1, np2] = bmpts(X(xk_in, :), -Mdir(1:n - np, :), Low, Upp, delta, Model.Par(3)); - for i = 1:n - np2 - D = Mdir1(i, :); - Res(i, 1) = D * (G + .5 * H * D'); - end - [a2, b] = min(Res(1:n - np2, 1)); - if a2 < a1 - Xsp = Mdir1(b, :); - end + Xsp = Model.model_improver(X(xk_in, :), Cres, Gres, Hres, Mdir, np, Low, Upp, delta, Model, combinemodels); + nf = nf + 1; X(nf, :) = min(Upp, max(Low, X(xk_in, :) + Xsp)); % Temp safeguard switch Model.Ffun_nargout From f0f0c800c6bea4473c44c9bb8c574b1ffff29bb3 Mon Sep 17 00:00:00 2001 From: Matt Menickelly Date: Wed, 22 May 2024 11:39:54 -0500 Subject: [PATCH 23/24] Temporary solution to different APIs required between using aux and not --- pounders/m/formquad.m | 2 +- pounders/m/pounders.m | 15 ++++++++------- 2 files changed, 9 insertions(+), 8 deletions(-) diff --git a/pounders/m/formquad.m b/pounders/m/formquad.m index cb6279f9..c27eb21b 100644 --- a/pounders/m/formquad.m +++ b/pounders/m/formquad.m @@ -32,7 +32,7 @@ % H [dbl] [n-by-n-by-m] Array of model Hessians centered at X(xk_in, :) % Mind [int] [np_max-by-1] Integer vector of model interpolation indices % -function [Mdir, np, valid, G, H, Mind] = formquad(X, F, delta, xk_in, np_max, Pars, vf) +function [Mdir, np, valid, G, H, Mind] = formquad(X, F, delta, xk_in, np_max, Pars, vf, ignored_arg) % --DEPENDS ON------------------------------------------------------------- % phi2eval : Evaluates the quadratic basis for vector inputs diff --git a/pounders/m/pounders.m b/pounders/m/pounders.m index f6a33f7a..a279619c 100644 --- a/pounders/m/pounders.m +++ b/pounders/m/pounders.m @@ -143,6 +143,7 @@ if nfs == 0 % Need to do the first evaluation X = [X_0; zeros(nf_max - 1, n)]; % Stores the point locations eval_data.F = zeros(nf_max, m); % Stores the function values + eval_data.aux = []; hF = zeros(nf_max, 1); % Stores the sum of squares of evaluated points nf = 1; switch Model.Ffun_nargout @@ -196,8 +197,8 @@ D = X(i, :) - X(xk_in, :); Res(i, :) = eval_data.F(i, :) - Cres - .5 * D * reshape(D * reshape(Hres, n, m * n), n, m); end - [Mdir, np, valid, Gres, Hresdel, Mind] = ... - Model.model_builder(X(1:nf, :), Res(1:nf, :), delta, xk_in, np_max, Model.Par, 0); + [Mdir, np, valid, Gres, Hresdel, Mind] = ... + Model.model_builder(X(1:nf, :), Res(1:nf, :), delta, xk_in, np_max, Model.Par, 0, eval_data.aux); if np < n % Must obtain and evaluate bounded geometry points [Mdir, np] = bmpts(X(xk_in, :), Mdir(1:n - np, :), Low, Upp, delta, Model.Par(3)); for i = 1:min(n - np, nf_max - nf) @@ -226,7 +227,7 @@ break end [~, np, valid, Gres, Hresdel, Mind] = ... - Model.model_builder(X(1:nf, :), Res(1:nf, :), delta, xk_in, np_max, Model.Par, 0); + Model.model_builder(X(1:nf, :), Res(1:nf, :), delta, xk_in, np_max, Model.Par, 0, eval_data.aux); if np < n [X, F, hF, flag] = prepare_outputs_before_return(X, eval_data.F, hF, nf, -5); return @@ -266,7 +267,7 @@ % Check to see if the model is valid within a region of size g_tol delta = max(g_tol, max(abs(X(xk_in, :))) * eps); % Safety for tiny g_tol values [Mdir, ~, valid] = ... - Model.model_builder(X(1:nf, :), eval_data.F(1:nf, :), delta, xk_in, np_max, Model.Par, 1); + Model.model_builder(X(1:nf, :), eval_data.F(1:nf, :), delta, xk_in, np_max, Model.Par, 1, eval_data.aux); if ~valid % Make model valid in this small region [Mdir, np] = bmpts(X(xk_in, :), Mdir, Low, Upp, delta, Model.Par(3)); for i = 1:min(n - np, nf_max - nf) @@ -294,7 +295,7 @@ end % Recalculate gradient based on a MFN model [~, ~, valid, Gres, Hres, Mind] = ... - Model.model_builder(X(1:nf, :), eval_data.F(1:nf, :), delta, xk_in, np_max, Model.Par, 0); + Model.model_builder(X(1:nf, :), eval_data.F(1:nf, :), delta, xk_in, np_max, Model.Par, 0, eval_data.aux); [G, H] = combinemodels(Cres, Gres, Hres); ind_Lnotbinding = and(X(xk_in, :) > Low, G' > 0); ind_Unotbinding = and(X(xk_in, :) < Upp, G' < 0); @@ -404,7 +405,7 @@ if ~(valid) && (nf < nf_max) && (rho < eta_1) % Implies xk_in & delta unchanged % Need to check because model may be valid after Xsp evaluation [Mdir, np, valid] = ... - Model.model_builder(X(1:nf, :), eval_data.F(1:nf, :), delta, xk_in, np_max, Model.Par, 1); + Model.model_builder(X(1:nf, :), eval_data.F(1:nf, :), delta, xk_in, np_max, Model.Par, 1, eval_data.aux); if ~(valid) % ! One strategy for choosing model-improving point: % Update model (exists because delta & xk_in unchanged) for i = 1:nf @@ -412,7 +413,7 @@ Res(i, :) = eval_data.F(i, :) - Cres - .5 * D * reshape(D * reshape(Hres, n, m * n), n, m); end [~, ~, valid, Gres, Hresdel, Mind] = ... - Model.model_builder(X(1:nf, :), Res(1:nf, :), delta, xk_in, np_max, Model.Par, 0); + Model.model_builder(X(1:nf, :), Res(1:nf, :), delta, xk_in, np_max, Model.Par, 0, eval_data.aux); Hres = Hres + Hresdel; Xsp = Model.model_improver(X(xk_in, :), Cres, Gres, Hres, Mdir, np, Low, Upp, delta, Model, combinemodels); From 4efa73fced4ea8170b920eea440526a63a4b85c0 Mon Sep 17 00:00:00 2001 From: Matt Menickelly Date: Wed, 22 May 2024 12:22:23 -0500 Subject: [PATCH 24/24] An example I need to demo using GN models with true Jacobians --- pounders/m/gaussnewton.m | 14 ++++++++++++++ pounders/m/pounders.m | 15 ++++++++++++--- 2 files changed, 26 insertions(+), 3 deletions(-) create mode 100644 pounders/m/gaussnewton.m diff --git a/pounders/m/gaussnewton.m b/pounders/m/gaussnewton.m new file mode 100644 index 00000000..8c583b86 --- /dev/null +++ b/pounders/m/gaussnewton.m @@ -0,0 +1,14 @@ +function [Mdir, np, valid, Gres, Hres, Mind] = gaussnewton(X, F, delta, xk_in, np_max, Pars, vf, aux) + + Gres = aux{xk_in}'; + [n, m] = size(Gres); + Hres = zeros(n, n, m); + + % silly things to make everything compatible in pounders + valid = true; + np = n; + + Mdir = []; % because valid is always true, pounders should never try to read this. + Mind = []; % ditto above + +end \ No newline at end of file diff --git a/pounders/m/pounders.m b/pounders/m/pounders.m index a279619c..24412aba 100644 --- a/pounders/m/pounders.m +++ b/pounders/m/pounders.m @@ -193,9 +193,18 @@ while nf < nf_max % 1a. Compute the interpolation set. - for i = 1:nf - D = X(i, :) - X(xk_in, :); - Res(i, :) = eval_data.F(i, :) - Cres - .5 * D * reshape(D * reshape(Hres, n, m * n), n, m); + % TEMPORARY HACK: We need a programmatic solution around this, in + % particular, something that indicates if we're doing minimum norm + % change models or just minimum norm models. This is a problem for me + % only because Gauss-Newton models are never going to be "minimum + % change" models. + if Model.Ffun_nargout == 2 + Res(xk_in, :) = eval_data.F(xk_in, :); + else + for i = 1:nf + D = X(i, :) - X(xk_in, :); + Res(i, :) = eval_data.F(i, :) - Cres - .5 * D * reshape(D * reshape(Hres, n, m * n), n, m); + end end [Mdir, np, valid, Gres, Hresdel, Mind] = ... Model.model_builder(X(1:nf, :), Res(1:nf, :), delta, xk_in, np_max, Model.Par, 0, eval_data.aux);