diff --git a/docs/code/stochastic_processes/karhunen_loeve/README.rst b/docs/code/stochastic_processes/karhunen_loeve_1d/README.rst similarity index 100% rename from docs/code/stochastic_processes/karhunen_loeve/README.rst rename to docs/code/stochastic_processes/karhunen_loeve_1d/README.rst diff --git a/docs/code/stochastic_processes/karhunen_loeve/plot_karhunen_loeve.py b/docs/code/stochastic_processes/karhunen_loeve_1d/plot_karhunen_loeve_1d.py similarity index 100% rename from docs/code/stochastic_processes/karhunen_loeve/plot_karhunen_loeve.py rename to docs/code/stochastic_processes/karhunen_loeve_1d/plot_karhunen_loeve_1d.py diff --git a/docs/code/stochastic_processes/karhunen_loeve_2d/README.rst b/docs/code/stochastic_processes/karhunen_loeve_2d/README.rst new file mode 100644 index 000000000..fe518bb33 --- /dev/null +++ b/docs/code/stochastic_processes/karhunen_loeve_2d/README.rst @@ -0,0 +1,4 @@ +Karhunen Loeve Expansion 2D Examples +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + + diff --git a/docs/code/stochastic_processes/karhunen_loeve_2d/plot_karhunen_loeve_2d.py b/docs/code/stochastic_processes/karhunen_loeve_2d/plot_karhunen_loeve_2d.py new file mode 100644 index 000000000..11ffd7a0a --- /dev/null +++ b/docs/code/stochastic_processes/karhunen_loeve_2d/plot_karhunen_loeve_2d.py @@ -0,0 +1,72 @@ +""" + +Karhunen Loeve Expansion 2 Dimesional +================================================================= + +In this example, the KL Expansion is used to generate 2 dimensional stochastic fields from a prescribed Autocorrelation +Function. This example illustrates how to use the :class:`.KarhunenLoeveExpansion2D` class for a 2 dimensional +random field and compare the statistics of the generated random field with the expected values. + +""" + +# %% md +# +# Import the necessary libraries. Here we import standard libraries such as numpy and matplotlib, but also need to +# import the :class:`.KarhunenLoeveExpansionTwoDimension` class from the :py:mod:`stochastic_processes` module of UQpy. + +# %% +from matplotlib import pyplot as plt +from UQpy.stochastic_process import KarhunenLoeveExpansion2D +import numpy as np + +plt.style.use('seaborn') + +# %% md +# +# The input parameters necessary for the generation of the stochastic processes are given below: + +# %% + +n_samples = 100 # Num of samples +nx, nt = 20, 10 +dx, dt = 0.05, 0.1 + +x = np.linspace(0, (nx - 1) * dx, nx) +t = np.linspace(0, (nt - 1) * dt, nt) +xt_list = np.meshgrid(x, x, t, t, indexing='ij') # R(t_1, t_2, x_1, x_2) + +# %% md + +# Defining the Autocorrelation Function. + +# %% + +R = np.exp(-(xt_list[0] - xt_list[1]) ** 2 - (xt_list[2] - xt_list[3]) ** 2) +# R(x_1, x_2, t_1, t_2) = exp(-(x_1 - x_2) ** 2 -(t_1 - t_2) ** 2) + +KLE_Object = KarhunenLoeveExpansion2D(n_samples=n_samples, correlation_function=R, time_intervals=np.array([dt, dx]), + thresholds=[4, 5], random_state=128) + +# %% md + +# Simulating the samples. + +# %% + +samples = KLE_Object.samples + +# %% md + +# Plotting a sample of the stochastic field. + +# %% + +fig = plt.figure() +plt.title('Realisation of the Karhunen Loeve Expansion for a 2D stochastic field') +plt.imshow(samples[0, 0]) +plt.ylabel('t (Time)') +plt.xlabel('x (Space)') +plt.show() + +print('The mean of the samples is ', np.mean(samples), 'whereas the expected mean is 0.000') +print('The variance of the samples is ', np.var(samples), 'whereas the expected variance is 1.000') diff --git a/docs/source/_static/architecture/adaptive_kriging_functions.png b/docs/source/_static/architecture/adaptive_kriging_functions.png new file mode 100644 index 000000000..a4e254fd2 Binary files /dev/null and b/docs/source/_static/architecture/adaptive_kriging_functions.png differ diff --git a/docs/source/_static/architecture/dimension_reduction.png b/docs/source/_static/architecture/dimension_reduction.png new file mode 100644 index 000000000..03db86e7e Binary files /dev/null and b/docs/source/_static/architecture/dimension_reduction.png differ diff --git a/docs/source/_static/architecture/distributions.png b/docs/source/_static/architecture/distributions.png new file mode 100644 index 000000000..4043d2fcb Binary files /dev/null and b/docs/source/_static/architecture/distributions.png differ diff --git a/docs/source/_static/architecture/inference.png b/docs/source/_static/architecture/inference.png new file mode 100644 index 000000000..1e529041e Binary files /dev/null and b/docs/source/_static/architecture/inference.png differ diff --git a/docs/source/_static/architecture/mcmc.png b/docs/source/_static/architecture/mcmc.png new file mode 100644 index 000000000..676665f12 Binary files /dev/null and b/docs/source/_static/architecture/mcmc.png differ diff --git a/docs/source/_static/architecture/reliability.png b/docs/source/_static/architecture/reliability.png new file mode 100644 index 000000000..abd9ac78f Binary files /dev/null and b/docs/source/_static/architecture/reliability.png differ diff --git a/docs/source/_static/architecture/run_model.png b/docs/source/_static/architecture/run_model.png new file mode 100644 index 000000000..f96855d4a Binary files /dev/null and b/docs/source/_static/architecture/run_model.png differ diff --git a/docs/source/_static/architecture/sensitivity.png b/docs/source/_static/architecture/sensitivity.png new file mode 100644 index 000000000..b89da0e12 Binary files /dev/null and b/docs/source/_static/architecture/sensitivity.png differ diff --git a/docs/source/_static/architecture/stochastic_process.png b/docs/source/_static/architecture/stochastic_process.png new file mode 100644 index 000000000..5966f7f16 Binary files /dev/null and b/docs/source/_static/architecture/stochastic_process.png differ diff --git a/docs/source/_static/architecture/stratified_sampling.png b/docs/source/_static/architecture/stratified_sampling.png new file mode 100644 index 000000000..d264ef53d Binary files /dev/null and b/docs/source/_static/architecture/stratified_sampling.png differ diff --git a/docs/source/_static/architecture/surrogates.png b/docs/source/_static/architecture/surrogates.png new file mode 100644 index 000000000..a60f2a351 Binary files /dev/null and b/docs/source/_static/architecture/surrogates.png differ diff --git a/docs/source/_static/architecture/transformations.png b/docs/source/_static/architecture/transformations.png new file mode 100644 index 000000000..e5b732795 Binary files /dev/null and b/docs/source/_static/architecture/transformations.png differ diff --git a/docs/source/architecture.rst b/docs/source/architecture.rst new file mode 100644 index 000000000..3e4a926ed --- /dev/null +++ b/docs/source/architecture.rst @@ -0,0 +1,169 @@ +UQpy architecture +================== + +Distributions +----------------------- + +The UML diagram below explains the class hierarchy of the distributions module. In case of distributions a series of +abstract base classes are defined. :class:`.Distribution` is the basis, with classes such as :class:`.Distribution1D` \ +and :class:`.DistributionND` abstract base classes refining this hierarchy according to the distribution dimensionality. +The 1D case if further refined with :class:`.DistributionContinuous1D` and :class:`.DistributionDiscrete1D` to take into +account different types of random variables. Finally, the :class:`.Copula` abstract base class serves as the common +interface for all copula implementations. Note that all the above are abstract classes and cannot be directly +instantiated. + +.. image:: _static/architecture/distributions.png + :scale: 30 % + :alt: UML diagram of distributions module. + :align: center + +Sampling +----------------------- + +The sampling module contains multiple methods and the UML diagrams below are focused to one method at a time to reduce +the complexity of the graphs. Before starting with hierarchy of each different method, we should mention that the +classes :class:`.MonteCarloSampling`, :class:`.SimplexSampling` and :class:`.ImportanceSampling` are independent object +and do not extend any baseclass. + +Starting with the Markov Chain Monte Carlo algorithms, the diagram below makes obvious the dependency between the +different implementations. the :class:`.MCMC` abstract base class includes the common functionality between the methods, +while the specific algorithms only override the abstract methods of that base class. + + +.. image:: _static/architecture/mcmc.png + :scale: 30 % + :alt: UML diagram of MCMC hierarchy. + :align: center + +The case of stratified sampling is more elaborate, as it includes a family of algorithms. The abstract base class +:class:`.StratifiedSampling` serves as the interface between the concrete implementations :class:`.LatinHypercubeSampling`, +:class:`.TrueStratifiedSampling` and :class:`.RefinedStratifiedSampling`. In the case of :class:`.LatinHypercubeSampling` +the Strategy Design Pattern was used to extract the Latin Hypercube criteria from the sampling. Here the +:class:`.Criterion` base class provides the interface that the specific criteria need to override. In a similar manner, +the geometric stratification in the case of :class:`.TrueStratifiedSampling` is extracted under the +:class:`.Strata` abstract base class. Last but not least, in the case of :class:`.RefinedStratifiedSampling` +the different strata refinement strategies are extracted using the :class:`.Refinement` baseclass as their common interface. + + +.. image:: _static/architecture/stratified_sampling.png + :scale: 30 % + :alt: UML diagram of Stratified sampling class hierarchy. + :align: center + +In the case of :class:`.AdaptiveKriging` sampling methods, again the different learning functions are extracted into +separate classes under the common :class:`.LearningFunction` class. + +.. image:: _static/architecture/adaptive_kriging_functions.png + :scale: 30 % + :alt: UML diagram of Adaptive Kriging Hierarchy. + :align: center + + +Transformations +----------------------- + +The transformations module is one of the most simple in :py:mod:`UQpy` with three independent classes available, namely +:class:`.Nataf`, :class:`.Correlate` and :class:`.Decorrelate`. + +.. image:: _static/architecture/transformations.png + :scale: 30 % + :alt: UML diagram of Transformations module. + :align: center + + +Stochastic Processes +----------------------- + +The stochastic process module is has again simple structure with five independent classes available. + +.. image:: _static/architecture/stochastic_process.png + :scale: 30 % + :alt: UML diagram of Stochastic Process module. + :align: center + +Run Model +----------------------- + +In case of the RunModel module, the final algorithm to run is constructed by object composition of two different inputs. +Initially, the type of the model to run, with :class:`.PythonModel` and :class:`.ThirdPartyModel` being the two +available options, while the execution part is delegated to either the :class:`.SerialExecution` or :class:`.ParallelExecution` +alternatives. + +.. image:: _static/architecture/run_model.png + :scale: 30 % + :alt: UML diagram of Run Model module. + :align: center + +Inference +----------------------- + +Compared to v3, the inference module has undergone a major refactoring towards v4. The initial :class:`.InferenceModel` +class that contained all cases of computing the posterior log-likelihood is now split into three independent cases. Given +the inference models, backward uncertainty propagation can be performed be choosing between :class:`.MLE`, +:class:`.BayesParameterEstimation` to infer the parameter distributions of a model, or :class:`.InformationModelSelection` +and :class:`.BayesModelSelection` to select the model that best describes the available data. In the case of +:class:`.InformationModelSelection` the selection criteria have been extracted into separate classes under the +:class:`.InformationCriterion` baseclass. Similarly, the evidence methods of :class:`.BayesModelSelection` are also +parameters that implement the abstract base class :class:`.EvidenceMethod`. + +.. image:: _static/architecture/inference.png + :scale: 30 % + :alt: UML diagram of Inference module. + :align: center + +Reliability +----------------------- + +The reliability module maintained the same class hierarachy as in v3, with :class:`.SubsetSimulation` being an +independent class and :class:`.FORM` and :class:`.SORM` methods providing concrete implementations to the +:class:`.TaylorSeries` abstract base class. + +.. image:: _static/architecture/reliability.png + :scale: 30 % + :alt: UML diagram of Reliability module. + :align: center + +Surrogates +----------------------- + +Another module that has extensively restructured in v4 is the surrogates. Apart from the :class:`.SROM` method which +was retained as an independent algorithm, the previous Kriging functionality was removed. It is now replaced with +:class:`.GaussianProcessRegression`. The functionality of the Gaussian is constructed using object composition, +and the specific implementation of :class:`.Regression` and :class:`.Kernel` abstract base classes. An additional +functionality of constrained surrogates is added by implementing the :class:`.ConstraintsGPR` abstract class. The +functionality of :class:`.PolynomialChaosExpansion` was rewritten from scratch to address some performance issues of v3. +The Strategy Design pattern was used here as well, with three abstract base classes :class:`.Polynomials`, +:class:`.PolynomialBasis` ans :class:`.Regression` serving as the interface for the concrete classes. + +.. image:: _static/architecture/surrogates.png + :scale: 30 % + :alt: UML diagram of Surrogates module. + :align: center + +Sensitivity +----------------------- + +The sensitivity module has significantly benefited from the enhanced of modularity of the code introduced in v4. +Apart from the existing independent :class:`.MorrisSensitivity` method, the :class:`.PceSensitivity` was added as an +independent class. Finally, based on the common :class:`.Sensitivity` abstract base class, a series of new algorithms +were introduced such as :class:`.SobolSensitivity`, :class:`.GeneralizedSobolSensitivity`, :class:`.ChatterjeeSensitivity` +and :class:`.CramerVonMisesSensitivity`. + +.. image:: _static/architecture/sensitivity.png + :scale: 30 % + :alt: UML diagram of Sensitivity module. + :align: center + +Dimension Reduction +----------------------- + +The final but one of the most import modules in :py:mod:`UQpy` is dimension reduction. The :class:`.SnapshotPOD` and +:class:`.DirectPOD` methods were retained under the :class:`.POD` abstract base class. :class:`.HigherOrderSVD` method +was introduced as independent class, while special attention was given to Grassmann Manifolds. +The abstract base class :class:`.GrassmannProjection` serves as an interface for different methods to project data on +the Grassmann Manifold, with :class:`.GrassmannOperations` and :class:`.GrassmannInterpolation` support all related operations. + +.. image:: _static/architecture/dimension_reduction.png + :scale: 30 % + :alt: UML diagram of Dimension Reduction module. + :align: center diff --git a/docs/source/bibliography.bib b/docs/source/bibliography.bib index f46fc199a..fd51140c8 100644 --- a/docs/source/bibliography.bib +++ b/docs/source/bibliography.bib @@ -764,7 +764,7 @@ @article{Chatterjee } @misc{gamboa2020global, - title={Global Sensitivity Analysis: a new generation of mighty estimators based on rank statistics}, + title={Global Sensitivity Analysis: a new generation of mighty estimators based on rank statistics}, author={Fabrice Gamboa and Pierre Gremaud and Thierry Klein and Agnès Lagnoux}, year={2020}, eprint={2003.01772}, @@ -827,4 +827,16 @@ @article{saltelli_2002 url = {https://www.sciencedirect.com/science/article/pii/S0010465502002801}, author = {Andrea Saltelli}, keywords = {Sensitivity analysis, Sensitivity measures, Sensitivity indices, Importance measures}, -} \ No newline at end of file +} +@article{Kle2D, +title = {Simulation of multi-dimensional random fields by Karhunen–Loève expansion}, +journal = {Computer Methods in Applied Mechanics and Engineering}, +volume = {324}, +pages = {221-247}, +year = {2017}, +issn = {0045-7825}, +doi = {https://doi.org/10.1016/j.cma.2017.05.022}, +url = {https://www.sciencedirect.com/science/article/pii/S0045782516318692}, +author = {Zhibao Zheng and Hongzhe Dai}, +keywords = {Multi-dimensional random field, Karhunen–Loève expansion, Random field simulation, Fredholm integral equation}, +} diff --git a/docs/source/conf.py b/docs/source/conf.py index 3c7c10316..0fe8f1b3c 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -97,7 +97,8 @@ "../code/sensitivity/generalised_sobol", "../code/sensitivity/comparison", "../code/stochastic_processes/bispectral", - "../code/stochastic_processes/karhunen_loeve", + '../code/stochastic_processes/karhunen_loeve_1d', + '../code/stochastic_processes/karhunen_loeve_2d', "../code/stochastic_processes/spectral", "../code/stochastic_processes/translation", "../code/reliability/form", @@ -136,7 +137,8 @@ "auto_examples/sensitivity/generalised_sobol", "auto_examples/sensitivity/comparison", "auto_examples/stochastic_processes/bispectral", - "auto_examples/stochastic_processes/karhunen_loeve", + 'auto_examples/stochastic_processes/karhunen_loeve_1d', + 'auto_examples/stochastic_processes/karhunen_loeve_2d', "auto_examples/stochastic_processes/spectral", "auto_examples/stochastic_processes/translation", "auto_examples/reliability/form", diff --git a/docs/source/index.rst b/docs/source/index.rst index d641b0695..c3a61ae63 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -106,6 +106,7 @@ The default logging level is set to **ERROR**. The user can change the logging s /surrogates/index /transformations/index /utilities/index + architecture.rst paper.rst bibliography.rst news_doc diff --git a/docs/source/stochastic_process/index.rst b/docs/source/stochastic_process/index.rst index 4dfe5ff5e..425812e59 100644 --- a/docs/source/stochastic_process/index.rst +++ b/docs/source/stochastic_process/index.rst @@ -47,5 +47,6 @@ the desired functionality. It does not require modification of any existing clas Spectral Representation Method Bispectral Representation Method - Karhunen Loeve Expansion + Karhunen Loeve Expansion + Karhunen Loeve Expansion 2D Translation Processes diff --git a/docs/source/stochastic_process/karhunen_loeve.rst b/docs/source/stochastic_process/karhunen_loeve_1d.rst similarity index 100% rename from docs/source/stochastic_process/karhunen_loeve.rst rename to docs/source/stochastic_process/karhunen_loeve_1d.rst diff --git a/docs/source/stochastic_process/karhunen_loeve_2d.rst b/docs/source/stochastic_process/karhunen_loeve_2d.rst new file mode 100644 index 000000000..5e23d00a6 --- /dev/null +++ b/docs/source/stochastic_process/karhunen_loeve_2d.rst @@ -0,0 +1,33 @@ +Karhunen Loève Expansion for Multi-Dimensional Fields +---------------------------- + +The Karhunen Loève Expansion expands the stochastic field as follows: + +.. math:: A(x, t) = \sum_{n=1}^{\infty} \sum_{k=1}^{\infty}\eta_{nk}(\theta) \sqrt{\lambda_n(x)} f_n(t, x) \sqrt{\mu_{nk}} g_{nk}(x) + +where :math:`\eta_{nk}(\theta)` are uncorrelated standardized normal random variables and :math:`\lambda_n(x)` and :math:`f_n(x, t)` are the eigenvalues and eigenvectors repsectively of the "quasi" one dimensional covariance function :math:`C(x, t_1, t_2)`. :math:`\mu_{nk}` and :math:`g_{nk}(x)` are the eigenvalues and eigenvectors of the derived "one" dimensional covariance function :math:`H(x_1, x_2)`. Additional details regarding the simulation formula can be found +at :cite:`Kle2D` + +KarhunenLoeve2D Class +^^^^^^^^^^^^^^^^^^^^ + +The :class:`.KarhunenLoeve2D` class is imported using the following command: + +>>> from UQpy.stochastic_process.KarhunenLoeveExpansionTwoDimension2D import KarhunenLoeveExpansion + +Methods +""""""" +.. autoclass:: UQpy.stochastic_process.KarhunenLoeveExpansion2D + :members: run + +Attributes +"""""""""" +.. autoattribute:: UQpy.stochastic_process.KarhunenLoeveExpansion2D.samples +.. autoattribute:: UQpy.stochastic_process.KarhunenLoeveExpansion2D.xi + +Examples +"""""""""" + +.. toctree:: + + Karhunen Loeve Examples <../auto_examples/stochastic_processes/karhunen_loeve_2d/index> \ No newline at end of file diff --git a/requirements.txt b/requirements.txt index 4dc4412a3..d11b73a6c 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,4 +1,4 @@ -numpy == 1.22.3 +numpy == 1.21.4 scipy == 1.8.0 matplotlib == 3.5.2 scikit-learn == 1.0.2 diff --git a/src/UQpy/distributions/baseclass/Copula.py b/src/UQpy/distributions/baseclass/Copula.py index f60a27f9a..a80690255 100644 --- a/src/UQpy/distributions/baseclass/Copula.py +++ b/src/UQpy/distributions/baseclass/Copula.py @@ -2,7 +2,7 @@ DistributionContinuous1D, ) from abc import ABC - +from typing import Union class Copula(ABC): @@ -44,7 +44,7 @@ def update_parameters(self, **kwargs: dict): self.parameters[key] = kwargs[key] @staticmethod - def check_marginals(marginals: list[DistributionContinuous1D]): + def check_marginals(marginals: Union[list, DistributionContinuous1D]): """ Perform some checks on the marginals, raise errors if necessary. diff --git a/src/UQpy/reliability/taylor_series/FORM.py b/src/UQpy/reliability/taylor_series/FORM.py index dee48a4fc..47eff2c67 100644 --- a/src/UQpy/reliability/taylor_series/FORM.py +++ b/src/UQpy/reliability/taylor_series/FORM.py @@ -185,7 +185,7 @@ def run(self, seed_x: Union[list, np.ndarray] = None, g_record.append(0.0) dg_u_record = np.zeros([self.n_iterations + 1, self.dimension]) - while not converged: + while not converged and k < self.n_iterations: self.logger.info("Number of iteration: %i", k) # FORM always starts from the standard normal space if k == 0: @@ -314,8 +314,6 @@ def run(self, seed_x: Union[list, np.ndarray] = None, self.logger.info("Error: %s", error_record[-1]) - if converged is True or k > self.n_iterations: - break if k > self.n_iterations: self.logger.info("UQpy: Maximum number of iterations {0} was reached before convergence." diff --git a/src/UQpy/stochastic_process/KarhunenLoeveExpansion.py b/src/UQpy/stochastic_process/KarhunenLoeveExpansion.py index 9b67496a2..2e941fc33 100644 --- a/src/UQpy/stochastic_process/KarhunenLoeveExpansion.py +++ b/src/UQpy/stochastic_process/KarhunenLoeveExpansion.py @@ -1,5 +1,4 @@ from scipy.linalg import sqrtm -import numpy as np from UQpy.utilities import * @@ -7,12 +6,13 @@ class KarhunenLoeveExpansion: # TODO: Test this for non-stationary processes. def __init__( - self, - n_samples: int, - correlation_function: np.ndarray, - time_interval: Union[np.ndarray, float], - threshold: int = None, - random_state: RandomStateType = None, + self, + n_samples: int, + correlation_function: np.ndarray, + time_interval: Union[np.ndarray, float], + threshold: int = None, + random_state: RandomStateType = None, + random_variables: np.ndarray = None, ): """ A class to simulate stochastic processes from a given auto-correlation function based on the Karhunen-Loeve @@ -23,10 +23,10 @@ def __init__( provided, then the :class:`.KarhunenLoeveExpansion` object is created but samples are not generated. :param correlation_function: The correlation function of the stochastic process of size :code:`(n_time_intervals, n_time_intervals)` - :param time_interval: The length of time discretization. + :param time_interval: The length of time interval. :param threshold: The threshold number of eigenvalues to be used in the expansion. :param random_state: Random seed used to initialize the pseudo-random number generator. Default is :any:`None`. - + :param random_variables: The random variables used to generate the stochastic process. Default is :any:`None`. """ self.correlation_function = correlation_function self.time_interval = time_interval @@ -46,7 +46,7 @@ def __init__( """The independent gaussian random variables used in the expansion.""" if self.n_samples is not None: - self.run(n_samples=self.n_samples) + self.run(n_samples=self.n_samples, random_variables=random_variables) def _simulate(self, xi): lam, phi = np.linalg.eig(self.correlation_function) @@ -60,7 +60,7 @@ def _simulate(self, xi): samples = samples[:, np.newaxis] return samples - def run(self, n_samples: int): + def run(self, n_samples: int, random_variables: np.ndarray = None): """ Execute the random sampling in the :class:`.KarhunenLoeveExpansion` class. @@ -85,14 +85,23 @@ def run(self, n_samples: int): self.logger.info("UQpy: Stochastic Process: Running Karhunen Loeve Expansion.") self.logger.info("UQpy: Stochastic Process: Starting simulation of Stochastic Processes.") - xi = np.random.normal(size=(self.n_eigenvalues, self.n_samples)) - samples = self._simulate(xi) + + if random_variables is not None and random_variables.shape == (self.n_eigenvalues, self.n_samples): + self.logger.info('UQpy: Stochastic Process: Using user defined random variables') + else: + self.logger.info('UQpy: Stochastic Process; Using computer generated random variables') + random_variables = np.random.normal(size=(self.n_eigenvalues, self.n_samples)) + + # xi = np.random.normal(size=(self.n_eigenvalues, self.n_samples)) + samples = self._simulate(random_variables) if self.samples is None: self.samples = samples - self.xi = xi + self.random_variables = random_variables else: self.samples = np.concatenate((self.samples, samples), axis=0) - self.xi = np.concatenate((self.xi, xi), axis=0) + self.random_variables = np.concatenate((self.random_variables, random_variables), axis=0) self.logger.info("UQpy: Stochastic Process: Karhunen-Loeve Expansion Complete.") + + diff --git a/src/UQpy/stochastic_process/KarhunenLoeveExpansion2D.py b/src/UQpy/stochastic_process/KarhunenLoeveExpansion2D.py new file mode 100644 index 000000000..f6fa146ad --- /dev/null +++ b/src/UQpy/stochastic_process/KarhunenLoeveExpansion2D.py @@ -0,0 +1,119 @@ +from typing import Union + +import numpy as np + +from UQpy.utilities.ValidationTypes import RandomStateType +from UQpy.stochastic_process.KarhunenLoeveExpansion import KarhunenLoeveExpansion + + +class KarhunenLoeveExpansion2D: + + def __init__( + self, + n_samples: int, + correlation_function: np.ndarray, + time_intervals: Union[np.ndarray, float], + thresholds: Union[list, int] = None, + random_state: RandomStateType = None, + random_variables=None + ): + """ + A class to simulate two dimensional stochastic fields from a given auto-correlation function based on the + Karhunen-Loeve Expansion + + :param n_samples: Number of samples of the stochastic field to be simulated. + The :meth:`run` method is automatically called if `n_samples` is provided. If `n_samples` is not + provided, then the :class:`.KarhunenLoeveExpansion2D` object is created but samples are not generated. + :param correlation_function: The correlation function of the stochastic process of size + :code:`(n_time_intervals_dim1, n_time_intervals_dim1, n_time_intervals_dim2, n_time_intervals_dim2)` + :param time_intervals: The length of time discretizations. + :param thresholds: The threshold number of eigenvalues to be used in the expansion for two dimensions. + :param random_state: Random seed used to initialize the pseudo-random number generator. Default is :any:`None`. + :param random_variables: The random variables used to generate the stochastic field. + """ + + self.n_samples = n_samples + self.correlation_function = correlation_function + assert (len(self.correlation_function.shape) == 4) + self.time_intervals = time_intervals + self.thresholds = thresholds + self.random_state = random_state + + if isinstance(self.random_state, int): + np.random.seed(self.random_state) + elif not isinstance(self.random_state, (type(None), np.random.RandomState)): + raise TypeError('UQpy: random_state must be None, an int or an np.random.RandomState object.') + + self.samples = None + """Array of generated samples.""" + self.random_variables = random_variables + + self._precompute_one_dimensional_correlation_function() + + if self.n_samples is not None: + self.run(n_samples=self.n_samples, random_variables=random_variables) + + def _precompute_one_dimensional_correlation_function(self): + self.quasi_correlation_function = np.zeros( + [self.correlation_function.shape[1], self.correlation_function.shape[2], + self.correlation_function.shape[3]]) + for i in range(self.correlation_function.shape[0]): + self.quasi_correlation_function[i] = self.correlation_function[i, i] + self.w, self.v = np.linalg.eig(self.quasi_correlation_function) + if np.linalg.norm(np.imag(self.w)) > 0: + print('Complex in the eigenvalues, check the positive definiteness') + self.w = np.real(self.w) + self.v = np.real(self.v) + if self.thresholds is not None: + self.w = self.w[:, :self.thresholds[1]] + self.v = self.v[:, :, :self.thresholds[1]] + self.one_dimensional_correlation_function = np.einsum('uvxy, uxn, vyn, un, vn -> nuv', + self.correlation_function, self.v, self.v, + 1 / np.sqrt(self.w), 1 / np.sqrt(self.w)) + + def run(self, n_samples, random_variables=None): + """ + Execute the random sampling in the :class:`.KarhunenLoeveExpansion` class. + + The :meth:`run` method is the function that performs random sampling in the :class:`.KarhunenLoeveExpansion2D` + class. If `n_samples` is provided when the :class:`.KarhunenLoeveExpansion2D` object is defined, the :meth:`run` + method is automatically called. The user may also call the :meth:`run` method directly to generate samples. + The :meth:`run`` method of the :class:`.KarhunenLoeveExpansion2D` class can be invoked many times and each time + the generated samples are appended to the existing samples. + + :param n_samples: Number of samples of the stochastic process to be simulated. + If the :meth:`run` method is invoked multiple times, the newly generated samples will be appended to the + existing samples. + + The :meth:`run` method has no returns, although it creates and/or appends the :py:attr:`samples` attribute of + the :class:`KarhunenLoeveExpansion2D` class. + """ + samples = np.zeros((n_samples, self.correlation_function.shape[0], self.correlation_function.shape[2])) + if random_variables is None: + random_variables = np.random.normal(size=[self.thresholds[1], self.thresholds[0], n_samples]) + else: + assert (random_variables.shape == (self.thresholds[1], self.thresholds[0], n_samples)) + for i in range(self.one_dimensional_correlation_function.shape[0]): + if self.thresholds is not None: + samples += np.einsum('x, xt, nx -> nxt', np.sqrt(self.w[:, i]), self.v[:, :, i], + KarhunenLoeveExpansion(n_samples=n_samples, + correlation_function= + self.one_dimensional_correlation_function[i], + time_interval=self.time_intervals, + threshold=self.thresholds[0], + random_variables=random_variables[i]).samples[:, 0, :]) + else: + samples += np.einsum('x, xt, nx -> nxt', np.sqrt(self.w[:, i]), self.v[:, :, i], + KarhunenLoeveExpansion(n_samples=n_samples, + correlation_function= + self.one_dimensional_correlation_function[i], + time_interval=self.time_intervals, + random_variables=random_variables[i]).samples[:, 0, :]) + samples = np.reshape(samples, [samples.shape[0], 1, samples.shape[1], samples.shape[2]]) + + if self.samples is None: + self.samples = samples + self.xi = random_variables + else: + self.samples = np.concatenate((self.samples, samples), axis=0) + self.xi = np.concatenate((self.xi, random_variables), axis=2) diff --git a/src/UQpy/stochastic_process/__init__.py b/src/UQpy/stochastic_process/__init__.py index 5a0320310..4c36e35f3 100644 --- a/src/UQpy/stochastic_process/__init__.py +++ b/src/UQpy/stochastic_process/__init__.py @@ -1,6 +1,7 @@ from UQpy.stochastic_process.BispectralRepresentation import BispectralRepresentation from UQpy.stochastic_process.InverseTranslation import InverseTranslation from UQpy.stochastic_process.KarhunenLoeveExpansion import KarhunenLoeveExpansion +from UQpy.stochastic_process.KarhunenLoeveExpansion2D import KarhunenLoeveExpansion2D from UQpy.stochastic_process.SpectralRepresentation import SpectralRepresentation from UQpy.stochastic_process.Translation import Translation diff --git a/tests/unit_tests/stochastic_process/test_karhunen_loeve.py b/tests/unit_tests/stochastic_process/test_karhunen_loeve_1d.py similarity index 100% rename from tests/unit_tests/stochastic_process/test_karhunen_loeve.py rename to tests/unit_tests/stochastic_process/test_karhunen_loeve_1d.py diff --git a/tests/unit_tests/stochastic_process/test_karhunen_loeve_2d.py b/tests/unit_tests/stochastic_process/test_karhunen_loeve_2d.py new file mode 100644 index 000000000..1bd2fcd21 --- /dev/null +++ b/tests/unit_tests/stochastic_process/test_karhunen_loeve_2d.py @@ -0,0 +1,37 @@ +from UQpy.stochastic_process.KarhunenLoeveExpansion2D import KarhunenLoeveExpansion2D +import numpy as np + +n_samples = 100_000 # Num of samples +nx, nt = 20, 10 +dx, dt = 0.05, 0.1 + +x = np.linspace(0, (nx - 1) * dx, nx) +t = np.linspace(0, (nt - 1) * dt, nt) +xt_list = np.meshgrid(x, x, t, t, indexing='ij') # R(t_1, t_2, x_1, x_2) + +R = np.exp(-(xt_list[0] - xt_list[1]) ** 2 - (xt_list[2] - xt_list[3]) ** 2) +# R(x_1, x_2, t_1, t_2) = exp(-(x_1 - x_2) ** 2 -(t_1 - t_2) ** 2) + +KLE_Object = KarhunenLoeveExpansion2D(n_samples=n_samples, correlation_function=R, + time_intervals=np.array([dt, dx]), thresholds=[4, 5], random_state=128) +samples = KLE_Object.samples + + +def test_samples_shape(): + assert samples.shape == (n_samples, 1, len(x), len(t)) + + +def test_samples_values(): + assert np.isclose(samples[13, 0, 13, 6], -1.8616264088843137, rtol=0.01) + + +def test_run_method(): + nsamples_second_run = 100_000 + KLE_Object.run(n_samples=nsamples_second_run) + samples = KLE_Object.samples + assert samples.shape == (n_samples + nsamples_second_run, 1, len(x), len(t)) + + +def test_empirical_correlation(): + assert np.isclose(np.mean(samples[:, 0, 5, 6] * samples[:, 0, 3, 2]), R[5, 3, 6, 2], rtol=0.01) + assert np.isclose(np.mean(samples[:, 0, 13, 3] * samples[:, 0, 8, 4]), R[13, 8, 3, 4], rtol=0.01)