Release version 0.0.10

AUTHORS

@@ -9,6 +9,7 @@ Contributors:
 - Dominik Bongartz                  (KU Leuven)
   Clara Witte                       (AVT.SVT, RWTH Aachen University)
   Aron Zingler                      (AVT.SVT, RWTH Aachen University)
+  Jannik Lüthje                     (AVT.SVT, RWTH Aachen University)
     Maintenance
 
 - Jaromil Najman                    (AVT.SVT, RWTH Aachen University)

ReleaseNotes.txt

@@ -2,6 +2,11 @@
 Online available: https://git.rwth-aachen.de/avt-svt/public/MeLOn/
 --------------------------------------------------------------------------------
 
+Release version 0.0.10 (July 26th, 2023):
+    - Fixes in Python training script for GPs
+    - Fix in Python bindings
+    - Updates in submodules
+
 Release version 0.0.9 (May 26th, 2023):
     - Minor updates in documentation
     - Updates in submodules (including switch to relative paths)

common/src/melonpy.cpp

@@ -40,7 +40,7 @@ PYBIND11_MODULE(melonpy, m) {
 	// enums
 
 	py::enum_<melon::SCALER_TYPE>(m, "SCALER_TYPE")
-		.value("IDENTIY", melon::SCALER_TYPE::IDENTITY)
+		.value("IDENTITY", melon::SCALER_TYPE::IDENTITY)
 		.value("MINMAX", melon::SCALER_TYPE::MINMAX)
 		.value("STANDARD", melon::SCALER_TYPE::STANDARD)
 		.export_values();

gaussian process/training/python/example_training_gp.py

 ##
-#  @file example_training_of_ANN_with_pruning.py
+#  @file example_training_gp.py
 #
-#  @brief Training of artificial neural network in Keras with pruning and export to file that is readable by MeLOn.
+#  @brief Training of gaussian process in GPyTorch and export to file that is readable by MeLOn.
 #
 # ==============================================================================\n
 #   Aachener Verfahrenstechnik-Systemverfahrenstechnik, RWTH Aachen University  \n
@@ -41,7 +41,7 @@ MATERN = 5
 N_TRAINING_ITER = 250 
 
 # Set output location
-FILE_PATH = "./output""
+FILE_PATH = "./output"
 FILE_NAME = "testGP.json"
 
 
@@ -81,10 +81,15 @@ class ExactGPModel(gpytorch.models.ExactGP): #definition of class
         #super() returns a proxy class which inherits gpytorch.models.ExactGP, its parent
         #in short, we create a subclass by extending the parent
         self.mean_module = gpytorch.means.ConstantMean() # constant mean
-        self.covar_module = gpytorch.kernels.ScaleKernel(
-            gpytorch.kernels.MaternKernel(nu=matern/2, ard_num_dims = 2)
-            )
-        #self.covar_module = gpytorch.kernels.MaternKernel(nu=matern/2, ard_num_dims = 2) #matern class
+        # self.mean_module = gpytorch.means.ZeroMean() # zero mean
+        if matern == 999:
+            self.covar_module = gpytorch.kernels.ScaleKernel(
+                gpytorch.kernels.RBFKernel(ard_num_dims = DX)
+                )
+        else:
+            self.covar_module = gpytorch.kernels.ScaleKernel(
+                gpytorch.kernels.MaternKernel(nu=matern/2, ard_num_dims = DX)
+                )
         #we have flexible options for mean and covar module definition (i.e. with/without white noise, here or below
         
     def forward(self, x): #computes 
@@ -128,5 +133,5 @@ for i in range(N_TRAINING_ITER):
 training_time = training_time - time.time()
 
 ############################  SAVE MODEL IN JSON FILE ############################ 
-utils.save_model_to_json(FILE_PATH, FILE_NAME, model, X_train, y_train, MATERN, scaler)
+utils.save_model_to_json(FILE_PATH, FILE_NAME, model, likelihood, X_train, y_train, MATERN, scaler)
 

gaussian process/training/python/example_training_gp_melonpy.py

+##
+#  @file example_training_of_ANN_with_pruning.py
+#
+#  @brief Training of artificial neural network in Keras with pruning and export to file that is readable by MeLOn.
+#
+# ==============================================================================\n
+#   Aachener Verfahrenstechnik-Systemverfahrenstechnik, RWTH Aachen University  \n
+# ==============================================================================\n
+#
+#  @author Artur M. Schweidtmann, Linus Netze, and Alexander Mitsos
+#  @date 18. February 2021
+##
+
+import utils
+import numpy as np
+import torch
+import gpytorch
+from pyDOE import lhs
+from sklearn.preprocessing import MinMaxScaler, StandardScaler
+import time
+
+from maingopy import *
+# Import MeLOn
+from maingopy.melonpy import *
+# If you are not getting MeLOn through MAiNGO, use:
+# from melonpy import *
+
+############################ SET PARAMETERS ############################ 
+
+# Set parameters for generation of training data 
+DX = 2 #input dims of data/GP
+DY = 1
+PROBLEM_LB = np.array([-3,-3]) #lower bound of inputs
+PROBLEM_UB = np.array([3,3]) #upper bounds of inputs
+N_TRAINING_POINTS = 20
+test_func = lambda x: (
+    3*(1-x[:,0])**2*np.exp(-(x[:,0]**2) - (x[:,1]+1)**2)
+    - 10*(x[:,0]/5 - x[:,0]**3 - x[:,1]**5)*np.exp(-x[:,0]**2-x[:,1]**2)
+    - 1/3*np.exp(-(x[:,0]+1)**2 - x[:,1]**2)
+    )
+
+# Set prameter for model
+MATERN = 5
+
+# Set parameter for training
+N_TRAINING_ITER = 250 
+#gpytorch.settings.fast_computations(False, False, False)
+#gpytorch.settings.fast_pred_samples(False)
+#gpytorch.settings.lazily_evaluate_kernels(False)
+#gpytorch.settings.cholesky_jitter(1e-4)
+
+############################ GENERATE DATA ############################ 
+
+# Generate sampling points using latin hypercube sampling
+X = lhs(DX, samples=N_TRAINING_POINTS)
+
+# Scale to bounds (lhs only yields values from 0 to 1)
+X = PROBLEM_LB + (PROBLEM_UB - PROBLEM_LB)*X
+
+y = test_func(X)
+
+############################  SCALE DATA ############################ 
+
+scaler = dict()
+
+# scale inputs to [0 , 1]
+scaler['input'] = MinMaxScaler(feature_range=(0,1 )) 
+X_scaled = scaler['input'].fit_transform(X)
+X_train = torch.from_numpy(X_scaled)
+
+# scale outputs to zero mean and unit variance
+scaler['output'] = StandardScaler()
+y_scaled = scaler['output'].fit_transform(y.reshape(-1, 1)).squeeze()
+y_train = torch.from_numpy(y_scaled)
+
+
+############################  Define MODEL ############################ 
+
+# A basic exact GP regression model (xompare gpytorch documentation: https://docs.gpytorch.ai/en/stable/examples/01_Exact_GPs/Simple_GP_Regression.html)
+# Example from documentaion is extended by using a Matern instead of a radial basis function kernel.
+
+class ExactGPModel(gpytorch.models.ExactGP): #definition of class
+    def __init__(self, train_x, train_y, likelihood, matern): # def creates objects needed for forward
+        super(ExactGPModel, self).__init__(train_x, train_y, likelihood)
+        #super() returns a proxy class which inherits gpytorch.models.ExactGP, its parent
+        #in short, we create a subclass by extending the parent
+        self.mean_module = gpytorch.means.ConstantMean() # constant mean
+        # self.mean_module = gpytorch.means.ZeroMean() # zero mean
+        if matern == 999:
+            self.covar_module = gpytorch.kernels.ScaleKernel(
+                gpytorch.kernels.RBFKernel(ard_num_dims = DX)
+                )
+        else:
+            self.covar_module = gpytorch.kernels.ScaleKernel(
+                gpytorch.kernels.MaternKernel(nu=matern/2, ard_num_dims = DX)
+                )
+        #we have flexible options for mean and covar module definition (i.e. with/without white noise, here or below
+        
+    def forward(self, x): #computes 
+        mean_x = self.mean_module(x) #prior mean
+        covar_x = self.covar_module(x)# + self.white_noise_module(x) #prior covariance matrix from kernels
+        return gpytorch.distributions.MultivariateNormal(mean_x, covar_x) #multivar normal
+
+# initialize likelihood and model
+likelihood = gpytorch.likelihoods.GaussianLikelihood()
+model = ExactGPModel(X_train, y_train, likelihood, matern=MATERN)
+
+############################  TRAINING ############################ 
+
+training_time = time.time()
+
+# Find optimal model hyperparameters
+model.train()
+likelihood.train()
+
+# Use the adam optimizer
+optimizer = torch.optim.Adam(model.parameters(), lr=0.1)  # Includes GaussianLikelihood parameters
+
+# "Loss" for GPs - the marginal log likelihood
+mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model)
+
+for i in range(N_TRAINING_ITER):
+    # Zero gradients from previous iteration
+    optimizer.zero_grad()
+    # Output from model
+    output = model(*model.train_inputs)
+    # Calc loss and backprop gradients
+    loss = -mll(output, model.train_targets)
+    loss.backward()
+    print('Iter %d/%d \n\tLoss:\t%s\n\tlengthscale:\t%s\n\tnoise:\t%s' % (
+        i + 1, N_TRAINING_ITER, loss,
+        model.covar_module.base_kernel.lengthscale,
+        model.likelihood.noise
+    ))
+    optimizer.step()
+
+training_time = training_time - time.time()
+
+############################  GET MODEL OBJECT ############################ 
+gp_model = utils.generate_melon_gp_object(model, likelihood, X_train, y_train, MATERN, scaler)
+utils.save_model_to_json("", "testGP.json", model, likelihood, X_train, y_train, MATERN, scaler)
+
+############################  DEFINE OPTIMIZATION PROBLEM ############################ 
+
+class Model(MAiNGOmodel):
+
+    #####################################
+    # This function defines the optimization variables in the problem. In the above (MI)NLP, this corresponds to the x \in X part.
+    # A variable typically needs bounds to define the host set X and can optionally have a variable type (VT_CONTINUOUS, VT_INTEGER, VT_BINARY)
+    # and a name (and a branching priority, which is an algorithmic parameter and usually not relevant).
+    # Note that the name is used *only* in the screen or log output of MAiNGO. In particular, it can not be used for modelling...
+    def get_variables(self):
+        return [OptimizationVariable(Bounds(-3, 3), VT_CONTINUOUS, "x"),
+                OptimizationVariable(Bounds(-3, 3), VT_CONTINUOUS, "y")  ]
+
+    #####################################
+    # This function essentially needs to implement the functions f(x), h(x), and g(x) in the aove (MI)NLP.
+    # Unfortunately, right now we cannot use the variable objects defined in the get_variables function above
+    # directly for modeling (nor other variable types defined elsewhere) for technical reasons.
+    # Instead, the "vars" list that MAiNGO hands to this evaluate function contains the same number of elements as the
+    # list that we returned in the get_variables function, and it is only through their position in this list that we can
+    # map the variables used herein to the ones we defined in get_variables.
+    def evaluate(self,vars):
+        
+        # First read in GP parameters from file "fileName" at "filePath"
+        gp = GaussianProcess()
+        gp.load_model(gp_model)
+        
+        # Evaluate the Gaussian process
+        # Input of the GP are the optimiation variables "vars", as defined in the "get_variables" function (cf. discussion above)
+        mu = gp.calculate_prediction_reduced_space(vars)
+        variance = gp.calculate_variance_reduced_space(vars)
+        sigma = sqrt(variance)
+
+        result = EvaluationContainer()
+        result.objective = mu
+        
+        # Optionally, we can define OutputVariables. These are things that we would like to have
+        # evaluated at the optimal solution but that do not form part of the (MI)NLP itself.
+        result.output = [ OutputVariable("mu ", mu),
+                          OutputVariable("sigma: ", sigma) ]
+        return result
+
+
+############################  SOLVE OPTIMIZATION PROBLEM ############################ 
+
+# To actually solve the problem, we instantiate a model, hand it to a MAiNGO instance.
+myModel = Model()
+myMAiNGO = MAiNGO(myModel)
+
+# There are lots of settings regarding algorithmic details, subsolvers, and output.
+# E.g., to disable screen output:
+# myMAiNGO.set_option("outstreamVerbosity", 0)
+
+# Actually solve the problem
+myMAiNGO.solve()
+
+# Query results
+optimalObjective = myMAiNGO.get_objective_value()
+print("\nOptimal objective value: {}".format(optimalObjective))
+optimalSolutionPoint = myMAiNGO.get_solution_point()
\ No newline at end of file

gaussian process/training/python/utils.py

-import matplotlib.pyplot as plt
 import numpy as np
 from sklearn.preprocessing import MinMaxScaler, StandardScaler
-import os
+from pathlib import Path
 import json
-import pymelon
 import torch
+from gpytorch.means import ConstantMean, ZeroMean
+# Import MeLOn
+import maingopy.melonpy as melonpy
+# If you are not getting MeLOn through MAiNGO, use:
+# import melonpy 
 
 def generate_melon_scaler_object(scaler):
-        scaler_data = pymelon.ScalerData()
+        scaler_data = melonpy.ScalerData()
         
         if scaler is None:
-            scaler_data = pymelon.SCALER_TYPE.IDENTITY
+            scaler_data = melonpy.SCALER_TYPE.IDENTITY
             scaler_data.parameters = {}
         elif isinstance(scaler, MinMaxScaler):
-            scaler_data.type = pymelon.SCALER_TYPE.MINMAX
+            scaler_data.type = melonpy.SCALER_TYPE.MINMAX
             
             scaled_bounds = scaler.get_params()['feature_range']
             scaler_data.parameters = {
-                pymelon.SCALER_PARAMETER.UPPER_BOUNDS : scaler.data_max_.tolist(),
-                pymelon.SCALER_PARAMETER.LOWER_BOUNDS : scaler.data_min_.tolist(),
-                pymelon.SCALER_PARAMETER.SCALED_LOWER_BOUNDS : [scaled_bounds[0]]*len(scaler.data_max_.tolist()), 
-                pymelon.SCALER_PARAMETER.SCALED_UPPER_BOUNDS : [scaled_bounds[1]]*len(scaler.data_max_.tolist())
+                melonpy.SCALER_PARAMETER.UPPER_BOUNDS : scaler.data_max_.tolist(),
+                melonpy.SCALER_PARAMETER.LOWER_BOUNDS : scaler.data_min_.tolist(),
+                melonpy.SCALER_PARAMETER.SCALED_LOWER_BOUNDS : [scaled_bounds[0]]*len(scaler.data_max_.tolist()), 
+                melonpy.SCALER_PARAMETER.SCALED_UPPER_BOUNDS : [scaled_bounds[1]]*len(scaler.data_max_.tolist())
             }
         elif isinstance(scaler, StandardScaler):
-            scaler_data.type = pymelon.SCALER_TYPE.STANDARD          
+            scaler_data.type = melonpy.SCALER_TYPE.STANDARD          
             scaler_data.parameters = {
-                pymelon.SCALER_PARAMETER.STD_DEV : np.sqrt(scaler.var_).tolist(),
-                pymelon.SCALER_PARAMETER.MEAN : scaler.mean_.tolist()
+                melonpy.SCALER_PARAMETER.STD_DEV : np.sqrt(scaler.var_).tolist(),
+                melonpy.SCALER_PARAMETER.MEAN : scaler.mean_.tolist()
             }
         else:
             raise Exception("Unsupported scaler type. Scaler has to be either a scikit-learn MinMaxScaler or StandardScaler instance or None (=identity(no scaling))")
         
         return scaler_data 
 
-def generate_melon_gp_object(GP_model, X, y, matern, scaler):
-    gp_data = pymelon.GPData()
+def generate_melon_gp_object(GP_model, GP_likelihood, X, y, matern, scaler):
+    gp_data = melonpy.GPData()
     
     gp_data.X = X.numpy()
     gp_data.Y = y.numpy()
@@ -42,17 +45,17 @@ def generate_melon_gp_object(GP_model, X, y, matern, scaler):
     gp_data.nX, gp_data.DX = X.shape
 
     if len(y.shape) == 1:
-        gp_data.DY = y.shape[0]
+        gp_data.DY = 1
     else:
         gp_data.DY = y.shape[1]
 
+    noise = GP_likelihood.noise.detach().numpy()
     cov_mat = GP_model.covar_module(X)
-    gp_data.K = cov_mat.numpy() 
-    inv_cov_mat = torch.inverse(cov_mat.evaluate())
-    gp_data.invK = inv_cov_mat.detach().numpy()
+    gp_data.K = cov_mat.numpy() + noise * np.eye(N=gp_data.nX)
+    gp_data.invK = np.linalg.inv(gp_data.K)
 
     gp_data.matern = matern
-    kernel_data = pymelon.KernelData()
+    kernel_data = melonpy.KernelData()
     kernel_data.sf2 = GP_model.covar_module.outputscale.detach().numpy().astype(float).item()               #outputscale sigma*K
     kernel_data.ell = GP_model.covar_module.base_kernel.lengthscale.detach().numpy().squeeze().tolist()     #lenghtscales kernel
     gp_data.kernelData = kernel_data
@@ -68,29 +71,43 @@ def generate_melon_gp_object(GP_model, X, y, matern, scaler):
     gp_data.predictionScalerData = generate_melon_scaler_object(scaler['output'])
 
     gp_data.stdOfOutput = np.sqrt(scaler['output'].var_)[0]
-    gp_data.meanFunction = 0
+    if isinstance(GP_model.mean_module, ConstantMean):
+        gp_data.meanFunction = GP_model.mean_module.constant.detach().numpy().tolist()
+    elif isinstance(GP_model.mean_module, ZeroMean):
+        gp_data.meanFunction = 0
+    else: 
+        raise Exception(f'GP uses {type(GP_model.mean_module)} as a mean function. Currently only ConstantMean or ZeroMean are supported as mean modules.') 
 
     return gp_data
 
 
-def save_model_to_json(filepath, filename, GP_model, X, y, matern, scalers=dict()):
+def save_model_to_json(filepath, filename, GP_model, GP_likelihood, X, y, matern, scalers=dict()):
     prediction_parameters = dict()
     prediction_parameters["nX"] = X.shape[0]
     prediction_parameters["DX"] = X.shape[1]
 
     if len(y.shape) == 1:
-        prediction_parameters["DY"] = y.shape[0]
+        prediction_parameters["DY"] = 1
     else:
         prediction_parameters["DY"] = y.shape[1]
 
     prediction_parameters["matern"] = matern
-    prediction_parameters["meanfunction"] = 0
+    
+    if isinstance(GP_model.mean_module, ConstantMean):
+        prediction_parameters["meanfunction"] = GP_model.mean_module.constant.detach().numpy().tolist()
+    elif isinstance(GP_model.mean_module, ZeroMean):
+        prediction_parameters["meanfunction"] = 0
+    else: 
+        raise Exception(f'GP uses {type(GP_model.mean_module)} as a mean function. Currently only ConstantMean or ZeroMean are supported as mean modules.') 
+
     prediction_parameters["X"] = X.numpy().tolist()
     prediction_parameters["Y"] = y.numpy().tolist()
 
-    cov_mat = GP_model.covar_module(X)
-    prediction_parameters["K"] = cov_mat.numpy().tolist()
-    prediction_parameters["invK"] = np.linalg.inv(prediction_parameters["K"]).tolist()
+    noise = GP_likelihood.noise.detach().numpy()
+
+    cov_mat = GP_model.covar_module(X).numpy()
+    prediction_parameters["K"] = (cov_mat + noise * np.eye(N=prediction_parameters["nX"])).tolist()
+    prediction_parameters["invK"] = np.linalg.inv(cov_mat + noise * np.eye(N=prediction_parameters["nX"])).tolist()
 
     if not 'input' in scalers or not isinstance(scalers['input'], MinMaxScaler):
         raise Exception("There has to be an inputscaler which is a scikit-learn MinMaxScaler instance")
@@ -105,10 +122,9 @@ def save_model_to_json(filepath, filename, GP_model, X, y, matern, scalers=dict(
     prediction_parameters["meanOfOutput"] = scalers['output'].mean_.item()
 
     prediction_parameters["sf2"] = GP_model.covar_module.outputscale.detach().numpy().astype(float).item()           #outputscale sigma*K
-    prediction_parameters["ell"] = GP_model.covar_module.base_kernel.lengthscale.detach().numpy().squeeze().tolist()     #lenghtscales kernel
+    prediction_parameters["ell"] = GP_model.covar_module.base_kernel.lengthscale.detach().numpy().flatten().tolist()     #lenghtscales kernel
 
-    if not os.path.exists(filepath):
-        os.makedirs(filepath)
+    Path(filepath).mkdir(parents=True, exist_ok=True)
     
-    with open(os.path.join(filepath,filename), 'w') as outfile:
+    with Path(filepath).joinpath(filename).open('w') as outfile:
         json.dump(prediction_parameters, outfile)