Display method and some the implementations

pysurrogate/experimental/run.py

@@ -0,0 +1,27 @@
+
+import numpy as np
+import matplotlib.pyplot as plt
+from pysurrogate.optimize import fit, predict
+
+
+if __name__ == '__main__':
+
+    # number of samples we will use for this example
+    n_samples = 100
+
+    # ---------------------------------------------------------
+    # Example 1: One input variable and one target
+    # ---------------------------------------------------------
+
+    X = np.random.rand(n_samples, 20) * 4 * np.pi
+    Y = np.cos(X)
+
+    # fit the model and predict the data
+    model = fit(X, Y, methods='sklearn_dacefit', disp=True, debug=True)
+
+    _X = np.linspace(0, 4 * np.pi, 1000)
+    _Y = predict(model, _X)
+
+    plt.scatter(X, Y, label="Observations")
+    plt.plot(_X, _Y, label="True")
+    plt.show()

pysurrogate/impl/my_dacefit/corr.py

@@ -0,0 +1,16 @@
+
+import numpy as np
+
+# function to calculate the correlation matrix all in one
+def calc_kernel_matrix(A, B, func, theta):
+    D = np.repeat(A, B.shape[0], axis=0) - np.tile(B, (A.shape[0], 1))
+    K = func(D, theta)
+    return np.reshape(K, (A.shape[0], B.shape[0]))
+
+
+# -------------------------------
+# Correlation Functions
+# -------------------------------
+
+def corr_gauss(D, theta):
+    return np.exp(np.sum(np.square(D) * -theta, axis=1))

pysurrogate/impl/my_dacefit/fit.py

@@ -0,0 +1,45 @@
+import warnings
+
+import numpy as np
+from numpy.linalg import LinAlgError
+
+from pysurrogate.impl.my_dacefit.corr import corr_gauss, calc_kernel_matrix
+from pysurrogate.impl.my_dacefit.regr import regr_constant
+
+
+def fit(X, Y, theta, regr=regr_constant, kernel=corr_gauss):
+    # attributes used for convenience
+    n_sample, n_var, n_target = X.shape[0], X.shape[1], Y.shape[1]
+
+    # calculate the kernel matrix R
+    R = calc_kernel_matrix(X, X, kernel, theta)
+    R += np.eye(n_sample) * (10 + n_sample) * 2.220446049250313e-16
+
+    # do the cholesky decomposition
+    try:
+        C = np.linalg.cholesky(R)
+    except LinAlgError:
+        warnings.warn("Error while doing Cholesky Decomposition.")
+        return {'obj': np.inf}
+
+    # fit the least squares for regression
+    F = regr(X)
+    Ft = np.linalg.lstsq(C, F, rcond=None)[0]
+    Q, G = np.linalg.qr(Ft)
+    rcond = 1.0 / np.linalg.cond(G)
+    if rcond > 1e15:
+        raise Exception('F is too ill conditioned: Poor combination of regression model and design sites')
+    Yt = np.linalg.solve(C, Y)
+    beta = np.linalg.lstsq(G, Q.T @ Yt, rcond=None)[0]
+
+    # calculate the residual to fit with gaussian process and calculate objective function
+    rho = Yt - Ft @ beta
+    sigma2 = np.sum(np.square(rho), axis=0) / n_sample
+    detR = np.prod(np.power(np.diag(C), (2 / n_sample)))
+    obj = np.sum(sigma2) * detR
+
+    # finally gamma to predict values
+    gamma = np.linalg.solve(C.T, rho)
+
+    return {'R': R, 'C': C, 'F': F, 'Ft': Ft, 'Q': Q, 'G': G, 'Yt': Yt, 'beta': beta, 'rho': rho,
+            '_sigma2': sigma2, 'obj': obj, 'gamma': gamma}

pysurrogate/impl/my_dacefit/regr.py

@@ -0,0 +1,4 @@
+import numpy as np
+
+def regr_constant(X):
+    return np.ones((X.shape[0], 1))

pysurrogate/optimize.py

@@ -1,3 +1,4 @@
+import time
 import warnings
 
 import numpy as np
@@ -10,23 +11,38 @@
 
 def get_method(name):
     if name == 'scipy_rbf':
-        from pysurrogate.models.scipy_rbf import RBF
+        from pysurrogate.surrogates.scipy_rbf import RBF
         return RBF
     elif name == 'matlab_dacefit':
-        from pysurrogate.models.matlab_dacefit import Dacefit
+        from pysurrogate.surrogates.matlab_dacefit import Dacefit
         return Dacefit
     elif name == 'sklearn_polyregr':
-        from pysurrogate.models.sklearn_polyregr import PolynomialRegression
+        from pysurrogate.surrogates.sklearn_polyregr import PolynomialRegression
         return PolynomialRegression
     elif name == 'sklearn_dacefit':
-        from pysurrogate.models.sklearn_dacefit import Dacefit
+        from pysurrogate.surrogates.sklearn_dacefit import Dacefit
         return Dacefit
+    elif name == 'my_dacefit':
+        from pysurrogate.surrogates.my_dacefit import MyDacefit
+        return MyDacefit
     elif name == 'myrbf':
-        from pysurrogate.models.my_rbf import MyRBF
+        from pysurrogate.surrogates.my_rbf import MyRBF
         return MyRBF
-    elif name == 'nn':
-        from pysurrogate.models.torch_nn import Torch
+    elif name == 'george_gp':
+        from pysurrogate.surrogates.gpgeorge_gp import GPGeorge
+        return GPGeorge
+    elif name == 'gpy_gp':
+        from pysurrogate.surrogates.gpy_gp import GPyMetamodel
+        return GPyMetamodel
+    elif name == 'sklearn_gradient_boosting':
+        from pysurrogate.surrogates.sklearn_gradient_boosting import GradientBoosting
+        return GradientBoosting
+    elif name == 'torch_nn':
+        from pysurrogate.surrogates.torch_nn import Torch
         return Torch
+    elif name == 'active_subspace_gp':
+        from pysurrogate.surrogates.active_subspaces import ActiveSubspacesSurrogate
+        return ActiveSubspacesSurrogate
     else:
         raise Exception('Surrogate is unknown: %s' % name)
 
@@ -44,9 +60,13 @@ def get_method_and_params(entry):
             "Either a list of strings (each is a model), or a list of tuples (1) model (2) params.")
 
 
-def fit(X, Y, methods=['nn', 'scipy_rbf', 'sklearn_polyregr', 'sklearn_dacefit'], func_error=calc_mse, disp=False,
+
+
+
+def fit(X, Y, methods=['george_gp', 'gpy_gp', 'sklearn_gradient_boosting', 'my_dacefit', 'torch_nn', 'scipy_rbf', 'sklearn_polyregr',
+                       'sklearn_dacefit'], func_error=calc_mse, disp=False,
         normalize_X=False, normalize_Y=False,
-        do_crossvalidation=True, n_folds=10, crossvalidation_sets=None):
+        do_crossvalidation=True, n_folds=5, crossvalidation_sets=None, debug=False):
     """
 
     The is the public interface which fits a surrogate res that is able to predict more than one target value.
@@ -65,7 +85,8 @@ def fit(X, Y, methods=['nn', 'scipy_rbf', 'sklearn_polyregr', 'sklearn_dacefit']
         prediction F_hat of the res with the true values F.
     disp : bool
         Print output during the fitting of the surrogate archive with information about the error.
-
+    debug : bool
+        If true warnings and exceptions are shown. Otherwise they are suppressed.
     Returns
     -------
 
@@ -107,8 +128,11 @@ def fit(X, Y, methods=['nn', 'scipy_rbf', 'sklearn_polyregr', 'sklearn_dacefit']
 
         try:
             method, params = get_method_and_params(entry)
-        except:
-            warnings.warn("Not able to load model %s. Will be skipped." % entry)
+        except Exception as e:
+            if debug:
+                raise e
+                warnings.warn(str(e))
+                warnings.warn("Not able to load model %s. Will be skipped." % entry)
             continue
 
         for param in params:
@@ -132,29 +156,38 @@ def fit(X, Y, methods=['nn', 'scipy_rbf', 'sklearn_polyregr', 'sklearn_dacefit']
             result = []
 
             # for each method validate
-            for entry in surrogates:
-                name, method, param = entry['name'], entry['method'], entry['param']
+            for k, entry in enumerate(surrogates):
+
+                try:
+                    name, method, param = entry['name'], entry['method'], entry['param']
 
-                error = np.full(n_folds, np.inf)
+                    error = np.full(n_folds, np.inf)
+                    duration = np.full(n_folds, np.nan)
 
-                # on each validation set
-                for i, (training, test) in enumerate(crossvalidation_sets):
-                    impl = method(**param)
-                    impl.fit(X[training, :], Y[training, [m]])
+                    # on each validation set
+                    for i, (training, test) in enumerate(crossvalidation_sets):
+                        impl = method(**param)
 
-                    Y_hat = impl.predict(X[test, :], return_std=False)
-                    error[i] = func_error(Y[test, [m]], Y_hat)
+                        start_time = time.time()
+                        warnings.filterwarnings("ignore")
+                        impl.fit(X[training, :], Y[training, [m]])
+                        duration[i] = time.time() - start_time
 
-                result.append({'name': name, 'method': method, 'param': param, 'error': error})
+                        Y_hat = impl.predict(X[test, :], return_std=False)
+                        error[i] = func_error(Y[test, [m]], Y_hat)
+
+                except Exception as e:
+                    if debug:
+                        print(e)
+                        warnings.warn("Error while using fitting: %s %s %s" % (name, method, param))
+
+                result.append(
+                    {'name': name, 'method': method, 'param': param, 'error': error, 'duration': np.mean(duration)})
 
             result = sorted(result, key=lambda e: np.mean(e['error']))
 
             if disp:
-                print("Target %s" % (m + 1))
-                for i in range(len(result)):
-                    entry = result[i]
-                    print(entry['name'], entry['param'], entry['error'], np.mean(entry['error']))
-                print("=" * 40)
+                __display(result, str(m+1))
 
             crossvalidation.append(result)
 
@@ -177,29 +210,59 @@ def fit(X, Y, methods=['nn', 'scipy_rbf', 'sklearn_polyregr', 'sklearn_dacefit']
         impl.fit(X, Y[:, m])
         models.append(impl)
 
-    res['models'] = models
+    res['surrogates'] = models
 
     return res
 
 
 def predict(res, X):
-
     # if it is only one dimensional convert it
     if X.ndim == 1:
         X = X[:, None]
 
-    Y = np.full((X.shape[0], len(res['models'])), np.inf)
+    Y = np.full((X.shape[0], len(res['surrogates'])), np.inf)
 
     # denormalize if normalized before
     if res['normalize_X']:
         X = normalize(X, res['X_min'], res['X_max'])
 
     # for each target value to predict there exists a model
-    for m, model in enumerate(res['models']):
+    for m, model in enumerate(res['surrogates']):
         Y[:, m] = model.predict(X)
 
     # denormalize target if done while fitting
     if res['normalize_Y']:
         Y = denormalize(Y, res['Y_min'], res['Y_max'])
 
     return Y
+
+
+
+def display(model):
+    for m, result in enumerate(model['crossvalidation']):
+        __display(result,(m+1))
+
+
+def __display(result, lbl_obj):
+
+    for i in range(len(result)):
+
+        entry = result[i]
+
+        attrs = [('name', entry['name'], 27),
+                 ('error', "%.5f" % np.mean(entry['error']), 7),
+                 ('duration (sec)', "%.5f" % np.mean(entry['duration']), 7),
+                 ('param', entry['param'], 200),
+                 ]
+
+        regex = " | ".join(["{}"] * len(attrs))
+
+        if i == 0:
+            print("=" * 50)
+            print("Target %s" % lbl_obj)
+            print("=" * 50)
+            print(regex.format(*[name.ljust(width) for name, _, width in attrs]))
+            print("-" * 50)
+
+        print(regex.format(*[str(val).ljust(width) for _, val, width in attrs]))
+

pysurrogate/surrogates/gpgeorge_gp.py

@@ -0,0 +1,74 @@
+import george
+import numpy as np
+import scipy.optimize as op
+
+from pysurrogate.surrogate import Surrogate
+
+
+class GPGeorge(Surrogate):
+    def __init__(self, kernel):
+        Surrogate.__init__(self)
+        self.kernel = kernel
+        self.model = None
+        self.F = None
+
+    def _predict(self, X):
+        return self.model.predict(self.F, X, return_var=True)
+
+    def _fit(self, X, F):
+
+        self.F = F
+        n_var = X.shape[1]
+
+        if self.kernel == "linear":
+            kernel = george.kernels.LinearKernel(order=2, log_gamma2=0.2, ndim=n_var)
+        elif self.kernel == "expsqrt":
+            kernel = george.kernels.ExpSquaredKernel(metric=np.ones(n_var), ndim=n_var)
+        elif self.kernel == "rational_quad":
+            kernel = george.kernels.RationalQuadraticKernel(log_alpha=0.2, metric=np.ones(n_var), ndim=n_var)
+        elif self.kernel == "exp":
+            kernel = george.kernels.ExpKernel(metric=np.ones(n_var), ndim=n_var)
+        elif self.kernel == "polynomial":
+            kernel = george.kernels.PolynomialKernel(metric=np.ones(n_var))
+        else:
+            raise ValueError("Parameter %s for kernel unknown." % self.kernel)
+
+        gp = george.GP(kernel, fit_mean=True)
+
+        t = 0.1 * np.ones((n_var, 1))
+
+        # Define the objective function (negative log-likelihood in this case).
+        def nll(p):
+            gp.set_parameter_vector(p)
+            ll = gp.log_likelihood(y, quiet=True)
+            return -ll if np.isfinite(ll) else 1e25
+
+        # And the gradient of the objective function.
+        def grad_nll(p):
+            gp.set_parameter_vector(p)
+            return -gp.grad_log_likelihood(y, quiet=True)
+
+        # You need to compute the GP once before starting the optimization.
+        gp.compute(t)
+
+        # Print the initial ln-likelihood.
+        print(gp.log_likelihood(F))
+
+        # Run the optimization routine.
+        p0 = gp.get_parameter_vector()
+        results = op.minimize(nll, p0, jac=grad_nll, method="L-BFGS-B")
+
+        # Update the kernel and print the final log-likelihood.
+        gp.set_parameter_vector(results.x)
+        print(gp.log_likelihood(F))
+
+        gp.optimize(X, F)
+        self.model = gp
+
+
+    @staticmethod
+    def get_params():
+        val = []
+        for kernel in ['linear', 'expsqrt', 'rational_quad', 'exp']:  # , , 'exp', , 'polynomial']:
+            val.append({'kernel': kernel})
+        return val