Module bayesvalidrox.surrogate_models.apoly_construction
Expand source code
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
import numpy as np
def apoly_construction(Data, degree):
"""
Construction of Data-driven Orthonormal Polynomial Basis
Author: Dr.-Ing. habil. Sergey Oladyshkin
Department of Stochastic Simulation and Safety Research for Hydrosystems
Institute for Modelling Hydraulic and Environmental Systems
Universitaet Stuttgart, Pfaffenwaldring 5a, 70569 Stuttgart
E-mail: Sergey.Oladyshkin@iws.uni-stuttgart.de
http://www.iws-ls3.uni-stuttgart.de
The current script is based on definition of arbitrary polynomial chaos
expansion (aPC), which is presented in the following manuscript:
Oladyshkin, S. and W. Nowak. Data-driven uncertainty quantification using
the arbitrary polynomial chaos expansion. Reliability Engineering & System
Safety, Elsevier, V. 106, P. 179-190, 2012.
DOI: 10.1016/j.ress.2012.05.002.
Parameters
----------
Data : array
Raw data.
degree : int
Maximum polynomial degree.
Returns
-------
Polynomial : array
The coefficients of the univariate orthonormal polynomials.
"""
# Initialization
dd = degree + 1
nsamples = len(Data)
# Forward linear transformation (Avoiding numerical issues)
MeanOfData = np.mean(Data)
Data = Data/MeanOfData
# Compute raw moments of input data
raw_moments = [np.sum(np.power(Data, p))/nsamples for p in range(2*dd+2)]
# Main Loop for Polynomial with degree up to dd
PolyCoeff_NonNorm = np.empty((0, 1))
Polynomial = np.zeros((dd+1, dd+1))
for degree in range(dd+1):
Mm = np.zeros((degree+1, degree+1))
Vc = np.zeros((degree+1))
# Define Moments Matrix Mm
for i in range(degree+1):
for j in range(degree+1):
if (i < degree):
Mm[i, j] = raw_moments[i+j]
elif (i == degree) and (j == degree):
Mm[i, j] = 1
# Numerical Optimization for Matrix Solver
Mm[i] = Mm[i] / max(abs(Mm[i]))
# Defenition of Right Hand side ortogonality conditions: Vc
for i in range(degree+1):
Vc[i] = 1 if i == degree else 0
# Solution: Coefficients of Non-Normal Orthogonal Polynomial: Vp Eq.(4)
try:
Vp = np.linalg.solve(Mm, Vc)
except:
inv_Mm = np.linalg.pinv(Mm)
Vp = np.dot(inv_Mm, Vc.T)
if degree == 0:
PolyCoeff_NonNorm = np.append(PolyCoeff_NonNorm, Vp)
if degree != 0:
if degree == 1:
zero = [0]
else:
zero = np.zeros((degree, 1))
PolyCoeff_NonNorm = np.hstack((PolyCoeff_NonNorm, zero))
PolyCoeff_NonNorm = np.vstack((PolyCoeff_NonNorm, Vp))
if 100*abs(sum(abs(np.dot(Mm, Vp)) - abs(Vc))) > 0.5:
print('\n---> Attention: Computational Error too high !')
print('\n---> Problem: Convergence of Linear Solver')
# Original Numerical Normalization of Coefficients with Norm and
# orthonormal Basis computation Matrix Storrage
# Note: Polynomial(i,j) correspont to coefficient number "j-1"
# of polynomial degree "i-1"
P_norm = 0
for i in range(nsamples):
Poly = 0
for k in range(degree+1):
if degree == 0:
Poly += PolyCoeff_NonNorm[k] * (Data[i]**k)
else:
Poly += PolyCoeff_NonNorm[degree, k] * (Data[i]**k)
P_norm += Poly**2 / nsamples
P_norm = np.sqrt(P_norm)
for k in range(degree+1):
if degree == 0:
Polynomial[degree, k] = PolyCoeff_NonNorm[k]/P_norm
else:
Polynomial[degree, k] = PolyCoeff_NonNorm[degree, k]/P_norm
# Backward linear transformation to the real data space
Data *= MeanOfData
for k in range(len(Polynomial)):
Polynomial[:, k] = Polynomial[:, k] / (MeanOfData**(k))
return Polynomial
Functions
def apoly_construction(Data, degree)
-
Construction of Data-driven Orthonormal Polynomial Basis Author: Dr.-Ing. habil. Sergey Oladyshkin Department of Stochastic Simulation and Safety Research for Hydrosystems Institute for Modelling Hydraulic and Environmental Systems Universitaet Stuttgart, Pfaffenwaldring 5a, 70569 Stuttgart E-mail: Sergey.Oladyshkin@iws.uni-stuttgart.de http://www.iws-ls3.uni-stuttgart.de The current script is based on definition of arbitrary polynomial chaos expansion (aPC), which is presented in the following manuscript: Oladyshkin, S. and W. Nowak. Data-driven uncertainty quantification using the arbitrary polynomial chaos expansion. Reliability Engineering & System Safety, Elsevier, V. 106, P. 179-190, 2012. DOI: 10.1016/j.ress.2012.05.002.
Parameters
Data
:array
- Raw data.
degree
:int
- Maximum polynomial degree.
Returns
Polynomial
:array
- The coefficients of the univariate orthonormal polynomials.
Expand source code
def apoly_construction(Data, degree): """ Construction of Data-driven Orthonormal Polynomial Basis Author: Dr.-Ing. habil. Sergey Oladyshkin Department of Stochastic Simulation and Safety Research for Hydrosystems Institute for Modelling Hydraulic and Environmental Systems Universitaet Stuttgart, Pfaffenwaldring 5a, 70569 Stuttgart E-mail: Sergey.Oladyshkin@iws.uni-stuttgart.de http://www.iws-ls3.uni-stuttgart.de The current script is based on definition of arbitrary polynomial chaos expansion (aPC), which is presented in the following manuscript: Oladyshkin, S. and W. Nowak. Data-driven uncertainty quantification using the arbitrary polynomial chaos expansion. Reliability Engineering & System Safety, Elsevier, V. 106, P. 179-190, 2012. DOI: 10.1016/j.ress.2012.05.002. Parameters ---------- Data : array Raw data. degree : int Maximum polynomial degree. Returns ------- Polynomial : array The coefficients of the univariate orthonormal polynomials. """ # Initialization dd = degree + 1 nsamples = len(Data) # Forward linear transformation (Avoiding numerical issues) MeanOfData = np.mean(Data) Data = Data/MeanOfData # Compute raw moments of input data raw_moments = [np.sum(np.power(Data, p))/nsamples for p in range(2*dd+2)] # Main Loop for Polynomial with degree up to dd PolyCoeff_NonNorm = np.empty((0, 1)) Polynomial = np.zeros((dd+1, dd+1)) for degree in range(dd+1): Mm = np.zeros((degree+1, degree+1)) Vc = np.zeros((degree+1)) # Define Moments Matrix Mm for i in range(degree+1): for j in range(degree+1): if (i < degree): Mm[i, j] = raw_moments[i+j] elif (i == degree) and (j == degree): Mm[i, j] = 1 # Numerical Optimization for Matrix Solver Mm[i] = Mm[i] / max(abs(Mm[i])) # Defenition of Right Hand side ortogonality conditions: Vc for i in range(degree+1): Vc[i] = 1 if i == degree else 0 # Solution: Coefficients of Non-Normal Orthogonal Polynomial: Vp Eq.(4) try: Vp = np.linalg.solve(Mm, Vc) except: inv_Mm = np.linalg.pinv(Mm) Vp = np.dot(inv_Mm, Vc.T) if degree == 0: PolyCoeff_NonNorm = np.append(PolyCoeff_NonNorm, Vp) if degree != 0: if degree == 1: zero = [0] else: zero = np.zeros((degree, 1)) PolyCoeff_NonNorm = np.hstack((PolyCoeff_NonNorm, zero)) PolyCoeff_NonNorm = np.vstack((PolyCoeff_NonNorm, Vp)) if 100*abs(sum(abs(np.dot(Mm, Vp)) - abs(Vc))) > 0.5: print('\n---> Attention: Computational Error too high !') print('\n---> Problem: Convergence of Linear Solver') # Original Numerical Normalization of Coefficients with Norm and # orthonormal Basis computation Matrix Storrage # Note: Polynomial(i,j) correspont to coefficient number "j-1" # of polynomial degree "i-1" P_norm = 0 for i in range(nsamples): Poly = 0 for k in range(degree+1): if degree == 0: Poly += PolyCoeff_NonNorm[k] * (Data[i]**k) else: Poly += PolyCoeff_NonNorm[degree, k] * (Data[i]**k) P_norm += Poly**2 / nsamples P_norm = np.sqrt(P_norm) for k in range(degree+1): if degree == 0: Polynomial[degree, k] = PolyCoeff_NonNorm[k]/P_norm else: Polynomial[degree, k] = PolyCoeff_NonNorm[degree, k]/P_norm # Backward linear transformation to the real data space Data *= MeanOfData for k in range(len(Polynomial)): Polynomial[:, k] = Polynomial[:, k] / (MeanOfData**(k)) return Polynomial