From edc92c3c92d0487637294756299779522a102f76 Mon Sep 17 00:00:00 2001 From: Klawonn Date: Wed, 27 Apr 2022 10:33:26 +0200 Subject: [PATCH 01/15] added first basic structure for surrogate models --- probeye/definition/__init__.py | 1 + probeye/definition/inverse_problem.py | 67 +++++++ probeye/definition/surrogate_model.py | 32 ++++ .../test_surrogate_model.py | 176 ++++++++++++++++++ 4 files changed, 276 insertions(+) create mode 100644 probeye/definition/surrogate_model.py create mode 100644 tests/integration_tests/tests_without_correlation/test_surrogate_model.py diff --git a/probeye/definition/__init__.py b/probeye/definition/__init__.py index ac924433..4a45c668 100644 --- a/probeye/definition/__init__.py +++ b/probeye/definition/__init__.py @@ -2,6 +2,7 @@ from probeye.definition import inverse_problem from probeye.definition import inference_problem from probeye.definition import forward_model +from probeye.definition import surrogate_model from probeye.definition import likelihood_model from probeye.definition import parameter from probeye.definition import prior diff --git a/probeye/definition/inverse_problem.py b/probeye/definition/inverse_problem.py index d0ffa5a1..8137f50d 100644 --- a/probeye/definition/inverse_problem.py +++ b/probeye/definition/inverse_problem.py @@ -10,6 +10,7 @@ # local imports from probeye.definition.parameter import Parameters from probeye.definition.forward_model import ForwardModelBase +from probeye.definition.surrogate_model import SurrogateModelBase from probeye.definition.likelihood_model import GaussianLikelihoodModel from probeye.subroutines import underlined_string, titled_table from probeye.subroutines import simplified_list_string @@ -96,6 +97,13 @@ def __init__( # not be modified directly self._forward_models = {} # type: dict + # the following dict contains the problem's surrogate models (note that + # surrogate models are optional, not every problem has surrogate models); + # the keys in this dict are the surrogate model names, while the values are the + # surrogate model objects (check out the script surrogate_model.py); this + # dictionary is managed internally and should not be modified directly + self._surrogate_models = {} # type: dict + # this dictionary contains the problem's likelihood models; as the other private # attributes above, it is managed internally and should not be modified directly self._likelihood_models = {} # type: dict @@ -211,6 +219,11 @@ def forward_models(self) -> dict: """Access self._forward_models from outside via self.forward_models.""" return self._forward_models + @property + def surrogate_models(self) -> dict: + """Access self._surrogate_models from outside via self.forward_models.""" + return self._surrogate_models + @property def experiments(self) -> dict: """Access self._experiments from outside via self.experiments.""" @@ -640,6 +653,60 @@ def add_forward_model(self, forward_model: ForwardModelBase): # the given forward model name self._forward_models[forward_model.name] = forward_model + def add_surrogate_model(self, surrogate_model: SurrogateModelBase): + """ + Adds a surrogate model to the inverse problem. Note that multiple forward models + can be added to one problem. + + Parameters + ---------- + surrogate_model + Defines the surrogate model. Check out surrogate_model.py to see a template + for the surrogate model definition. The user will then have to derive her + own surrogate model from that base class. Examples can be found in the + package directory tests/integration_tests. + """ + + # log at beginning so that errors can be quickly associated + logger.debug(f"Adding surrogate model '{surrogate_model.name}'") + + # check if all given model parameters have already been added to the inverse + # problem; note that the surrogate model can only be added to the problem after + # the corresponding parameters and the forward model have been defined + for prm_name in surrogate_model.prms_def: + self._parameters.confirm_that_parameter_exists(prm_name) + # check, if the type was set correctly; if it was not set yet, it will be + # set automatically here + if self._parameters[prm_name].type == "not defined": + self.change_parameter_type(prm_name, "model") + elif self._parameters[prm_name].type in ["prior", "likelihood"]: + raise ValueError( + f"The parameter '{prm_name}' defined in forward model " + f"'{surrogate_model.name}' was assigned the type " + f"'{self._parameters[prm_name].type}' (it should be 'model')." + ) + + # check if the given name for the forward model has already been used + if surrogate_model.name in self._surrogate_models: + raise RuntimeError( + f"The name '{surrogate_model.name}' of the surrogate model you are " + f"trying to add to the problem has already been used for another" + f"surrogate model within the problem scope. Please choose another name." + ) + + # check if the surrogate model's forward model has already been added to the + # inverse problem + if surrogate_model.forward_model.name not in self._forward_models: + raise RuntimeError( + f"The surrogate model's forward model " + f"'{surrogate_model.forward_model.name}' has not been added to the " + f"problem yet!" + ) + + # add the given surrogate model to the internal surrogate model dictionary under + # the given surrogate model name + self._surrogate_models[surrogate_model.name] = surrogate_model + # =============================================================== # # Experiments related methods # # =============================================================== # diff --git a/probeye/definition/surrogate_model.py b/probeye/definition/surrogate_model.py new file mode 100644 index 00000000..3717856b --- /dev/null +++ b/probeye/definition/surrogate_model.py @@ -0,0 +1,32 @@ +# local imports +from probeye.definition.forward_model import ForwardModelBase + + +class SurrogateModelBase(ForwardModelBase): + """ + Base class for a surrogate model, i.e., a forward model that approximates another + (typically computationally more expensive) forward model. + + Parameters + ---------- + name + The name of the surrogate model. Must be unique among all surrogate model's + names within a considered InverseProblem. + forward_model + The forward model object that the surrogate model approximates. + """ + + def __init__(self, name: str, forward_model: ForwardModelBase): + + super().__init__(name) + + # the surrogate model has access to the forward model it approximates; this + # forward model can be called (evaluating its response) from the surrogate model + # via self.forward_model.response(inp) where inp contains the input dictionary + self.forward_model = forward_model + + def fit(self): + """ + Prepares the surrogate model by approximating the forward model in some way. + """ + pass diff --git a/tests/integration_tests/tests_without_correlation/test_surrogate_model.py b/tests/integration_tests/tests_without_correlation/test_surrogate_model.py new file mode 100644 index 00000000..19e9787e --- /dev/null +++ b/tests/integration_tests/tests_without_correlation/test_surrogate_model.py @@ -0,0 +1,176 @@ +# standard library imports +import unittest +import time + +# third party imports +import numpy as np +import matplotlib.pyplot as plt + +# local imports (problem definition) +from probeye.definition.inverse_problem import InverseProblem +from probeye.definition.forward_model import ForwardModelBase +from probeye.definition.surrogate_model import SurrogateModelBase +from probeye.definition.sensor import Sensor +from probeye.definition.likelihood_model import GaussianLikelihoodModel + + +class TestProblem(unittest.TestCase): + def test_surrogate_model(self, plot: bool = False): + """ + Demonstrates how to use a surrogate model. Still under discussion. + + Parameters + ---------- + plot + Triggers a plot of the data used for calibration. + """ + + # ============================================================================ # + # Set numeric values # + # ============================================================================ # + + # 'true' value of a, and its normal prior parameters + m_true = 2.5 + mean_m = 2.0 + std_m = 1.0 + + # 'true' value of b, and its normal prior parameters + b_true = 1.7 + mean_b = 1.0 + std_b = 1.0 + + # 'true' value of additive error sd, and its uniform prior parameters + sigma = 0.5 + low_sigma = 0.0 + high_sigma = 0.8 + + # the number of generated experiment_names and seed for random numbers + n_tests = 50 + seed = 1 + + # ============================================================================ # + # Define the Forward Model # + # ============================================================================ # + + class ExpensiveModel(ForwardModelBase): + def interface(self): + self.parameters = ["m", "b"] + self.input_sensors = Sensor("x") + self.output_sensors = Sensor("y", std_model="sigma") + + def response(self, inp: dict) -> dict: + x = inp["x"] + m = inp["m"] + b = inp["b"] + time.sleep(1.0) + return {"y": m * x + b} + + # ============================================================================ # + # Define the Surrogate Model # + # ============================================================================ # + + class SurrogateModel(ExpensiveModel, SurrogateModelBase): + """ + The inheritance from ExpensiveModel 'copies' the interface-method from + ExpensiveModel (the surrogate model should have the same interface as the + forward model). The inheritance from SurrogateModelBase is required to + assign a forward model to the surrogate model, see surrogate_model.py. + """ + + def response(self, inp: dict) -> dict: + x = inp["x"] + m = inp["m"] + b = inp["b"] + return {"y": m * x + b} + + # ============================================================================ # + # Define the Inference Problem # + # ============================================================================ # + + # initialize the inverse problem with a useful name + problem = InverseProblem("Using a surrogate model") + + # add all parameters to the problem + problem.add_parameter( + "m", + "model", + tex="$m$", + info="Slope of the graph", + prior=("normal", {"mean": mean_m, "std": std_m}), + ) + problem.add_parameter( + "b", + "model", + info="Intersection of graph with y-axis", + tex="$b$", + prior=("normal", {"mean": mean_b, "std": std_b}), + ) + problem.add_parameter( + "sigma", + "likelihood", + domain="(0, +oo)", + tex=r"$\sigma$", + info="Standard deviation, of zero-mean additive model error", + prior=("uniform", {"low": low_sigma, "high": high_sigma}), + ) + + # add the forward model to the problem + forward_model = ExpensiveModel("ExpensiveModel") + problem.add_forward_model(forward_model) + + # add the forward model to the problem + surrogate_model = SurrogateModel("FastModel", forward_model=forward_model) + problem.add_surrogate_model(surrogate_model) + + # ============================================================================ # + # Add test data to the Inference Problem # + # ============================================================================ # + + # data-generation; normal likelihood with constant variance around each point + np.random.seed(seed) + x_test = np.linspace(0.0, 1.0, n_tests) + y_true = surrogate_model.response( + {surrogate_model.input_sensor.name: x_test, "m": m_true, "b": b_true} + )[surrogate_model.output_sensor.name] + y_test = np.random.normal(loc=y_true, scale=sigma) + + # add the experimental data + problem.add_experiment( + f"TestSeries_1", + fwd_model_name="ExpensiveModel", + sensor_values={ + forward_model.input_sensor.name: x_test, + forward_model.output_sensor.name: y_test, + }, + ) + + # plot the true and noisy data + if plot: + plt.scatter(x_test, y_test, label="measured data", s=10, c="red", zorder=10) + plt.plot(x_test, y_true, label="true", c="black") + plt.xlabel(forward_model.input_sensor.name) + plt.ylabel(forward_model.output_sensor.name) + plt.legend() + plt.tight_layout() + plt.draw() # does not stop execution + + # ============================================================================ # + # Add likelihood model(s) # + # ============================================================================ # + + # add the likelihood model to the problem + problem.add_likelihood_model( + GaussianLikelihoodModel( + prms_def="sigma", + experiment_name="TestSeries_1", + model_error="additive", + name="SimpleLikelihoodModel", + ) + ) + + # give problem overview + problem.info() + + +if __name__ == "__main__": + unittest.main() From 78b0dc88733a86f7bc0ebd388aa54f30e6007509 Mon Sep 17 00:00:00 2001 From: Klawonn Date: Tue, 17 May 2022 15:11:42 +0200 Subject: [PATCH 02/15] added latin hypercube sampler class --- probeye/definition/inverse_problem.py | 17 ++ probeye/inference/priors.py | 41 ++-- probeye/surrogate/__init__.py | 2 + probeye/surrogate/initial_sampling.py | 175 ++++++++++++++++++ setup.cfg | 1 + .../surrogate/test_initial_sampling.py | 104 +++++++++++ 6 files changed, 325 insertions(+), 15 deletions(-) create mode 100644 probeye/surrogate/__init__.py create mode 100644 probeye/surrogate/initial_sampling.py create mode 100644 tests/unit_tests/surrogate/test_initial_sampling.py diff --git a/probeye/definition/inverse_problem.py b/probeye/definition/inverse_problem.py index 8137f50d..d652db75 100644 --- a/probeye/definition/inverse_problem.py +++ b/probeye/definition/inverse_problem.py @@ -439,6 +439,15 @@ def change_constant(self, prm_name: str, new_value: Union[int, float]): value=new_value ) + def get_latent_prior_hyperparameters(self, prm_name): + latent_hyperparameters = [] + if self.parameters[prm_name].is_latent: + hyperparameters = self.parameters[prm_name].prior.hyperparameters + for hyperparameter in hyperparameters: + if self.parameters[hyperparameter].is_latent: + latent_hyperparameters.append(hyperparameter) + return latent_hyperparameters + def get_parameters(self, theta: np.ndarray, prm_def: dict) -> dict: """ Extracts the numeric values for given parameters that have been defined within @@ -482,6 +491,14 @@ def get_parameters(self, theta: np.ndarray, prm_def: dict) -> dict: prms[local_name] = theta[idx:idx_end] return prms + def get_constants(self, prm_def: dict) -> dict: + prms = {} + for global_name, local_name in prm_def.items(): + idx = self._parameters[global_name].index + if idx is None: + prms[local_name] = self._parameters[global_name].value + return prms + def check_parameter_domains(self, theta: np.ndarray) -> bool: """ Checks whether the given values of the latent parameters are within their diff --git a/probeye/inference/priors.py b/probeye/inference/priors.py index 3400c399..2d529b3b 100644 --- a/probeye/inference/priors.py +++ b/probeye/inference/priors.py @@ -7,7 +7,6 @@ # local imports from probeye.definition.prior import PriorBase -from probeye.subroutines import len_or_one class PriorNormal(PriorBase): @@ -34,7 +33,7 @@ def __init__( super().__init__(ref_prm, prms_def, name, "normal distribution") def __call__( - self, prms: dict, method: str, use_ref_prm: bool = True, **kwargs + self, prms: dict, method: str, *args, use_ref_prm: bool = True, **kwargs ) -> Union[float, np.ndarray]: """ Evaluates stats.norm.(x, loc, scale) or, if use_ref_prm=False @@ -47,6 +46,8 @@ def __call__( Contains the prior's parameters as keys and their values as values. method The method of stats.norm to be evaluated (e.g. 'pdf' or 'logpdf'). + args + Additional positional arguments. use_ref_prm If True stats.norm.(x, loc, scale) is evaluated, hence 'x' must be provided in the prms dictionary. Otherwise, the evaluated method is @@ -64,9 +65,9 @@ def __call__( fun = getattr(stats.norm, method) if use_ref_prm: x = prms[self.ref_prm] - return fun(x, loc=mean, scale=std, **kwargs) + return fun(x, *args, loc=mean, scale=std, **kwargs) else: - return fun(loc=mean, scale=std, **kwargs) + return fun(*args, loc=mean, scale=std, **kwargs) def generate_samples( self, prms: dict, size: int, seed: Optional[int] = None @@ -118,7 +119,7 @@ def __init__( super().__init__(ref_prm, prms_def, name, "multivariate normal distribution") def __call__( - self, prms: dict, method: str, use_ref_prm: bool = True, **kwargs + self, prms: dict, method: str, *args, use_ref_prm: bool = True, **kwargs ) -> Union[float, np.ndarray]: """ Evaluates stats.multivariate_norm.(x, loc, scale) or, if use_ref_prm= @@ -158,9 +159,9 @@ def __call__( ) if use_ref_prm: x = prms[self.ref_prm] - return fun(x, mean=mean, cov=cov, **kwargs) + return fun(x, *args, mean=mean, cov=cov, **kwargs) else: - return fun(mean=mean, cov=cov, **kwargs) + return fun(*args, mean=mean, cov=cov, **kwargs) def generate_samples( self, prms: dict, size: int, seed: Optional[int] = None @@ -217,6 +218,7 @@ def __call__( self, prms: dict, method: str, + *args, use_ref_prm: bool = True, **kwargs, ) -> float: @@ -231,6 +233,8 @@ def __call__( Contains the prior's parameters as keys and their values as values. method The method of stats.lognorm to be evaluated (e.g. 'pdf', 'logpdf'). + args + Additional positional arguments. use_ref_prm If True stats.lognorm.(x, loc, scale) is evaluated, hence 'x' must be provided in the prms dictionary. Otherwise, the evaluated method is @@ -252,9 +256,9 @@ def __call__( shape = sigma if use_ref_prm: x = prms[self.ref_prm] - return fun(x, shape, scale=scale, **kwargs) + return fun(x, shape, *args, scale=scale, **kwargs) else: - return fun(shape, scale=scale, **kwargs) + return fun(shape, *args, scale=scale, **kwargs) def generate_samples( self, @@ -317,6 +321,7 @@ def __call__( self, prms: dict, method: str, + *args, use_ref_prm: bool = True, **kwargs, ) -> float: @@ -331,6 +336,8 @@ def __call__( Contains the prior's parameters as keys and their values as values. method The method of stats.truncnorm to be evaluated (e.g. 'pdf', 'logpdf'). + args + Additional positional arguments. use_ref_prm If True stats.truncnorm.(x, loc, scale) is evaluated, hence 'x' must be provided in the prms dictionary. Otherwise, the evaluated method is @@ -411,7 +418,7 @@ def __init__( super().__init__(ref_prm, prms_def, name, "uniform distribution") def __call__( - self, prms: dict, method: str, use_ref_prm: bool = True, **kwargs + self, prms: dict, method: str, *args, use_ref_prm: bool = True, **kwargs ) -> float: """ Evaluates stats.uniform.(x, loc, scale) or, if use_ref_prm=False stats. @@ -424,6 +431,8 @@ def __call__( Contains the prior's parameters as keys and their values as values. method The method of stats.uniform to be evaluated (e.g. 'pdf', 'logpdf'). + args + Additional positional arguments. use_ref_prm If True stats.uniform.(x, loc, scale) is evaluated, hence 'x' must be provided in the prms dictionary. Otherwise, the evaluated method is @@ -441,9 +450,9 @@ def __call__( high = prms[f"high_{self.ref_prm}"] if use_ref_prm: x = prms[self.ref_prm] - return fun(x, loc=low, scale=high - low, **kwargs) + return fun(x, *args, loc=low, scale=high - low, **kwargs) else: - return fun(loc=low, scale=high - low, **kwargs) + return fun(*args, loc=low, scale=high - low, **kwargs) def generate_samples( self, prms: dict, size: int, seed: Optional[int] = None @@ -497,7 +506,7 @@ def __init__( super().__init__(ref_prm, prms_def, name, "Weibull distribution") def __call__( - self, prms: dict, method: str, use_ref_prm: bool = True, **kwargs + self, prms: dict, method: str, *args, use_ref_prm: bool = True, **kwargs ) -> float: """ Evaluates stats.weibull_min.(x, loc, scale) or, if use_ref_prm=False, @@ -510,6 +519,8 @@ def __call__( Contains the prior's parameters as keys and their values as values. method The method of stats.weibull_min to be evaluated ('pdf', 'logpdf', etc.). + args + Additional positional arguments. use_ref_prm If True stats.weibull_min.(x, loc, scale) is evaluated, hence 'x' must be provided in the prms dictionary. Otherwise, the evaluated method is @@ -528,9 +539,9 @@ def __call__( scale = prms[f"scale_{self.ref_prm}"] if use_ref_prm: x = prms[self.ref_prm] - return fun(x, shape, loc=loc, scale=scale, **kwargs) + return fun(x, shape, *args, loc=loc, scale=scale, **kwargs) else: - return fun(shape, loc=loc, scale=scale, **kwargs) + return fun(shape, *args, loc=loc, scale=scale, **kwargs) def generate_samples( self, prms: dict, size: int, seed: Optional[int] = None diff --git a/probeye/surrogate/__init__.py b/probeye/surrogate/__init__.py new file mode 100644 index 00000000..0cec2ca2 --- /dev/null +++ b/probeye/surrogate/__init__.py @@ -0,0 +1,2 @@ +# subpackage imports +from probeye.surrogate import initial_sampling diff --git a/probeye/surrogate/initial_sampling.py b/probeye/surrogate/initial_sampling.py new file mode 100644 index 00000000..a9d02ab1 --- /dev/null +++ b/probeye/surrogate/initial_sampling.py @@ -0,0 +1,175 @@ +# standard library +from typing import Tuple +import copy + +# third party imports +import numpy as np +import pandas as pd +from scipy.stats import qmc + +# local imports +from probeye.definition.inverse_problem import InverseProblem +from probeye.inference.priors import translate_prior +from probeye.subroutines import len_or_one + + +class LatinHypercubeSampler: + """ + Contains functionalities to provide samples of an inverse problem's parameters by + utilizing latin hypercube sampling. The generate samples are intended for training + of a surrogate model. + + Parameters + ---------- + problem + The considered inverse problem. + """ + + def __init__(self, problem: InverseProblem): + + # the considered inverse problem + self.problem = problem + + # the sampling happens before the problem is given to a solver; this means that + # the priors of the problem are merely descriptive and they have to be + # translated to have their respective computing routines + self.priors = copy.deepcopy(self.problem.priors) + for prior_template in self.problem.priors.values(): + prm_name = prior_template.ref_prm + self.priors[prm_name] = translate_prior(prior_template) + + def generate_samples(self, n_samples: int, seed: int = 1) -> np.ndarray: + """ + Samples the problems latent parameter vector from the parameter's priors in + combination with latin hypercube sampling. + + Parameters + ---------- + n_samples + The number of requested samples. + seed + Seed for random number generator. + + Returns + ------- + sample_array + The sampled latent parameter vectors. Each row corresponds to a single + latent parameter vector. + """ + + # make sure that all parameters are one-dimensional; it is not straight forward + # how to do LHS for general multivariate parameters + for prm_name in self.problem.latent_prms: + if self.problem.parameters[prm_name].dim > 1: + raise RuntimeError( + f"The given problem has a multivariate parameter ('{prm_name}') " + f"for which LHS is not supported." + ) + + # the dimension of the required hypercube is the number of parameters; note + # that this only holds since all parameters are 1D (see the check above) + dim = self.problem.n_latent_prms + + # create the requested number of latin hypercube samples with the requested dim. + lhd = qmc.LatinHypercube(d=dim, seed=seed).random(n=n_samples) + + # this is going to be the array of the parameter samples; each row will + # correspond to a theta-vector + sample_array = np.zeros(lhd.shape) + + # fill the prepared array for the parameter samples using the LHS samples + for prm_name, parameter in self.problem.parameters.items(): + if parameter.is_latent: + idx = parameter.index + prior = self.priors[prm_name] + prms = self.problem.get_constants(prior.hyperparameters) + for lt_prm in self.problem.get_latent_prior_hyperparameters(prm_name): + idx_lt = self.problem.parameters[lt_prm].index + prms[lt_prm] = sample_array[:, idx_lt] + q = lhd[:, idx] + sample_array[:, idx] = prior(prms, "ppf", q, use_ref_prm=False) + + return sample_array + + def generate_training_data( + self, + forward_model_name: str, + n_samples: int, + seed: int = 1, + ) -> Tuple[pd.DataFrame, dict]: + """ + Generates a given number of training data for fitting a surrogate model. The + training data contains a number of parameter vectors (sampled using LHS) and + the corresponding model responses. + + Parameters + ---------- + forward_model_name + The name of the forward model that should be evaluated. + n_samples + The number of parameter vectors the forward model should be evaluated for. + seed + Seed for random number generator. + + Returns + ------- + prm_samples_pd + The parameter samples the forward model was evaluated at. + responses_over_experiments + The keys are the names of the experiment the forward model is associated + with, while the values are 3D-arrays containing the forward model's + response. responses_over_experiments[i][j] will contain the forward model's + response with the ith parameter-vector for the jth output sensor. + """ + + # get the forward model object with the given name and prepare the corresponding + # experimental in- and output dictionaries + forward_model = self.problem.forward_models[forward_model_name] + forward_model.prepare_experimental_inputs_and_outputs() + + # generate the latent parameter samples and convert it to a data frame to have + # the association between columns and parameter names + prm_samples = self.generate_samples(n_samples, seed=seed) + prm_samples_pd = pd.DataFrame( + prm_samples, columns=self.problem.get_theta_names() + ) + + # this dictionary will contain the forward model responses for each of the + # experiments associated with the forward model; so the keys will be experiment + # names while the values will be 3D-arrays with the forward model's responses; + # responses_over_experiments[i] will correspond to the response of the ith + # parameter vector; responses_over_experiments[i][j] will contain the forward + # model's response with the ith parameter vector for the jth output sensor + responses_over_experiments = {} + + # here, the length of the vectors of the forward model's output sensors is + # determined; to that end, the forward model is evaluated once + first_exp_name = forward_model.experiment_names[0] + exp_inp = forward_model.input_from_experiments[first_exp_name] + first_theta = prm_samples[0] + prms_model = self.problem.get_parameters(first_theta, forward_model.prms_def) + inp = {**exp_inp, **prms_model} + response_dict = forward_model.response(inp) + # make sure that the vectors returned by each of the forward model's output + # sensors has the same length; otherwise an AssertionError is raised + length_set = set() + for value in response_dict.values(): + length_set.add(len_or_one(value)) + assert len(length_set) == 1 + n_out_values = list(length_set)[0] + + # evaluate the forward model for each experiment/parameter vector + for exp_name in forward_model.experiment_names: + exp_inp = forward_model.input_from_experiments[exp_name] + response_array = np.zeros( + (n_samples, forward_model.n_output_sensors, n_out_values) + ) + for i, theta in enumerate(prm_samples): + prms_model = self.problem.get_parameters(theta, forward_model.prms_def) + inp = {**exp_inp, **prms_model} # adds the two dictionaries + response_dict = forward_model.response(inp) + for j, response_vector in enumerate(response_dict.values()): + response_array[i, j, :] = response_vector + responses_over_experiments[exp_name] = response_array + + return prm_samples_pd, responses_over_experiments diff --git a/setup.cfg b/setup.cfg index c72ec079..e4c53338 100644 --- a/setup.cfg +++ b/setup.cfg @@ -33,6 +33,7 @@ install_requires = owlready2<1 dynesty tri-py<1 + pandas<2 [options.package_data] probeye = diff --git a/tests/unit_tests/surrogate/test_initial_sampling.py b/tests/unit_tests/surrogate/test_initial_sampling.py new file mode 100644 index 00000000..0fc68742 --- /dev/null +++ b/tests/unit_tests/surrogate/test_initial_sampling.py @@ -0,0 +1,104 @@ +# standard library +import unittest + +# third party imports +import numpy as np + +# local imports +from probeye.definition.inverse_problem import InverseProblem +from probeye.definition.forward_model import ForwardModelBase +from probeye.definition.sensor import Sensor +from probeye.surrogate.initial_sampling import LatinHypercubeSampler + + +class LinearModel(ForwardModelBase): + def interface(self): + self.parameters = ["a", "b"] + self.input_sensors = Sensor("x") + self.output_sensors = Sensor("y", std_model="sigma") + + def response(self, inp: dict) -> dict: + x = inp["x"] + m = inp["a"] + b = inp["b"] + return {"y": m * x + b} + + +class TestProblem(unittest.TestCase): + def test_LatinHypercubeSampler_generate_samples_normal_case(self): + + # prepare a dummy problem as far as necessary for this test + problem = InverseProblem("Test") + problem.add_parameter("a", prior=("normal", {"mean": 2.0, "std": 1.0})) + problem.add_parameter("b", prior=("normal", {"mean": 1.0, "std": 0.5})) + + # create the samples from the priors + lhs = LatinHypercubeSampler(problem) + sample_array = lhs.generate_samples(10, seed=1) + + # check if everything is as expected + self.assertEqual(sample_array.shape, (10, 2)) + self.assertAlmostEqual(float(np.mean(sample_array[:, 0])), 2.0, delta=0.3) + self.assertAlmostEqual(float(np.std(sample_array[:, 0])), 1.0, delta=0.3) + self.assertAlmostEqual(float(np.mean(sample_array[:, 1])), 1.0, delta=0.3) + self.assertAlmostEqual(float(np.std(sample_array[:, 1])), 0.5, delta=0.2) + + def test_LatinHypercubeSampler_generate_samples_mv_prior(self): + + # prepare a dummy problem as far as necessary for this test + problem = InverseProblem("Test") + problem.add_parameter( + "ab", + dim=2, + prior=( + "multivariate-normal", + { + "mean": np.array([2.0, 1.0]), + "cov": np.array([[1.0, 0], [0, 0.5**2]]), + }, + ), + ) + + # try creating the samples from the priors; this should not work since the + # problem contains a multivariate prior + lhs = LatinHypercubeSampler(problem) + with self.assertRaises(RuntimeError): + _ = lhs.generate_samples(10, seed=1) + + def test_LatinHypercubeSampler_generate_training_data(self): + + # prepare a dummy problem as far as necessary for this test + problem = InverseProblem("Test") + problem.add_parameter("a", prior=("normal", {"mean": 2.0, "std": 1.0})) + problem.add_parameter("b", prior=("normal", {"mean": 1.0, "std": 0.5})) + problem.add_forward_model(LinearModel("LinearModel")) + problem.add_experiment( + exp_name="Exp1", + fwd_model_name="LinearModel", + sensor_values={ + "x": np.array([0.0, 1.0, 2.0]), + "y": np.array([1.1, 2.9, 5.3]), + }, + ) + + # create the training data + lhs = LatinHypercubeSampler(problem) + prm_samples_pd, responses_over_experiments = lhs.generate_training_data( + "LinearModel", 10, seed=1 + ) + + # check if everything is as expected (with prm_samples_pd) + sample_array = prm_samples_pd.values + self.assertEqual(sample_array.shape, (10, 2)) + self.assertTrue((prm_samples_pd.columns == ["a", "b"]).all()) + self.assertAlmostEqual(float(np.mean(sample_array[:, 0])), 2.0, delta=0.3) + self.assertAlmostEqual(float(np.std(sample_array[:, 0])), 1.0, delta=0.3) + self.assertAlmostEqual(float(np.mean(sample_array[:, 1])), 1.0, delta=0.3) + self.assertAlmostEqual(float(np.std(sample_array[:, 1])), 0.5, delta=0.2) + + # check if everything is as expected (with responses_over_experiments) + self.assertEqual(responses_over_experiments["Exp1"].shape, (10, 1, 3)) + + +if __name__ == "__main__": + unittest.main() From b0c0653af4272f0a21bf73c25767e45e11a67a3a Mon Sep 17 00:00:00 2001 From: Klawonn Date: Tue, 17 May 2022 15:21:45 +0200 Subject: [PATCH 03/15] removed compatability with Python 3.6 --- setup.cfg | 1 - 1 file changed, 1 deletion(-) diff --git a/setup.cfg b/setup.cfg index e4c53338..c4f865dc 100644 --- a/setup.cfg +++ b/setup.cfg @@ -7,7 +7,6 @@ description = A general framework for setting up parameter estimation problems. long_description = file: README.md long_description_content_type = text/markdown classifiers = - Programming Language :: Python :: 3.6 Programming Language :: Python :: 3.7 Programming Language :: Python :: 3.8 Programming Language :: Python :: 3.9 From 6598753aeeacc271432e371b150e7c84cad49642 Mon Sep 17 00:00:00 2001 From: Klawonn Date: Tue, 17 May 2022 15:26:28 +0200 Subject: [PATCH 04/15] had to do more changes to remove the Python 3.6 flag --- .github/workflows/push.yaml | 2 +- setup.cfg | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/push.yaml b/.github/workflows/push.yaml index 9bba1720..7af3e42e 100644 --- a/.github/workflows/push.yaml +++ b/.github/workflows/push.yaml @@ -34,7 +34,7 @@ jobs: needs: [ lint_and_type_check ] strategy: matrix: - python-version: [ 3.6, 3.7, 3.8, 3.9 ] + python-version: [ 3.7, 3.8, 3.9 ] steps: - uses: actions/checkout@v2 diff --git a/setup.cfg b/setup.cfg index c4f865dc..a55f1109 100644 --- a/setup.cfg +++ b/setup.cfg @@ -17,7 +17,7 @@ classifiers = license_files = LICENSE [options] -python_requires = >= 3.6 +python_requires = >= 3.7 packages = find: include_package_data = True install_requires = From c88e33e8f75cadbbfa027f9ac894e898072cb0f8 Mon Sep 17 00:00:00 2001 From: Klawonn Date: Thu, 19 May 2022 10:49:09 +0200 Subject: [PATCH 05/15] removed add_surrogate_model; some other small changes --- probeye/definition/inverse_problem.py | 130 ++++++++---------- probeye/surrogate/initial_sampling.py | 8 +- .../test_surrogate_model.py | 9 +- .../surrogate/test_initial_sampling.py | 5 +- 4 files changed, 66 insertions(+), 86 deletions(-) diff --git a/probeye/definition/inverse_problem.py b/probeye/definition/inverse_problem.py index d652db75..d846d3cd 100644 --- a/probeye/definition/inverse_problem.py +++ b/probeye/definition/inverse_problem.py @@ -10,15 +10,18 @@ # local imports from probeye.definition.parameter import Parameters from probeye.definition.forward_model import ForwardModelBase -from probeye.definition.surrogate_model import SurrogateModelBase from probeye.definition.likelihood_model import GaussianLikelihoodModel -from probeye.subroutines import underlined_string, titled_table -from probeye.subroutines import simplified_list_string -from probeye.subroutines import make_list, len_or_one -from probeye.subroutines import print_probeye_header -from probeye.subroutines import logging_setup -from probeye.subroutines import add_index_to_tex_prm_name -from probeye.subroutines import translate_simple_correlation +from probeye.subroutines import ( + underlined_string, + titled_table, + simplified_list_string, + make_list, + len_or_one, + print_probeye_header, + logging_setup, + add_index_to_tex_prm_name, + translate_simple_correlation, +) class InverseProblem: @@ -97,13 +100,6 @@ def __init__( # not be modified directly self._forward_models = {} # type: dict - # the following dict contains the problem's surrogate models (note that - # surrogate models are optional, not every problem has surrogate models); - # the keys in this dict are the surrogate model names, while the values are the - # surrogate model objects (check out the script surrogate_model.py); this - # dictionary is managed internally and should not be modified directly - self._surrogate_models = {} # type: dict - # this dictionary contains the problem's likelihood models; as the other private # attributes above, it is managed internally and should not be modified directly self._likelihood_models = {} # type: dict @@ -219,11 +215,6 @@ def forward_models(self) -> dict: """Access self._forward_models from outside via self.forward_models.""" return self._forward_models - @property - def surrogate_models(self) -> dict: - """Access self._surrogate_models from outside via self.forward_models.""" - return self._surrogate_models - @property def experiments(self) -> dict: """Access self._experiments from outside via self.experiments.""" @@ -439,13 +430,35 @@ def change_constant(self, prm_name: str, new_value: Union[int, float]): value=new_value ) - def get_latent_prior_hyperparameters(self, prm_name): + def get_latent_prior_hyperparameters(self, prm_name: str) -> list: + """ + Returns a list of the latent hyperparameters of a parameter's prior. In most + cases there will be none, so an empty list will be returned. + + Parameters + ---------- + prm_name + The name of the parameter the prior of which should be checked for latent + hyperparameters. + + Returns + ------- + latent_hyperparameters + Contains the global names of latent hyperparameters in the prior of the + parameter 'prm_name'. + """ + + # first, make sure that the given parameter exists + self._parameters.confirm_that_parameter_exists(prm_name) + + # now look for possible latent hyperparameters latent_hyperparameters = [] if self.parameters[prm_name].is_latent: hyperparameters = self.parameters[prm_name].prior.hyperparameters for hyperparameter in hyperparameters: if self.parameters[hyperparameter].is_latent: latent_hyperparameters.append(hyperparameter) + return latent_hyperparameters def get_parameters(self, theta: np.ndarray, prm_def: dict) -> dict: @@ -492,6 +505,27 @@ def get_parameters(self, theta: np.ndarray, prm_def: dict) -> dict: return prms def get_constants(self, prm_def: dict) -> dict: + """ + Similar to 'get_parameters' with the difference that this method only extracts + the numeric values of constants that have been defined within the problem scope. + For that reason it does not need the 'theta' argument as in 'get_parameters'. If + 'prm_def' also contains latent parameters, they will simply be ignored. + + Parameters + ---------- + prm_def + Defines which constants to extract. The keys of this dictionary are the + global parameter names, while the values are the local parameter names. In + most cases global and local names will be identical, but sometimes it is + convenient to define a local parameter name, e.g. in the forward model. + + Returns + ------- + prms + Contains : <(global) parameter value> pairs. If a + parameter is scalar, its value will be returned as a float. In case of a + vector-valued parameter, its value will be returned as a np.ndarray. + """ prms = {} for global_name, local_name in prm_def.items(): idx = self._parameters[global_name].index @@ -670,60 +704,6 @@ def add_forward_model(self, forward_model: ForwardModelBase): # the given forward model name self._forward_models[forward_model.name] = forward_model - def add_surrogate_model(self, surrogate_model: SurrogateModelBase): - """ - Adds a surrogate model to the inverse problem. Note that multiple forward models - can be added to one problem. - - Parameters - ---------- - surrogate_model - Defines the surrogate model. Check out surrogate_model.py to see a template - for the surrogate model definition. The user will then have to derive her - own surrogate model from that base class. Examples can be found in the - package directory tests/integration_tests. - """ - - # log at beginning so that errors can be quickly associated - logger.debug(f"Adding surrogate model '{surrogate_model.name}'") - - # check if all given model parameters have already been added to the inverse - # problem; note that the surrogate model can only be added to the problem after - # the corresponding parameters and the forward model have been defined - for prm_name in surrogate_model.prms_def: - self._parameters.confirm_that_parameter_exists(prm_name) - # check, if the type was set correctly; if it was not set yet, it will be - # set automatically here - if self._parameters[prm_name].type == "not defined": - self.change_parameter_type(prm_name, "model") - elif self._parameters[prm_name].type in ["prior", "likelihood"]: - raise ValueError( - f"The parameter '{prm_name}' defined in forward model " - f"'{surrogate_model.name}' was assigned the type " - f"'{self._parameters[prm_name].type}' (it should be 'model')." - ) - - # check if the given name for the forward model has already been used - if surrogate_model.name in self._surrogate_models: - raise RuntimeError( - f"The name '{surrogate_model.name}' of the surrogate model you are " - f"trying to add to the problem has already been used for another" - f"surrogate model within the problem scope. Please choose another name." - ) - - # check if the surrogate model's forward model has already been added to the - # inverse problem - if surrogate_model.forward_model.name not in self._forward_models: - raise RuntimeError( - f"The surrogate model's forward model " - f"'{surrogate_model.forward_model.name}' has not been added to the " - f"problem yet!" - ) - - # add the given surrogate model to the internal surrogate model dictionary under - # the given surrogate model name - self._surrogate_models[surrogate_model.name] = surrogate_model - # =============================================================== # # Experiments related methods # # =============================================================== # diff --git a/probeye/surrogate/initial_sampling.py b/probeye/surrogate/initial_sampling.py index a9d02ab1..a4924f10 100644 --- a/probeye/surrogate/initial_sampling.py +++ b/probeye/surrogate/initial_sampling.py @@ -9,6 +9,7 @@ # local imports from probeye.definition.inverse_problem import InverseProblem +from probeye.definition.forward_model import ForwardModelBase from probeye.inference.priors import translate_prior from probeye.subroutines import len_or_one @@ -93,7 +94,7 @@ def generate_samples(self, n_samples: int, seed: int = 1) -> np.ndarray: def generate_training_data( self, - forward_model_name: str, + forward_model: ForwardModelBase, n_samples: int, seed: int = 1, ) -> Tuple[pd.DataFrame, dict]: @@ -104,8 +105,8 @@ def generate_training_data( Parameters ---------- - forward_model_name - The name of the forward model that should be evaluated. + forward_model + The forward model that should be evaluated. n_samples The number of parameter vectors the forward model should be evaluated for. seed @@ -124,7 +125,6 @@ def generate_training_data( # get the forward model object with the given name and prepare the corresponding # experimental in- and output dictionaries - forward_model = self.problem.forward_models[forward_model_name] forward_model.prepare_experimental_inputs_and_outputs() # generate the latent parameter samples and convert it to a data frame to have diff --git a/tests/integration_tests/tests_without_correlation/test_surrogate_model.py b/tests/integration_tests/tests_without_correlation/test_surrogate_model.py index 19e9787e..d39754df 100644 --- a/tests/integration_tests/tests_without_correlation/test_surrogate_model.py +++ b/tests/integration_tests/tests_without_correlation/test_surrogate_model.py @@ -114,13 +114,12 @@ def response(self, inp: dict) -> dict: prior=("uniform", {"low": low_sigma, "high": high_sigma}), ) - # add the forward model to the problem + # add the surrogate model (which is essentially a forward model) to the problem; + # whether or not the surrogate model gets the original forward model is not + # finally decided on, but it might not be necessary forward_model = ExpensiveModel("ExpensiveModel") - problem.add_forward_model(forward_model) - - # add the forward model to the problem surrogate_model = SurrogateModel("FastModel", forward_model=forward_model) - problem.add_surrogate_model(surrogate_model) + problem.add_forward_model(surrogate_model) # ============================================================================ # # Add test data to the Inference Problem # diff --git a/tests/unit_tests/surrogate/test_initial_sampling.py b/tests/unit_tests/surrogate/test_initial_sampling.py index 0fc68742..eb79d73d 100644 --- a/tests/unit_tests/surrogate/test_initial_sampling.py +++ b/tests/unit_tests/surrogate/test_initial_sampling.py @@ -71,7 +71,8 @@ def test_LatinHypercubeSampler_generate_training_data(self): problem = InverseProblem("Test") problem.add_parameter("a", prior=("normal", {"mean": 2.0, "std": 1.0})) problem.add_parameter("b", prior=("normal", {"mean": 1.0, "std": 0.5})) - problem.add_forward_model(LinearModel("LinearModel")) + forward_model = LinearModel("LinearModel") + problem.add_forward_model(forward_model) problem.add_experiment( exp_name="Exp1", fwd_model_name="LinearModel", @@ -84,7 +85,7 @@ def test_LatinHypercubeSampler_generate_training_data(self): # create the training data lhs = LatinHypercubeSampler(problem) prm_samples_pd, responses_over_experiments = lhs.generate_training_data( - "LinearModel", 10, seed=1 + forward_model, 10, seed=1 ) # check if everything is as expected (with prm_samples_pd) From 650f621c2649fa1a393c0ebabe0c8c27a1ce91f7 Mon Sep 17 00:00:00 2001 From: Klawonn Date: Thu, 19 May 2022 10:56:31 +0200 Subject: [PATCH 06/15] fixed bug in test_surrogate_model.py --- .../tests_without_correlation/test_surrogate_model.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/integration_tests/tests_without_correlation/test_surrogate_model.py b/tests/integration_tests/tests_without_correlation/test_surrogate_model.py index d39754df..a6cd4e6e 100644 --- a/tests/integration_tests/tests_without_correlation/test_surrogate_model.py +++ b/tests/integration_tests/tests_without_correlation/test_surrogate_model.py @@ -136,7 +136,7 @@ def response(self, inp: dict) -> dict: # add the experimental data problem.add_experiment( f"TestSeries_1", - fwd_model_name="ExpensiveModel", + fwd_model_name="FastModel", sensor_values={ forward_model.input_sensor.name: x_test, forward_model.output_sensor.name: y_test, From e899df08465d933e09eb226a8d54eb2c6f21ae72 Mon Sep 17 00:00:00 2001 From: Klawonn Date: Thu, 19 May 2022 11:36:48 +0200 Subject: [PATCH 07/15] added a missing import statement in inverse_problem.py --- probeye/definition/inverse_problem.py | 1 + 1 file changed, 1 insertion(+) diff --git a/probeye/definition/inverse_problem.py b/probeye/definition/inverse_problem.py index 9740337a..9b1ee6e5 100644 --- a/probeye/definition/inverse_problem.py +++ b/probeye/definition/inverse_problem.py @@ -21,6 +21,7 @@ logging_setup, add_index_to_tex_prm_name, translate_simple_correlation, + safe_string, ) From 3ea64be2a42899f5728e3b3d727abb4c64dc981c Mon Sep 17 00:00:00 2001 From: Koune Date: Thu, 9 Jun 2022 13:48:16 +0200 Subject: [PATCH 08/15] added example of inference with surrogate of hartman 6D function --- examples/example_surrogate_hartman6D.py | 424 ++++++++++++++++++++++++ 1 file changed, 424 insertions(+) create mode 100644 examples/example_surrogate_hartman6D.py diff --git a/examples/example_surrogate_hartman6D.py b/examples/example_surrogate_hartman6D.py new file mode 100644 index 00000000..5868e487 --- /dev/null +++ b/examples/example_surrogate_hartman6D.py @@ -0,0 +1,424 @@ +""" +Example of a Bayesian parameter estimation problem using surrogate +modeling with probeye. + +The Hartmann test function f:[0, 1]^6 -> R^1 is used to simulate a +physical model. The last two dimensions are considered as space and +time coordinates, while the first four dimensions are taken as +latent variables to be inferred. Measurements are generated by +adding I.i.d. Gaussian noise to samples from this function. + +Notes: + * torch, botorch and harlow are required to run this example + * the sampling and surrogating package `harlow` will be made + public soon +""" + +# ========================================================================= +# Imports +# ========================================================================= + +# third party imports +import numpy as np + +# local imports (problem definition) +from probeye.definition.inverse_problem import InverseProblem +from probeye.definition.forward_model import ForwardModelBase +from probeye.definition.sensor import Sensor +from probeye.definition.likelihood_model import GaussianLikelihoodModel +from probeye.surrogate.initial_sampling import LatinHypercubeSampler +from probeye.inference.emcee.solver import EmceeSolver + +# local imports (inference data post-processing) +from probeye.postprocessing.sampling import create_pair_plot +from probeye.postprocessing.sampling import create_posterior_plot + +# Surrogate model imports +from harlow.surrogate_model import ModelListGaussianProcess + +import torch +from botorch.test_functions.synthetic import Hartmann +from matplotlib import pyplot as plt + +# ========================================================================= +# General settings +# ========================================================================= + +plot = True + +# Emcee settings +n_steps = 1_000 +n_init_steps = 200 +n_walkers = 20 + +# Surrogate settings +N_train = 50 +N_train_iter = 200 +show_progress = True + +# ========================================================================= +# Define parameters +# ========================================================================= + +# Ground truth +X_true = np.array([0.5, 0.5, 0.5, 0.5]) + +# Bounds for function defined on unit hypercube +X_low = 0.0 +X_high = 1.0 + +# Ground truth and prior for measurement uncertainty std. dev. +std_true = 0.05 +std_low = 0.0 +std_high = 1.0 + +# ========================================================================= +# Define physical model +# ========================================================================= + +# Number of sensors and number of points in timeseries +Nx = 3 +Nt = 5 + +# Sensor names and positions +sensor_names = ["S" + str(i + 1) for i in range(Nx)] +x_vec = np.linspace(0, 1, Nx) +t_vec = np.linspace(0, 1, Nt) + +isensor = Sensor("t") +osensor_list = [ + Sensor(sensor_names[i], x=float(x_vec[i]), std_model="sigma") for i in range(Nx) +] + +# ========================================================================= +# Define forward model +# ========================================================================= + +# Initialize model +expensive_model = Hartmann(noise_std=0.00001) + + +class SyntheticModel(ForwardModelBase): + def interface(self): + self.parameters = ["X" + str(i + 1) for i in range(4)] + self.input_sensors = isensor + self.output_sensors = osensor_list + + def response(self, inp: dict) -> dict: + + # Arange input vector + params = np.tile([inp["X" + str(i + 1)] for i in range(4)], (Nx * Nt, 1)) + xt = np.array(np.meshgrid(x_vec, t_vec)).T.reshape(-1, 2) + X = torch.tensor(np.hstack((params, xt))) + + # Evaluate function and arange output on grid + f = np.array(expensive_model(X)).reshape(Nx, Nt) + + # Store prediction as dict + response = dict() + for idx_x, os in enumerate(self.output_sensors): + response[os.name] = f[idx_x, :] + return response + + +# ========================================================================= +# Define surrogate model +# ========================================================================= + + +class SurrogateModel(SyntheticModel): + """ + The inheritance from ExpensiveModel 'copies' the interface-method from + ExpensiveModel (the surrogate model should have the same interface as the + forward model). The inheritance from SurrogateModelBase is required to + assign a forward model to the surrogate model, see surrogate_model.py. + """ + + def __init__( + self, + name: str, + forward_model: ForwardModelBase, + surrogate_model, + surrogate_kwargs, + input_transform=None, + output_transform=None, + ): + + super().__init__(name, forward_model) + self.forward_model = forward_model + self.surrogate_model = surrogate_model + self.surrogate_kwargs = surrogate_kwargs + self.input_transform = input_transform + self.output_transform = output_transform + + # List of model names for ModelListGP + self.model_names = [] + for idx_x, os in enumerate(self.output_sensors): + self.model_names.append(os.name) + + if not self.input_transform: + self.input_transform = lambda x: x + if not self.output_transform: + self.output_transform = lambda y: y + + def fit(self, train_X, train_y): + + # Initialize model + self.model = self.surrogate_model( + self.input_transform(train_X), + self.output_transform(train_y), + self.model_names, + **self.surrogate_kwargs, + ) + + # Fit model + self.model.fit() + + def response(self, inp: dict) -> dict: + + params = np.tile([inp["X" + str(i + 1)] for i in range(4)], (Nt, 1)) + X = input_transform(torch.tensor(np.hstack((params, t_vec.reshape(-1, 1))))) + + # # ========================================== + # # If surrogate is regular GP + # # ========================================== + response = dict() + + # Evaluate function and arange output on grid + _ = self.model.predict(X, return_std=False) + f = np.array(self.model.mean).T + + for idx_x, os in enumerate(self.output_sensors): + response[os.name] = f[idx_x, :] + + return response + + +# ========================================================================= +# Define inference problem +# ========================================================================= +problem = InverseProblem("Parameter estimation using surrogate model") + +# Parameters of the Hartmann function +for i in range(4): + problem.add_parameter( + "X" + str(i + 1), + "model", + prior=("uniform", {"low": X_low, "high": X_high}), + info="Parameter of the 6D Hartmann function", + tex=r"$X_{{{}}}$".format(i + 1), + ) + +# Noise std. dev. +problem.add_parameter( + "sigma", + "likelihood", + prior=("uniform", {"low": std_low, "high": std_high}), + info="Std. dev. of zero-mean noise model", + tex=r"$\sigma$", +) + +# add the forward model to the problem +forward_model = SyntheticModel("ExpensiveModel") + +# Add Expensive FE model to forward models +problem.add_forward_model(forward_model) + +# ========================================================================= +# Create surrogate model +# ========================================================================= + +list_params = [[0, 1, 2, 3, 4]] * len(sensor_names) + +surrogate_kwargs = { + "training_max_iter": N_train_iter, + "list_params": list_params, + "show_progress": True, + "silence_warnings": True, + "fast_pred_var": True, +} + + +def input_transform(X): + return torch.tensor(X).float() + + +def output_transform(y): + return torch.tensor(y).float() + + +surrogate_model = SurrogateModel( + "FastModel", + forward_model, + ModelListGaussianProcess, + surrogate_kwargs=surrogate_kwargs, + input_transform=input_transform, + output_transform=output_transform, +) + +problem.add_forward_model(surrogate_model) + +# ========================================================================= +# Add test data to the inference problem +# ========================================================================= + + +def generate_data(): + inp = {"X" + str(idx + 1): X_i for idx, X_i in enumerate(X_true)} + sensors = forward_model(inp) + for sname, svals in sensors.items(): + sensors[sname] = list( + np.array(svals) + np.random.normal(loc=0.0, scale=std_true, size=Nt) + ) + sensors[isensor.name] = t_vec + problem.add_experiment( + "TestSeriesSurrogate", sensor_values=sensors, fwd_model_name="FastModel" + ) + problem.add_experiment( + "TestSeriesFull", sensor_values=sensors, fwd_model_name="ExpensiveModel" + ) + + +generate_data() + +# ========================================================================= +# Create training data +# ========================================================================= +sampler = LatinHypercubeSampler(problem=problem) +test_samples = sampler.generate_samples(N_train) +train_samples, train_data = sampler.generate_training_data(forward_model, N_train) +train_X = train_samples.to_numpy()[:, :-1] +train_y = train_data["TestSeriesFull"] + +# Reshape `train_X` and `train_y +t_arr = np.tile(t_vec, N_train).reshape(-1, 1) +train_X_arr = np.repeat(train_X, Nt, axis=0) +train_X_arr = np.hstack((train_X_arr, t_arr)) + +train_y_arr = np.transpose(train_y, axes=[0, 2, 1]) +train_y_arr = train_y_arr.reshape(N_train * Nt, Nx) + +# Train surrogate +model_names = sensor_names +surrogate_model.fit(train_X_arr, train_y_arr) + + +# ==================================================================== +# Plot surrogate vs. FE model prediction +# ==================================================================== + +# Physical model prediction +inp = {"X" + str(idx + 1): X_i for idx, X_i in enumerate(X_true)} +y_true = forward_model.response(inp) +y_pred = surrogate_model.response(inp) +upper = list(surrogate_model.model.cr_u.squeeze().detach().numpy().T) +lower = list(surrogate_model.model.cr_l.squeeze().detach().numpy().T) + +# Initialize plots +nrows = 3 +ncols = int(np.ceil(len(sensor_names) / 3)) + +# Initialize plot +f, axes = plt.subplots(nrows, ncols, sharex=True, figsize=(3 * ncols, 3 * nrows)) + +for j in range(len(sensor_names)): + + ax_i = axes.ravel()[j] + grid_idx = np.unravel_index(j, (nrows, ncols)) + + # Plot training data as black stars + ax_i.plot(t_vec, train_y[:, j, :].T, "k*", label="Observations") + + # Predictive mean as blue line + ax_i.plot(t_vec, y_pred[sensor_names[j]], "b", label="Mean") + + # Shade in confidence + ax_i.fill_between( + t_vec, + lower[j], + upper[j], + alpha=0.5, + label="Confidence", + ) + ax_i.plot( + t_vec, + y_true[sensor_names[j]], + color="red", + linestyle="dashed", + label="Model", + ) + ax_i.set_title(f"Sensor: {sensor_names[j]}") + +axes = np.atleast_2d(axes) +# axes[0, 0].legend() +[ax_i.set_xlabel(r"$X_6$") for ax_i in axes[-1, :]] +[ax_i.set_ylabel("Stress [MPa]") for ax_i in axes[:, 0]] +plt.show() + + +# ========================================================================= +# Add noise models +# ========================================================================= + +# add the likelihood models to the problem +for osensor in osensor_list: + problem.add_likelihood_model( + GaussianLikelihoodModel( + experiment_name="TestSeriesSurrogate", + prms_def=["sigma"], + model_error="additive", + additive_measurement_error=False, + correlation_variables=None, + ) + ) + + +# Problem overview +problem.info() + +true_values = { + "X1": X_true[0], + "X2": X_true[1], + "X3": X_true[2], + "X4": X_true[3], + "sigma": std_true, +} + +# ========================================================================= +# Initialize and run solver +# ========================================================================= + +emcee_solver = EmceeSolver( + problem, + show_progress=True, +) + + +inference_data = emcee_solver.run_mcmc( + n_walkers=n_walkers, n_steps=n_steps, n_initial_steps=n_init_steps, vectorize=False +) + +# ========================================================================= +# Plotting +# ========================================================================= +create_pair_plot( + inference_data, + problem, + show=False, + true_values=true_values, + title="Joint posterior", +) + +create_posterior_plot( + inference_data, + problem, + show=False, + true_values=true_values, + title="Marginal posteriors", +) + + +if plot: + plt.show() # shows all plots at once due to 'show=False' above +else: + plt.close("all") From 72eb7823ee704bf8228679da012bf78017922b09 Mon Sep 17 00:00:00 2001 From: Koune Date: Sun, 17 Jul 2022 23:54:22 +0200 Subject: [PATCH 09/15] added surrogating example and temporarily disabled mypy from pre-commit hooks --- .pre-commit-config.yaml | 10 +- examples/example_surrogate_hartman6D.py | 177 +++++++++++++++++------- 2 files changed, 135 insertions(+), 52 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 21c96cde..df79e386 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -6,11 +6,11 @@ repos: args: ['--maxkb=100'] - id: check-yaml -- repo: https://github.com/pre-commit/mirrors-mypy - rev: v0.910 - hooks: - - id: mypy - args: ['--install-types', '--non-interactive', '--ignore-missing-imports'] +#- repo: https://github.com/pre-commit/mirrors-mypy +# rev: v0.910 +# hooks: +# - id: mypy +# args: ['--install-types', '--non-interactive', '--ignore-missing-imports'] - repo: https://github.com/psf/black rev: 22.3.0 diff --git a/examples/example_surrogate_hartman6D.py b/examples/example_surrogate_hartman6D.py index 5868e487..51c46d00 100644 --- a/examples/example_surrogate_hartman6D.py +++ b/examples/example_surrogate_hartman6D.py @@ -34,7 +34,11 @@ from probeye.postprocessing.sampling import create_posterior_plot # Surrogate model imports -from harlow.surrogate_model import ModelListGaussianProcess +from harlow.sampling import Sampler +from harlow.sampling import LatinHypercube +from harlow.surrogating import Surrogate +from harlow.surrogating import ModelListGaussianProcess +from harlow.utils.transforms import ExpandDims import torch from botorch.test_functions.synthetic import Hartmann @@ -126,17 +130,16 @@ def response(self, inp: dict) -> dict: # ========================================================================= -class SurrogateModel(SyntheticModel): +class HarlowSurrogate(ForwardModelBase): """ The inheritance from ExpensiveModel 'copies' the interface-method from ExpensiveModel (the surrogate model should have the same interface as the - forward model). The inheritance from SurrogateModelBase is required to - assign a forward model to the surrogate model, see surrogate_model.py. + forward model). """ def __init__( self, - name: str, + name, forward_model: ForwardModelBase, surrogate_model, surrogate_kwargs, @@ -144,49 +147,61 @@ def __init__( output_transform=None, ): - super().__init__(name, forward_model) + # Get interface from existing forward model self.forward_model = forward_model - self.surrogate_model = surrogate_model - self.surrogate_kwargs = surrogate_kwargs + self.interface = self.forward_model.interface + + super().__init__(name, forward_model) + self.input_transform = input_transform self.output_transform = output_transform + self.surrogate_kwargs = surrogate_kwargs + self._surrogate = surrogate_model + self.lower_bounds = [] + self.upper_bounds = [] - # List of model names for ModelListGP - self.model_names = [] - for idx_x, os in enumerate(self.output_sensors): - self.model_names.append(os.name) + # Common input arguments for surrogate models + if self.surrogate_kwargs["model_names"] == None: + self.model_names = [] + for idx_x, os in enumerate(self.output_sensors): + self.model_names.append(os.name) + self.surrogate_kwargs["model_names"] = self.model_names - if not self.input_transform: - self.input_transform = lambda x: x - if not self.output_transform: - self.output_transform = lambda y: y + if self.surrogate_kwargs["num_features"] == None: + self.num_features = len(self.output_sensors) + self.surrogate_kwargs["num_features"] = self.num_features - def fit(self, train_X, train_y): + # Get bounds + # TODO: - # Initialize model - self.model = self.surrogate_model( - self.input_transform(train_X), - self.output_transform(train_y), - self.model_names, - **self.surrogate_kwargs, - ) + # Get param names + + # Initialize surrogate + self.surrogate_model = self._surrogate(**self.surrogate_kwargs) - # Fit model - self.model.fit() + def _cast_to_numpy(self, inp: dict) -> np.ndarray: + raise NotImplementedError + + def _cast_to_dict(self, X: np.ndarray) -> dict: + raise NotImplementedError + + def target(self, X: np.ndarray) -> np.ndarray: + inp = self._cast_to_dict(X) + response = self.forward_model.response(inp) + y = [[]] * len(self.output_sensors) + for idx, os in enumerate(self.output_sensors): + y[idx] = response[os.name] + return np.array(y).T def response(self, inp: dict) -> dict: params = np.tile([inp["X" + str(i + 1)] for i in range(4)], (Nt, 1)) - X = input_transform(torch.tensor(np.hstack((params, t_vec.reshape(-1, 1))))) + X = torch.tensor(np.hstack((params, t_vec.reshape(-1, 1)))) - # # ========================================== - # # If surrogate is regular GP - # # ========================================== response = dict() # Evaluate function and arange output on grid - _ = self.model.predict(X, return_std=False) - f = np.array(self.model.mean).T + f = self.model.predict(X, return_std=False) for idx_x, os in enumerate(self.output_sensors): response[os.name] = f[idx_x, :] @@ -194,6 +209,85 @@ def response(self, inp: dict) -> dict: return response +# class HarlowSurrogate(ForwardModelBase): +# """ +# The inheritance from ExpensiveModel 'copies' the interface-method from +# ExpensiveModel (the surrogate model should have the same interface as the +# forward model). The inheritance from SurrogateModelBase is required to +# assign a forward model to the surrogate model, see surrogate_model.py. +# """ +# +# def __init__( +# self, +# name: str, +# forward_model: ForwardModelBase, +# surrogate_model, +# surrogate_kwargs, +# input_transform=None, +# output_transform=None, +# ): +# +# super().__init__(name, forward_model) +# self.forward_model = forward_model +# self.surrogate_model = surrogate_model +# self.surrogate_kwargs = surrogate_kwargs +# self.input_transform = input_transform +# self.output_transform = output_transform +# +# # List of model names for ModelListGP +# self.model_names = [] +# for idx_x, os in enumerate(self.output_sensors): +# self.model_names.append(os.name) +# +# if not self.input_transform: +# self.input_transform = lambda x: x +# if not self.output_transform: +# self.output_transform = lambda y: y +# +# def fit(self, train_X, train_y): +# +# # Initialize model +# self.model = self.surrogate_model( +# self.input_transform(train_X), +# self.output_transform(train_y), +# **self.surrogate_kwargs, +# ) +# +# # Fit model +# self.model.fit() +# +# +# def update(self, train_X, train_y): +# +# # Initialize model +# self.model = self.surrogate_model( +# self.input_transform(train_X), +# self.output_transform(train_y), +# **self.surrogate_kwargs, +# ) +# +# # Fit model +# self.model.fit() +# +# def response(self, inp: dict) -> dict: +# +# params = np.tile([inp["X" + str(i + 1)] for i in range(4)], (Nt, 1)) +# X = input_transform(torch.tensor(np.hstack((params, t_vec.reshape(-1, 1))))) +# +# # # ========================================== +# # # If surrogate is regular GP +# # # ========================================== +# response = dict() +# +# # Evaluate function and arange output on grid +# _ = self.model.predict(X, return_std=False) +# f = np.array(self.model.mean).T +# +# for idx_x, os in enumerate(self.output_sensors): +# response[os.name] = f[idx_x, :] +# +# return response + # ========================================================================= # Define inference problem # ========================================================================= @@ -238,22 +332,11 @@ def response(self, inp: dict) -> dict: "fast_pred_var": True, } - -def input_transform(X): - return torch.tensor(X).float() - - -def output_transform(y): - return torch.tensor(y).float() - - -surrogate_model = SurrogateModel( - "FastModel", - forward_model, - ModelListGaussianProcess, +surrogate_model = HarlowSurrogate( + name="FastModel", + forward_model=forward_model, + surrogate_model=ModelListGaussianProcess, surrogate_kwargs=surrogate_kwargs, - input_transform=input_transform, - output_transform=output_transform, ) problem.add_forward_model(surrogate_model) From 6226ca0269e4d38c2addc1b5aa3d5d3b22dce9ab Mon Sep 17 00:00:00 2001 From: Koune Date: Mon, 18 Jul 2022 15:21:40 +0200 Subject: [PATCH 10/15] draft example of combining probeye with harlow for surrogating and sampling --- examples/example_surrogate_hartman6D.py | 316 ++++++++++-------------- 1 file changed, 137 insertions(+), 179 deletions(-) diff --git a/examples/example_surrogate_hartman6D.py b/examples/example_surrogate_hartman6D.py index 51c46d00..9f1edd38 100644 --- a/examples/example_surrogate_hartman6D.py +++ b/examples/example_surrogate_hartman6D.py @@ -26,16 +26,16 @@ from probeye.definition.forward_model import ForwardModelBase from probeye.definition.sensor import Sensor from probeye.definition.likelihood_model import GaussianLikelihoodModel -from probeye.surrogate.initial_sampling import LatinHypercubeSampler from probeye.inference.emcee.solver import EmceeSolver +from probeye.definition.distribution import Uniform # local imports (inference data post-processing) -from probeye.postprocessing.sampling import create_pair_plot -from probeye.postprocessing.sampling import create_posterior_plot +from probeye.postprocessing.sampling_plots import create_pair_plot +from probeye.postprocessing.sampling_plots import create_posterior_plot # Surrogate model imports from harlow.sampling import Sampler -from harlow.sampling import LatinHypercube +from harlow.sampling import LatinHypercube, FuzzyLolaVoronoi from harlow.surrogating import Surrogate from harlow.surrogating import ModelListGaussianProcess from harlow.utils.transforms import ExpandDims @@ -81,8 +81,8 @@ # ========================================================================= # Number of sensors and number of points in timeseries -Nx = 3 -Nt = 5 +Nx = 1 +Nt = 10 # Sensor names and positions sensor_names = ["S" + str(i + 1) for i in range(Nx)] @@ -130,163 +130,117 @@ def response(self, inp: dict) -> dict: # ========================================================================= -class HarlowSurrogate(ForwardModelBase): - """ - The inheritance from ExpensiveModel 'copies' the interface-method from - ExpensiveModel (the surrogate model should have the same interface as the - forward model). - """ +def harlow_model_factory( + problem: InverseProblem, + forward_model: ForwardModelBase, + surrogate_model: Surrogate, + **kwargs, +) -> ForwardModelBase: + class HarlowSurrogate(ForwardModelBase): + """ + The inheritance from ExpensiveModel 'copies' the interface-method from + ExpensiveModel (the surrogate model should have the same interface as the + forward model). + """ + + def __init__( + self, + name, + ): + + self.name = name + self.kwargs = kwargs + self.lower_bounds = [] + self.upper_bounds = [] + self.bounds = {} + + self.forward_model = forward_model + self.problem = problem + self._surrogate = surrogate_model - def __init__( - self, - name, - forward_model: ForwardModelBase, - surrogate_model, - surrogate_kwargs, - input_transform=None, - output_transform=None, - ): + self.model_names = [] + for idx_x, os in enumerate(forward_model.output_sensors): + self.model_names.append(os.name) + self.kwargs["model_names"] = self.model_names - # Get interface from existing forward model - self.forward_model = forward_model - self.interface = self.forward_model.interface + self.num_features = len(forward_model.output_sensors) + self.kwargs["num_features"] = self.num_features - super().__init__(name, forward_model) + # Check inputs + if "input_transform" not in self.kwargs.keys(): + self.input_transform = None + else: + self.input_transform = self.kwargs["input_transform"] - self.input_transform = input_transform - self.output_transform = output_transform - self.surrogate_kwargs = surrogate_kwargs - self._surrogate = surrogate_model - self.lower_bounds = [] - self.upper_bounds = [] + if "output_transform" not in self.kwargs.keys(): + self.output_transform = None + else: + self.output_transform = self.kwargs["output_transform"] - # Common input arguments for surrogate models - if self.surrogate_kwargs["model_names"] == None: - self.model_names = [] - for idx_x, os in enumerate(self.output_sensors): - self.model_names.append(os.name) - self.surrogate_kwargs["model_names"] = self.model_names + # Initialize surrogate + self.model = self._surrogate(**kwargs) + self.initialize_forward_model() - if self.surrogate_kwargs["num_features"] == None: - self.num_features = len(self.output_sensors) - self.surrogate_kwargs["num_features"] = self.num_features + # Get interface from existing forward model + def interface(self): + self.input_sensors = self.forward_model.input_sensors + self.output_sensors = self.forward_model.output_sensors + self.parameters = self.forward_model.parameters - # Get bounds - # TODO: + def initialize_forward_model(self): + super().__init__(self.name) + self._get_bounds() - # Get param names + # TODO: This will likely fail for multidimensional parameters. FIX! + def _get_bounds(self): - # Initialize surrogate - self.surrogate_model = self._surrogate(**self.surrogate_kwargs) + # Find the priors corresponding to model parameters + for key, value in self.priors.items(): - def _cast_to_numpy(self, inp: dict) -> np.ndarray: - raise NotImplementedError + # Error if any of the model priors are not uniform + if value.prior_type != "uniform": + raise ValueError( + f"Non-uniform distribution of type {value.prior_type}" + f" found for prior {key}. The `HarlowSurrogate` currently" + f" only supports uniform priors." + ) - def _cast_to_dict(self, X: np.ndarray) -> dict: - raise NotImplementedError + for param in self.problem.model_prms: + lb = self.const_prms_dict["low_" + param] + ub = self.const_prms_dict["high_" + param] + self.bounds[param] = {"low": lb, "high": ub} - def target(self, X: np.ndarray) -> np.ndarray: - inp = self._cast_to_dict(X) - response = self.forward_model.response(inp) - y = [[]] * len(self.output_sensors) - for idx, os in enumerate(self.output_sensors): - y[idx] = response[os.name] - return np.array(y).T + def _cast_to_numpy(self, inp: dict) -> np.ndarray: + raise NotImplementedError - def response(self, inp: dict) -> dict: + def _cast_to_dict(self, X: np.ndarray) -> dict: + raise NotImplementedError - params = np.tile([inp["X" + str(i + 1)] for i in range(4)], (Nt, 1)) - X = torch.tensor(np.hstack((params, t_vec.reshape(-1, 1)))) + def target(self, X: np.ndarray) -> np.ndarray: + inp = self._cast_to_dict(X) + response = self.forward_model.response(inp) + y = [[]] * len(self.output_sensors) + for idx, os in enumerate(self.output_sensors): + y[idx] = response[os.name] + return np.array(y).T - response = dict() + def response(self, inp: dict) -> dict: - # Evaluate function and arange output on grid - f = self.model.predict(X, return_std=False) + params = np.tile([inp["X" + str(i + 1)] for i in range(4)], (Nt, 1)) + X = np.hstack((params, t_vec.reshape(-1, 1))) - for idx_x, os in enumerate(self.output_sensors): - response[os.name] = f[idx_x, :] + response = dict() - return response + # Evaluate function and arange output on grid + f = self.model.predict(X, return_std=False) + + for idx_x, os in enumerate(self.output_sensors): + response[os.name] = f[idx_x, :] + return response + + return HarlowSurrogate -# class HarlowSurrogate(ForwardModelBase): -# """ -# The inheritance from ExpensiveModel 'copies' the interface-method from -# ExpensiveModel (the surrogate model should have the same interface as the -# forward model). The inheritance from SurrogateModelBase is required to -# assign a forward model to the surrogate model, see surrogate_model.py. -# """ -# -# def __init__( -# self, -# name: str, -# forward_model: ForwardModelBase, -# surrogate_model, -# surrogate_kwargs, -# input_transform=None, -# output_transform=None, -# ): -# -# super().__init__(name, forward_model) -# self.forward_model = forward_model -# self.surrogate_model = surrogate_model -# self.surrogate_kwargs = surrogate_kwargs -# self.input_transform = input_transform -# self.output_transform = output_transform -# -# # List of model names for ModelListGP -# self.model_names = [] -# for idx_x, os in enumerate(self.output_sensors): -# self.model_names.append(os.name) -# -# if not self.input_transform: -# self.input_transform = lambda x: x -# if not self.output_transform: -# self.output_transform = lambda y: y -# -# def fit(self, train_X, train_y): -# -# # Initialize model -# self.model = self.surrogate_model( -# self.input_transform(train_X), -# self.output_transform(train_y), -# **self.surrogate_kwargs, -# ) -# -# # Fit model -# self.model.fit() -# -# -# def update(self, train_X, train_y): -# -# # Initialize model -# self.model = self.surrogate_model( -# self.input_transform(train_X), -# self.output_transform(train_y), -# **self.surrogate_kwargs, -# ) -# -# # Fit model -# self.model.fit() -# -# def response(self, inp: dict) -> dict: -# -# params = np.tile([inp["X" + str(i + 1)] for i in range(4)], (Nt, 1)) -# X = input_transform(torch.tensor(np.hstack((params, t_vec.reshape(-1, 1))))) -# -# # # ========================================== -# # # If surrogate is regular GP -# # # ========================================== -# response = dict() -# -# # Evaluate function and arange output on grid -# _ = self.model.predict(X, return_std=False) -# f = np.array(self.model.mean).T -# -# for idx_x, os in enumerate(self.output_sensors): -# response[os.name] = f[idx_x, :] -# -# return response # ========================================================================= # Define inference problem @@ -298,7 +252,7 @@ def response(self, inp: dict) -> dict: problem.add_parameter( "X" + str(i + 1), "model", - prior=("uniform", {"low": X_low, "high": X_high}), + prior=Uniform(low=X_low, high=X_high), info="Parameter of the 6D Hartmann function", tex=r"$X_{{{}}}$".format(i + 1), ) @@ -307,7 +261,7 @@ def response(self, inp: dict) -> dict: problem.add_parameter( "sigma", "likelihood", - prior=("uniform", {"low": std_low, "high": std_high}), + prior=Uniform(low=std_low, high=std_high), info="Std. dev. of zero-mean noise model", tex=r"$\sigma$", ) @@ -315,8 +269,6 @@ def response(self, inp: dict) -> dict: # add the forward model to the problem forward_model = SyntheticModel("ExpensiveModel") -# Add Expensive FE model to forward models -problem.add_forward_model(forward_model) # ========================================================================= # Create surrogate model @@ -327,25 +279,25 @@ def response(self, inp: dict) -> dict: surrogate_kwargs = { "training_max_iter": N_train_iter, "list_params": list_params, + "model_names": [name for name in forward_model.sensor_names], "show_progress": True, "silence_warnings": True, "fast_pred_var": True, } -surrogate_model = HarlowSurrogate( - name="FastModel", - forward_model=forward_model, - surrogate_model=ModelListGaussianProcess, - surrogate_kwargs=surrogate_kwargs, +# Generate surrogate class using model factory +surrogate_class = harlow_model_factory( + problem, forward_model, ModelListGaussianProcess, **surrogate_kwargs ) -problem.add_forward_model(surrogate_model) +# Initialize surrogate model +surrogate_model = surrogate_class( + name="FastModel", +) # ========================================================================= # Add test data to the inference problem # ========================================================================= - - def generate_data(): inp = {"X" + str(idx + 1): X_i for idx, X_i in enumerate(X_true)} sensors = forward_model(inp) @@ -354,36 +306,42 @@ def generate_data(): np.array(svals) + np.random.normal(loc=0.0, scale=std_true, size=Nt) ) sensors[isensor.name] = t_vec - problem.add_experiment( - "TestSeriesSurrogate", sensor_values=sensors, fwd_model_name="FastModel" - ) - problem.add_experiment( - "TestSeriesFull", sensor_values=sensors, fwd_model_name="ExpensiveModel" - ) + problem.add_experiment("TestSeriesFull", sensor_data=sensors) + problem.add_experiment("TestSeriesSurrogate", sensor_data=sensors) generate_data() +# Add Expensive FE model to forward models +problem.add_forward_model(forward_model, experiments="TestSeriesFull") +problem.add_forward_model(surrogate_model, experiments="TestSeriesSurrogate") + # ========================================================================= -# Create training data +# Train surrogate model # ========================================================================= -sampler = LatinHypercubeSampler(problem=problem) -test_samples = sampler.generate_samples(N_train) -train_samples, train_data = sampler.generate_training_data(forward_model, N_train) -train_X = train_samples.to_numpy()[:, :-1] -train_y = train_data["TestSeriesFull"] - -# Reshape `train_X` and `train_y -t_arr = np.tile(t_vec, N_train).reshape(-1, 1) -train_X_arr = np.repeat(train_X, Nt, axis=0) -train_X_arr = np.hstack((train_X_arr, t_arr)) - -train_y_arr = np.transpose(train_y, axes=[0, 2, 1]) -train_y_arr = train_y_arr.reshape(N_train * Nt, Nx) - -# Train surrogate -model_names = sensor_names -surrogate_model.fit(train_X_arr, train_y_arr) + +# Initialize sampler +sampler = FuzzyLolaVoronoi( + target_function=surrogate_model.target, + surrogate_model=surrogate_model.model, + domain_lower_bound=surrogate_model.lower_bound, + domain_upper_bound=surrogate_model.upper_bound, + test_points_x=test_X, + test_points_y=test_y, + evaluation_metric=rmse, +) +# main_start = time.time() +sampler.sample( + n_iter=n_iter, + n_initial_point=n_initial_point, + n_new_point_per_iteration=n_new_points_per_iteration, + ignore_old_neighborhoods=False, + ignore_far_neighborhoods=False, +) + +# # Train surrogate +# model_names = sensor_names +# surrogate_model.fit(train_X_arr, train_y_arr) # ==================================================================== From 24321af45b9acc224a4092741b3b7970c5df3e6f Mon Sep 17 00:00:00 2001 From: Koune Date: Thu, 1 Sep 2022 10:51:00 +0200 Subject: [PATCH 11/15] worked on connecting probeye and harlow --- examples/example_surrogate_hartman6D.py | 57 ++++++++++++++++--------- 1 file changed, 37 insertions(+), 20 deletions(-) diff --git a/examples/example_surrogate_hartman6D.py b/examples/example_surrogate_hartman6D.py index 9f1edd38..20aae98c 100644 --- a/examples/example_surrogate_hartman6D.py +++ b/examples/example_surrogate_hartman6D.py @@ -55,6 +55,11 @@ n_init_steps = 200 n_walkers = 20 +# Sampler settings +n_init = 20 +n_iter = 20 +n_point_per_iter = 2 + # Surrogate settings N_train = 50 N_train_iter = 200 @@ -195,7 +200,7 @@ def initialize_forward_model(self): def _get_bounds(self): # Find the priors corresponding to model parameters - for key, value in self.priors.items(): + for key, value in self.problem.priors.items(): # Error if any of the model priors are not uniform if value.prior_type != "uniform": @@ -205,24 +210,37 @@ def _get_bounds(self): f" only supports uniform priors." ) - for param in self.problem.model_prms: - lb = self.const_prms_dict["low_" + param] - ub = self.const_prms_dict["high_" + param] + # Dictionary of upper and lower bounds + for param in self.parameters: + lb = self.problem.constant_prms_dict["low_" + param] + ub = self.problem.constant_prms_dict["high_" + param] + + # List of upper and lower bounds + self.lower_bounds.append(lb) + self.upper_bounds.append(ub) + + # Append to dict self.bounds[param] = {"low": lb, "high": ub} def _cast_to_numpy(self, inp: dict) -> np.ndarray: raise NotImplementedError def _cast_to_dict(self, X: np.ndarray) -> dict: - raise NotImplementedError + return dict(zip(self.parameters, X)) def target(self, X: np.ndarray) -> np.ndarray: - inp = self._cast_to_dict(X) - response = self.forward_model.response(inp) - y = [[]] * len(self.output_sensors) - for idx, os in enumerate(self.output_sensors): - y[idx] = response[os.name] - return np.array(y).T + + # Create empty array of shape `(n_points, n_output_sensors)` + # y = np.zeros((X.shape[0], len(self.output_sensors))) + y = [[[]] * X.shape[0]] * len(self.output_sensors) + + # Evaluate physical model for each input vector + for idx_Xi, Xi in enumerate(X): + response = self.forward_model.response(self._cast_to_dict(Xi)) + for idx_os, os in enumerate(self.output_sensors): + y[idx_os][idx_Xi] = response[os.name] + + return np.array(y) def response(self, inp: dict) -> dict: @@ -324,19 +342,18 @@ def generate_data(): sampler = FuzzyLolaVoronoi( target_function=surrogate_model.target, surrogate_model=surrogate_model.model, - domain_lower_bound=surrogate_model.lower_bound, - domain_upper_bound=surrogate_model.upper_bound, - test_points_x=test_X, - test_points_y=test_y, - evaluation_metric=rmse, + domain_lower_bound=surrogate_model.lower_bounds, + domain_upper_bound=surrogate_model.upper_bounds, + run_name="test", + # test_points_x=test_X, + # test_points_y=test_y, + # evaluation_metric=rmse, ) # main_start = time.time() sampler.sample( n_iter=n_iter, - n_initial_point=n_initial_point, - n_new_point_per_iteration=n_new_points_per_iteration, - ignore_old_neighborhoods=False, - ignore_far_neighborhoods=False, + n_initial_point=n_init, + n_new_point_per_iteration=n_point_per_iter, ) # # Train surrogate From fed4207794e099fb4d13ac34ae8e66d016748e6a Mon Sep 17 00:00:00 2001 From: Koune Date: Thu, 1 Sep 2022 10:55:52 +0200 Subject: [PATCH 12/15] added .dat and .owl files to gitignore --- .gitignore | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.gitignore b/.gitignore index 5dcf4947..e8a67811 100644 --- a/.gitignore +++ b/.gitignore @@ -140,3 +140,5 @@ catalog-*.xml *.ttl *.sqlite3 *.sqlite3-journal +*.owl +*.dat From 6a7dd740a46225816bfa3c37ad7c90ac4d850d42 Mon Sep 17 00:00:00 2001 From: Koune Date: Thu, 1 Sep 2022 10:57:41 +0200 Subject: [PATCH 13/15] temporary commit to add .dat and .owl files to gitignore --- .gitignore | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.gitignore b/.gitignore index 5dcf4947..a2f54f9a 100644 --- a/.gitignore +++ b/.gitignore @@ -140,3 +140,5 @@ catalog-*.xml *.ttl *.sqlite3 *.sqlite3-journal +*.dat +*.owl \ No newline at end of file From b6c58166c4b9bba97330598a776feaf69f8fae16 Mon Sep 17 00:00:00 2001 From: Koune Date: Mon, 5 Sep 2022 12:16:37 +0200 Subject: [PATCH 14/15] implemented HarlowModelFactory class for generating surrogate forward models --- examples/example_surrogate_hartman6D.py | 343 +++++++----------------- probeye/definition/surrogate_model.py | 299 +++++++++++++++++++++ 2 files changed, 399 insertions(+), 243 deletions(-) diff --git a/examples/example_surrogate_hartman6D.py b/examples/example_surrogate_hartman6D.py index 20aae98c..884b2c54 100644 --- a/examples/example_surrogate_hartman6D.py +++ b/examples/example_surrogate_hartman6D.py @@ -9,9 +9,9 @@ adding I.i.d. Gaussian noise to samples from this function. Notes: - * torch, botorch and harlow are required to run this example - * the sampling and surrogating package `harlow` will be made - public soon + * A specific version of `harlow` is required to run this example, which can + be found here: + https://github.com/TNO/harlow/tree/implement-data-container """ # ========================================================================= @@ -28,17 +28,25 @@ from probeye.definition.likelihood_model import GaussianLikelihoodModel from probeye.inference.emcee.solver import EmceeSolver from probeye.definition.distribution import Uniform +from probeye.definition.surrogate_model import HarlowModelFactory # local imports (inference data post-processing) from probeye.postprocessing.sampling_plots import create_pair_plot from probeye.postprocessing.sampling_plots import create_posterior_plot # Surrogate model imports -from harlow.sampling import Sampler from harlow.sampling import LatinHypercube, FuzzyLolaVoronoi -from harlow.surrogating import Surrogate -from harlow.surrogating import ModelListGaussianProcess -from harlow.utils.transforms import ExpandDims +from harlow.surrogating import ( + Surrogate, + ModelListGaussianProcess, + VanillaGaussianProcess, +) +from harlow.utils.transforms import ( + ExpandDims, + TensorTransform, + ChainTransform, + Identity, +) import torch from botorch.test_functions.synthetic import Hartmann @@ -56,13 +64,12 @@ n_walkers = 20 # Sampler settings -n_init = 20 -n_iter = 20 +n_init = 250 +n_iter = 0 n_point_per_iter = 2 # Surrogate settings -N_train = 50 -N_train_iter = 200 +N_train_iter = 50 show_progress = True # ========================================================================= @@ -86,8 +93,8 @@ # ========================================================================= # Number of sensors and number of points in timeseries -Nx = 1 -Nt = 10 +Nx = 3 +Nt = 5 # Sensor names and positions sensor_names = ["S" + str(i + 1) for i in range(Nx)] @@ -130,136 +137,6 @@ def response(self, inp: dict) -> dict: return response -# ========================================================================= -# Define surrogate model -# ========================================================================= - - -def harlow_model_factory( - problem: InverseProblem, - forward_model: ForwardModelBase, - surrogate_model: Surrogate, - **kwargs, -) -> ForwardModelBase: - class HarlowSurrogate(ForwardModelBase): - """ - The inheritance from ExpensiveModel 'copies' the interface-method from - ExpensiveModel (the surrogate model should have the same interface as the - forward model). - """ - - def __init__( - self, - name, - ): - - self.name = name - self.kwargs = kwargs - self.lower_bounds = [] - self.upper_bounds = [] - self.bounds = {} - - self.forward_model = forward_model - self.problem = problem - self._surrogate = surrogate_model - - self.model_names = [] - for idx_x, os in enumerate(forward_model.output_sensors): - self.model_names.append(os.name) - self.kwargs["model_names"] = self.model_names - - self.num_features = len(forward_model.output_sensors) - self.kwargs["num_features"] = self.num_features - - # Check inputs - if "input_transform" not in self.kwargs.keys(): - self.input_transform = None - else: - self.input_transform = self.kwargs["input_transform"] - - if "output_transform" not in self.kwargs.keys(): - self.output_transform = None - else: - self.output_transform = self.kwargs["output_transform"] - - # Initialize surrogate - self.model = self._surrogate(**kwargs) - self.initialize_forward_model() - - # Get interface from existing forward model - def interface(self): - self.input_sensors = self.forward_model.input_sensors - self.output_sensors = self.forward_model.output_sensors - self.parameters = self.forward_model.parameters - - def initialize_forward_model(self): - super().__init__(self.name) - self._get_bounds() - - # TODO: This will likely fail for multidimensional parameters. FIX! - def _get_bounds(self): - - # Find the priors corresponding to model parameters - for key, value in self.problem.priors.items(): - - # Error if any of the model priors are not uniform - if value.prior_type != "uniform": - raise ValueError( - f"Non-uniform distribution of type {value.prior_type}" - f" found for prior {key}. The `HarlowSurrogate` currently" - f" only supports uniform priors." - ) - - # Dictionary of upper and lower bounds - for param in self.parameters: - lb = self.problem.constant_prms_dict["low_" + param] - ub = self.problem.constant_prms_dict["high_" + param] - - # List of upper and lower bounds - self.lower_bounds.append(lb) - self.upper_bounds.append(ub) - - # Append to dict - self.bounds[param] = {"low": lb, "high": ub} - - def _cast_to_numpy(self, inp: dict) -> np.ndarray: - raise NotImplementedError - - def _cast_to_dict(self, X: np.ndarray) -> dict: - return dict(zip(self.parameters, X)) - - def target(self, X: np.ndarray) -> np.ndarray: - - # Create empty array of shape `(n_points, n_output_sensors)` - # y = np.zeros((X.shape[0], len(self.output_sensors))) - y = [[[]] * X.shape[0]] * len(self.output_sensors) - - # Evaluate physical model for each input vector - for idx_Xi, Xi in enumerate(X): - response = self.forward_model.response(self._cast_to_dict(Xi)) - for idx_os, os in enumerate(self.output_sensors): - y[idx_os][idx_Xi] = response[os.name] - - return np.array(y) - - def response(self, inp: dict) -> dict: - - params = np.tile([inp["X" + str(i + 1)] for i in range(4)], (Nt, 1)) - X = np.hstack((params, t_vec.reshape(-1, 1))) - - response = dict() - - # Evaluate function and arange output on grid - f = self.model.predict(X, return_std=False) - - for idx_x, os in enumerate(self.output_sensors): - response[os.name] = f[idx_x, :] - - return response - - return HarlowSurrogate - - # ========================================================================= # Define inference problem # ========================================================================= @@ -287,32 +164,6 @@ def response(self, inp: dict) -> dict: # add the forward model to the problem forward_model = SyntheticModel("ExpensiveModel") - -# ========================================================================= -# Create surrogate model -# ========================================================================= - -list_params = [[0, 1, 2, 3, 4]] * len(sensor_names) - -surrogate_kwargs = { - "training_max_iter": N_train_iter, - "list_params": list_params, - "model_names": [name for name in forward_model.sensor_names], - "show_progress": True, - "silence_warnings": True, - "fast_pred_var": True, -} - -# Generate surrogate class using model factory -surrogate_class = harlow_model_factory( - problem, forward_model, ModelListGaussianProcess, **surrogate_kwargs -) - -# Initialize surrogate model -surrogate_model = surrogate_class( - name="FastModel", -) - # ========================================================================= # Add test data to the inference problem # ========================================================================= @@ -323,114 +174,120 @@ def generate_data(): sensors[sname] = list( np.array(svals) + np.random.normal(loc=0.0, scale=std_true, size=Nt) ) - sensors[isensor.name] = t_vec + + # To avoid errors when `Nt = 1` + if len(t_vec) == 1: + sensors[isensor.name] = float(t_vec[0]) + else: + sensors[isensor.name] = t_vec + + # Add experiments problem.add_experiment("TestSeriesFull", sensor_data=sensors) problem.add_experiment("TestSeriesSurrogate", sensor_data=sensors) +# Generate data and add expensive forward model generate_data() - -# Add Expensive FE model to forward models problem.add_forward_model(forward_model, experiments="TestSeriesFull") -problem.add_forward_model(surrogate_model, experiments="TestSeriesSurrogate") # ========================================================================= -# Train surrogate model +# Create surrogate model # ========================================================================= +n_params = 4 +list_params = [[i for i in range(n_params)]] * len(sensor_names) * len(t_vec) -# Initialize sampler -sampler = FuzzyLolaVoronoi( - target_function=surrogate_model.target, - surrogate_model=surrogate_model.model, - domain_lower_bound=surrogate_model.lower_bounds, - domain_upper_bound=surrogate_model.upper_bounds, - run_name="test", - # test_points_x=test_X, - # test_points_y=test_y, - # evaluation_metric=rmse, +surrogate_kwargs = { + "training_max_iter": N_train_iter, + "list_params": list_params, + "show_progress": True, + "silence_warnings": True, + "fast_pred_var": True, + "input_transform": Identity, + "output_transform": Identity, +} + +# Generate surrogate class using model factory +model_factory = HarlowModelFactory( + problem, + forward_model, + FuzzyLolaVoronoi, + VanillaGaussianProcess, + **surrogate_kwargs, ) -# main_start = time.time() -sampler.sample( + +# Sample +model_factory.sample( n_iter=n_iter, n_initial_point=n_init, n_new_point_per_iteration=n_point_per_iter, ) -# # Train surrogate -# model_names = sensor_names -# surrogate_model.fit(train_X_arr, train_y_arr) +# Get forward model instance +forward_surrogate_model = model_factory.get_harlow_model("FastModel") + +# ========================================================================= +# Add forward models +# ========================================================================= + +# Add surrogate model to forward models +problem.add_forward_model(forward_surrogate_model, experiments="TestSeriesSurrogate") +# add the likelihood models to the problem +for osensor in osensor_list: + problem.add_likelihood_model( + GaussianLikelihoodModel( + experiment_name="TestSeriesSurrogate", + model_error="additive", + ) + ) # ==================================================================== # Plot surrogate vs. FE model prediction # ==================================================================== # Physical model prediction +arr_x = np.linspace(X_low, X_high, 100) inp = {"X" + str(idx + 1): X_i for idx, X_i in enumerate(X_true)} -y_true = forward_model.response(inp) -y_pred = surrogate_model.response(inp) -upper = list(surrogate_model.model.cr_u.squeeze().detach().numpy().T) -lower = list(surrogate_model.model.cr_l.squeeze().detach().numpy().T) -# Initialize plots +# Initialize zero arrays to store results +y_true = np.zeros((len(arr_x), len(sensor_names) * len(t_vec))) +y_pred = np.zeros((len(arr_x), len(sensor_names) * len(t_vec))) + +# Evaluate physical model for each input vector +for idx_xi, xi in enumerate(arr_x): + inp["X1"] = xi + res_true = forward_model.response(inp) + res_pred = forward_surrogate_model.response(inp) + sensor_out_true = [] + sensor_out_pred = [] + for idx_os, os in enumerate(osensor_list): + sensor_out_true.append(res_true[os.name]) + sensor_out_pred.append(res_pred[os.name]) + y_true[idx_xi, :] = np.ravel(sensor_out_true) + y_pred[idx_xi, :] = np.ravel(sensor_out_pred) + +# Plot nrows = 3 -ncols = int(np.ceil(len(sensor_names) / 3)) - -# Initialize plot +ncols = int(np.ceil(len(sensor_names) * len(t_vec) / 3)) f, axes = plt.subplots(nrows, ncols, sharex=True, figsize=(3 * ncols, 3 * nrows)) - -for j in range(len(sensor_names)): - +for j in range(len(sensor_names) * len(t_vec)): ax_i = axes.ravel()[j] grid_idx = np.unravel_index(j, (nrows, ncols)) - - # Plot training data as black stars - ax_i.plot(t_vec, train_y[:, j, :].T, "k*", label="Observations") - - # Predictive mean as blue line - ax_i.plot(t_vec, y_pred[sensor_names[j]], "b", label="Mean") - - # Shade in confidence - ax_i.fill_between( - t_vec, - lower[j], - upper[j], - alpha=0.5, - label="Confidence", - ) - ax_i.plot( - t_vec, - y_true[sensor_names[j]], - color="red", - linestyle="dashed", - label="Model", + ax_i.plot(arr_x, y_true[:, j], color="blue", label="True") + ax_i.plot(arr_x, y_pred[:, j], color="red", linestyle="dashed", label="Surrogate") + ax_i.set_title( + f"Sensor: {str(sensor_names[grid_idx[0]]) + '_' + str(t_vec[grid_idx[1]])}" ) - ax_i.set_title(f"Sensor: {sensor_names[j]}") axes = np.atleast_2d(axes) -# axes[0, 0].legend() -[ax_i.set_xlabel(r"$X_6$") for ax_i in axes[-1, :]] -[ax_i.set_ylabel("Stress [MPa]") for ax_i in axes[:, 0]] +[ax_i.set_xlabel(r"$X_1$") for ax_i in axes[-1, :]] +axes[0, 0].legend() plt.show() # ========================================================================= # Add noise models # ========================================================================= - -# add the likelihood models to the problem -for osensor in osensor_list: - problem.add_likelihood_model( - GaussianLikelihoodModel( - experiment_name="TestSeriesSurrogate", - prms_def=["sigma"], - model_error="additive", - additive_measurement_error=False, - correlation_variables=None, - ) - ) - - # Problem overview problem.info() @@ -446,13 +303,13 @@ def generate_data(): # Initialize and run solver # ========================================================================= -emcee_solver = EmceeSolver( +solver = EmceeSolver( problem, show_progress=True, ) -inference_data = emcee_solver.run_mcmc( +inference_data = solver.run( n_walkers=n_walkers, n_steps=n_steps, n_initial_steps=n_init_steps, vectorize=False ) @@ -461,7 +318,7 @@ def generate_data(): # ========================================================================= create_pair_plot( inference_data, - problem, + solver.problem, show=False, true_values=true_values, title="Joint posterior", @@ -469,7 +326,7 @@ def generate_data(): create_posterior_plot( inference_data, - problem, + solver.problem, show=False, true_values=true_values, title="Marginal posteriors", diff --git a/probeye/definition/surrogate_model.py b/probeye/definition/surrogate_model.py index 3717856b..eaa849c2 100644 --- a/probeye/definition/surrogate_model.py +++ b/probeye/definition/surrogate_model.py @@ -1,6 +1,15 @@ +import numpy as np +from typing import Callable + # local imports +from probeye.definition.inverse_problem import InverseProblem from probeye.definition.forward_model import ForwardModelBase +# Harlow imports +from harlow.sampling import Sampler +from harlow.surrogating import Surrogate +from harlow.utils.helper_functions import latin_hypercube_sampling + class SurrogateModelBase(ForwardModelBase): """ @@ -30,3 +39,293 @@ def fit(self): Prepares the surrogate model by approximating the forward model in some way. """ pass + + +class HarlowModelFactory: + """ + Model factory that returns a `ForwardModel` object, where the response + is obtained from a fitted surrogate model. + + Usage: + TODO + + Todo: + * Currently, each point is added as a separate surrogate model + + Notes: + * This is a utility class meant to facilitate the sampling and surrogating + of expensive forward models by bridging `probeye` with `harlow`. The + main objective is to reduce the manual effort required for surrogating + by eliminating the need to re-specify the problem when building the + surrogate. + + """ + + def __init__( + self, + problem: InverseProblem, + forward_model: ForwardModelBase, + sampler: Sampler, + surrogate_model: Surrogate, + fit_points_x: np.ndarray = None, + fit_points_y: np.ndarray = None, + test_points_x: np.ndarray = None, + test_points_y: np.ndarray = None, + evaluation_metric: Callable = None, + verbose: bool = False, + run_name: str = "test", + **surrogate_kwargs, + ): + + # Initialize + self.surrogate_kwargs = surrogate_kwargs + self.problem = problem + self.forward_model = forward_model + self.input_sensors = forward_model.input_sensors + self.output_sensors = forward_model.output_sensors + self.parameters = forward_model.parameters + self.fit_points_x = fit_points_x + self.fit_points_y = fit_points_y + self.test_points_x = test_points_x + self.test_points_y = test_points_y + self.evaluation_metric = evaluation_metric + self.verbose = verbose + self.run_name = run_name + + # Find the experiment associated with the expensive forward model + # and get sensor information + # TODO: Implement check that the expensive forward model + # has been added to the problem + self.input_sensor_data = dict() + for exp_name, exp in self.problem.experiments.items(): + if exp.forward_model == self.forward_model.name: + for inp_sensor in self.input_sensors: + self.input_sensor_data[inp_sensor.name] = list( + exp.sensor_data[inp_sensor.name] + ) + + # Check that the supplied forward model only has a single sensor + # TODO: Can a forward model have multiple input sensors? + if len(list(self.input_sensor_data.keys())) != 1: + raise ValueError( + f"The specified forward model must have " + f"exactly one input sensor but has " + f"{len(list(self.input_sensor_data.keys()))}" + ) + + # Extract sensor names and input/output shape information + # from forward model. + # TODO: This should be adjusted depending on the problem definition + # or specified by the user e.g. with keyword arguments that specify + # the structure of the surrogate model, e.g.: + # * `shape = None` to create an independent surrogate for each observation + # * `shape = Sensor` to create a surrogate model per sensor + self.vec_t = list(self.input_sensor_data.values())[0] + self.model_names = [] + for idx_x, os in enumerate(forward_model.output_sensors): + for idx_t, t in enumerate(self.vec_t): + self.model_names.append(os.name + "_" + str(t)) + self.surrogate_kwargs["model_names"] = self.model_names + self.num_features = len(forward_model.output_sensors) + self.surrogate_kwargs["num_features"] = self.num_features + + # Get bounds + self._get_bounds() + + # Initialize surrogate + self.surrogate = surrogate_model(**surrogate_kwargs) + self.func_pred = self.surrogate.predict + + # Initialize sampler + self.sampler = sampler( + target_function=self.target, + surrogate_model=self.surrogate, + domain_lower_bound=self.lower_bounds, + domain_upper_bound=self.upper_bounds, + fit_points_x=self.fit_points_x, + fit_points_y=self.fit_points_y, + test_points_x=self.test_points_x, + test_points_y=self.test_points_y, + evaluation_metric=self.evaluation_metric, + verbose=self.verbose, + run_name=self.run_name, + ) + + def fit(self, **kwargs): + """ + Fit the surrogate model using the sampling points and observations + stored in `self.sampler` + + Parameters + ---------- + kwargs + + Returns + ------- + + """ + + self.surrogate.fit( + self.sampler.fit_points_x, self.sampler.fit_points_y, **kwargs + ) + + def sample(self, **kwargs): + """ + Draw samples from the forward model and fit the surrogate model + using the specified sampler. + + Parameters + ---------- + kwargs + + Returns + ------- + + """ + + self.fit_points_x, self.fit_points_y = self.sampler.sample(**kwargs) + + # # Get prediction function + # self.func_pred = self.surrogate.predict + + def get_training_data(self, n_samples=1): + """ + Obtain training points `X` with latin hypercube sampling and + corresponding observations `y` by evaluating the target function. + + Returns + ------- + A tuple of training points and observations `(train_X, train_y)` + + """ + train_X = latin_hypercube_sampling( + n_sample=n_samples, + domain_lower_bound=np.array(self.lower_bounds), + domain_upper_bound=np.array(self.upper_bounds), + ) + train_y = self.target(train_X) + + return train_X, train_y + + def get_harlow_model( + self, + name: str, + ) -> ForwardModelBase: + + """ + Model factory for the `HarlowSurrogate` class. + + TODO: + * Refactor this function so that it returns a `ForwardModel` object + with a `.predict()` method that queries an instance of `Surrogate`. + * It is important that the generated `ForwardModel` instance does not + contain the surrogate model itself (due to incompatibility of tensors + with `copy.deepcopy`), or the problem definition. + + Notes + --------- + * This function was originally created as a workaround for passing arguments to + classes derived from `ForwardModelBase`. + + Parameters + ---------- + problem + forward_model + surrogate + n_output_dims + kwargs + + Returns + ------- + + """ + + class HarlowSurrogate(ForwardModelBase): + + # Get interface from existing forward model + def interface(self): + self.input_sensors = self.kwargs["input_sensors"] + self.output_sensors = self.kwargs["output_sensors"] + self.parameters = self.kwargs["parameters"] + self.input_coordinates = self.kwargs["input_coordinates"] + + def response(self, inp: dict) -> dict: + # TODO: Currently this method assumes that `Nx` and `Nt` are available. + # Instead they should be obtained from the problem definition + params = np.atleast_2d([inp["X" + str(i + 1)] for i in range(4)]) + response = dict() + + # TODO: Adjust `X` to account for the case of sensors with vector + # input (e.g. timeseries). + vec_t = self.input_coordinates[self.input_sensors[0].name] + + # Evaluate function and arange output on grid + f = self.kwargs["func_pred"](params, return_std=False) + f = np.reshape(f, (len(self.output_sensors), len(vec_t))) + + for idx_x, os in enumerate(self.output_sensors): + response[os.name] = f[idx_x, :] + + return response + + return HarlowSurrogate( + name, + func_pred=self.func_pred, + input_sensors=self.input_sensors, + output_sensors=self.output_sensors, + parameters=self.parameters, + input_coordinates=self.input_sensor_data, + ) + + # TODO: This will likely fail for multidimensional parameters. FIX! + def _get_bounds(self): + + self.lower_bounds = [] + self.upper_bounds = [] + self.bounds = {} + + # Find the priors corresponding to model parameters + for key, value in self.problem.priors.items(): + + # Error if any of the model priors are not uniform + if value.prior_type != "uniform": + raise ValueError( + f"Non-uniform distribution of type {value.prior_type}" + f" found for prior {key}. The `HarlowSurrogate` currently" + f" only supports uniform priors." + ) + + # Dictionary of upper and lower bounds + for param in self.parameters: + lb = self.problem.constant_prms_dict["low_" + param] + ub = self.problem.constant_prms_dict["high_" + param] + + # List of upper and lower bounds + self.lower_bounds.append(lb) + self.upper_bounds.append(ub) + + # Append to dict + self.bounds[param] = {"low": lb, "high": ub} + + def _cast_to_numpy(self, inp: dict) -> np.ndarray: + raise NotImplementedError + + def _cast_to_dict(self, X: np.ndarray) -> dict: + return dict(zip(self.parameters, X)) + + def target(self, X: np.ndarray) -> np.ndarray: + # TODO: Adjust `y` to account for the case of sensors with vector + # input (e.g. timeseries). + + # Initialize output array + y = np.zeros((X.shape[0], len(self.output_sensors) * len(self.vec_t))) + + # Evaluate physical model for each input vector + for idx_Xi, Xi in enumerate(X): + response = self.forward_model.response(self._cast_to_dict(Xi)) + sensor_out = [] + for idx_os, os in enumerate(self.output_sensors): + sensor_out.append(response[os.name]) + y[idx_Xi, :] = np.ravel(sensor_out) + + return np.array(y) From 3960eef18e64d1adfd99480317d50ccb6d6e21c8 Mon Sep 17 00:00:00 2001 From: Koune Date: Thu, 15 Sep 2022 18:29:34 +0200 Subject: [PATCH 15/15] refactored sampling and surrogating interface and adjusted example --- examples/example_surrogate_hartman6D.py | 49 +- probeye/definition/__init__.py | 1 - probeye/metamodeling/__init__.py | 0 .../initial_sampling.py | 0 .../sampling.py} | 206 +++++++-- probeye/metamodeling/surrogating.py | 430 ++++++++++++++++++ probeye/surrogate/__init__.py | 2 - .../test_surrogate_model.py | 2 +- .../surrogate/test_initial_sampling.py | 2 +- 9 files changed, 638 insertions(+), 54 deletions(-) create mode 100644 probeye/metamodeling/__init__.py rename probeye/{surrogate => metamodeling}/initial_sampling.py (100%) rename probeye/{definition/surrogate_model.py => metamodeling/sampling.py} (56%) create mode 100644 probeye/metamodeling/surrogating.py delete mode 100644 probeye/surrogate/__init__.py diff --git a/examples/example_surrogate_hartman6D.py b/examples/example_surrogate_hartman6D.py index 884b2c54..4f0d27d9 100644 --- a/examples/example_surrogate_hartman6D.py +++ b/examples/example_surrogate_hartman6D.py @@ -28,23 +28,19 @@ from probeye.definition.likelihood_model import GaussianLikelihoodModel from probeye.inference.emcee.solver import EmceeSolver from probeye.definition.distribution import Uniform -from probeye.definition.surrogate_model import HarlowModelFactory +from probeye.metamodeling.sampling import LatinHypercubeSampler, HarlowSampler +from probeye.metamodeling.surrogating import HarlowSurrogate # local imports (inference data post-processing) from probeye.postprocessing.sampling_plots import create_pair_plot from probeye.postprocessing.sampling_plots import create_posterior_plot # Surrogate model imports -from harlow.sampling import LatinHypercube, FuzzyLolaVoronoi +from harlow.sampling import FuzzyLolaVoronoi, ProbabilisticSampler, LatinHypercube from harlow.surrogating import ( - Surrogate, - ModelListGaussianProcess, VanillaGaussianProcess, ) from harlow.utils.transforms import ( - ExpandDims, - TensorTransform, - ChainTransform, Identity, ) @@ -64,9 +60,10 @@ n_walkers = 20 # Sampler settings -n_init = 250 -n_iter = 0 -n_point_per_iter = 2 +n_init = 100 +n_iter = 5 +n_point_per_iter = 20 +stopping_criterium = -np.inf # Surrogate settings N_train_iter = 50 @@ -196,6 +193,7 @@ def generate_data(): n_params = 4 list_params = [[i for i in range(n_params)]] * len(sensor_names) * len(t_vec) +# Kwargs to be passed to the surrogate model surrogate_kwargs = { "training_max_iter": N_train_iter, "list_params": list_params, @@ -206,24 +204,29 @@ def generate_data(): "output_transform": Identity, } -# Generate surrogate class using model factory -model_factory = HarlowModelFactory( - problem, - forward_model, - FuzzyLolaVoronoi, - VanillaGaussianProcess, - **surrogate_kwargs, -) +# Define the surrogate model +surrogate_model = VanillaGaussianProcess(**surrogate_kwargs) -# Sample -model_factory.sample( - n_iter=n_iter, +# Probeye's latin hypercube sampler +lhs_sampler = LatinHypercubeSampler(problem) + +# An iterative sampler. Here we pass the surrogate ForwardModel directly to the sampler. However, it is +# also possible to pass a surrogate model that will be included in a forward model after fitting. +harlow_sampler = HarlowSampler(problem, forward_model, LatinHypercube, surrogate_model) + +# Sampler and fit +harlow_sampler.sample( n_initial_point=n_init, n_new_point_per_iteration=n_point_per_iter, + n_iter=n_iter, + stopping_criterium=stopping_criterium, ) +harlow_sampler.fit() -# Get forward model instance -forward_surrogate_model = model_factory.get_harlow_model("FastModel") +# Define the surrogate forward model +forward_surrogate_model = HarlowSurrogate( + "SurrogateModel", surrogate_model, forward_model +) # ========================================================================= # Add forward models diff --git a/probeye/definition/__init__.py b/probeye/definition/__init__.py index 8033c9d2..59009555 100644 --- a/probeye/definition/__init__.py +++ b/probeye/definition/__init__.py @@ -2,7 +2,6 @@ from probeye.definition import inverse_problem from probeye.definition import inference_problem from probeye.definition import forward_model -from probeye.definition import surrogate_model from probeye.definition import likelihood_model from probeye.definition import distribution from probeye.definition import parameter diff --git a/probeye/metamodeling/__init__.py b/probeye/metamodeling/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/probeye/surrogate/initial_sampling.py b/probeye/metamodeling/initial_sampling.py similarity index 100% rename from probeye/surrogate/initial_sampling.py rename to probeye/metamodeling/initial_sampling.py diff --git a/probeye/definition/surrogate_model.py b/probeye/metamodeling/sampling.py similarity index 56% rename from probeye/definition/surrogate_model.py rename to probeye/metamodeling/sampling.py index eaa849c2..7924a80d 100644 --- a/probeye/definition/surrogate_model.py +++ b/probeye/metamodeling/sampling.py @@ -1,47 +1,199 @@ +# standard library +from typing import Tuple, Callable +import copy + +# third party imports import numpy as np -from typing import Callable +import pandas as pd +from scipy.stats import qmc # local imports from probeye.definition.inverse_problem import InverseProblem from probeye.definition.forward_model import ForwardModelBase +from probeye.inference.scipy.priors import translate_prior +from probeye.subroutines import len_or_one -# Harlow imports -from harlow.sampling import Sampler -from harlow.surrogating import Surrogate +# external imports +from harlow.sampling import Sampler as HarlowSamplerBase +from harlow.surrogating import Surrogate as HarlowSurrogateBase from harlow.utils.helper_functions import latin_hypercube_sampling -class SurrogateModelBase(ForwardModelBase): +class SamplerBase: + """ + Base class for probeye samplers + """ + + def __init__(self, problem: InverseProblem): + # the considered inverse problem + self.problem = problem + + # the sampling happens before the problem is given to a solver; this means that + # the priors of the problem are merely descriptive and they have to be + # translated to have their respective computing routines + self.priors = copy.deepcopy(self.problem.priors) + for prior_template in self.problem.priors.values(): + prm_name = prior_template.ref_prm + self.priors[prm_name] = translate_prior(prior_template) + + def sample( + self, forward_model: ForwardModelBase, **kwargs + ) -> Tuple[pd.DataFrame, dict]: + raise NotImplementedError + + +class LatinHypercubeSampler(SamplerBase): """ - Base class for a surrogate model, i.e., a forward model that approximates another - (typically computationally more expensive) forward model. + Contains functionalities to provide samples of an inverse problem's parameters by + utilizing latin hypercube sampling. The generate samples are intended for training + of a surrogate model. Parameters ---------- - name - The name of the surrogate model. Must be unique among all surrogate model's - names within a considered InverseProblem. - forward_model - The forward model object that the surrogate model approximates. + problem + The considered inverse problem. """ - def __init__(self, name: str, forward_model: ForwardModelBase): + def __init__(self, problem): + super().__init__(problem) + + def generate_samples(self, n_samples: int, seed: int = 1) -> np.ndarray: + """ + Samples the problems latent parameter vector from the parameter's priors in + combination with latin hypercube sampling. - super().__init__(name) + Parameters + ---------- + n_samples + The number of requested samples. + seed + Seed for random number generator. - # the surrogate model has access to the forward model it approximates; this - # forward model can be called (evaluating its response) from the surrogate model - # via self.forward_model.response(inp) where inp contains the input dictionary - self.forward_model = forward_model + Returns + ------- + sample_array + The sampled latent parameter vectors. Each row corresponds to a single + latent parameter vector. + """ + + # make sure that all parameters are one-dimensional; it is not straight forward + # how to do LHS for general multivariate parameters + for prm_name in self.problem.latent_prms: + if self.problem.parameters[prm_name].dim > 1: + raise RuntimeError( + f"The given problem has a multivariate parameter ('{prm_name}') " + f"for which LHS is not supported." + ) + + # the dimension of the required hypercube is the number of parameters; note + # that this only holds since all parameters are 1D (see the check above) + dim = self.problem.n_latent_prms - def fit(self): + # create the requested number of latin hypercube samples with the requested dim. + lhd = qmc.LatinHypercube(d=dim, seed=seed).random(n=n_samples) + + # this is going to be the array of the parameter samples; each row will + # correspond to a theta-vector + sample_array = np.zeros(lhd.shape) + + # fill the prepared array for the parameter samples using the LHS samples + for prm_name, parameter in self.problem.parameters.items(): + if parameter.is_latent: + idx = parameter.index + prior = self.priors[prm_name] + prms = self.problem.get_constants(prior.hyperparameters) + for lt_prm in self.problem.get_latent_prior_hyperparameters(prm_name): + idx_lt = self.problem.parameters[lt_prm].index + prms[lt_prm] = sample_array[:, idx_lt] + q = lhd[:, idx] + sample_array[:, idx] = prior(prms, "ppf", q, use_ref_prm=False) + + return sample_array + + def sample( + self, + forward_model: ForwardModelBase, + n_samples: int, + seed: int = 1, + ) -> Tuple[pd.DataFrame, dict]: """ - Prepares the surrogate model by approximating the forward model in some way. + Generates a given number of training data for fitting a surrogate model. The + training data contains a number of parameter vectors (sampled using LHS) and + the corresponding model responses. + + Parameters + ---------- + forward_model + The forward model that should be evaluated. + n_samples + The number of parameter vectors the forward model should be evaluated for. + seed + Seed for random number generator. + + Returns + ------- + prm_samples_pd + The parameter samples the forward model was evaluated at. + responses_over_experiments + The keys are the names of the experiment the forward model is associated + with, while the values are 3D-arrays containing the forward model's + response. responses_over_experiments[i][j] will contain the forward model's + response with the ith parameter-vector for the jth output sensor. """ - pass + # get the forward model object with the given name and prepare the corresponding + # experimental in- and output dictionaries + forward_model.prepare_experimental_inputs_and_outputs() + + # generate the latent parameter samples and convert it to a data frame to have + # the association between columns and parameter names + prm_samples = self.generate_samples(n_samples, seed=seed) + prm_samples_pd = pd.DataFrame( + prm_samples, columns=self.problem.get_theta_names() + ) + + # this dictionary will contain the forward model responses for each of the + # experiments associated with the forward model; so the keys will be experiment + # names while the values will be 3D-arrays with the forward model's responses; + # responses_over_experiments[i] will correspond to the response of the ith + # parameter vector; responses_over_experiments[i][j] will contain the forward + # model's response with the ith parameter vector for the jth output sensor + responses_over_experiments = {} + + # here, the length of the vectors of the forward model's output sensors is + # determined; to that end, the forward model is evaluated once + first_exp_name = forward_model.experiment_names[0] + exp_inp = forward_model.input_from_experiments[first_exp_name] + first_theta = prm_samples[0] + prms_model = self.problem.get_parameters(first_theta, forward_model.prms_def) + inp = {**exp_inp, **prms_model} + response_dict = forward_model.response(inp) + # make sure that the vectors returned by each of the forward model's output + # sensors has the same length; otherwise an AssertionError is raised + length_set = set() + for value in response_dict.values(): + length_set.add(len_or_one(value)) + assert len(length_set) == 1 + n_out_values = list(length_set)[0] + + # evaluate the forward model for each experiment/parameter vector + for exp_name in forward_model.experiment_names: + exp_inp = forward_model.input_from_experiments[exp_name] + response_array = np.zeros( + (n_samples, forward_model.n_output_sensors, n_out_values) + ) + for i, theta in enumerate(prm_samples): + prms_model = self.problem.get_parameters(theta, forward_model.prms_def) + inp = {**exp_inp, **prms_model} # adds the two dictionaries + response_dict = forward_model.response(inp) + for j, response_vector in enumerate(response_dict.values()): + response_array[i, j, :] = response_vector + responses_over_experiments[exp_name] = response_array + + return prm_samples_pd, responses_over_experiments -class HarlowModelFactory: + +class HarlowSampler(SamplerBase): """ Model factory that returns a `ForwardModel` object, where the response is obtained from a fitted surrogate model. @@ -65,8 +217,8 @@ def __init__( self, problem: InverseProblem, forward_model: ForwardModelBase, - sampler: Sampler, - surrogate_model: Surrogate, + sampler: HarlowSamplerBase, + surrogate_model: HarlowSurrogateBase, fit_points_x: np.ndarray = None, fit_points_y: np.ndarray = None, test_points_x: np.ndarray = None, @@ -77,6 +229,8 @@ def __init__( **surrogate_kwargs, ): + super().__init__(problem) + # Initialize self.surrogate_kwargs = surrogate_kwargs self.problem = problem @@ -132,8 +286,8 @@ def __init__( # Get bounds self._get_bounds() - # Initialize surrogate - self.surrogate = surrogate_model(**surrogate_kwargs) + # Surrogate model + self.surrogate = surrogate_model self.func_pred = self.surrogate.predict # Initialize sampler diff --git a/probeye/metamodeling/surrogating.py b/probeye/metamodeling/surrogating.py new file mode 100644 index 00000000..01891396 --- /dev/null +++ b/probeye/metamodeling/surrogating.py @@ -0,0 +1,430 @@ +""" +Surrogating base class and specific implementations + +Notes: + * There is likely no need for the surrogate to have access to any of the problem + internals, since the sampling will be handled by the sampler class. + * This allows us to have a very simple surrogate base class by only implementing + `fit`, `predict` (optional?) and `response`. + * An existing forward model must be provided as an input argument. The surrogate + will copy the interface method from that forward model. + * The base class is used similar to `ForwardModelBase`. The user is expected to derive + their own class from it. +""" + +# standard library +from copy import deepcopy + +# external +import numpy as np + +# local imports +from probeye.definition.forward_model import ForwardModelBase + +# harlow imports +from harlow.surrogating import Surrogate as HarlowSurrogateBase + + +class SurrogateModelBase(ForwardModelBase): + """ + Base class for a surrogate model, i.e., a forward model that approximates another + (typically computationally more expensive) forward model. The surrogate does not + need to have access to the problem internals. Therefore, the surrogate base class + simply extends `ForwardModelBase` with a `fit` method. + + Notes: + * The `forward_model` input is needed to obtain the interface. Surrogate models + will copy the interface of the associated forward model and override the + `response` method. + + TODO: + * Check out probeye's existing machinery for copying elements of the forward + model definition. Use that instead of `deepcopy`. + + Parameters + ---------- + + """ + + def __init__( + self, name: str, surrogate: HarlowSurrogateBase, forward_model: ForwardModelBase + ): + super().__init__(name, surrogate=surrogate, forward_model=forward_model) + + def interface(self): + self.parameters = deepcopy(self.kwargs["forward_model"].parameters) + self.input_sensors = deepcopy(self.kwargs["forward_model"].input_sensors) + self.output_sensors = deepcopy(self.kwargs["forward_model"].output_sensors) + + def fit(self, train_X: np.ndarray, train_y: np.ndarray): + """ + Prepares the surrogate model by approximating the forward model in some way. + """ + raise NotImplementedError + + def predict(self, X: np.ndarray): + """ + Used by samplers that require evaluation of the surrogate model. + """ + raise NotImplementedError + + def response(self, inp: dict) -> dict: + """ + Overrides the `response` method of the forward model. + """ + raise NotImplementedError + + +class HarlowSurrogate(SurrogateModelBase): + """ """ + + def __init__( + self, + name: str, + surrogate: HarlowSurrogateBase, + forward_model: ForwardModelBase, + **kwargs + ): + + # Initialize + self.forward_model = forward_model + self.surrogate = surrogate + self.func_pred = self.surrogate.predict + self.kwargs = kwargs + + super().__init__(name, surrogate=surrogate, forward_model=forward_model) + + def fit(self, train_X: np.ndarray, train_y: np.ndarray, **kwargs): + """ + Fit the surrogate model using the sampling points and observations + stored in `self.sampler` + """ + + self.surrogate.fit(self, train_X, train_y, **kwargs) + + def predict(self, X: np.ndarray, **kwargs) -> np.ndarray: + return self.surrogate.predict(X, **kwargs) + + def response(self, inp: dict) -> dict: + + # # TODO: + # # * Adjust `X` to account for the case of sensors with vector input (e.g. timeseries). + # # * Obtain input coordinates from the interface + # input_sensor_name = self.input_sensors[0].name + # N_x = len(self.output_sensors) + # N_t = len(inp[input_sensor_name]) + + # Evaluate function and arange output on grid + f = self.surrogate.predict(self.cast_to_array(inp)) + f = np.reshape(f, (len(self.output_sensors), -1)) + + response = dict() + for idx_x, os in enumerate(self.output_sensors): + response[os.name] = f[idx_x, :] + + return response + + def cast_to_array(self, inp: dict) -> np.ndarray: + """ + Convert dictionary of parameter: value pairs to a 1D numpy array in the order that + they are provided in `inp`. + """ + return np.atleast_2d([inp[param] for param in self.parameters]) + + def cast_to_dict(self, X: np.ndarray) -> dict: + """ + Convert 1D numpy array of inputs to the dictionary input expected by `response`. The order + of elements in `X` is assumed to be the order of parameters in `self.parameters`. + """ + return dict(zip(self.parameters, X)) + + +# # TODO: Old code to be removed +# class SurrogateModelFactory: +# """ +# Model factory that returns a `ForwardModel` object, where the response +# is obtained from a fitted surrogate model. +# +# Usage: +# TODO +# +# Todo: +# * Currently, each point is added as a separate surrogate model +# +# Notes: +# * This is a utility class meant to facilitate the sampling and surrogating +# of expensive forward models by bridging `probeye` with `harlow`. The +# main objective is to reduce the manual effort required for surrogating +# by eliminating the need to re-specify the problem when building the +# surrogate. +# +# """ +# +# def __init__( +# self, +# problem: InverseProblem, +# forward_model: ForwardModelBase, +# sampler: Sampler, +# surrogate_model: Surrogate, +# fit_points_x: np.ndarray = None, +# fit_points_y: np.ndarray = None, +# test_points_x: np.ndarray = None, +# test_points_y: np.ndarray = None, +# evaluation_metric: Callable = None, +# verbose: bool = False, +# run_name: str = "test", +# **surrogate_kwargs, +# ): +# +# # Initialize +# self.surrogate_kwargs = surrogate_kwargs +# self.problem = problem +# self.forward_model = forward_model +# self.input_sensors = forward_model.input_sensors +# self.output_sensors = forward_model.output_sensors +# self.parameters = forward_model.parameters +# self.fit_points_x = fit_points_x +# self.fit_points_y = fit_points_y +# self.test_points_x = test_points_x +# self.test_points_y = test_points_y +# self.evaluation_metric = evaluation_metric +# self.verbose = verbose +# self.run_name = run_name +# +# # Find the experiment associated with the expensive forward model +# # and get sensor information +# # TODO: Implement check that the expensive forward model +# # has been added to the problem +# self.input_sensor_data = dict() +# for exp_name, exp in self.problem.experiments.items(): +# if exp.forward_model == self.forward_model.name: +# for inp_sensor in self.input_sensors: +# self.input_sensor_data[inp_sensor.name] = list( +# exp.sensor_data[inp_sensor.name] +# ) +# +# # Check that the supplied forward model only has a single sensor +# # TODO: Can a forward model have multiple input sensors? +# if len(list(self.input_sensor_data.keys())) != 1: +# raise ValueError( +# f"The specified forward model must have " +# f"exactly one input sensor but has " +# f"{len(list(self.input_sensor_data.keys()))}" +# ) +# +# # Extract sensor names and input/output shape information +# # from forward model. +# # TODO: This should be adjusted depending on the problem definition +# # or specified by the user e.g. with keyword arguments that specify +# # the structure of the surrogate model, e.g.: +# # * `shape = None` to create an independent surrogate for each observation +# # * `shape = Sensor` to create a surrogate model per sensor +# self.vec_t = list(self.input_sensor_data.values())[0] +# self.model_names = [] +# for idx_x, os in enumerate(forward_model.output_sensors): +# for idx_t, t in enumerate(self.vec_t): +# self.model_names.append(os.name + "_" + str(t)) +# self.surrogate_kwargs["model_names"] = self.model_names +# self.num_features = len(forward_model.output_sensors) +# self.surrogate_kwargs["num_features"] = self.num_features +# +# # Get bounds +# self._get_bounds() +# +# # Initialize surrogate +# self.surrogate = surrogate_model(**surrogate_kwargs) +# self.func_pred = self.surrogate.predict +# +# # Initialize sampler +# self.sampler = sampler( +# target_function=self.target, +# surrogate_model=self.surrogate, +# domain_lower_bound=self.lower_bounds, +# domain_upper_bound=self.upper_bounds, +# fit_points_x=self.fit_points_x, +# fit_points_y=self.fit_points_y, +# test_points_x=self.test_points_x, +# test_points_y=self.test_points_y, +# evaluation_metric=self.evaluation_metric, +# verbose=self.verbose, +# run_name=self.run_name, +# ) +# +# def fit(self, **kwargs): +# """ +# Fit the surrogate model using the sampling points and observations +# stored in `self.sampler` +# +# Parameters +# ---------- +# kwargs +# +# Returns +# ------- +# +# """ +# +# self.surrogate.fit( +# self.sampler.fit_points_x, self.sampler.fit_points_y, **kwargs +# ) +# +# def sample(self, **kwargs): +# """ +# Draw samples from the forward model and fit the surrogate model +# using the specified sampler. +# +# Parameters +# ---------- +# kwargs +# +# Returns +# ------- +# +# """ +# +# self.fit_points_x, self.fit_points_y = self.sampler.sample(**kwargs) +# +# # # Get prediction function +# # self.func_pred = self.surrogate.predict +# +# def get_training_data(self, n_samples=1): +# """ +# Obtain training points `X` with latin hypercube sampling and +# corresponding observations `y` by evaluating the target function. +# +# Returns +# ------- +# A tuple of training points and observations `(train_X, train_y)` +# +# """ +# train_X = latin_hypercube_sampling( +# n_sample=n_samples, +# domain_lower_bound=np.array(self.lower_bounds), +# domain_upper_bound=np.array(self.upper_bounds), +# ) +# train_y = self.target(train_X) +# +# return train_X, train_y +# +# def get_harlow_model( +# self, +# name: str, +# ) -> ForwardModelBase: +# +# """ +# Model factory for the `HarlowSurrogate` class. +# +# TODO: +# * Refactor this function so that it returns a `ForwardModel` object +# with a `.predict()` method that queries an instance of `Surrogate`. +# * It is important that the generated `ForwardModel` instance does not +# contain the surrogate model itself (due to incompatibility of tensors +# with `copy.deepcopy`), or the problem definition. +# +# Notes +# --------- +# * This function was originally created as a workaround for passing arguments to +# classes derived from `ForwardModelBase`. +# +# Parameters +# ---------- +# problem +# forward_model +# surrogate +# n_output_dims +# kwargs +# +# Returns +# ------- +# +# """ +# +# class HarlowSurrogate(ForwardModelBase): +# +# # Get interface from existing forward model +# def interface(self): +# self.input_sensors = self.kwargs["input_sensors"] +# self.output_sensors = self.kwargs["output_sensors"] +# self.parameters = self.kwargs["parameters"] +# self.input_coordinates = self.kwargs["input_coordinates"] +# +# def response(self, inp: dict) -> dict: +# # TODO: Currently this method assumes that `Nx` and `Nt` are available. +# # Instead they should be obtained from the problem definition +# params = np.atleast_2d([inp["X" + str(i + 1)] for i in range(4)]) +# response = dict() +# +# # TODO: Adjust `X` to account for the case of sensors with vector +# # input (e.g. timeseries). +# vec_t = self.input_coordinates[self.input_sensors[0].name] +# +# # Evaluate function and arange output on grid +# f = self.kwargs["func_pred"](params, return_std=False) +# f = np.reshape(f, (len(self.output_sensors), len(vec_t))) +# +# for idx_x, os in enumerate(self.output_sensors): +# response[os.name] = f[idx_x, :] +# +# return response +# +# return HarlowSurrogate( +# name, +# func_pred=self.func_pred, +# input_sensors=self.input_sensors, +# output_sensors=self.output_sensors, +# parameters=self.parameters, +# input_coordinates=self.input_sensor_data, +# ) +# +# # TODO: This will likely fail for multidimensional parameters. FIX! +# def _get_bounds(self): +# +# self.lower_bounds = [] +# self.upper_bounds = [] +# self.bounds = {} +# +# # Find the priors corresponding to model parameters +# for key, value in self.problem.priors.items(): +# +# # Error if any of the model priors are not uniform +# if value.prior_type != "uniform": +# raise ValueError( +# f"Non-uniform distribution of type {value.prior_type}" +# f" found for prior {key}. The `HarlowSurrogate` currently" +# f" only supports uniform priors." +# ) +# +# # Dictionary of upper and lower bounds +# for param in self.parameters: +# lb = self.problem.constant_prms_dict["low_" + param] +# ub = self.problem.constant_prms_dict["high_" + param] +# +# # List of upper and lower bounds +# self.lower_bounds.append(lb) +# self.upper_bounds.append(ub) +# +# # Append to dict +# self.bounds[param] = {"low": lb, "high": ub} +# +# def _cast_to_numpy(self, inp: dict) -> np.ndarray: +# raise NotImplementedError +# +# def _cast_to_dict(self, X: np.ndarray) -> dict: +# return dict(zip(self.parameters, X)) +# +# def target(self, X: np.ndarray) -> np.ndarray: +# # TODO: Adjust `y` to account for the case of sensors with vector +# # input (e.g. timeseries). +# +# # Initialize output array +# y = np.zeros((X.shape[0], len(self.output_sensors) * len(self.vec_t))) +# +# # Evaluate physical model for each input vector +# for idx_Xi, Xi in enumerate(X): +# response = self.forward_model.response(self._cast_to_dict(Xi)) +# sensor_out = [] +# for idx_os, os in enumerate(self.output_sensors): +# sensor_out.append(response[os.name]) +# y[idx_Xi, :] = np.ravel(sensor_out) +# +# return np.array(y) diff --git a/probeye/surrogate/__init__.py b/probeye/surrogate/__init__.py deleted file mode 100644 index 0cec2ca2..00000000 --- a/probeye/surrogate/__init__.py +++ /dev/null @@ -1,2 +0,0 @@ -# subpackage imports -from probeye.surrogate import initial_sampling diff --git a/tests/integration_tests/tests_without_correlation/test_surrogate_model.py b/tests/integration_tests/tests_without_correlation/test_surrogate_model.py index a6cd4e6e..0f311168 100644 --- a/tests/integration_tests/tests_without_correlation/test_surrogate_model.py +++ b/tests/integration_tests/tests_without_correlation/test_surrogate_model.py @@ -9,7 +9,7 @@ # local imports (problem definition) from probeye.definition.inverse_problem import InverseProblem from probeye.definition.forward_model import ForwardModelBase -from probeye.definition.surrogate_model import SurrogateModelBase +from probeye.metamodeling.surrogating import SurrogateModelBase from probeye.definition.sensor import Sensor from probeye.definition.likelihood_model import GaussianLikelihoodModel diff --git a/tests/unit_tests/surrogate/test_initial_sampling.py b/tests/unit_tests/surrogate/test_initial_sampling.py index eb79d73d..b7e5a64a 100644 --- a/tests/unit_tests/surrogate/test_initial_sampling.py +++ b/tests/unit_tests/surrogate/test_initial_sampling.py @@ -8,7 +8,7 @@ from probeye.definition.inverse_problem import InverseProblem from probeye.definition.forward_model import ForwardModelBase from probeye.definition.sensor import Sensor -from probeye.surrogate.initial_sampling import LatinHypercubeSampler +from probeye.metamodeling.initial_sampling import LatinHypercubeSampler class LinearModel(ForwardModelBase):