diff --git a/docs/code/RunModel/ClusterScript_Example/add_numbers.py b/docs/code/RunModel/ClusterScript_Example/add_numbers.py new file mode 100644 index 000000000..89e104291 --- /dev/null +++ b/docs/code/RunModel/ClusterScript_Example/add_numbers.py @@ -0,0 +1,25 @@ +import sys +import os +import json +import numpy as np + +def addNumbers(): + inputPath = sys.argv[1] + outputPath = sys.argv[2] + + # Open JSON file + with open(inputPath, "r") as jsonFile: + data = json.load(jsonFile) + + # Read generated numbers + number1 = data["number1"] + number2 = data["number2"] + + randomAddition = number1 + number2 + + # Write addition to file + with open(outputPath, 'w') as outputFile: + outputFile.write('{}\n'.format(randomAddition)) + +if __name__ == '__main__': + addNumbers() diff --git a/docs/code/RunModel/ClusterScript_Example/addition_run.py b/docs/code/RunModel/ClusterScript_Example/addition_run.py new file mode 100644 index 000000000..3adfa62a0 --- /dev/null +++ b/docs/code/RunModel/ClusterScript_Example/addition_run.py @@ -0,0 +1,18 @@ +import os +import shutil +import fire + +def runAddition(index): + index = int(index) + + inputRealizationPath = os.path.join(os.getcwd(), 'run_' + str(index), 'InputFiles', 'inputRealization_' \ + + str(index) + ".json") + outputPath = os.path.join(os.getcwd(), 'OutputFiles') + + # This is where pre-processing commands would be executed prior to running the cluster script. + command1 = ("echo \"This is where pre-processing would be happening\"") + + os.system(command1) + +if __name__ == '__main__': + fire.Fire(runAddition) diff --git a/docs/code/RunModel/ClusterScript_Example/inputRealization.json b/docs/code/RunModel/ClusterScript_Example/inputRealization.json new file mode 100644 index 000000000..8de6601f6 --- /dev/null +++ b/docs/code/RunModel/ClusterScript_Example/inputRealization.json @@ -0,0 +1,4 @@ +{ + "number1" : , + "number2" : +} diff --git a/docs/code/RunModel/ClusterScript_Example/process_addition_output.py b/docs/code/RunModel/ClusterScript_Example/process_addition_output.py new file mode 100644 index 000000000..393c66fe3 --- /dev/null +++ b/docs/code/RunModel/ClusterScript_Example/process_addition_output.py @@ -0,0 +1,25 @@ +import numpy as np +from pathlib import Path + +class OutputProcessor: + + def __init__(self, index): + filePath = Path("./OutputFiles/qoiFile_" + str(index) + ".txt") + self.numberOfColumns = 0 + self.numberOfLines = 0 + addedNumbers = [] + + # Check if file exists + if filePath.is_file(): + # Now, open and read data + with open(filePath) as f: + for line in f: + currentLine = line.split() + + if len(currentLine) != 0: + addedNumbers.append(currentLine[:]) + + if len(addedNumbers) != 0: + self.qoi = np.vstack(addedNumbers) + else: + self.qoi = np.empty(shape=(0,0)) diff --git a/docs/code/RunModel/ClusterScript_Example/run_script.sh b/docs/code/RunModel/ClusterScript_Example/run_script.sh new file mode 100755 index 000000000..c46b62ae3 --- /dev/null +++ b/docs/code/RunModel/ClusterScript_Example/run_script.sh @@ -0,0 +1,57 @@ +#!/bin/bash + +# NOTE: The job configuration etc. would be in the batch script that launches +# your python script that uses UQpy. This script would then utilize those +# resources by using the appropriate commands here to launch parallel jobs. For +# example, TACC uses slurm and ibrun, so you would launch your python script in +# the slurm batch script and then use ibrun here to tile parallel runs. + +# This function is where you can define all the parts of a single +taskFunction(){ + coresPerProc=$1 + runNumber=$2 + host=$3 + + let offset=$coresPerProc*$runNumber # Sometimes, this might be necessary to pass as an argument to launch jobs. Not used here. + + cd run_$runNumber + # Here, we launch a parallel job. The example uses multiple cores to add numbers, + # which is somewhat pointless. This is just to illustrate the process for how tiled + # parallel jobs are launched and where MPI-capable applications would be initiated + mkdir -p ./OutputFiles + mpirun -n $coresPerProc --host $host:$coresPerProc python3 ../add_numbers.py ./InputFiles/inputRealization_$runNumber.json ./OutputFiles/qoiFile_$runNumber.txt + cd .. +} + +# Get list of hosts +echo $SLURM_NODELIST > hostfile + +# Split by comma +IFS="," read -ra HOSTS < hostfile + +# This is the loop that launches taskFunction in parallel +coresPerProcess=$1 +numberOfJobs=$2 +# This number will vary depending on the number of cores per node. In this case, it is 32. +N=32 + +echo +echo "Starting parallel job launch" + +declare -i index=0 + +for i in $(seq 0 $((numberOfJobs-1))) +do + # Launch task function and put into the background + echo "Launching job number ${i} on ${HOSTS[$index]}" + taskFunction $coresPerProcess $i ${HOSTS[$index]}& + + # Increment host when all nodes allocated on current node + if !((${i}%N)) && [ $i -ne 0 ] + then + index=${index}+1 + fi +done + +wait # This wait call is necessary so that loop above completes before script returns +echo "Analyses done!" diff --git a/docs/code/RunModel/cluster_script_example.py b/docs/code/RunModel/cluster_script_example.py new file mode 100644 index 000000000..8abecd6f1 --- /dev/null +++ b/docs/code/RunModel/cluster_script_example.py @@ -0,0 +1,78 @@ +""" + +Cluster Script Example for Third-party +====================================== +""" + +# %% md +# +# In this case, we're just running a simple addition of random numbers, but +# the process is exactly the same for more complicated workflows. The pre- +# and post-processing is done through `model_script` and `output_script` +# respectively, while the computationally intensive portion of the workflow +# is launched in `cluster_script. The example below provides a minimal framework +# from which more complex cases can be constructed. +# +# Import the necessary libraries + +# %% +from UQpy.sampling import LatinHypercubeSampling +from UQpy.run_model.RunModel import RunModel +from UQpy.run_model.model_execution.ThirdPartyModel import ThirdPartyModel +from UQpy.distributions import Uniform +import numpy as np +import time +import csv + +# %% md +# +# Define the distribution objects. + +# %% + +var_names=["var_1", "var_2"] +distributions = [Uniform(250.0, 40.0), Uniform(66.0, 24.0)] + +# %% md +# +# Draw the samples using Latin Hypercube Sampling. + +# %% + +x_lhs = LatinHypercubeSampling(distributions, nsamples=64) + +# %% md +# +# Run the model. + +# %% + +model = ThirdPartyModel(var_names=var_names, input_template='inputRealization.json', model_script='addition_run.py', + output_script='process_addition_output.py', output_object_name='OutputProcessor', + model_dir='AdditionRuns') + +t = time.time() +modelRunner = RunModel(model=model, samples=x_lhs.samples, ntasks=1, + cores_per_task=2, nodes=1, resume=False, + run_type='CLUSTER', cluster_script='./run_script.sh') + +t_total = time.time() - t +print("\nTotal time for all experiments:") +print(t_total, "\n") + +# %% md +# +# Print model results--this is just for illustration + +# %% +for index, experiment in enumerate(modelRunner.qoi_list, 0): + if len(experiment.qoi) != 0: + for item in experiment.qoi: + print("These are the random numbers for sample {}:".format(index)) + for sample in x_lhs.samples[index]: + print("{}\t".format(sample)) + + print("This is their sum:") + for result in item: + print("{}\t".format(result)) + print() diff --git a/docs/code/sampling/monte_carlo/plot_monte_carlo.py b/docs/code/sampling/monte_carlo/monte_carlo.py similarity index 100% rename from docs/code/sampling/monte_carlo/plot_monte_carlo.py rename to docs/code/sampling/monte_carlo/monte_carlo.py diff --git a/docs/source/runmodel_doc.rst b/docs/source/runmodel_doc.rst index dd094dcb1..9c8ffaf72 100644 --- a/docs/source/runmodel_doc.rst +++ b/docs/source/runmodel_doc.rst @@ -106,6 +106,21 @@ at `https://www.open-mpi.org/faq/?category=building self.n_iterations: break diff --git a/src/UQpy/run_model/RunModel.py b/src/UQpy/run_model/RunModel.py index db426bb5f..8ddfaf548 100755 --- a/src/UQpy/run_model/RunModel.py +++ b/src/UQpy/run_model/RunModel.py @@ -22,16 +22,21 @@ import numpy as np from beartype import beartype +from enum import Enum, auto from UQpy.utilities.ValidationTypes import NumpyFloatArray +class RunType(Enum): + LOCAL = auto() + CLUSTER = auto() + class RunModel: # Authors: - # B.S. Aakash, Lohit Vandanapu, Michael D.Shields + # B.S. Aakash, Lohit Vandanapu, Michael D.Shields, Michael H. Gardner # # Last - # modified: 5 / 8 / 2020 by Michael D.Shields + # modified: 8 / 31 / 2022 by Michael H. Gardner @beartype def __init__( self, @@ -41,6 +46,8 @@ def __init__( cores_per_task: int = 1, nodes: int = 1, resume: bool = False, + run_type: str = 'LOCAL', + cluster_script: str = None ): """ Run a computational model at specified sample points. @@ -85,6 +92,9 @@ def __init__( self.model = model # Save option for resuming parallel execution self.resume = resume + # Set location for model runs + self.run_type = RunType[run_type] + self.cluster_script = cluster_script self.nodes = nodes self.ntasks = ntasks @@ -189,9 +199,20 @@ def parallel_execution(self): pickle.dump(self.model, filehandle) with open('samples.pkl', 'wb') as filehandle: pickle.dump(self.samples, filehandle) - os.system(f"mpirun python -m " - f"UQpy.run_model.model_execution.ParallelExecution {self.n_existing_simulations} " - f"{self.n_new_simulations}") + + if self.run_type is RunType.LOCAL: + os.system(f"mpirun python -m " + f"UQpy.run_model.model_execution.ParallelExecution {self.n_existing_simulations} " + f"{self.n_new_simulations}") + + elif self.run_type is RunType.CLUSTER: + if self.cluster_script is None: + raise ValueError("\nUQpy: User-provided slurm script not input, please provide this input\n") + os.system(f"python -m UQpy.run_model.model_execution.ClusterExecution {self.cores_per_task} " + f"{self.n_new_simulations} {self.n_existing_simulations} {self.cluster_script}") + else: + raise ValueError("\nUQpy: RunType is not in currently supported list of cluster types\n") + with open('qoi.pkl', 'rb') as filehandle: results = pickle.load(filehandle) diff --git a/src/UQpy/run_model/model_execution/ClusterExecution.py b/src/UQpy/run_model/model_execution/ClusterExecution.py new file mode 100644 index 000000000..b094e2841 --- /dev/null +++ b/src/UQpy/run_model/model_execution/ClusterExecution.py @@ -0,0 +1,69 @@ +# pragma: no cover +from __future__ import print_function + +import math +import sys + +import numpy as np +import os +import pickle + +try: + model = None + samples = None + samples_per_process = 0 + samples_shape = None + samples_list = None + ranges_list = None + local_ranges = None + local_samples = None + + cores_per_task = int(sys.argv[1]) + n_new_simulations = int(sys.argv[2]) + n_existing_simulations = int(sys.argv[3]) + cluster_script = str(sys.argv[4]) + + with open('model.pkl', 'rb') as filehandle: + model = pickle.load(filehandle) + + with open('samples.pkl', 'rb') as filehandle: + samples = pickle.load(filehandle) + + # Loop over the number of samples and create input files in a folder in current directory + for i in range(len(samples)): + new_text = model._find_and_replace_var_names_with_values(samples[i]) + folder_to_write = 'run_' + str(i+n_existing_simulations) + '/InputFiles' + # Write the new text to the input file + model._create_input_files(file_name=model.input_template, num=i+n_existing_simulations, + text=new_text, new_folder=folder_to_write) + + # Use model script to perform necessary preprocessing prior to model execution + for i in range(len(samples)): + sample = 'sample' # Sample input in original third-party model, though doesn't seem to use it + model.execute_single_sample(i, sample) + + # Run user-provided cluster script--for now, it is assumed the user knows how to + # tile jobs in the script + os.system(f"{cluster_script} {cores_per_task} {n_new_simulations} {n_existing_simulations}") + + results = [] + + for i in range(len(samples)): + # Change current working directory to model run directory + work_dir = os.path.join(model.model_dir, "run_" + str(i)) + # if model.verbose: + # print('\nUQpy: Changing to the following directory for output processing:\n' + work_dir) + os.chdir(work_dir) + + output = model._output_serial(i) + results.append(output) + + # Change back to model directory + os.chdir(model.model_dir) + + with open('qoi.pkl', 'wb') as filehandle: + pickle.dump(results, filehandle) + +except Exception as e: + print(e) + diff --git a/src/UQpy/run_model/model_execution/ThirdPartyModel.py b/src/UQpy/run_model/model_execution/ThirdPartyModel.py index 548ed5a83..0d09334fc 100644 --- a/src/UQpy/run_model/model_execution/ThirdPartyModel.py +++ b/src/UQpy/run_model/model_execution/ThirdPartyModel.py @@ -310,7 +310,7 @@ def _find_and_replace_var_names_with_values(self, sample): print("\nUQpy: Index Error: {0}\n".format(err)) raise IndexError("{0}".format(err)) - if isinstance(temp, collections.Iterable): + if isinstance(temp, collections.abc.Iterable): # If it is iterable, flatten and write as text file with designated separator temp = np.array(temp).flatten() to_add = "" diff --git a/src/UQpy/sampling/SimplexSampling.py b/src/UQpy/sampling/SimplexSampling.py index 2b7426635..613bba256 100644 --- a/src/UQpy/sampling/SimplexSampling.py +++ b/src/UQpy/sampling/SimplexSampling.py @@ -10,7 +10,7 @@ def __init__( self, nodes: Union[list, Numpy2DFloatArray] = None, nsamples: PositiveInteger = None, - random_state: RandomStateType = None, + random_state: Union[RandomStateType, np.random.Generator] = None, ): """ Generate uniform random samples inside an n-dimensional simplex. diff --git a/src/UQpy/sampling/stratified_sampling/RefinedStratifiedSampling.py b/src/UQpy/sampling/stratified_sampling/RefinedStratifiedSampling.py index 6a233ea1e..e891651c1 100644 --- a/src/UQpy/sampling/stratified_sampling/RefinedStratifiedSampling.py +++ b/src/UQpy/sampling/stratified_sampling/RefinedStratifiedSampling.py @@ -45,7 +45,7 @@ def __init__( self.random_state = random_state if isinstance(self.random_state, int): - self.random_state = np.random.RandomState(self.random_state) + self.random_state = np.random.default_rng(self.random_state) elif not isinstance(self.random_state, (type(None), np.random.RandomState)): raise TypeError('UQpy: random_state must be None, an int or an np.random.Generator object.') if self.random_state is None: diff --git a/src/UQpy/sampling/stratified_sampling/refinement/baseclass/Refinement.py b/src/UQpy/sampling/stratified_sampling/refinement/baseclass/Refinement.py index 9e3f37f78..0ca80e523 100644 --- a/src/UQpy/sampling/stratified_sampling/refinement/baseclass/Refinement.py +++ b/src/UQpy/sampling/stratified_sampling/refinement/baseclass/Refinement.py @@ -46,7 +46,7 @@ def finalize(self, samples, samples_per_iteration): def identify_bins(strata_metrics, points_to_add, random_state): bins2break = np.array([]) points_left = points_to_add - while (np.where(strata_metrics == strata_metrics.max())[0].shape[0] < points_left): + while np.where(strata_metrics == strata_metrics.max())[0].shape[0] < points_left: bin = np.where(strata_metrics == strata_metrics.max())[0] bins2break = np.hstack([bins2break, bin]) strata_metrics[bin] = 0 diff --git a/src/UQpy/surrogates/__init__.py b/src/UQpy/surrogates/__init__.py index 168ebdbfa..fe76e50e2 100644 --- a/src/UQpy/surrogates/__init__.py +++ b/src/UQpy/surrogates/__init__.py @@ -1,6 +1,5 @@ from UQpy.surrogates.polynomial_chaos import * from UQpy.surrogates.stochastic_reduced_order_models import * -from UQpy.surrogates.kriging import * from UQpy.surrogates.gaussian_process import * from UQpy.surrogates.baseclass import * diff --git a/src/UQpy/surrogates/kriging/Kriging.py b/src/UQpy/surrogates/kriging/Kriging.py deleted file mode 100755 index 19cbfaf1d..000000000 --- a/src/UQpy/surrogates/kriging/Kriging.py +++ /dev/null @@ -1,354 +0,0 @@ -import logging -from typing import Callable - -import numpy as np -from scipy.linalg import cholesky -import scipy.stats as stats -from beartype import beartype - -from UQpy.utilities import MinimizeOptimizer -from UQpy.utilities.Utilities import process_random_state -from UQpy.surrogates.baseclass.Surrogate import Surrogate -from UQpy.utilities.ValidationTypes import RandomStateType -from UQpy.surrogates.kriging.correlation_models.baseclass.Correlation import Correlation -from UQpy.surrogates.kriging.regression_models.baseclass.Regression import Regression - - -class Kriging(Surrogate): - @beartype - def __init__( - self, - regression_model: Regression, - correlation_model: Correlation, - correlation_model_parameters: list, - optimizer, - bounds: list = None, - optimize: bool = True, - optimizations_number: int = 1, - normalize: bool = True, - random_state: RandomStateType = None, - ): - """ - Κriging generates an Gaussian process regression-based surrogate model to predict the model output at new sample - points. - - :param regression_model: `regression_model` specifies and evaluates the basis functions and their coefficients, - which defines the trend of the model. Built-in options: :class:`Constant`, :class:`Linear`, :class:`Quadratic` - :param correlation_model: `correlation_model` specifies and evaluates the correlation function. - Built-in options: :class:`Exponential`, :class:`Gaussian`, :class:`Linear`, :class:`Spherical`, - :class:`Cubic`, :class:`Spline` - :param correlation_model_parameters: List or array of initial values for the correlation model - hyperparameters/scale parameters. - :param bounds: Bounds on the hyperparameters used to solve optimization problem to estimate maximum likelihood - estimator. This should be a closed bound. - Default: :math:`[0.001, 10^7]` for each hyperparameter. - :param optimize: Indicator to solve MLE problem or not. If :any:'True' corr_model_params will be used as initial - solution for optimization problem. Otherwise, correlation_model_parameters will be directly use as the - hyperparamters. - Default: :any:`True`. - :param optimizations_number: Number of times MLE optimization problem is to be solved with a random starting - point. Default: :math:`1`. - :param normalize: Boolean flag used in case data normalization is required. - :param optimizer: Object of the :class:`Optimizer` optimizer used during the Kriging surrogate. - Default: :class:`.MinimizeOptimizer`. - :param random_state: Random seed used to initialize the pseudo-random number generator. If an :any:`int` is - provided, this sets the seed for an object of :class:`numpy.random.RandomState`. Otherwise, the - object itself can be passed directly. - """ - self.regression_model = regression_model - self.correlation_model = correlation_model - self.correlation_model_parameters = np.array(correlation_model_parameters) - self.bounds = bounds - self.optimizer = optimizer - self.optimizations_number = optimizations_number - self.optimize = optimize - self.normalize = normalize - self.logger = logging.getLogger(__name__) - self.random_state = random_state - - # Variables are used outside the __init__ - self.samples = None - self.values = None - self.sample_mean, self.sample_std = None, None - self.value_mean, self.value_std = None, None - self.rmodel, self.cmodel = None, None - self.beta: list = None - """Regression coefficients.""" - self.gamma = None - self.err_var: float = None - """Variance of the Gaussian random process.""" - self.F_dash = None - self.C_inv = None - self.G = None - self.F, self.R = None, None - - if isinstance(self.optimizer, str): - raise ValueError("The optimization function provided a input parameter cannot be None.") - - if optimizer._bounds is None: - optimizer.update_bounds([[0.001, 10 ** 7]] * self.correlation_model_parameters.shape[0]) - - self.jac = optimizer.supports_jacobian() - self.random_state = process_random_state(random_state) - - def fit( - self, - samples, - values, - optimizations_number: int = None, - correlation_model_parameters: list = None, - ): - """ - Fit the surrogate model using the training samples and the corresponding model values. - - The user can run this method multiple time after initiating the :class:`.Kriging` class object. - - This method updates the samples and parameters of the :class:`.Kriging` object. This method uses - `correlation_model_parameters` from previous run as the starting point for MLE problem unless user provides a - new starting point. - - :param samples: :class:`numpy.ndarray` containing the training points. - :param values: :class:`numpy.ndarray` containing the model evaluations at the training points. - :param optimizations_number: number of optimization iterations - :param correlation_model_parameters: List or array of initial values for the correlation model - hyperparameters/scale parameters. - - The :meth:`fit` method has no returns, although it creates the :py:attr:`beta`, :py:attr:`err_var` and - :py:attr:`C_inv` attributes of the :class:`.Kriging` class. - """ - self.logger.info("UQpy: Running kriging.fit") - - if optimizations_number is not None: - self.optimizations_number = optimizations_number - if correlation_model_parameters is not None: - self.correlation_model_parameters = np.array(correlation_model_parameters) - self.samples = np.array(samples) - - # Number of samples and dimensions of samples and values - nsamples, input_dim = self.samples.shape - output_dim = int(np.size(values) / nsamples) - - self.values = np.array(values).reshape(nsamples, output_dim) - - # Normalizing the data - if self.normalize: - self.sample_mean, self.sample_std = np.mean(self.samples, 0), np.std(self.samples, 0) - self.value_mean, self.value_std = np.mean(self.values, 0), np.std(self.values, 0) - s_ = (self.samples - self.sample_mean) / self.sample_std - y_ = (self.values - self.value_mean) / self.value_std - else: - s_ = self.samples - y_ = self.values - - self.F, jf_ = self.regression_model.r(s_) - - # Maximum Likelihood Estimation : Solving optimization problem to calculate hyperparameters - if self.optimize: - starting_point = self.correlation_model_parameters - - minimizer, fun_value = np.zeros([self.optimizations_number, input_dim]),\ - np.zeros([self.optimizations_number, 1]) - for i__ in range(self.optimizations_number): - p_ = self.optimizer.optimize(function=Kriging.log_likelihood, - initial_guess=starting_point, - args=(self.correlation_model, s_, self.F, y_, self.jac), - jac=self.jac) - print(p_.success) - # print(self.kwargs_optimizer) - minimizer[i__, :] = p_.x - fun_value[i__, 0] = p_.fun - # Generating new starting points using log-uniform distribution - if i__ != self.optimizations_number - 1: - starting_point = stats.reciprocal.rvs([j[0] for j in self.optimizer._bounds], - [j[1] for j in self.optimizer._bounds], 1, - random_state=self.random_state) - print(starting_point) - - if min(fun_value) == np.inf: - raise NotImplementedError("Maximum likelihood estimator failed: Choose different starting point or " - "increase nopt") - t = np.argmin(fun_value) - self.correlation_model_parameters = minimizer[t, :] - - # Updated Correlation matrix corresponding to MLE estimates of hyperparameters - self.R = self.correlation_model.c(x=s_, s=s_, params=self.correlation_model_parameters) - - self.beta, self.gamma, tmp = self._compute_additional_parameters(self.R) - self.C_inv, self.F_dash, self.G, self.err_var = tmp[1], tmp[3], tmp[2], tmp[5] - - self.logger.info("UQpy: kriging fit complete.") - - def _compute_additional_parameters(self, correlation_matrix): - if self.normalize: - y_ = (self.values - self.value_mean) / self.value_std - else: - y_ = self.values - # Compute the regression coefficient (solving this linear equation: F * beta = Y) - # Eq: 3.8, DACE - c = cholesky(correlation_matrix + (10 + self.samples.shape[0]) * 2 ** (-52) * np.eye(self.samples.shape[0]), - lower=True, check_finite=False) - c_inv = np.linalg.inv(c) - f_dash = np.linalg.solve(c, self.F) - y_dash = np.linalg.solve(c, y_) - q_, g_ = np.linalg.qr(f_dash) # Eq: 3.11, DACE - # Check if F is a full rank matrix - if np.linalg.matrix_rank(g_) != min(np.size(self.F, 0), np.size(self.F, 1)): - raise NotImplementedError("Chosen regression functions are not sufficiently linearly independent") - # Design parameters (beta: regression coefficient) - beta = np.linalg.solve(g_, np.matmul(np.transpose(q_), y_dash)) - - # Design parameter (R * gamma = Y - F * beta = residual) - gamma = np.linalg.solve(c.T, (y_dash - np.matmul(f_dash, beta))) - - # Computing the process variance (Eq: 3.13, DACE) - err_var = np.zeros(self.values.shape[1]) - for i in range(self.values.shape[1]): - err_var[i] = (1 / self.samples.shape[0]) * (np.linalg.norm(y_dash[:, i] - - np.matmul(f_dash, beta[:, i])) ** 2) - - return beta, gamma, (c, c_inv, g_, f_dash, y_dash, err_var) - - def predict(self, points: np.ndarray, return_std: bool = False, correlation_model_parameters: list = None): - """ - Predict the model response at new points. - - This method evaluates the regression and correlation model at new sample points. Then, it predicts the function - value and standard deviation. - - :param points: Points at which to predict the model response. - :param return_std: Indicator to estimate standard deviation. - :param correlation_model_parameters: Hyperparameters for correlation model. - :return: Predicted values at the new points, Standard deviation of predicted values at the new points - """ - x_ = np.atleast_2d(points) - if self.normalize: - x_ = (x_ - self.sample_mean) / self.sample_std - s_ = (self.samples - self.sample_mean) / self.sample_std - else: - s_ = self.samples - fx, jf = self.regression_model.r(x_) - if correlation_model_parameters is None: - correlation_model_parameters = self.correlation_model_parameters - rx = self.correlation_model.c( - x=x_, s=s_, params=correlation_model_parameters - ) - if correlation_model_parameters is None: - beta, gamma = self.beta, self.gamma - c_inv, f_dash, g_, err_var = self.C_inv, self.F_dash, self.G, self.err_var - else: - beta, gamma, tmp = self._compute_additional_parameters( - self.correlation_model.c(x=s_, s=s_, params=correlation_model_parameters)) - c_inv, f_dash, g_, err_var = tmp[1], tmp[3], tmp[2], tmp[5] - y = np.einsum("ij,jk->ik", fx, beta) + np.einsum( - "ij,jk->ik", rx, gamma - ) - if self.normalize: - y = self.value_mean + y * self.value_std - if x_.shape[1] == 1: - y = y.flatten() - if return_std: - r_dash = np.matmul(c_inv, rx.T) - u = np.matmul(f_dash.T, r_dash) - fx.T - norm1 = np.linalg.norm(r_dash, 2, 0) - norm2 = np.linalg.norm(np.linalg.solve(g_, u), 2, 0) - mse = np.sqrt(err_var * np.atleast_2d(1 + norm2 - norm1).T) - if self.normalize: - mse = self.value_std * mse - if x_.shape[1] == 1: - mse = mse.flatten() - return y, mse - else: - return y - - def jacobian(self, points: np.ndarray): - """ - Predict the gradient of the model at new points. - - This method evaluates the regression and correlation model at new sample point. Then, it predicts the gradient - using the regression coefficients and the training second_order_tensor. - - :param points: Points at which to evaluate the gradient. - :return: Gradient of the surrogate model evaluated at the new points. - """ - x_ = np.atleast_2d(points) - if self.normalize: - x_ = (x_ - self.sample_mean) / self.sample_std - s_ = (self.samples - self.sample_mean) / self.sample_std - else: - s_ = self.samples - - fx, jf = self.regression_model.r(x_) - rx, drdx = self.correlation_model.c( - x=x_, s=s_, params=self.correlation_model_parameters, dx=True - ) - y_grad = np.einsum("ikj,jm->ik", jf, self.beta) + np.einsum( - "ijk,jm->ki", drdx.T, self.gamma - ) - if self.normalize: - y_grad = y_grad * self.value_std / self.sample_std - if x_.shape[1] == 1: - y_grad = y_grad.flatten() - return y_grad - - @staticmethod - def log_likelihood(p0, cm, s, f, y, return_grad): - # Return the log-likelihood function and it's gradient. Gradient is calculate using Central Difference - m = s.shape[0] - n = s.shape[1] - r__, dr_ = cm.c(x=s, s=s, params=p0, dt=True) - try: - cc = cholesky(r__ + 2 ** (-52) * np.eye(m), lower=True) - except np.linalg.LinAlgError: - return np.inf, np.zeros(n) - - # Product of diagonal terms is negligible sometimes, even when cc exists. - if np.prod(np.diagonal(cc)) == 0: - return np.inf, np.zeros(n) - - cc_inv = np.linalg.inv(cc) - r_inv = np.matmul(cc_inv.T, cc_inv) - f__ = cc_inv.dot(f) - y__ = cc_inv.dot(y) - - q__, g__ = np.linalg.qr(f__) # Eq: 3.11, DACE - - # Check if F is a full rank matrix - if np.linalg.matrix_rank(g__) != min(np.size(f__, 0), np.size(f__, 1)): - raise NotImplementedError( - "Chosen regression functions are not sufficiently linearly independent" - ) - - # Design parameters - beta_ = np.linalg.solve(g__, np.matmul(np.transpose(q__), y__)) - - # Computing the process variance (Eq: 3.13, DACE) - sigma_ = np.zeros(y.shape[1]) - - ll = 0 - for out_dim in range(y.shape[1]): - sigma_[out_dim] = (1 / m) * ( - np.linalg.norm(y__[:, out_dim] - np.matmul(f__, beta_[:, out_dim])) ** 2) - # Objective function:= log(det(sigma**2 * R)) + constant - ll = (ll + ( np.log(np.linalg.det(sigma_[out_dim] * r__)) + m * (np.log(2 * np.pi) + 1)) / 2) - - # Gradient of loglikelihood - # Reference: C. E. Rasmussen & C. K. I. Williams, Gaussian Processes for Machine Learning, the MIT Press, - # 2006, ISBN 026218253X. (Page 114, Eq.(5.9)) - residual = y - np.matmul(f, beta_) - gamma = np.matmul(r_inv, residual) - grad_mle = np.zeros(n) - for in_dim in range(n): - r_inv_derivative = np.matmul(r_inv, np.matmul(dr_[:, :, in_dim], r_inv)) - tmp = np.matmul(residual.T, np.matmul(r_inv_derivative, residual)) - for out_dim in range(y.shape[1]): - alpha = gamma / sigma_[out_dim] - tmp1 = np.matmul(alpha, alpha.T) - r_inv / sigma_[out_dim] - cov_der = sigma_[out_dim] * dr_[:, :, in_dim] + tmp * r__ / m - grad_mle[in_dim] = grad_mle[in_dim] - 0.5 * np.trace( - np.matmul(tmp1, cov_der) - ) - - if return_grad: - return ll, grad_mle - else: - return ll diff --git a/src/UQpy/surrogates/kriging/__init__.py b/src/UQpy/surrogates/kriging/__init__.py deleted file mode 100644 index 55a50199d..000000000 --- a/src/UQpy/surrogates/kriging/__init__.py +++ /dev/null @@ -1,4 +0,0 @@ -from UQpy.surrogates.kriging.Kriging import Kriging - -from UQpy.surrogates.kriging.regression_models import * -from UQpy.surrogates.kriging.correlation_models import * diff --git a/src/UQpy/surrogates/kriging/correlation_models/CubicCorrelation.py b/src/UQpy/surrogates/kriging/correlation_models/CubicCorrelation.py deleted file mode 100644 index 3e909507f..000000000 --- a/src/UQpy/surrogates/kriging/correlation_models/CubicCorrelation.py +++ /dev/null @@ -1,28 +0,0 @@ -from UQpy.surrogates.kriging.correlation_models.baseclass.Correlation import * - - -class CubicCorrelation(Correlation): - def c(self, x, s, params, dt=False, dx=False): - zeta_matrix, dtheta_derivs, dx_derivs = Correlation.derivatives( - x_=x, s_=s, params=params - ) - # Initial matrices containing derivates for all values in array. Note since - # dtheta_s and dx_s already accounted for where derivative should be zero, all - # that must be done is multiplying the |dij| or thetaj matrix on top of a - # matrix of derivates w.r.t zeta (in this case, dzeta = -6zeta+6zeta**2) - drdt = (-6 * zeta_matrix + 6 * zeta_matrix ** 2) * dtheta_derivs - drdx = (-6 * zeta_matrix + 6 * zeta_matrix ** 2) * dx_derivs - # Also, create matrix for values of equation, 1 - 3zeta**2 + 2zeta**3, for loop - zeta_function_cubic = 1 - 3 * zeta_matrix ** 2 + 2 * zeta_matrix ** 3 - rx = np.prod(zeta_function_cubic, 2) - if dt: - # Same as previous example, loop over zeta matrix by shifting index - for i in range(len(params) - 1): - drdt = drdt * np.roll(zeta_function_cubic, i + 1, axis=2) - return rx, drdt - if dx: - # Same as previous example, loop over zeta matrix by shifting index - for i in range(len(params) - 1): - drdx = drdx * np.roll(zeta_function_cubic, i + 1, axis=2) - return rx, drdx - return rx diff --git a/src/UQpy/surrogates/kriging/correlation_models/ExponentialCorrelation.py b/src/UQpy/surrogates/kriging/correlation_models/ExponentialCorrelation.py deleted file mode 100644 index 94702760c..000000000 --- a/src/UQpy/surrogates/kriging/correlation_models/ExponentialCorrelation.py +++ /dev/null @@ -1,20 +0,0 @@ -from UQpy.surrogates.kriging.correlation_models.baseclass.Correlation import * - - -class ExponentialCorrelation(Correlation): - def c(self, x, s, params, dt=False, dx=False): - stack = Correlation.check_samples_and_return_stack(x, s) - rx = np.exp(np.sum(-params * abs(stack), axis=2)) - if dt: - drdt = -abs(stack) * np.transpose( - np.tile(rx, (np.size(x, 1), 1, 1)), (1, 2, 0) - ) - return rx, drdt - if dx: - drdx = ( - -params - * np.sign(stack) - * np.transpose(np.tile(rx, (np.size(x, 1), 1, 1)), (1, 2, 0)) - ) - return rx, drdx - return rx diff --git a/src/UQpy/surrogates/kriging/correlation_models/GaussianCorrelation.py b/src/UQpy/surrogates/kriging/correlation_models/GaussianCorrelation.py deleted file mode 100644 index 05ce09830..000000000 --- a/src/UQpy/surrogates/kriging/correlation_models/GaussianCorrelation.py +++ /dev/null @@ -1,21 +0,0 @@ -from UQpy.surrogates.kriging.correlation_models.baseclass.Correlation import * - - -class GaussianCorrelation(Correlation): - def c(self, x, s, params, dt=False, dx=False): - stack = Correlation.check_samples_and_return_stack(x, s) - rx = np.exp(np.sum(-params * (stack ** 2), axis=2)) - if dt: - drdt = -(stack ** 2) * np.transpose( - np.tile(rx, (np.size(x, 1), 1, 1)), (1, 2, 0) - ) - return rx, drdt - if dx: - drdx = ( - -2 - * params - * stack - * np.transpose(np.tile(rx, (np.size(x, 1), 1, 1)), (1, 2, 0)) - ) - return rx, drdx - return rx diff --git a/src/UQpy/surrogates/kriging/correlation_models/LinearCorrelation.py b/src/UQpy/surrogates/kriging/correlation_models/LinearCorrelation.py deleted file mode 100644 index 69d7f1506..000000000 --- a/src/UQpy/surrogates/kriging/correlation_models/LinearCorrelation.py +++ /dev/null @@ -1,35 +0,0 @@ -from UQpy.surrogates.kriging.correlation_models.baseclass.Correlation import * - - -class LinearCorrelation(Correlation): - def c(self, x, s, params, dt=False, dx=False): - stack = Correlation.check_samples_and_return_stack(x, s) - # Taking stack and turning each d value into 1-theta*dij - after_parameters = 1 - params * abs(stack) - # Define matrix of zeros to compare against (not necessary to be defined separately, - # but the line is bulky if this isn't defined first, and it is used more than once) - comp_zero = np.zeros((np.size(x, 0), np.size(s, 0), np.size(s, 1))) - # Compute matrix of max{0,1-theta*d} - max_matrix = np.maximum(after_parameters, comp_zero) - rx = np.prod(max_matrix, 2) - # Create matrix that has 1s where max_matrix is nonzero - # -Essentially, this acts as a way to store the indices of where the values are nonzero - ones_and_zeros = max_matrix.astype(bool).astype(int) - # Set initial derivatives as if all were positive - first_dtheta = -abs(stack) - first_dx = np.negative(params) * np.sign(stack) - # Multiply derivs by ones_and_zeros...this will set the values where the - # derivative should be zero to zero, and keep all other values the same - drdt = np.multiply(first_dtheta, ones_and_zeros) - drdx = np.multiply(first_dx, ones_and_zeros) - if dt: - # Loop over parameters, shifting max_matrix and multiplying over derivative matrix with each iter - for i in range(len(params) - 1): - drdt = drdt * np.roll(max_matrix, i + 1, axis=2) - return rx, drdt - if dx: - # Loop over parameters, shifting max_matrix and multiplying over derivative matrix with each iter - for i in range(len(params) - 1): - drdx = drdx * np.roll(max_matrix, i + 1, axis=2) - return rx, drdx - return rx diff --git a/src/UQpy/surrogates/kriging/correlation_models/SphericalCorrelation.py b/src/UQpy/surrogates/kriging/correlation_models/SphericalCorrelation.py deleted file mode 100644 index 1f6b8173d..000000000 --- a/src/UQpy/surrogates/kriging/correlation_models/SphericalCorrelation.py +++ /dev/null @@ -1,28 +0,0 @@ -from UQpy.surrogates.kriging.correlation_models.baseclass.Correlation import * - - -class SphericalCorrelation(Correlation): - def c(self, x, s, params, dt=False, dx=False): - zeta_matrix, dtheta_derivs, dx_derivs = Correlation.derivatives( - x_=x, s_=s, params=params - ) - # Initial matrices containing derivates for all values in array. Note since - # dtheta_s and dx_s already accounted for where derivative should be zero, all - # that must be done is multiplying the |dij| or thetaj matrix on top of a - # matrix of derivates w.r.t zeta (in this case, dzeta = -1.5+1.5zeta**2) - drdt = (-1.5 + 1.5 * zeta_matrix ** 2) * dtheta_derivs - drdx = (-1.5 + 1.5 * zeta_matrix ** 2) * dx_derivs - # Also, create matrix for values of equation, 1 - 1.5zeta + 0.5zeta**3, for loop - zeta_function = 1 - 1.5 * zeta_matrix + 0.5 * zeta_matrix ** 3 - rx = np.prod(zeta_function, 2) - if dt: - # Same as previous example, loop over zeta matrix by shifting index - for i in range(len(params) - 1): - drdt = drdt * np.roll(zeta_function, i + 1, axis=2) - return rx, drdt - if dx: - # Same as previous example, loop over zeta matrix by shifting index - for i in range(len(params) - 1): - drdx = drdx * np.roll(zeta_function, i + 1, axis=2) - return rx, drdx - return rx diff --git a/src/UQpy/surrogates/kriging/correlation_models/SplineCorrelation.py b/src/UQpy/surrogates/kriging/correlation_models/SplineCorrelation.py deleted file mode 100644 index 0aa6282d1..000000000 --- a/src/UQpy/surrogates/kriging/correlation_models/SplineCorrelation.py +++ /dev/null @@ -1,58 +0,0 @@ -from UQpy.surrogates.kriging.correlation_models.baseclass.Correlation import * - - -class SplineCorrelation(Correlation): - def c(self, x, s, params, dt=False, dx=False): - # x_, s_ = np.atleast_2d(x_), np.atleast_2d(s_) - # # Create stack matrix, where each block is x_i with all s - # stack = np.tile(np.swapaxes(np.atleast_3d(x_), 1, 2), (1, np.size(s_, 0), 1)) - np.tile(s_, ( - # np.size(x_, 0), - # 1, 1)) - stack = Correlation.check_samples_and_return_stack(x, s) - # In this case, the zeta value is just abs(stack)*parameters, no comparison - zeta_matrix = abs(stack) * params - # So, dtheta and dx are just |dj| and theta*sgn(dj), respectively - dtheta_derivs = abs(stack) - # dx_derivs = np.ones((np.size(x,0),np.size(s,0),np.size(s,1)))*parameters - dx_derivs = np.sign(stack) * params - - # Initialize empty sigma and dsigma matrices - sigma = np.ones( - (zeta_matrix.shape[0], zeta_matrix.shape[1], zeta_matrix.shape[2]) - ) - dsigma = np.ones( - (zeta_matrix.shape[0], zeta_matrix.shape[1], zeta_matrix.shape[2]) - ) - - # Loop over cases to create zeta_matrix and subsequent dR matrices - for i in range(zeta_matrix.shape[0]): - for j in range(zeta_matrix.shape[1]): - for k in range(zeta_matrix.shape[2]): - y = zeta_matrix[i, j, k] - if 0 <= y <= 0.2: - sigma[i, j, k] = 1 - 15 * y ** 2 + 30 * y ** 3 - dsigma[i, j, k] = -30 * y + 90 * y ** 2 - elif 0.2 < y < 1.0: - sigma[i, j, k] = 1.25 * (1 - y) ** 3 - dsigma[i, j, k] = 3.75 * (1 - y) ** 2 * -1 - elif y >= 1: - sigma[i, j, k] = 0 - dsigma[i, j, k] = 0 - - rx = np.prod(sigma, 2) - - if dt: - # Initialize derivative matrices incorporating chain rule - drdt = dsigma * dtheta_derivs - # Loop over to create proper matrices - for i in range(len(params) - 1): - drdt = drdt * np.roll(sigma, i + 1, axis=2) - return rx, drdt - if dx: - # Initialize derivative matrices incorporating chain rule - drdx = dsigma * dx_derivs - # Loop over to create proper matrices - for i in range(len(params) - 1): - drdx = drdx * np.roll(sigma, i + 1, axis=2) - return rx, drdx - return rx diff --git a/src/UQpy/surrogates/kriging/correlation_models/__init__.py b/src/UQpy/surrogates/kriging/correlation_models/__init__.py deleted file mode 100644 index 10f39dafc..000000000 --- a/src/UQpy/surrogates/kriging/correlation_models/__init__.py +++ /dev/null @@ -1,7 +0,0 @@ -from UQpy.surrogates.kriging.correlation_models.baseclass import * -from UQpy.surrogates.kriging.correlation_models.CubicCorrelation import CubicCorrelation -from UQpy.surrogates.kriging.correlation_models.ExponentialCorrelation import ExponentialCorrelation -from UQpy.surrogates.kriging.correlation_models.GaussianCorrelation import GaussianCorrelation -from UQpy.surrogates.kriging.correlation_models.LinearCorrelation import LinearCorrelation -from UQpy.surrogates.kriging.correlation_models.SphericalCorrelation import SphericalCorrelation -from UQpy.surrogates.kriging.correlation_models.SplineCorrelation import SplineCorrelation diff --git a/src/UQpy/surrogates/kriging/correlation_models/baseclass/Correlation.py b/src/UQpy/surrogates/kriging/correlation_models/baseclass/Correlation.py deleted file mode 100644 index 703461b5f..000000000 --- a/src/UQpy/surrogates/kriging/correlation_models/baseclass/Correlation.py +++ /dev/null @@ -1,47 +0,0 @@ -from abc import ABC, abstractmethod -import numpy as np - - -class Correlation(ABC): - """ - Abstract base class of all Correlations. Serves as a template for creating new Kriging correlation - functions. - """ - - @abstractmethod - def c(self, x, s, params, dt=False, dx=False): - """ - Abstract method that needs to be implemented by the user when creating a new Correlation function. - """ - pass - - @staticmethod - def check_samples_and_return_stack(x, s): - x_, s_ = np.atleast_2d(x), np.atleast_2d(s) - # Create stack matrix, where each block is x_i with all s - stack = np.tile( - np.swapaxes(np.atleast_3d(x_), 1, 2), (1, np.size(s_, 0), 1) - ) - np.tile(s_, (np.size(x_, 0), 1, 1)) - return stack - - @staticmethod - def derivatives(x_, s_, params): - stack = Correlation.check_samples_and_return_stack(x_, s_) - # Taking stack and creating array of all thetaj*dij - after_parameters = params * abs(stack) - # Create matrix of all ones to compare - comp_ones = np.ones((np.size(x_, 0), np.size(s_, 0), np.size(s_, 1))) - # zeta_matrix has all values min{1,theta*dij} - zeta_matrix_ = np.minimum(after_parameters, comp_ones) - # Copy zeta_matrix to another matrix that will used to find where derivative should be zero - indices = zeta_matrix_.copy() - # If value of min{1,theta*dij} is 1, the derivative should be 0. - # So, replace all values of 1 with 0, then perform the .astype(bool).astype(int) - # operation like in the linear example, so you end up with an array of 1's where - # the derivative should be caluclated and 0 where it should be zero - indices[indices == 1] = 0 - # Create matrix of all |dij| (where non zero) to be used in calculation of dR/dtheta - dtheta_derivs_ = indices.astype(bool).astype(int) * abs(stack) - # Same as above, but for matrix of all thetaj where non-zero - dx_derivs_ = indices.astype(bool).astype(int) * params * np.sign(stack) - return zeta_matrix_, dtheta_derivs_, dx_derivs_ diff --git a/src/UQpy/surrogates/kriging/correlation_models/baseclass/__init__.py b/src/UQpy/surrogates/kriging/correlation_models/baseclass/__init__.py deleted file mode 100644 index e8cf1815d..000000000 --- a/src/UQpy/surrogates/kriging/correlation_models/baseclass/__init__.py +++ /dev/null @@ -1 +0,0 @@ -from UQpy.surrogates.kriging.correlation_models.baseclass.Correlation import Correlation diff --git a/src/UQpy/surrogates/kriging/regression_models/ConstantRegression.py b/src/UQpy/surrogates/kriging/regression_models/ConstantRegression.py deleted file mode 100644 index 0e4f9e984..000000000 --- a/src/UQpy/surrogates/kriging/regression_models/ConstantRegression.py +++ /dev/null @@ -1,10 +0,0 @@ -import numpy as np -from UQpy.surrogates.kriging.regression_models.baseclass.Regression import Regression - - -class ConstantRegression(Regression): - def r(self, s): - s = np.atleast_2d(s) - fx = np.ones([np.size(s, 0), 1]) - jf = np.zeros([np.size(s, 0), np.size(s, 1), 1]) - return fx, jf diff --git a/src/UQpy/surrogates/kriging/regression_models/LinearRegression.py b/src/UQpy/surrogates/kriging/regression_models/LinearRegression.py deleted file mode 100644 index 118d8d73c..000000000 --- a/src/UQpy/surrogates/kriging/regression_models/LinearRegression.py +++ /dev/null @@ -1,12 +0,0 @@ -import numpy as np -from UQpy.surrogates.kriging.regression_models.baseclass.Regression import Regression - - -class LinearRegression(Regression): - def r(self, s): - s = np.atleast_2d(s) - fx = np.concatenate((np.ones([np.size(s, 0), 1]), s), 1) - jf_b = np.zeros([np.size(s, 0), np.size(s, 1), np.size(s, 1)]) - np.einsum("jii->ji", jf_b)[:] = 1 - jf = np.concatenate((np.zeros([np.size(s, 0), np.size(s, 1), 1]), jf_b), 2) - return fx, jf diff --git a/src/UQpy/surrogates/kriging/regression_models/QuadraticRegression.py b/src/UQpy/surrogates/kriging/regression_models/QuadraticRegression.py deleted file mode 100644 index fdddefbb5..000000000 --- a/src/UQpy/surrogates/kriging/regression_models/QuadraticRegression.py +++ /dev/null @@ -1,37 +0,0 @@ -import numpy as np -from UQpy.surrogates.kriging.regression_models.baseclass.Regression import Regression - - -class QuadraticRegression(Regression): - def r(self, s): - s = np.atleast_2d(s) - fx = np.zeros( - [np.size(s, 0), int((np.size(s, 1) + 1) * (np.size(s, 1) + 2) / 2)] - ) - jf = np.zeros( - [ - np.size(s, 0), - np.size(s, 1), - int((np.size(s, 1) + 1) * (np.size(s, 1) + 2) / 2), - ] - ) - for i in range(np.size(s, 0)): - temp = np.hstack((1, s[i, :])) - for j in range(np.size(s, 1)): - temp = np.hstack((temp, s[i, j] * s[i, j::])) - fx[i, :] = temp - # definie H matrix - h_ = 0 - for j in range(np.size(s, 1)): - tmp_ = s[i, j] * np.eye(np.size(s, 1)) - t1 = np.zeros([np.size(s, 1), np.size(s, 1)]) - t1[j, :] = s[i, :] - tmp = tmp_ + t1 - if j == 0: - h_ = tmp[:, j::] - else: - h_ = np.hstack((h_, tmp[:, j::])) - jf[i, :, :] = np.hstack( - (np.zeros([np.size(s, 1), 1]), np.eye(np.size(s, 1)), h_) - ) - return fx, jf diff --git a/src/UQpy/surrogates/kriging/regression_models/__init__.py b/src/UQpy/surrogates/kriging/regression_models/__init__.py deleted file mode 100644 index e6da265b3..000000000 --- a/src/UQpy/surrogates/kriging/regression_models/__init__.py +++ /dev/null @@ -1,4 +0,0 @@ -from UQpy.surrogates.kriging.regression_models.baseclass import * -from UQpy.surrogates.kriging.regression_models.ConstantRegression import ConstantRegression -from UQpy.surrogates.kriging.regression_models.LinearRegression import LinearRegression -from UQpy.surrogates.kriging.regression_models.QuadraticRegression import QuadraticRegression diff --git a/src/UQpy/surrogates/kriging/regression_models/baseclass/Regression.py b/src/UQpy/surrogates/kriging/regression_models/baseclass/Regression.py deleted file mode 100644 index 3d435e51d..000000000 --- a/src/UQpy/surrogates/kriging/regression_models/baseclass/Regression.py +++ /dev/null @@ -1,14 +0,0 @@ -from abc import ABC, abstractmethod - - -class Regression(ABC): - """ - Abstract base class of all Regressions. Serves as a template for creating new Kriging regression - functions. - """ - @abstractmethod - def r(self, s): - """ - Abstract method that needs to be implemented by the user when creating a new Regression function. - """ - pass diff --git a/src/UQpy/surrogates/kriging/regression_models/baseclass/__init__.py b/src/UQpy/surrogates/kriging/regression_models/baseclass/__init__.py deleted file mode 100644 index 004e95e23..000000000 --- a/src/UQpy/surrogates/kriging/regression_models/baseclass/__init__.py +++ /dev/null @@ -1 +0,0 @@ -from UQpy.surrogates.kriging.regression_models.baseclass.Regression import Regression diff --git a/src/UQpy/transformations/Nataf.py b/src/UQpy/transformations/Nataf.py index 90c3adf2c..66bd16c71 100644 --- a/src/UQpy/transformations/Nataf.py +++ b/src/UQpy/transformations/Nataf.py @@ -32,6 +32,7 @@ def __init__( itam_threshold1: Union[float, int] = 0.001, itam_threshold2: Union[float, int] = 0.1, itam_max_iter: int = 100, + n_gauss_points: int = 128 ): """ Transform random variables using the Nataf or Inverse Nataf transformation @@ -63,6 +64,8 @@ def __init__( for iteration :math:`i`. Default: :math:`0.01` :param itam_max_iter: Maximum number of iterations for the ITAM method. Default: :math:`100` + :param n_gauss_points: The number of integration points used for the numerical integration of the + correlation matrix (:math:`\mathbf{C_Z}`) of the standard normal random vector **Z** """ self.n_dimensions = 0 if isinstance(distributions, list): @@ -108,7 +111,7 @@ def __init__( elif all(isinstance(x, Normal) for x in distributions): self.corr_x = self.corr_z else: - self.corr_x = self.distortion_z2x(self.dist_object, self.corr_z, n_gauss_points=128) + self.corr_x = self.distortion_z2x(self.dist_object, self.corr_z, n_gauss_points=n_gauss_points) self.H: NumpyFloatArray = cholesky(self.corr_z, lower=True) """The lower triangular matrix resulting from the Cholesky decomposition of the correlation matrix diff --git a/src/UQpy/utilities/MinimizeOptimizer.py b/src/UQpy/utilities/MinimizeOptimizer.py index aa76a99ec..d15904b01 100644 --- a/src/UQpy/utilities/MinimizeOptimizer.py +++ b/src/UQpy/utilities/MinimizeOptimizer.py @@ -24,11 +24,11 @@ def optimize(self, function, initial_guess, args=(), jac=False): return minimize(function, initial_guess, args=args, method=self.method, bounds=self._bounds, constraints=self.constraints, jac=jac, - options={'disp': True, 'maxiter': 10000, 'catol': 0.002}) + options={'disp': False, 'maxiter': 10000, 'catol': 0.002}) else: return minimize(function, initial_guess, args=args, method=self.method, bounds=self._bounds, jac=jac, - options={'disp': True, 'maxiter': 10000, 'catol': 0.002}) + options={'disp': False, 'maxiter': 10000, 'catol': 0.002}) def apply_constraints(self, constraints): if self.method.lower() in ['cobyla', 'slsqp', 'trust-constr']: diff --git a/tests/unit_tests/sampling/test_adaptive_kriging.py b/tests/unit_tests/sampling/test_adaptive_kriging.py index ef932fb7f..4fcf3b1e6 100644 --- a/tests/unit_tests/sampling/test_adaptive_kriging.py +++ b/tests/unit_tests/sampling/test_adaptive_kriging.py @@ -1,179 +1,166 @@ import pytest +from UQpy import GaussianProcessRegression, RBF, LinearRegression from UQpy.run_model.model_execution.PythonModel import PythonModel from UQpy.utilities.MinimizeOptimizer import MinimizeOptimizer -from UQpy.surrogates.kriging.Kriging import Kriging from UQpy.sampling import MonteCarloSampling, AdaptiveKriging from UQpy.run_model.RunModel import RunModel from UQpy.distributions.collection import Normal from UQpy.sampling.adaptive_kriging_functions import * -import shutil def test_akmcs_weighted_u(): - from UQpy.surrogates.kriging.regression_models.LinearRegression import LinearRegression - from UQpy.surrogates.kriging.correlation_models.ExponentialCorrelation import ExponentialCorrelation - marginals = [Normal(loc=0., scale=4.), Normal(loc=0., scale=4.)] x = MonteCarloSampling(distributions=marginals, nsamples=20, random_state=0) model = PythonModel(model_script='series.py', model_object_name="series") rmodel = RunModel(model=model) - regression_model = LinearRegression() - correlation_model = ExponentialCorrelation() - K = Kriging(regression_model=regression_model, correlation_model=correlation_model, - optimizer=MinimizeOptimizer('l-bfgs-b'), - optimizations_number=10, correlation_model_parameters=[1, 1], random_state=1) + + kernel1 = RBF() + bounds_1 = [[10 ** (-4), 10 ** 3], [10 ** (-3), 10 ** 2], [10 ** (-3), 10 ** 2]] + optimizer1 = MinimizeOptimizer(method='L-BFGS-B', bounds=bounds_1) + gpr = GaussianProcessRegression(kernel=kernel1, hyperparameters=[1, 10 ** (-3), 10 ** (-2)], optimizer=optimizer1, + optimizations_number=10, noise=False, regression_model=LinearRegression(), + random_state=1) # OPTIONS: 'U', 'EFF', 'Weighted-U' learning_function = WeightedUFunction(weighted_u_stop=2) - a = AdaptiveKriging(distributions=marginals, runmodel_object=rmodel, surrogate=K, - learning_nsamples=10**3, n_add=1, learning_function=learning_function, + a = AdaptiveKriging(distributions=marginals, runmodel_object=rmodel, surrogate=gpr, + learning_nsamples=10 ** 3, n_add=1, learning_function=learning_function, random_state=2) a.run(nsamples=25, samples=x.samples) - assert a.samples[23, 0] == 1.083176685073489 - assert a.samples[20, 1] == 0.20293978126855253 - + assert a.samples[23, 0] == -0.48297825309989356 + assert a.samples[20, 1] == 0.39006110248010434 def test_akmcs_u(): - from UQpy.surrogates.kriging.regression_models.LinearRegression import LinearRegression - from UQpy.surrogates.kriging.correlation_models.ExponentialCorrelation import ExponentialCorrelation marginals = [Normal(loc=0., scale=4.), Normal(loc=0., scale=4.)] x = MonteCarloSampling(distributions=marginals, nsamples=20, random_state=1) model = PythonModel(model_script='series.py', model_object_name="series") rmodel = RunModel(model=model) - regression_model = LinearRegression() - correlation_model = ExponentialCorrelation() - K = Kriging(regression_model=regression_model, correlation_model=correlation_model, - optimizer=MinimizeOptimizer('l-bfgs-b'), - optimizations_number=10, correlation_model_parameters=[1, 1], random_state=0) + kernel1 = RBF() + bounds_1 = [[10 ** (-4), 10 ** 3], [10 ** (-3), 10 ** 2], [10 ** (-3), 10 ** 2]] + optimizer1 = MinimizeOptimizer(method='L-BFGS-B', bounds=bounds_1) + gpr = GaussianProcessRegression(kernel=kernel1, hyperparameters=[1, 10 ** (-3), 10 ** (-2)], optimizer=optimizer1, + optimizations_number=100, noise=False, regression_model=LinearRegression(), + random_state=0) # OPTIONS: 'U', 'EFF', 'Weighted-U' learning_function = UFunction(u_stop=2) - a = AdaptiveKriging(distributions=marginals, runmodel_object=rmodel, surrogate=K, - learning_nsamples=10**3, n_add=1, learning_function=learning_function, + a = AdaptiveKriging(distributions=marginals, runmodel_object=rmodel, surrogate=gpr, + learning_nsamples=10 ** 3, n_add=1, learning_function=learning_function, random_state=2) a.run(nsamples=25, samples=x.samples) - assert a.samples[23, 0] == -4.141979058326188 - assert a.samples[20, 1] == -1.6476534435429009 - + assert a.samples[23, 0] == -3.781937137406927 + assert a.samples[20, 1] == 0.17610325620498946 def test_akmcs_expected_feasibility(): - from UQpy.surrogates.kriging.regression_models.LinearRegression import LinearRegression - from UQpy.surrogates.kriging.correlation_models.ExponentialCorrelation import ExponentialCorrelation marginals = [Normal(loc=0., scale=4.), Normal(loc=0., scale=4.)] x = MonteCarloSampling(distributions=marginals, nsamples=20, random_state=1) model = PythonModel(model_script='series.py', model_object_name="series") rmodel = RunModel(model=model) - regression_model = LinearRegression() - correlation_model = ExponentialCorrelation() - K = Kriging(regression_model=regression_model, correlation_model=correlation_model, - optimizations_number=10, correlation_model_parameters=[1, 1], - optimizer=MinimizeOptimizer('l-bfgs-b'),) + kernel1 = RBF() + bounds_1 = [[10 ** (-4), 10 ** 3], [10 ** (-3), 10 ** 2], [10 ** (-3), 10 ** 2]] + optimizer1 = MinimizeOptimizer(method='L-BFGS-B', bounds=bounds_1) + gpr = GaussianProcessRegression(kernel=kernel1, hyperparameters=[1, 10 ** (-3), 10 ** (-2)], optimizer=optimizer1, + optimizations_number=100, noise=False, regression_model=LinearRegression(), + random_state=0) # OPTIONS: 'U', 'EFF', 'Weighted-U' learning_function = ExpectedFeasibility(eff_a=0, eff_epsilon=2, eff_stop=0.001) - a = AdaptiveKriging(distributions=marginals, runmodel_object=rmodel, surrogate=K, - learning_nsamples=10**3, n_add=1, learning_function=learning_function, + a = AdaptiveKriging(distributions=marginals, runmodel_object=rmodel, surrogate=gpr, + learning_nsamples=10 ** 3, n_add=1, learning_function=learning_function, random_state=2) a.run(nsamples=25, samples=x.samples) - assert a.samples[23, 0] == 1.366058523912817 - assert a.samples[20, 1] == -12.914668932772358 - + assert a.samples[23, 0] == 5.423754197908594 + assert a.samples[20, 1] == 2.0355505295053384 def test_akmcs_expected_improvement(): - from UQpy.surrogates.kriging.regression_models.LinearRegression import LinearRegression - from UQpy.surrogates.kriging.correlation_models.ExponentialCorrelation import ExponentialCorrelation marginals = [Normal(loc=0., scale=4.), Normal(loc=0., scale=4.)] x = MonteCarloSampling(distributions=marginals, nsamples=20, random_state=1) model = PythonModel(model_script='series.py', model_object_name="series") rmodel = RunModel(model=model) - regression_model = LinearRegression() - correlation_model = ExponentialCorrelation() - K = Kriging(regression_model=regression_model, correlation_model=correlation_model, - optimizations_number=10, correlation_model_parameters=[1, 1], - optimizer=MinimizeOptimizer('l-bfgs-b'),) + kernel1 = RBF() + bounds_1 = [[10 ** (-4), 10 ** 3], [10 ** (-3), 10 ** 2], [10 ** (-3), 10 ** 2]] + optimizer1 = MinimizeOptimizer(method='L-BFGS-B', bounds=bounds_1) + gpr = GaussianProcessRegression(kernel=kernel1, hyperparameters=[1, 10 ** (-3), 10 ** (-2)], optimizer=optimizer1, + optimizations_number=50, noise=False, regression_model=LinearRegression(), + random_state=0) # OPTIONS: 'U', 'EFF', 'Weighted-U' learning_function = ExpectedImprovement() - a = AdaptiveKriging(distributions=marginals, runmodel_object=rmodel, surrogate=K, - learning_nsamples=10**3, n_add=1, learning_function=learning_function, + a = AdaptiveKriging(distributions=marginals, runmodel_object=rmodel, surrogate=gpr, + learning_nsamples=10 ** 3, n_add=1, learning_function=learning_function, random_state=2) a.run(nsamples=25, samples=x.samples) - assert a.samples[23, 0] == 4.553078100499578 - assert a.samples[20, 1] == -3.508949564718469 - + assert a.samples[21, 0] == 6.878734574049913 + assert a.samples[20, 1] == -6.3410533857909215 def test_akmcs_expected_improvement_global_fit(): - from UQpy.surrogates.kriging.regression_models.LinearRegression import LinearRegression - from UQpy.surrogates.kriging.correlation_models.ExponentialCorrelation import ExponentialCorrelation marginals = [Normal(loc=0., scale=4.), Normal(loc=0., scale=4.)] x = MonteCarloSampling(distributions=marginals, nsamples=20, random_state=1) model = PythonModel(model_script='series.py', model_object_name="series") rmodel = RunModel(model=model) - regression_model = LinearRegression() - correlation_model = ExponentialCorrelation() - K = Kriging(regression_model=regression_model, correlation_model=correlation_model, - optimizations_number=10, correlation_model_parameters=[1, 1], - optimizer=MinimizeOptimizer('l-bfgs-b'),) + kernel1 = RBF() + bounds_1 = [[10 ** (-4), 10 ** 3], [10 ** (-3), 10 ** 2], [10 ** (-3), 10 ** 2]] + optimizer1 = MinimizeOptimizer(method='L-BFGS-B', bounds=bounds_1) + gpr = GaussianProcessRegression(kernel=kernel1, hyperparameters=[1, 10 ** (-3), 10 ** (-2)], optimizer=optimizer1, + optimizations_number=50, noise=False, regression_model=LinearRegression(), + random_state=0) # OPTIONS: 'U', 'EFF', 'Weighted-U' learning_function = ExpectedImprovementGlobalFit() - a = AdaptiveKriging(distributions=marginals, runmodel_object=rmodel, surrogate=K, - learning_nsamples=10**3, n_add=1, learning_function=learning_function, + a = AdaptiveKriging(distributions=marginals, runmodel_object=rmodel, surrogate=gpr, + learning_nsamples=10 ** 3, n_add=1, learning_function=learning_function, random_state=2) a.run(nsamples=25, samples=x.samples) - assert a.samples[23, 0] == 11.939859785098493 - assert a.samples[20, 1] == -8.429899469300118 + assert a.samples[23, 0] == -10.24267076486663 + assert a.samples[20, 1] == -11.419510366469687 def test_akmcs_samples_error(): - from UQpy.surrogates.kriging.regression_models.LinearRegression import LinearRegression - from UQpy.surrogates.kriging.correlation_models.ExponentialCorrelation import ExponentialCorrelation marginals = [Normal(loc=0., scale=4.), Normal(loc=0., scale=4.)] x = MonteCarloSampling(distributions=marginals, nsamples=20, random_state=0) model = PythonModel(model_script='series.py', model_object_name="series") rmodel = RunModel(model=model) - regression_model = LinearRegression() - correlation_model = ExponentialCorrelation() - K = Kriging(regression_model=regression_model, correlation_model=correlation_model, - optimizer=MinimizeOptimizer('l-bfgs-b'), - optimizations_number=10, correlation_model_parameters=[1, 1], random_state=1) + kernel1 = RBF() + bounds_1 = [[10 ** (-4), 10 ** 3], [10 ** (-3), 10 ** 2], [10 ** (-3), 10 ** 2]] + optimizer1 = MinimizeOptimizer(method='L-BFGS-B', bounds=bounds_1) + gpr = GaussianProcessRegression(kernel=kernel1, hyperparameters=[1, 10 ** (-3), 10 ** (-2)], optimizer=optimizer1, + optimizations_number=50, noise=False, regression_model=LinearRegression(), + random_state=0) # OPTIONS: 'U', 'EFF', 'Weighted-U' learning_function = WeightedUFunction(weighted_u_stop=2) with pytest.raises(NotImplementedError): - a = AdaptiveKriging(distributions=[Normal(loc=0., scale=4.)]*3, runmodel_object=rmodel, surrogate=K, - learning_nsamples=10**3, n_add=1, learning_function=learning_function, + a = AdaptiveKriging(distributions=[Normal(loc=0., scale=4.)] * 3, runmodel_object=rmodel, surrogate=gpr, + learning_nsamples=10 ** 3, n_add=1, learning_function=learning_function, random_state=2, samples=x.samples) def test_akmcs_u_run_from_init(): - from UQpy.surrogates.kriging.regression_models.LinearRegression import LinearRegression - from UQpy.surrogates.kriging.correlation_models.ExponentialCorrelation import ExponentialCorrelation - marginals = [Normal(loc=0., scale=4.), Normal(loc=0., scale=4.)] x = MonteCarloSampling(distributions=marginals, nsamples=20, random_state=1) model = PythonModel(model_script='series.py', model_object_name="series") rmodel = RunModel(model=model) - regression_model = LinearRegression() - correlation_model = ExponentialCorrelation() - K = Kriging(regression_model=regression_model, correlation_model=correlation_model, - optimizer=MinimizeOptimizer('l-bfgs-b'), - optimizations_number=10, correlation_model_parameters=[1, 1], random_state=0) + kernel1 = RBF() + bounds_1 = [[10 ** (-4), 10 ** 3], [10 ** (-3), 10 ** 2], [10 ** (-3), 10 ** 2]] + optimizer1 = MinimizeOptimizer(method='L-BFGS-B', bounds=bounds_1) + gpr = GaussianProcessRegression(kernel=kernel1, hyperparameters=[1, 10 ** (-3), 10 ** (-2)], optimizer=optimizer1, + optimizations_number=100, noise=False, regression_model=LinearRegression(), + random_state=0) # OPTIONS: 'U', 'EFF', 'Weighted-U' learning_function = UFunction(u_stop=2) - a = AdaptiveKriging(distributions=marginals, runmodel_object=rmodel, surrogate=K, - learning_nsamples=10**3, n_add=1, learning_function=learning_function, + a = AdaptiveKriging(distributions=marginals, runmodel_object=rmodel, surrogate=gpr, + learning_nsamples=10 ** 3, n_add=1, learning_function=learning_function, random_state=2, nsamples=25, samples=x.samples) - assert a.samples[23, 0] == -4.141979058326188 - assert a.samples[20, 1] == -1.6476534435429009 + assert a.samples[23, 0] == -3.781937137406927 + assert a.samples[20, 1] == 0.17610325620498946 diff --git a/tests/unit_tests/sampling/test_refined_stratified.py b/tests/unit_tests/sampling/test_refined_stratified.py index 15a3ccd12..13afea15f 100644 --- a/tests/unit_tests/sampling/test_refined_stratified.py +++ b/tests/unit_tests/sampling/test_refined_stratified.py @@ -1,6 +1,7 @@ import pytest from beartype.roar import BeartypeCallHintPepParamException +from UQpy import GaussianProcessRegression, RBF, LinearRegression from UQpy.run_model.model_execution.PythonModel import PythonModel from UQpy.utilities.MinimizeOptimizer import MinimizeOptimizer from UQpy.sampling.stratified_sampling.refinement.GradientEnhancedRefinement import GradientEnhancedRefinement @@ -9,7 +10,6 @@ from UQpy.sampling.stratified_sampling.refinement.RandomRefinement import * from UQpy.sampling.stratified_sampling.strata.VoronoiStrata import * from UQpy.run_model.RunModel import * -from UQpy.surrogates.kriging.Kriging import Kriging def test_rss_simple_rectangular(): @@ -23,10 +23,10 @@ def test_rss_simple_rectangular(): samples_per_iteration=2, refinement_algorithm=algorithm, random_state=2) - assert y.samples[16, 0] == 0.06614276178462988 - assert y.samples[16, 1] == 0.7836449863362334 - assert y.samples[17, 0] == 0.1891972651582183 - assert y.samples[17, 1] == 0.2961099664117288 + assert y.samples[16, 0] == 0.22677821757428504 + assert y.samples[16, 1] == 0.2729789855337742 + assert y.samples[17, 0] == 0.07501256574570675 + assert y.samples[17, 1] == 0.9321401317029486 def test_rss_simple_voronoi(): @@ -40,10 +40,10 @@ def test_rss_simple_voronoi(): samples_per_iteration=2, refinement_algorithm=algorithm, random_state=2) - assert np.round(y.samples[16, 0], 6) == 0.363793 - assert np.round(y.samples[16, 1], 6) == 0.467625 - assert np.round(y.samples[17, 0], 6) == 0.424586 - assert np.round(y.samples[17, 1], 6) == 0.217301 + assert np.round(y.samples[16, 0], 6) == 0.324738 + assert np.round(y.samples[16, 1], 6) == 0.488029 + assert np.round(y.samples[17, 0], 6) == 0.349367 + assert np.round(y.samples[17, 1], 6) == 0.132426 def test_rect_rss(): @@ -57,10 +57,10 @@ def test_rect_rss(): refinement_algorithm=RandomRefinement(strata=strata)) assert np.allclose(y.samples, np.array([[0.417022, 0.36016225], [1.00011437, 0.15116629], [0.14675589, 0.5461693], [1.18626021, 0.67278036], - [0.77483124, 0.7176612], [1.7101839, 0.66516741]])) + [1.90711287, 0.04595797], [0.80005026, 0.86428026]])) assert np.allclose(np.array(y.samplesU01), np.array([[0.208511, 0.36016225], [0.50005719, 0.15116629], [0.07337795, 0.5461693], [0.59313011, 0.67278036], - [0.38741562, 0.7176612], [0.85509195, 0.66516741]])) + [0.95355644, 0.04595797], [0.40002513, 0.86428026]])) def test_rect_gerss(): @@ -72,31 +72,22 @@ def test_rect_gerss(): x = TrueStratifiedSampling(distributions=marginals, strata_object=strata, nsamples_per_stratum=1) model = PythonModel(model_script='python_model_function.py', model_object_name="y_func") rmodel = RunModel(model=model) - from UQpy.surrogates.kriging.regression_models import LinearRegression - from UQpy.surrogates.kriging.correlation_models import ExponentialCorrelation - K = Kriging(regression_model=LinearRegression(), correlation_model=ExponentialCorrelation(), optimizations_number=20, random_state=0, - correlation_model_parameters=[1, 1], optimizer=MinimizeOptimizer('l-bfgs-b'), ) - K.fit(samples=x.samples, values=rmodel.qoi_list) + kernel1 = RBF() + bounds_1 = [[10 ** (-4), 10 ** 3], [10 ** (-3), 10 ** 2], [10 ** (-3), 10 ** 2]] + optimizer1 = MinimizeOptimizer(method='L-BFGS-B', bounds=bounds_1) + gpr = GaussianProcessRegression(kernel=kernel1, hyperparameters=[1, 10 ** (-3), 10 ** (-2)], optimizer=optimizer1, + optimizations_number=100, noise=False, regression_model=LinearRegression(), + random_state=0) + # gpr.fit(samples=x.samples, values=rmodel.qoi_list) refinement = GradientEnhancedRefinement(strata=x.strata_object, runmodel_object=rmodel, - surrogate=K, nearest_points_number=4) + surrogate=gpr, nearest_points_number=4) z = RefinedStratifiedSampling(stratified_sampling=x, random_state=2, refinement_algorithm=refinement) z.run(nsamples=6) assert np.allclose(z.samples, np.array([[0.417022, 0.36016225], [1.00011437, 0.15116629], [0.14675589, 0.5461693], [1.18626021, 0.67278036], - [1.51296312, 0.77483124], [0.74237455, 0.66026822]])) - # assert np.allclose(z.samples, np.array([[0.417022, 0.36016225], [1.00011437, 0.15116629], - # [0.14675589, 0.5461693], [1.18626021, 0.67278036], - # [1.59254104, 0.96577043], [1.97386531, 0.24237455]])) - # assert np.allclose(z.samples, np.array([[0.417022, 0.36016225], [1.00011437, 0.15116629], - # [0.14675589, 0.5461693], [1.18626021, 0.67278036], - # [1.59254104, 0.96577043], [1.7176612, 0.2101839]])) - # assert np.allclose(z.samplesU01, np.array([[0.208511, 0.36016225], [0.50005719, 0.15116629], - # [0.07337795, 0.5461693], [0.59313011, 0.67278036], - # [0.79627052, 0.96577043], [0.98693265, 0.24237455]])) - # assert np.allclose(z.samplesU01, np.array([[0.208511, 0.36016225], [0.50005719, 0.15116629], - # [0.07337795, 0.5461693], [0.59313011, 0.67278036], - # [0.79627052, 0.96577043], [0.8588306 , 0.2101839]])) + [1.64924557, 0.90711287], + [0.54595797, 0.30005026]])) def test_vor_rss(): @@ -124,26 +115,25 @@ def test_vor_gerss(): marginals = [Uniform(loc=0., scale=2.), Uniform(loc=0., scale=1.)] strata_vor = VoronoiStrata(seeds_number=4, dimension=2, random_state=10) x_vor = TrueStratifiedSampling(distributions=marginals, strata_object=strata_vor, nsamples_per_stratum=1, ) - from UQpy.surrogates.kriging.regression_models.LinearRegression import LinearRegression - from UQpy.surrogates.kriging.correlation_models.ExponentialCorrelation import ExponentialCorrelation model = PythonModel(model_script='python_model_function.py', model_object_name="y_func") rmodel = RunModel(model=model) - K_ = Kriging(regression_model=LinearRegression(), correlation_model=ExponentialCorrelation(), optimizations_number=20, - optimizer=MinimizeOptimizer('l-bfgs-b'), random_state=0, - correlation_model_parameters=[1, 1]) - - K_.fit(samples=x_vor.samples, values=rmodel.qoi_list) - z_vor = RefinedStratifiedSampling(stratified_sampling=x_vor, nsamples=6, random_state=x_vor.random_state, + kernel1 = RBF() + bounds_1 = [[10 ** (-4), 10 ** 3], [10 ** (-3), 10 ** 2], [10 ** (-3), 10 ** 2]] + optimizer1 = MinimizeOptimizer(method='L-BFGS-B', bounds=bounds_1) + gpr = GaussianProcessRegression(kernel=kernel1, hyperparameters=[1, 10 ** (-3), 10 ** (-2)], optimizer=optimizer1, + optimizations_number=100, noise=False, regression_model=LinearRegression(), + random_state=0) + z_vor = RefinedStratifiedSampling(stratified_sampling=x_vor, nsamples=6, random_state=0, refinement_algorithm=GradientEnhancedRefinement(strata=x_vor.strata_object, runmodel_object=rmodel, - surrogate=K_, + surrogate=gpr, nearest_points_number=4)) assert np.allclose(z_vor.samples, np.array([[1.78345908, 0.01640854], [1.46201137, 0.70862104], [0.4021338, 0.05290083], [0.1062376, 0.88958226], - [0.61246269, 0.47160095], [1.16609055, 0.30832536]])) + [0.66730342, 0.46988084], [1.5015904 , 0.97050966]])) assert np.allclose(z_vor.samplesU01, np.array([[0.89172954, 0.01640854], [0.73100569, 0.70862104], [0.2010669, 0.05290083], [0.0531188, 0.88958226], - [0.30623134, 0.47160095], [0.58304527, 0.30832536]])) + [0.33365171, 0.46988084], [0.7507952 , 0.97050966]])) def test_rss_random_state(): @@ -165,17 +155,17 @@ def test_rss_runmodel_object(): marginals = [Uniform(loc=0., scale=2.), Uniform(loc=0., scale=1.)] strata = RectangularStrata(strata_number=[2, 2]) x = TrueStratifiedSampling(distributions=marginals, strata_object=strata, nsamples_per_stratum=1, random_state=1) - from UQpy.surrogates.kriging.regression_models import LinearRegression - from UQpy.surrogates.kriging.correlation_models import ExponentialCorrelation - - K = Kriging(regression_model=LinearRegression(), correlation_model=ExponentialCorrelation(), optimizations_number=20, - correlation_model_parameters=[1, 1], optimizer=MinimizeOptimizer('l-bfgs-b'), ) + kernel1 = RBF() + bounds_1 = [[10 ** (-4), 10 ** 3], [10 ** (-3), 10 ** 2], [10 ** (-3), 10 ** 2]] + optimizer1 = MinimizeOptimizer(method='L-BFGS-B', bounds=bounds_1) + gpr = GaussianProcessRegression(kernel=kernel1, hyperparameters=[1, 10 ** (-3), 10 ** (-2)], optimizer=optimizer1, + optimizations_number=100, noise=False, regression_model=LinearRegression(), + random_state=0) model = PythonModel(model_script='python_model_function.py', model_object_name="y_func") rmodel = RunModel(model=model) - K.fit(samples=x.samples, values=rmodel.qoi_list) with pytest.raises(BeartypeCallHintPepParamException): refinement = GradientEnhancedRefinement(strata=x.strata_object, runmodel_object='abc', - surrogate=K) + surrogate=gpr) RefinedStratifiedSampling(stratified_sampling=x, samples_number=6, samples_per_iteration=2, refinement_algorithm=refinement) diff --git a/tests/unit_tests/surrogates/test_kriging.py b/tests/unit_tests/surrogates/test_kriging.py deleted file mode 100644 index d25084b02..000000000 --- a/tests/unit_tests/surrogates/test_kriging.py +++ /dev/null @@ -1,198 +0,0 @@ -import pytest -from UQpy.utilities.MinimizeOptimizer import MinimizeOptimizer -from beartype.roar import BeartypeCallHintPepParamException - -from UQpy.surrogates.kriging.Kriging import Kriging -import numpy as np -from UQpy.surrogates.kriging.regression_models import LinearRegression, ConstantRegression -from UQpy.surrogates.kriging.correlation_models import GaussianCorrelation - - -samples = np.linspace(0, 5, 20).reshape(-1, 1) -values = np.cos(samples) -optimizer = MinimizeOptimizer(method="L-BFGS-B") -krig = Kriging(regression_model=LinearRegression(), correlation_model=GaussianCorrelation(), optimizer=optimizer, - correlation_model_parameters=[0.14], optimize=False, random_state=1) -krig.fit(samples=samples, values=values, correlation_model_parameters=[0.3]) - -optimizer = MinimizeOptimizer(method="L-BFGS-B") -krig2 = Kriging(regression_model=ConstantRegression(), correlation_model=GaussianCorrelation(), optimizer=optimizer, - correlation_model_parameters=[0.3], bounds=[[0.01, 5]], - optimize=False, optimizations_number=100, normalize=False, - random_state=2) -krig2.fit(samples=samples, values=values) - - -# Using the in-built linear regression model as a function -linear_regression_model = Kriging(regression_model=LinearRegression(), correlation_model=GaussianCorrelation(), optimizer=optimizer, - correlation_model_parameters=[1]).regression_model -optimizer = MinimizeOptimizer(method="L-BFGS-B") -gaussian_corrleation_model = Kriging(regression_model=LinearRegression(), correlation_model=GaussianCorrelation(), optimizer=optimizer, - correlation_model_parameters=[1]).correlation_model - -optimizer = MinimizeOptimizer(method="L-BFGS-B") -krig3 = Kriging(regression_model=linear_regression_model, correlation_model=gaussian_corrleation_model, - optimizer=optimizer, - correlation_model_parameters=[1], optimize=False, normalize=False, random_state=0) -krig3.fit(samples=samples, values=values) - - -def test_predict(): - prediction = np.round(krig.predict([[1], [np.pi/2], [np.pi]], True), 3) - expected_prediction = np.array([[0.54, 0., -1.], [0., 0., 0.]]) - assert (expected_prediction == prediction).all() - - -def test_predict1(): - prediction = np.round(krig2.predict([[1], [2*np.pi], [np.pi]], True), 3) - expected_prediction = np.array([[0.54, 1.009, -1.], [0., 0.031, 0.]]) - assert (expected_prediction == prediction).all() - - -def test_predict2(): - prediction = np.round(krig3.predict([[1], [np.pi/2], [np.pi]]), 3) - expected_prediction = np.array([[0.54, -0., -1.]]) - assert (expected_prediction == prediction).all() - - -def test_jacobian(): - jacobian = np.round(krig.jacobian([[np.pi], [np.pi/2]]), 3) - expected_jacobian = np.array([-0., -1.]) - assert (expected_jacobian == jacobian).all() - - -def test_jacobian1(): - jacobian = np.round(krig3.jacobian([[np.pi], [np.pi/2]]), 3) - expected_jacobian = np.array([0., -1.]) - assert (expected_jacobian == jacobian).all() - - -def test_regression_models(): - from UQpy.surrogates.kriging.regression_models import ConstantRegression, LinearRegression, QuadraticRegression - krig.regression_model = ConstantRegression() - tmp = krig.regression_model.r([[0], [1]]) - tmp_test1 = (tmp[0] == np.array([[1.], [1.]])).all() and (tmp[1] == np.array([[[0.]], [[0.]]])).all() - - krig.regression_model = LinearRegression() - tmp = krig.regression_model.r([[0], [1]]) - tmp_test2 = (tmp[0] == (np.array([[1., 0.], [1., 1.]]))).all() and \ - (tmp[1] == np.array([[[0., 1.]], [[0., 1.]]])).all() - - krig.regression_model = QuadraticRegression() - tmp = krig.regression_model.r([[-1, 1], [2, -0.5]]) - tmp_test3 = (tmp[0] == np.array([[1., -1., 1., 1., -1., 1.], [1., 2., -0.5, 4., -1., 0.25]])).all() and \ - (tmp[1] == np.array([[[0., 1., 0., -2., 1., 0.], - [0., 0., 1., 0., -1., 2.]], - [[0., 1., 0., 4., -0.5, 0.], - [0., 0., 1., 0., 2., -1.]]])).all() - - assert tmp_test1 and tmp_test2 and tmp_test3 - - -def test_correlation_models(): - from UQpy.surrogates.kriging.correlation_models import ExponentialCorrelation, LinearCorrelation, SphericalCorrelation, CubicCorrelation, SplineCorrelation - krig.correlation_model = ExponentialCorrelation() - rx_exponential = (np.round(krig.correlation_model.c([[0], [1], [2]], [[2]], np.array([1])), 3) == - np.array([[0.135], [0.368], [1.]])).all() - drdt_exponential = (np.round(krig.correlation_model.c([[0], [1], [2]], [[2]], np.array([1]), dt=True)[1], 3) == - np.array([[[-0.271]], [[-0.368]], [[0.]]])).all() - drdx_exponential = (np.round(krig.correlation_model.c([[0], [1], [2]], [[2]], np.array([1]), dx=True)[1], 3) == - np.array([[[0.135]], [[0.368]], [[0.]]])).all() - expon = rx_exponential and drdt_exponential and drdx_exponential - - krig.correlation_model = LinearCorrelation() - rx_linear = (np.round(krig.correlation_model.c([[0], [1], [2]], [[2]], np.array([1])), 3) == - np.array([[0.], [0.], [1.]])).all() - drdt_linear = (np.round(krig.correlation_model.c([[0.4], [0.5], [0.6]], [[0.5]], np.array([1]), dt=True)[1], 3) == - np.array([[[-0.1]], [[-0.]], [[-0.1]]])).all() - drdx_linear = (np.round(krig.correlation_model.c([[0.4], [0.5], [0.6]], [[0.5]], np.array([1]), dx=True)[1], 3) == - np.array([[[1.]], [[-0.]], [[-1.]]])).all() - linear = rx_linear and drdt_linear and drdx_linear - - krig.correlation_model = SphericalCorrelation() - rx_spherical = (np.round(krig.correlation_model.c([[0], [1], [2]], [[2]], np.array([1])), 3) == - np.array([[0.], [0.], [1.]])).all() - drdt_spherical = (np.round(krig.correlation_model.c([[0.4], [0.5], [0.6]], [[0.5]], np.array([1]), dt=True)[1], 3) - == np.array([[[-0.148]], [[-0.]], [[-0.148]]])).all() - drdx_spherical = (np.round(krig.correlation_model.c([[0.4], [0.5], [0.6]], [[0.5]], np.array([1]), dx=True)[1], 3) - == np.array([[[1.485]], [[-0.]], [[-1.485]]])).all() - spherical = rx_spherical and drdt_spherical and drdx_spherical - - krig.correlation_model = CubicCorrelation() - rx_cubic = (np.round(krig.correlation_model.c([[0.2], [0.5], [1]], [[0.5]], np.array([1])), 3) == - np.array([[0.784], [1.], [0.5]])).all() - drdt_cubic = (np.round(krig.correlation_model.c([[0.4], [0.5], [0.6]], [[0.5]], np.array([1]), dt=True)[1], 3) == - np.array([[[-0.054]], [[0.]], [[-0.054]]])).all() - drdx_cubic = (np.round(krig.correlation_model.c([[0.4], [0.5], [0.6]], [[0.5]], np.array([1]), dx=True)[1], 3) == - np.array([[[0.54]], [[0.]], [[-0.54]]])).all() - cubic = rx_cubic and drdt_cubic and drdx_cubic - - krig.correlation_model = SplineCorrelation() - rx_spline = (np.round(krig.correlation_model.c([[0.2], [0.5], [1]], [[0.5]], np.array([1])), 3) == - np.array([[0.429], [1.], [0.156]])).all() - drdt_spline = (np.round(krig.correlation_model.c([[0.4], [0.5], [0.6]], [[0.5]], np.array([1]), dt=True)[1], 3) == - np.array([[[-0.21]], [[0.]], [[-0.21]]])).all() - drdx_spline = (np.round(krig.correlation_model.c([[0.4], [0.5], [0.6]], [[0.5]], np.array([1]), dx=True)[1], 3) == - np.array([[[2.1]], [[0.]], [[-2.1]]])).all() - spline = rx_spline and drdt_spline and drdx_spline - - assert expon and linear and spherical and cubic and spline - - -def test_wrong_regression_model(): - """ - Raises an error if reg_model is not callable or a string of an in-built model. - """ - with pytest.raises(BeartypeCallHintPepParamException): - Kriging(regression_model='A', correlation_model=GaussianCorrelation(), correlation_model_parameters=[1]) - - -def test_wrong_correlation_model(): - """ - Raises an error if corr_model is not callable or a string of an in-built model. - """ - with pytest.raises(BeartypeCallHintPepParamException): - Kriging(regression_model=LinearRegression(), correlation_model='A', correlation_model_parameters=[1]) - - -def test_missing_correlation_model_parameters(): - """ - Raises an error if corr_model_params is not defined. - """ - with pytest.raises(TypeError): - Kriging(regression_model=LinearRegression(), correlation_model=GaussianCorrelation(), bounds=[[0.01, 5]], - optimizations_number=100, random_state=1) - - -def test_optimizer(): - """ - Raises an error if corr_model_params is not defined. - """ - with pytest.raises(ValueError): - Kriging(regression_model=LinearRegression(), correlation_model=GaussianCorrelation(), - correlation_model_parameters=[1], optimizer='A') - - -def test_random_state(): - """ - Raises an error if type of random_state is not correct. - """ - with pytest.raises(BeartypeCallHintPepParamException): - Kriging(regression_model=LinearRegression(), correlation_model=GaussianCorrelation(), - correlation_model_parameters=[1], random_state='A') - - -def test_log_likelihood(): - prediction = np.round(krig3.log_likelihood(np.array([1.5]), - krig3.correlation_model, np.array([[1], [2]]), - np.array([[1], [1]]), np.array([[1], [2]]), return_grad=False), 3) - expected_prediction = 1.679 - assert (expected_prediction == prediction).all() - - -def test_log_likelihood_derivative(): - prediction = np.round(krig3.log_likelihood(np.array([1.5]), krig3.correlation_model, np.array([[1], [2]]), - np.array([[1], [1]]), np.array([[1], [2]]), return_grad=True)[1], 3) - expected_prediction = np.array([-0.235]) - assert (expected_prediction == prediction).all() -