Merge pull request #362 from yisangriB/surrogate_fix

sy - surrogate fix

modules/performFEM/surrogateGP/gpPredict.py

@@ -10,6 +10,9 @@
 from scipy.stats import lognorm, norm
 from sklearn.linear_model import LinearRegression
 
+errFileName = os.path.join(os.getcwd(), 'workflow.err')  # noqa: N816, PTH109, PTH118
+sys.stderr = open(errFileName, 'a')  # noqa: SIM115, PTH123
+
 try:
     moduleName = 'GPy'  # noqa: N816
     import GPy as GPy  # noqa: PLC0414

modules/performUQ/SimCenterUQ/surrogateBuild.py

@@ -89,6 +89,7 @@
 print('Initializing error log file..')  # noqa: T201
 print(f'Current working dir (getcwd): {os.getcwd()}')  # noqa: T201, PTH109
 
+# errFileName = os.path.join(os.getcwd(),'dakota.err')
 work_dir_tmp = sys.argv[1].replace(os.sep, '/')
 errFileName = os.path.join(work_dir_tmp, 'dakota.err')  # noqa: N816, PTH118
 
@@ -105,7 +106,6 @@
     sys.stderr = open(errFileName, 'w')  # noqa: SIM115, PTH123
     print(f'Error file created at: {errFileName}')  # noqa: T201
 
-
 #
 # Modify GPy package
 #
@@ -254,7 +254,10 @@ def readJson(self):  # noqa: C901, N802, D102, PLR0912, PLR0915
                     self.isEEUQ = True
                 elif Evt[0]['EventClassification'] == 'Wind':
                     self.isWEUQ = True
-                elif Evt[0]['EventClassification'] == 'Hydro' or Evt[0]['EventClassification'] == 'Water':
+                elif (
+                    Evt[0]['EventClassification'] == 'Hydro'
+                    or Evt[0]['EventClassification'] == 'Water'
+                ):
                     self.isHydroUQ = True
 
         self.rv_name_ee = []
@@ -939,10 +942,16 @@ def predictStoVars(self, X_repl, Y_var_repl, X_new, Y_mean, counts):  # noqa: AR
             if parname.endswith('lengthscale'):
                 for nx in range(X_repl.shape[1]):
                     myrange = np.max(X_repl, axis=0) - np.min(X_repl, axis=0)
+                    lb = myrange[nx] / X_repl.shape[0]
+                    ub = myrange[nx] * 5
+                    if lb >= ub:
+                        lb = 0
+
                     m_var.Mat52.lengthscale[[nx]].constrain_bounded(
-                        myrange[nx] / X_repl.shape[0], myrange[nx] * 5, warning=False
+                        lb, ub, warning=False
                     )
                     m_var.Mat52.lengthscale[[nx]] = myrange[nx]  # initial points
+
                     # m_var.Gaussian_noise.value = 0.05
                     # m_var.Gaussian_noise.constrain_bounded(0.1/np.var(log_vars), 0.8/np.var(log_vars), warning=False)
                     # m_var.Mat52.lengthscale[[nx]].constrain_bounded(
@@ -978,7 +987,6 @@ def predictStoVars(self, X_repl, Y_var_repl, X_new, Y_mean, counts):  # noqa: AR
         else:
             norm_var_str = var_pred.T[0]  # if normalization was not used..
 
-
         # norm_var_str = (X_new+2)**2/max((X_new+2)**2)
         Y_metadata = {'variance_structure': norm_var_str / counts}  # noqa: N806
 
@@ -1122,8 +1130,8 @@ def calibrate(self):  # noqa: C901, D102, RUF100
         nugget_opt_tmp = self.nugget_opt
         nopt = self.nopt
 
-        parallel_calib = True
-        # parallel_calib = self.do_parallel
+        # parallel_calib = True
+        parallel_calib = self.do_parallel
 
         if parallel_calib:
             iterables = (
@@ -3436,12 +3444,17 @@ def calibrating(  # noqa: C901, D103
             m_tmp[variance_keyword].constrain_bounded(0.05, 2, warning=False)
             for parname in m_tmp.parameter_names():
                 if parname.endswith('lengthscale'):
-                    for nx in range(X.shape[1]):  # noqa: B007
-                        myrange = np.max(X, axis=0) - np.min(X, axis=0)  # noqa: F841, RUF100
+                    for nx in range(X.shape[1]):
+                        myrange = np.max(X, axis=0) - np.min(X, axis=0)
+                        lb = myrange[nx]
+                        ub = myrange[nx] / X.shape[0] * 10
+                        if lb >= ub:
+                            lb = 0
+
                         exec(  # noqa: S102
                             'm_tmp.'
                             + parname
-                            + '[[nx]].constrain_bounded(myrange[nx] / X.shape[0]*10, myrange[nx],warning=False)'
+                            + '[[nx]].constrain_bounded(lb, ub,warning=False)'
                         )
                         # m_tmp[parname][nx].constrain_bounded(myrange[nx] / X.shape[0], myrange[nx]*100)
         elif nugget_opt_tmp == 'Fixed Values':
@@ -3470,12 +3483,17 @@ def calibrating(  # noqa: C901, D103
 
             for parname in m_tmp.parameter_names():
                 if parname.endswith('lengthscale'):
-                    for nx in range(X.shape[1]):  # noqa: B007
-                        myrange = np.max(X, axis=0) - np.min(X, axis=0)  # noqa: F841
+                    for nx in range(X.shape[1]):
+                        myrange = np.max(X, axis=0) - np.min(X, axis=0)
+                        lb = myrange[nx] / X.shape[0] * 10
+                        ub = myrange[nx]
+                        if lb >= ub:
+                            lb = 0
+
                         exec(  # noqa: S102
                             'm_tmp.'
                             + parname
-                            + '[[nx]].constrain_bounded(myrange[nx] / X.shape[0]*10, myrange[nx],warning=False)'
+                            + '[[nx]].constrain_bounded(lb, ub, warning=False)'
                         )
                         exec(  # noqa: S102
                             'm_tmp.' + parname + '[[nx]] = myrange[nx]*1'
@@ -3527,7 +3545,7 @@ def calibrating(  # noqa: C901, D103
 
             print('Calibrating final surrogate')  # noqa: T201
             m_tmp = my_optimize_restart(m_tmp, nopt)
-  # noqa: RUF100, W293
+        # noqa: RUF100, W293
         # if develop_mode:
         #     print(m_tmp)
         #     #print(m_tmp.rbf.lengthscale)