diff --git a/AUTHORS b/AUTHORS index b32f48a0d6a67ec07c9c78a09049e8dd868daba4..6a504a173e3533ebe242acaeafb09bdfa74bc0da 100644 --- a/AUTHORS +++ b/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) diff --git a/ReleaseNotes.txt b/ReleaseNotes.txt index fa0f91406e84a4e09ff8e55501fc2fda35cb3bb2..7636db704b46288b08a9dd6f64e99ae42de8cdac 100644 --- a/ReleaseNotes.txt +++ b/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) diff --git a/common/src/melonpy.cpp b/common/src/melonpy.cpp index cf778778c038df80eb04988e7c3c26bbe9b74e88..02e2d5f0e1b8078ddf9329b14988641be8a4c944 100644 --- a/common/src/melonpy.cpp +++ b/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(); diff --git a/gaussian process/training/python/example_training_gp.py b/gaussian process/training/python/example_training_gp.py index 4d8966f7092c7fa6bdff704d5597b06bb726f2a6..482e71c8a7e529493cc24eff9424221409a9c7c4 100644 --- a/gaussian process/training/python/example_training_gp.py +++ b/gaussian process/training/python/example_training_gp.py @@ -1,7 +1,7 @@ ## -# @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) diff --git a/gaussian process/training/python/example_training_gp_melonpy.py b/gaussian process/training/python/example_training_gp_melonpy.py new file mode 100644 index 0000000000000000000000000000000000000000..b3f4dde447e34ce297922231c6b8c86eef4aa5d3 --- /dev/null +++ b/gaussian process/training/python/example_training_gp_melonpy.py @@ -0,0 +1,201 @@ +## +# @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 diff --git a/gaussian process/training/python/utils.py b/gaussian process/training/python/utils.py index c19d7f6f624eda4529e666761a1b298e246f561e..7831b302ba5085ed0619389b41435ff450a69330 100644 --- a/gaussian process/training/python/utils.py +++ b/gaussian process/training/python/utils.py @@ -1,40 +1,43 @@ -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)