get penalty actually working

examples/plot_penalty.py

@@ -0,0 +1,93 @@
+import numpy as np
+import sbopt
+import matplotlib.pyplot as plt
+
+
+def my_fun(x):
+    # define the rosenbrock function to minmize
+    A = 100.0*((x[1] - (x[0]**2))**2)
+    B = (1.0 - x[0])**2
+    return A + B
+
+
+bounds = np.zeros((2, 2))
+bounds[:, 0] = 0.5
+bounds[:, 1] = 1.5
+
+# set random seed for reproducibility
+np.random.seed(1234124)
+
+# initialize the RbfOpt object
+my_opt = sbopt.RbfOpt(my_fun,  # your objective function to minimize
+                      bounds,  # bounds for your design variables
+                      initial_design='latin',  # initial design type
+                      # 'latin' default, or 'random'
+                      initial_design_ndata=20,  # number of initial points
+                      n_local_optimze=5,  # number of local BFGS optimizers
+                      # scipy radial basis function parameters see
+                      # https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.Rbf.html
+                      rbf_function='linear',  # default for SbOpt is 'linear'
+                      # while the default for scipy.interpolate.Rbf is
+                      # 'multi-quadratic'
+                      epsilon=None,  # (default)
+                      smooth=0.0,  # (default)
+                      norm='euclidean'  # (default)
+                      )
+# run the optimizer
+result = my_opt.minimize(max_iter=1,  # maximum number of iterations
+                         # (default)
+                         n_same_best=20,  # number of iterations to run
+                         # without improving best function value (default)
+                         eps=1e-1,  # minimum distance a new design point
+                         # may be from an existing design point (default)
+                         verbose=1,  # number of iterations to go for
+                         # printing the status (default)
+                         initialize=True,  # boolean, wether or not to
+                         # perform the initial sampling (default)
+                         strategy='local_best',  # str, which minimize strategy
+                         # to use. strategy='local_best' (default) adds only
+                         # one design point per iteration, where this selection
+                         # first looks at the best local optimizer result. 
+                         # strategy='all_local' adds the results from each of
+                         # the local optima in a single iteration.
+                         )
+print('Best design variables:', result[0])
+print('Best function value:', result[1])
+print('Convergence by max iteration:', result[2])
+print('Convergence by n_same_best:', result[3])
+
+
+nval = 1000
+x = np.linspace(0.9, 1.1, nval)
+xx, yy = np.meshgrid(x, x)
+
+X = np.vstack((xx.flatten(), yy.flatten())).T
+Z = np.zeros(nval*nval)
+Z_rbf = np.zeros_like(Z)
+Z_pen = np.zeros_like(Z)
+for i, j in enumerate(X):
+    Z[i] = my_fun(j)
+    Z_rbf[i] = my_opt.rbf_eval(j)
+    Z_pen[i] = my_opt.my_obj(j)
+Z = Z.reshape((nval, nval))
+Z_rbf = Z_rbf.reshape((nval, nval))
+Z_pen = Z_pen.reshape((nval, nval))
+
+plt.figure()
+plt.title('True function')
+plt.contourf(xx, yy, Z)
+plt.colorbar()
+
+plt.figure()
+plt.title('RBF with ' + str(my_opt.X.shape[0]) + ' number of samples')
+plt.contourf(xx, yy, Z_rbf)
+plt.plot(my_opt.X[:, 0], my_opt.X[:, 1], 'xk')
+
+plt.colorbar()
+
+plt.figure()
+plt.title('RBF and Pen with ' + str(my_opt.X.shape[0]) + ' number of samples')
+plt.contourf(xx, yy, Z_pen)
+plt.plot(my_opt.X[:, 0], my_opt.X[:, 1], 'xk')
+plt.colorbar()
+plt.show()

sbopt/sbopt.py

@@ -95,7 +95,7 @@ def minimize(self, max_iter=100, n_same_best=20, eps=1e-6, verbose=1,
             # set one of the local optimizers to start from the best location
             y_best_ind = np.nanargmin(self.X[:, -1])
             x_samp[0] = self.X[y_best_ind, :self.n_dim]
-
+            self.y_mean = np.nanmean(self.X[:, -1])
             for j in range(self.n_local_optimze):
                 # print(x_samp[j], x_samp[j].shape)
                 res = fmin_l_bfgs_b(self.my_obj, x_samp[j],
@@ -190,7 +190,9 @@ def my_obj(self, x):
         dist -= self.eps
         d = np.where(dist >= 0., 0., dist)
         pen = np.dot(d.T, d)[0, 0]
-        return f + pen
+        p = np.sum(np.abs(d))
+        # print(f, self.y_mean*p, self.y_mean*pen)
+        return f + self.y_mean*p + self.y_mean*pen
 
     def transfrom_bounds(self, x):
         """