diff --git a/UQit/__init__.py b/UQit/__init__.py index cfa111f..dabf4c8 100644 --- a/UQit/__init__.py +++ b/UQit/__init__.py @@ -1,17 +1,17 @@ -__all__ = ['analyticTestFuncs','gpr_torch','lagInt','linAlg', - 'nodes','pce','ppce','reshaper','sampling', - 'sobol','stats','surr2surr','write'] +__all__ = ['analyticTestFuncs', 'gpr_torch', 'lagInt', 'linAlg', + 'nodes', 'pce', 'ppce', 'reshaper', 'sampling', + 'sobol', 'stats', 'surr2surr', 'write'] -from UQit.analyticTestFuncs import fEx1D,fEx2D,fEx3D +from UQit.analyticTestFuncs import fEx1D, fEx2D, fEx3D from UQit.gpr_torch import gpr, gprPost, gprPlot from UQit.lagInt import lagInt, lagInt_Quads2Line from UQit.linAlg import myLinearRegress -from UQit.nodes import Clenshaw_pts,ClenshawCurtis_pts,gllPts +from UQit.nodes import Clenshaw_pts, ClenshawCurtis_pts, gllPts from UQit.pce import pce, pceEval, convPlot from UQit.ppce import ppce -from UQit.reshaper import lengthVector,vecs2grid,vecsGlue -from UQit.sampling import trainSample,testSample,LHS_sampling +from UQit.reshaper import lengthVector, vecs2grid, vecsGlue +from UQit.sampling import trainSample, testSample, LHS_sampling from UQit.sobol import sobol -from UQit.stats import pdfFit_uniVar,pdfPredict_uniVar +from UQit.stats import pdfFit_uniVar, pdfPredict_uniVar from UQit.surr2surr import lagIntAtGQs from UQit.write import printRepeated diff --git a/UQit/analyticTestFuncs.py b/UQit/analyticTestFuncs.py index a06cc73..c026a2a 100644 --- a/UQit/analyticTestFuncs.py +++ b/UQit/analyticTestFuncs.py @@ -1,12 +1,12 @@ ################################################### -# Analytical model functions to test implementation +# Analytical model functions to test implementation # of different UQ techniques ################################################### -#-------------------------------------------------- +# -------------------------------------------------- # Saleh Rezaeiravesh, salehr@kth.se -#-------------------------------------------------- -#TODO: Generalize Sobol of Ishigami to [ai,bi] -#-------------------------------------------- +# -------------------------------------------------- +# TODO: Generalize Sobol of Ishigami to [ai,bi] +# -------------------------------------------- # import os import sys @@ -15,6 +15,8 @@ import UQit.sampling as sampling # # + + class fEx1D: """ Analytical test functions and their exact moments for 1D parameter @@ -44,62 +46,66 @@ class fEx1D: `var`: float V[f(q)] for `q` """ - def __init__(self,q,typ,qInfo=[]): - self.q=q - self.typ=typ - self.qInfo=qInfo + + def __init__(self, q, typ, qInfo=[]): + self.q = q + self.typ = typ + self.qInfo = qInfo self._check() self.eval() def _check(self): - if self.typ not in ['type1','type2']: - raise ValueError("Invalid 'typ'. Choose 'type1' (for q~Uniform) or 'type2' (for q~Normal)") - + if self.typ not in ['type1', 'type2']: + raise ValueError( + "Invalid 'typ'. Choose 'type1' (for q~Uniform) or 'type2' (for q~Normal)") + def eval(self): """ Value of f(q) at q """ z = np.array(self.q, copy=False, ndmin=1) - if self.typ=='type1': - val=(10.0+.7*np.sin(5.0*z)+3.*np.cos(z)) - val=np.array(val) - elif self.typ=='type2': - m=self.qInfo[0] - val=np.cos(m*(z-m)) - self.val=val - - def _mean1(self,q_): + if self.typ == 'type1': + val = (10.0+.7*np.sin(5.0*z)+3.*np.cos(z)) + val = np.array(val) + elif self.typ == 'type2': + m = self.qInfo[0] + val = np.cos(m*(z-m)) + self.val = val + + def _mean1(self, q_): return (10.*q_-0.7*mt.cos(5.*q_)/5.0+3.*mt.sin(q_)) - def _var1(self,q_): - tmp=100*q_+0.245*(q_-0.1*mt.sin(10*q_))+4.5*(q_+0.5*mt.sin(2*q_)) + def _var1(self, q_): + tmp = 100*q_+0.245*(q_-0.1*mt.sin(10*q_))+4.5*(q_+0.5*mt.sin(2*q_)) return (tmp-2.8*mt.cos(5*q_)+60*mt.sin(q_)-2.1*(mt.cos(6*q_)/6.+0.25*mt.cos(4*q_))) - def _mean2(self,q_): - m=q_[0] - v=q_[1] - return(mt.exp(-0.5*(m*v)**2)) + def _mean2(self, q_): + m = q_[0] + v = q_[1] + return(mt.exp(-0.5*(m*v)**2)) - def _var2(self,q_): - m=q_[0] - v=q_[1] - t_=(m*v)**2 + def _var2(self, q_): + m = q_[0] + v = q_[1] + t_ = (m*v)**2 return(0.5*(1+mt.exp(-2*t_))-mt.exp(-t_)) - def moments(self,qInfo): + def moments(self, qInfo): """ Analytical Mean and variance of f(q) """ - Q=qInfo - if self.typ=='type1': - fMean=(self._mean1(Q[1])-self._mean1(Q[0]))/(Q[1]-Q[0]) - fVar =(self._var1(Q[1])-self._var1(Q[0]))/(Q[1]-Q[0])-fMean**2. - elif self.typ=='type2': - fMean=self._mean2(Q) - fVar=self._var2(Q) - self.mean=fMean - self.var=fVar + Q = qInfo + if self.typ == 'type1': + fMean = (self._mean1(Q[1])-self._mean1(Q[0]))/(Q[1]-Q[0]) + fVar = (self._var1(Q[1])-self._var1(Q[0]))/(Q[1]-Q[0])-fMean**2. + elif self.typ == 'type2': + fMean = self._mean2(Q) + fVar = self._var2(Q) + self.mean = fMean + self.var = fVar # + + class fEx2D: """ Analytical test functions for 2D parameter @@ -133,65 +139,67 @@ class fEx2D: `Sij`: [S12], Dual interaction """ - def __init__(self,q1,q2,typ,method): - self.q1=q1 - self.q2=q2 - self.typ=typ - self.method=method + + def __init__(self, q1, q2, typ, method): + self.q1 = q1 + self.q2 = q2 + self.typ = typ + self.method = method self._check() self.eval() def _check(self): - method_valid=['comp','tensorProd'] + method_valid = ['comp', 'tensorProd'] if self.method not in method_valid: - raise ValueError("Invalid 'method'. Choose from ",method_valid) - typ_valid=['type1','type2','type3','Rosenbrock'] + raise ValueError("Invalid 'method'. Choose from ", method_valid) + typ_valid = ['type1', 'type2', 'type3', 'Rosenbrock'] if self.typ not in typ_valid: - raise ValueError("Invalid 'typ'. Choose from ",typ_valid) - - def _funVal(self,z1_,z2_,typ): - if typ=='type1': # from: https://se.mathworks.com/help/symbolic/graphics.html?s_tid=CRUX_lftnav - tmp1=3.*mt.exp(-(z2_+2)**2.-z1_**2) * (z1_-1.0)**2. - tmp2=-(mt.exp(-(z1_+2)**2.-z1_**2))/3. - tmp3=mt.exp(-(z1_**2.+z2_**2.))*(10.*z1_**3.-2.*z1_+10.*z2_**5.) - tmp=tmp1+tmp2+tmp3 - elif typ=='type2': - tmp1=3.*mt.exp(-(z2_+1)**2.-z1_**2) * (z1_-1.0)**2. - tmp2=-1.*mt.exp(-(z2_-1.)**2.-z1_**2) * (z1_-1.)**2. - tmp=tmp1+tmp2+0.001 - elif typ=='type3': #simple enough for analytical derivation of Sobol indices - tmp=z1_**2.+z1_*z2_ - elif typ=='Rosenbrock': - tmp=100*(z2_-z1_**2.)**2.+(1-z1_)**2. + raise ValueError("Invalid 'typ'. Choose from ", typ_valid) + + def _funVal(self, z1_, z2_, typ): + if typ == 'type1': # from: https://se.mathworks.com/help/symbolic/graphics.html?s_tid=CRUX_lftnav + tmp1 = 3.*mt.exp(-(z2_+2)**2.-z1_**2) * (z1_-1.0)**2. + tmp2 = -(mt.exp(-(z1_+2)**2.-z1_**2))/3. + tmp3 = mt.exp(-(z1_**2.+z2_**2.))*(10.*z1_**3.-2.*z1_+10.*z2_**5.) + tmp = tmp1+tmp2+tmp3 + elif typ == 'type2': + tmp1 = 3.*mt.exp(-(z2_+1)**2.-z1_**2) * (z1_-1.0)**2. + tmp2 = -1.*mt.exp(-(z2_-1.)**2.-z1_**2) * (z1_-1.)**2. + tmp = tmp1+tmp2+0.001 + elif typ == 'type3': # simple enough for analytical derivation of Sobol indices + tmp = z1_**2.+z1_*z2_ + elif typ == 'Rosenbrock': + tmp = 100*(z2_-z1_**2.)**2.+(1-z1_)**2. return tmp - + def eval(self): """ Evaluates f(q) at given q1, q2 """ z1 = np.array(self.q1, copy=False, ndmin=1) z2 = np.array(self.q2, copy=False, ndmin=1) - n1=z1.shape[0] - n2=z2.shape[0] - f=[] - if (self.method=='tensorProd'): - for i2 in range(n2): - z2_=z2[i2] - for i1 in range(n1): - z1_=z1[i1] - tmp=self._funVal(z1_,z2_,self.typ) - f.append(tmp) - elif (self.method=='comp'): - if (n1!=n2): - raise ValueError('Pairs of a paramater sample vector should have the same size') - for i in range(n1): - z1_=z1[i] - z2_=z2[i] - tmp=self._funVal(z1_,z2_,self.typ) - f.append(tmp) - self.val=np.asarray(f) - - def moments(self,distType,qInfo): + n1 = z1.shape[0] + n2 = z2.shape[0] + f = [] + if (self.method == 'tensorProd'): + for i2 in range(n2): + z2_ = z2[i2] + for i1 in range(n1): + z1_ = z1[i1] + tmp = self._funVal(z1_, z2_, self.typ) + f.append(tmp) + elif (self.method == 'comp'): + if (n1 != n2): + raise ValueError( + 'Pairs of a paramater sample vector should have the same size') + for i in range(n1): + z1_ = z1[i] + z2_ = z2[i] + tmp = self._funVal(z1_, z2_, self.typ) + f.append(tmp) + self.val = np.asarray(f) + + def moments(self, distType, qInfo): """ Mean and variance of f(q) estimated by the Monte-Carlo approach (These can be used as reference values instead of the analytical values) @@ -203,37 +211,37 @@ def moments(self,distType,qInfo): Information about the parameter range or distribution. * If `q` is Gaussian ('Norm' or 'normRand') => qInfo=[mean,sdev] * Otherwise, qInfo=[min(q),max(q)]=admissible range of q - + Returns: `mean`: float Expected value of f(q) estimated by the Monte-Carlo method `var`: float Variance of f(q) estimated by the Monte-Carlo method """ - nMC=100000 #number of MC samples - print('... Reference moments are calculated by the Monte-Carlo method with %d samples' %nMC) - qMC=[] - p=len(distType) - if p!=2: - raise ValueError("distType should have length 2") + nMC = 100000 # number of MC samples + print('... Reference moments are calculated by the Monte-Carlo method with %d samples' % nMC) + qMC = [] + p = len(distType) + if p != 2: + raise ValueError("distType should have length 2") for i in range(p): - if distType[i]=='Unif': - sampleType_='unifRand' - elif distType[i]=='Norm': - sampleType_='normRand' + if distType[i] == 'Unif': + sampleType_ = 'unifRand' + elif distType[i] == 'Norm': + sampleType_ = 'normRand' else: - raise ValueError("Invalid distType for parameter %d" %i) - samps=sampling.trainSample(sampleType=sampleType_,GQdistType=distType[i], - qInfo=qInfo[i],nSamp=nMC) + raise ValueError("Invalid distType for parameter %d" % i) + samps = sampling.trainSample(sampleType=sampleType_, GQdistType=distType[i], + qInfo=qInfo[i], nSamp=nMC) qMC.append(samps.q) - self.fVal_mc=fEx2D(qMC[0],qMC[1],self.typ,'comp').val - self.mean=np.mean(self.fVal_mc) - self.var=np.mean(self.fVal_mc**2.)-self.mean**2. + self.fVal_mc = fEx2D(qMC[0], qMC[1], self.typ, 'comp').val + self.mean = np.mean(self.fVal_mc) + self.var = np.mean(self.fVal_mc**2.)-self.mean**2. - def sobol(self,qBound): + def sobol(self, qBound): """ Sobol sensitivity indices of f(q) with respect to q1 and q2 - + Args: `qBound`: =[qBound1,qBound2] admissible range of q1, q2 @@ -245,64 +253,68 @@ def sobol(self,qBound): `Sij`: =[S12], dual interaction """ - if self.typ=='type3': - #Exact sobol indices for f(q1,q2)=q1^2+q1*q2 - #->Take a,b according to eval() for 'type3' - def _D1_Ex(a,b,q1): - """ - a,b: terms in derived expression f1(q1)=q1^2+a*q1+b - q1: sample for q1 at which the expression is evaluated - """ - return(0.2*q1**5.+0.5*a*q1**4.+(a**2.+2.*b)/3.*q1**3.+a*b*q1**2.+b**2.*q1) - def _D2_Ex(a,b,q2): - """ - a,b: terms in derived expression f2(q2)=a*q2+b - q2: sample for q2 at which the expression is evaluated - """ - return((a**2./3.)*q2**3.+a*b*q2**2.+b**2.*q2) - def _D12_Ex(a,b,c,q1Bound,q2Bound): - """ - a,b,c: terms in derived expression f12(q1,q2)=q1*q2+a*q1+b*q2+c - """ - q1_1=q1Bound[1]-q1Bound[0] - q1_2=q1Bound[1]**2.-q1Bound[0]**2. - q1_3=q1Bound[1]**3.-q1Bound[0]**3. - q2_1=q2Bound[1]-q2Bound[0] - q2_2=q2Bound[1]**2.-q2Bound[0]**2. - q2_3=q2Bound[1]**3.-q2Bound[0]**3. - return(q1_3*q2_3/9.+a/3.*q1_3*q2_2+b*q1_2*q2_3/3.+ - 0.5*(c+a*b)*q1_2*q2_2+a**2.*q1_3*q2_1/3.+b**2.*q1_1*q2_3/3.+ - a*c*q1_2*q2_1+b*c*q1_1*q2_2+c**2.*q1_1*q2_1) - #Variances in Sobol decomposition - a1=qBound[0][0] - b1=qBound[0][1] - a2=qBound[1][0] - b2=qBound[1][1] - f0=(a1**2.+b1**2.+a1*b1)/3.+0.25*(a1+b1)*(a2+b2) - a=(a2+b2)/2.0 - b=-f0 - D1_ex=(_D1_Ex(a,b,b1)-_D1_Ex(a,b,a1))/(b1-a1) - a=(a1+b1)/2.0 - b=(a1**2.+a1*b1+b1**2.)/3.-f0 - D2_ex=(_D2_Ex(a,b,b2)-_D2_Ex(a,b,a2))/(b2-a2) - a=-(a2+b2)/2. - b=-(a1+b1)/2. - c=f0-(a1**2.+a1*b1+b1**2.)/3. - D12_ex=(_D12_Ex(a,b,c,qBound[0],qBound[1]))/((b1-a1)*(b2-a2)) - #Sensitivity indices - D_ex=D1_ex+D2_ex+D12_ex - Si=[D1_ex/D_ex,D2_ex/D_ex] - Sij=[D12_ex/D_ex] - STi=[Si[0]+Sij[0],Si[1]+Sij[0]] - self.Si=Si - self.Sij=Sij - self.STi=STi + if self.typ == 'type3': + # Exact sobol indices for f(q1,q2)=q1^2+q1*q2 + # ->Take a,b according to eval() for 'type3' + def _D1_Ex(a, b, q1): + """ + a,b: terms in derived expression f1(q1)=q1^2+a*q1+b + q1: sample for q1 at which the expression is evaluated + """ + return(0.2*q1**5.+0.5*a*q1**4.+(a**2.+2.*b)/3.*q1**3.+a*b*q1**2.+b**2.*q1) + + def _D2_Ex(a, b, q2): + """ + a,b: terms in derived expression f2(q2)=a*q2+b + q2: sample for q2 at which the expression is evaluated + """ + return((a**2./3.)*q2**3.+a*b*q2**2.+b**2.*q2) + + def _D12_Ex(a, b, c, q1Bound, q2Bound): + """ + a,b,c: terms in derived expression f12(q1,q2)=q1*q2+a*q1+b*q2+c + """ + q1_1 = q1Bound[1]-q1Bound[0] + q1_2 = q1Bound[1]**2.-q1Bound[0]**2. + q1_3 = q1Bound[1]**3.-q1Bound[0]**3. + q2_1 = q2Bound[1]-q2Bound[0] + q2_2 = q2Bound[1]**2.-q2Bound[0]**2. + q2_3 = q2Bound[1]**3.-q2Bound[0]**3. + return(q1_3*q2_3/9.+a/3.*q1_3*q2_2+b*q1_2*q2_3/3. + + 0.5*(c+a*b)*q1_2*q2_2+a**2.*q1_3*q2_1/3.+b**2.*q1_1*q2_3/3. + + a*c*q1_2*q2_1+b*c*q1_1*q2_2+c**2.*q1_1*q2_1) + # Variances in Sobol decomposition + a1 = qBound[0][0] + b1 = qBound[0][1] + a2 = qBound[1][0] + b2 = qBound[1][1] + f0 = (a1**2.+b1**2.+a1*b1)/3.+0.25*(a1+b1)*(a2+b2) + a = (a2+b2)/2.0 + b = -f0 + D1_ex = (_D1_Ex(a, b, b1)-_D1_Ex(a, b, a1))/(b1-a1) + a = (a1+b1)/2.0 + b = (a1**2.+a1*b1+b1**2.)/3.-f0 + D2_ex = (_D2_Ex(a, b, b2)-_D2_Ex(a, b, a2))/(b2-a2) + a = -(a2+b2)/2. + b = -(a1+b1)/2. + c = f0-(a1**2.+a1*b1+b1**2.)/3. + D12_ex = (_D12_Ex(a, b, c, qBound[0], qBound[1]))/((b1-a1)*(b2-a2)) + # Sensitivity indices + D_ex = D1_ex+D2_ex+D12_ex + Si = [D1_ex/D_ex, D2_ex/D_ex] + Sij = [D12_ex/D_ex] + STi = [Si[0]+Sij[0], Si[1]+Sij[0]] + self.Si = Si + self.Sij = Sij + self.STi = STi else: - print("No Exact Sobol indices for 'typ' else than 'type3'") - self.Si=[] - self.Sij=[] - self.STi=[] + print("No Exact Sobol indices for 'typ' else than 'type3'") + self.Si = [] + self.Sij = [] + self.STi = [] # + + class fEx3D: """ Analytical test functions for 3D parameter @@ -326,67 +338,69 @@ class fEx3D: Analytical values of mean and variance of f(q) `sobol(qBound)`: Analytical Sobol indices - + Returns: `val`: 1D numpy array of size n Values of f(q) at q * If 'comp': n=n1=n2=n3 * If 'tensorProd': n=n1*n2*n3 """ - def __init__(self,q1,q2,q3,typ,method,opts): - self.q1=q1 - self.q2=q2 - self.q3=q3 - self.typ=typ - self.method=method - self.opts=opts + + def __init__(self, q1, q2, q3, typ, method, opts): + self.q1 = q1 + self.q2 = q2 + self.q3 = q3 + self.typ = typ + self.method = method + self.opts = opts self._check() self.eval() def _check(self): - method_valid=['comp','tensorProd'] + method_valid = ['comp', 'tensorProd'] if self.method not in method_valid: - raise ValueError("Invalid 'method'. Choose from ",method_valid) - typ_valid=['Ishigami'] + raise ValueError("Invalid 'method'. Choose from ", method_valid) + typ_valid = ['Ishigami'] if self.typ not in typ_valid: - raise ValueError("Invalid 'typ'. Choose from ",typ_valid) + raise ValueError("Invalid 'typ'. Choose from ", typ_valid) - def _funVal(self,z1_,z2_,z3_,typ,opts): - if typ=='Ishigami': # Ishigami Function - a=opts['a'] - b=opts['b'] - tmp=mt.sin(z1_)+a*(mt.sin(z2_))**2.+b*z3_**4.*mt.sin(z1_) + def _funVal(self, z1_, z2_, z3_, typ, opts): + if typ == 'Ishigami': # Ishigami Function + a = opts['a'] + b = opts['b'] + tmp = mt.sin(z1_)+a*(mt.sin(z2_))**2.+b*z3_**4.*mt.sin(z1_) return tmp def eval(self): z1 = np.array(self.q1, copy=False, ndmin=1) z2 = np.array(self.q2, copy=False, ndmin=1) z3 = np.array(self.q3, copy=False, ndmin=1) - n1=z1.shape[0] - n2=z2.shape[0] - n3=z3.shape[0] - f=[] - if (self.method=='tensorProd'): - for i3 in range(n3): - z3_=z3[i3] - for i2 in range(n2): - z2_=z2[i2] - for i1 in range(n1): - z1_=z1[i1] - tmp=self._funVal(z1_,z2_,z3_,self.typ,self.opts) - f.append(tmp) - elif (self.method=='comp'): - if (n1!=n2 or n1!=n3 or n2!=n3): - raise ValueError("q1,q2,q3 should have the same length for 'method':'comp'") - for i in range(n1): - z1_=z1[i] - z2_=z2[i] - z3_=z3[i] - tmp=self._funVal(z1_,z2_,z3_,self.typ,self.opts) - f.append(tmp) - self.val=np.asarray(f) - - def moments(self,qInfo): + n1 = z1.shape[0] + n2 = z2.shape[0] + n3 = z3.shape[0] + f = [] + if (self.method == 'tensorProd'): + for i3 in range(n3): + z3_ = z3[i3] + for i2 in range(n2): + z2_ = z2[i2] + for i1 in range(n1): + z1_ = z1[i1] + tmp = self._funVal(z1_, z2_, z3_, self.typ, self.opts) + f.append(tmp) + elif (self.method == 'comp'): + if (n1 != n2 or n1 != n3 or n2 != n3): + raise ValueError( + "q1,q2,q3 should have the same length for 'method':'comp'") + for i in range(n1): + z1_ = z1[i] + z2_ = z2[i] + z3_ = z3[i] + tmp = self._funVal(z1_, z2_, z3_, self.typ, self.opts) + f.append(tmp) + self.val = np.asarray(f) + + def moments(self, qInfo): """ Analytical mean and variance of f(q) @@ -400,37 +414,37 @@ def moments(self,qInfo): `var`: Variance of f(q) """ - if self.typ=='Ishigami': - #Analytical values of mean and variance of Ishigami function f(q1,q2,q3) - #Assuming q1~U[a11,a12], q3~U[q21,q22], q3~U[a31,a32], - # and q1,q2,q3 are mutually independent - a1=qInfo[0] - a2=qInfo[1] - a3=qInfo[2] - a1_=a1[1]-a1[0] #par range - a2_=a2[1]-a2[0] - a3_=a3[1]-a3[0] - m=-a2_*a3_*(mt.cos(a1[1])-mt.cos(a1[0]))+0.5*self.opts['a']*a1_*a3_*(a2_-0.5*(mt.sin(2.*a2[1])- - mt.sin(2*a2[0])))-0.2*self.opts['b']*a2_*(a3[1]**5.-a3[0]**5.)*(mt.cos(a1[1])-mt.cos(a1[0])) - m=(m/(a1_*a2_*a3_)) - - SIN2_1=(mt.sin(2*a1[1])-mt.sin(2*a1[0])) - SIN2_2=(mt.sin(2*a2[1])-mt.sin(2*a2[0])) - SIN4_2=(mt.sin(4*a2[1])-mt.sin(4*a2[0])) - v1=0.5*(a1_-0.5*SIN2_1)*(a3_+self.opts['b']**2.*(a3[1]**9.-a3[0]**9.)/9.0 + - 0.4*self.opts['b']*(a3[1]**5.-a3[0]**5.))*a2_ - v2=-self.opts['a']*(mt.cos(a1[1])-mt.cos(a1[0]))*(a2_-0.5*SIN2_2)*(a3_+0.2*self.opts['b']* - (a3[1]**5.-a3[0]**5.)) - v3=0.25*self.opts['a']**2.*a1_*a3_*(1.5*a2_-SIN2_2+0.125*SIN4_2) - v=v1+v2+v3 - v=v/(a1_*a2_*a3_) - m**2. - self.mean=m - self.var=v - - def sobol(self,qBound): + if self.typ == 'Ishigami': + # Analytical values of mean and variance of Ishigami function f(q1,q2,q3) + # Assuming q1~U[a11,a12], q3~U[q21,q22], q3~U[a31,a32], + # and q1,q2,q3 are mutually independent + a1 = qInfo[0] + a2 = qInfo[1] + a3 = qInfo[2] + a1_ = a1[1]-a1[0] # par range + a2_ = a2[1]-a2[0] + a3_ = a3[1]-a3[0] + m = -a2_*a3_*(mt.cos(a1[1])-mt.cos(a1[0]))+0.5*self.opts['a']*a1_*a3_*(a2_-0.5*(mt.sin(2.*a2[1]) - + mt.sin(2*a2[0])))-0.2*self.opts['b']*a2_*(a3[1]**5.-a3[0]**5.)*(mt.cos(a1[1])-mt.cos(a1[0])) + m = (m/(a1_*a2_*a3_)) + + SIN2_1 = (mt.sin(2*a1[1])-mt.sin(2*a1[0])) + SIN2_2 = (mt.sin(2*a2[1])-mt.sin(2*a2[0])) + SIN4_2 = (mt.sin(4*a2[1])-mt.sin(4*a2[0])) + v1 = 0.5*(a1_-0.5*SIN2_1)*(a3_+self.opts['b']**2.*(a3[1]**9.-a3[0]**9.)/9.0 + + 0.4*self.opts['b']*(a3[1]**5.-a3[0]**5.))*a2_ + v2 = -self.opts['a']*(mt.cos(a1[1])-mt.cos(a1[0]))*(a2_-0.5*SIN2_2)*(a3_+0.2*self.opts['b'] * + (a3[1]**5.-a3[0]**5.)) + v3 = 0.25*self.opts['a']**2.*a1_*a3_*(1.5*a2_-SIN2_2+0.125*SIN4_2) + v = v1+v2+v3 + v = v/(a1_*a2_*a3_) - m**2. + self.mean = m + self.var = v + + def sobol(self, qBound): """ Sobol sensitivity indices of f(q) with respect to q1, q2, and q3 - + Args: `qBound`: List of length 3 =[qBound1,qBound2,qBound3] admissible range of q1, q2, q3 @@ -441,31 +455,31 @@ def sobol(self,qBound): `Sij`: =[S12,S13,S23] dual interactions """ - pi=mt.pi - iFac=0 + pi = mt.pi + iFac = 0 for i in range(3): - if qBound[i][0]!=-pi or qBound[i][1]!=pi: - iFac=1 + if qBound[i][0] != -pi or qBound[i][1] != pi: + iFac = 1 + + if iFac == 1: + raise ValueError("qBound should be [-pi,pi] in all 3-directions.") - if iFac==1: - raise ValueError("qBound should be [-pi,pi] in all 3-directions.") - # Exact Sobol indices for Ishigami function for q_i\in[-pi,pi], i=1,2,3 - a=self.opts['a'] - b=self.opts['b'] - D1=b*pi**4./5.+b**2.*pi**8./50. + 0.5 - D2=a**2./8.0 - D3=0.0 - D13=b**2.*pi**8./50.0+7.*b**2.*pi**8./450 - D12=0.0 - D23=0.0 - D123=0.0 - - D=D1+D2+D3+D12+D13+D23+D123 - Si=[D1/D,D2/D,D3/D] - Sij=[D12/D,D13/D,D23/D] - self.Si=Si - self.Sij=Sij - STi=[Si[0]+Sij[0]+Sij[1],Si[1]+Sij[0]+Sij[2],Si[2]+Sij[1]+Sij[2]] - self.STi=STi + a = self.opts['a'] + b = self.opts['b'] + D1 = b*pi**4./5.+b**2.*pi**8./50. + 0.5 + D2 = a**2./8.0 + D3 = 0.0 + D13 = b**2.*pi**8./50.0+7.*b**2.*pi**8./450 + D12 = 0.0 + D23 = 0.0 + D123 = 0.0 + + D = D1+D2+D3+D12+D13+D23+D123 + Si = [D1/D, D2/D, D3/D] + Sij = [D12/D, D13/D, D23/D] + self.Si = Si + self.Sij = Sij + STi = [Si[0]+Sij[0]+Sij[1], Si[1]+Sij[0]+Sij[2], Si[2]+Sij[1]+Sij[2]] + self.STi = STi # diff --git a/UQit/gpr_torch.py b/UQit/gpr_torch.py index c231b4d..6e4edda 100644 --- a/UQit/gpr_torch.py +++ b/UQit/gpr_torch.py @@ -1,9 +1,9 @@ ################################################################ # Gaussian Process Regression (GPR) ################################################################# -#---------------------------------------------------------------- +# ---------------------------------------------------------------- # Saleh Rezaeiravesh, salehr@kth.se -#---------------------------------------------------------------- +# ---------------------------------------------------------------- """ Notes: 1. GPyTorch: Size of the training data cannot exceed 128. @@ -29,73 +29,84 @@ One can change the default kernel in `self.covar_module(...)` (see below) manually by modifying the source code of :code:`UQit`. -""" -#---------------------------------------------------------------- +""" +# ---------------------------------------------------------------- # import os import sys import numpy as np import math as mt from matplotlib import pyplot as plt -from mpl_toolkits.mplot3d import Axes3D #for 3d plot +from mpl_toolkits.mplot3d import Axes3D # for 3d plot import torch import gpytorch from gpytorch.likelihoods import ( - _MultitaskGaussianLikelihoodBase, - MultitaskGaussianLikelihood, - GaussianLikelihood, - _GaussianLikelihoodBase, - FixedNoiseGaussianLikelihood, - HeteroskedasticNoise, - ) + _MultitaskGaussianLikelihoodBase, + MultitaskGaussianLikelihood, + GaussianLikelihood, + _GaussianLikelihoodBase, + FixedNoiseGaussianLikelihood, + HeteroskedasticNoise, +) # + + class SingletaskGPModel(gpytorch.models.ExactGP): """ GPR for single-task output using GPyTorch, 1D input: y=f(x) in R, x in R """ + def __init__(self, train_x, train_y, likelihood): super(SingletaskGPModel, self).__init__(train_x, train_y, likelihood) self.mean_module = gpytorch.means.ConstantMean() self.covar_module = gpytorch.kernels.ScaleKernel( gpytorch.kernels.RBFKernel() ) + def forward(self, x): mean_x = self.mean_module(x) covar_x = self.covar_module(x) return gpytorch.distributions.MultivariateNormal(mean_x, covar_x) # + + class SingletaskGPModel_mIn(gpytorch.models.ExactGP): """ GPR for single-task output using GPyTorch, p-D input: y=f(x) in R, x in R^p, p>1 """ + def __init__(self, train_x, train_y, likelihood): - super(SingletaskGPModel_mIn, self).__init__(train_x, train_y, likelihood) + super(SingletaskGPModel_mIn, self).__init__( + train_x, train_y, likelihood) num_dims = train_x.size(-1) self.mean_module = gpytorch.means.ConstantMean() self.covar_module = gpytorch.kernels.ScaleKernel( - ##different length scales in different dimentions - #gpytorch.kernels.RBFKernel(ard_num_dims=num_dims) + # different length scales in different dimentions + # gpytorch.kernels.RBFKernel(ard_num_dims=num_dims) gpytorch.kernels.RBFKernel(ard_num_dims=num_dims, - lengthscale_constraint=None, - lengthscale_prior=gpytorch.priors.NormalPrior(10, 10)) #prior (optional), to make a better fit - #gpytorch.kernels.MaternKernel(nu=2.5,ard_num_dims=num_dims) - ##different length scales in different dimentions, Matern nu - #gpytorch.kernels.RBFKernel() #equal length scales in all input dimensions + lengthscale_constraint=None, + lengthscale_prior=gpytorch.priors.NormalPrior(10, 10)) # prior (optional), to make a better fit + # gpytorch.kernels.MaternKernel(nu=2.5,ard_num_dims=num_dims) + # different length scales in different dimentions, Matern nu + # gpytorch.kernels.RBFKernel() #equal length scales in all input dimensions ) + def forward(self, x): mean_x = self.mean_module(x) covar_x = self.covar_module(x) return gpytorch.distributions.MultivariateNormal(mean_x, covar_x) # + + class gpr: """ Constructs and evaluates a GPR over the space of uncertain parameters/inputs. The GPR model is - + .. math:: y=f(x)+\epsilon - + where :math:`\epsilon\sim N(M,V)` * The parameter (input) space is p-dimensional, where p=1,2,... @@ -130,121 +141,129 @@ class gpr: Attributes: `post_f`: Posterior of f(x) at `xTest`. - + `post_obs`: Predictive posterior (likelihood) at `xTest`. - + """ - def __init__(self,xTrain,yTrain,noiseV,xTest,gprOpts): - self.xTrain=xTrain - self.yTrain=yTrain.copy() - self.noiseV=noiseV - self.xTest=xTest - self.gprOpts=gprOpts + + def __init__(self, xTrain, yTrain, noiseV, xTest, gprOpts): + self.xTrain = xTrain + self.yTrain = yTrain.copy() + self.noiseV = noiseV + self.xTest = xTest + self.gprOpts = gprOpts self._info() self.train_pred() - + def _info(self): - self.p=self.xTrain.shape[-1] - self.nResp=self.yTrain.shape[-1] - self.shift=np.zeros(self.nResp) - self.scale=np.ones(self.nResp) - if 'standardizeYTrain' in self.gprOpts.keys(): - if self.gprOpts['standardizeYTrain']: - for i in range(self.nResp): - self.shift[i]=np.mean(self.yTrain[:,i]) - self.scale[i]=np.std(self.yTrain[:,i]) - self.yTrain[:,i]=(self.yTrain[:,i]-self.shift[i])/self.scale[i] - print('...... GPR: y-Training data is standardized by shifting and scaling.') - - def train_pred(self): + self.p = self.xTrain.shape[-1] + self.nResp = self.yTrain.shape[-1] + self.shift = np.zeros(self.nResp) + self.scale = np.ones(self.nResp) + if 'standardizeYTrain' in self.gprOpts.keys(): + if self.gprOpts['standardizeYTrain']: + for i in range(self.nResp): + self.shift[i] = np.mean(self.yTrain[:, i]) + self.scale[i] = np.std(self.yTrain[:, i]) + self.yTrain[:, i] = ( + self.yTrain[:, i]-self.shift[i])/self.scale[i] + print( + '...... GPR: y-Training data is standardized by shifting and scaling.') + + def train_pred(self): """ Constructor of the GPR (training and predicting at test samples) """ - if self.p==1: - self.gprTorch_1d() - elif self.p>1: - self.gprTorch_pd() + if self.p == 1: + self.gprTorch_1d() + elif self.p > 1: + self.gprTorch_pd() def gprTorch_1d(self): """ GPR for 1D uncertain parameter. Observations :math:`{(x_i,y_i)}_{i=1}^n` are assumed to be independent but their noise variance can be either the same (iid=homoscedastic) or non-identical (heteroscedastic). - """ - if self.nResp==1: #single-task (=single-response) - self.gprTorch_1d_singleTask() - else: #multi-task (=multi-response) - raise ValueError('Multitask version of GPR is not available yet!') + """ + if self.nResp == 1: # single-task (=single-response) + self.gprTorch_1d_singleTask() + else: # multi-task (=multi-response) + raise ValueError('Multitask version of GPR is not available yet!') def gprTorch_1d_singleTask(self): """ GPR for 1D uncertain parameter and single-variate response y. """ - xTrain=self.xTrain - yTrain=self.yTrain[:,0] - noiseSdev=self.noiseV - xTest=self.xTest - gprOpts=self.gprOpts - #(0) Assignments - nIter=gprOpts['nIter'] #number of iterations in optimization of hyperparameters - lr_ =gprOpts['lr'] #learning rate for the optimizaer of the hyperparameters - torch.set_printoptions(precision=8) #to avoid losing accuracy in print after converting to torch - #(1) convert numpy arrays to torch tensors - xTrain=torch.from_numpy(xTrain) - yTrain=torch.from_numpy(yTrain) - yLogVar=torch.from_numpy(np.log(noiseSdev**2.)) - - #(2) Construct GPR for noise + xTrain = self.xTrain + yTrain = self.yTrain[:, 0] + noiseSdev = self.noiseV + xTest = self.xTest + gprOpts = self.gprOpts + # (0) Assignments + # number of iterations in optimization of hyperparameters + nIter = gprOpts['nIter'] + # learning rate for the optimizaer of the hyperparameters + lr_ = gprOpts['lr'] + # to avoid losing accuracy in print after converting to torch + torch.set_printoptions(precision=8) + # (1) convert numpy arrays to torch tensors + xTrain = torch.from_numpy(xTrain) + yTrain = torch.from_numpy(yTrain) + yLogVar = torch.from_numpy(np.log(noiseSdev**2.)) + + # (2) Construct GPR for noise log_noise_model = SingletaskGPModel( - xTrain, - yLogVar, - GaussianLikelihood(), - ) - #(3) Construct GPR for f(q) + xTrain, + yLogVar, + GaussianLikelihood(), + ) + # (3) Construct GPR for f(q) # (a) Likelihood likelihood = _GaussianLikelihoodBase( - noise_covar=HeteroskedasticNoise(log_noise_model), - ) + noise_covar=HeteroskedasticNoise(log_noise_model), + ) # (b) prior GPR model model = SingletaskGPModel(xTrain, yTrain, likelihood) - #(4) Train the model + # (4) Train the model model.train().double() likelihood.train().double() - #(5) Optimize the model hyperparameters - optimizer = torch.optim.Adam([ #Adam optimizaer: https://arxiv.org/abs/1412.6980 - {'params': model.parameters()}, # Includes GaussianLikelihood parameters - ], lr=lr_) + # (5) Optimize the model hyperparameters + optimizer = torch.optim.Adam([ # Adam optimizaer: https://arxiv.org/abs/1412.6980 + # Includes GaussianLikelihood parameters + {'params': model.parameters()}, + ], lr=lr_) # "Loss" for GPs - mll: marginal log likelihood mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model) - losses=[] - lengthSc=[] + losses = [] + lengthSc = [] for i in range(nIter): - optimizer.zero_grad() - output = model(xTrain) - loss = -mll(output, yTrain, xTrain) - loss.backward() - optimizer.step() - losses.append(loss.item()) - lengthSc.append(model.covar_module.base_kernel.lengthscale.item()) - if i==0 or (i+1) % 100== 0: - print('...... GPR-hyperparameters Optimization, iter %d/%d - loss: %.3f - lengthsc: %.3f' % (i + 1, nIter, losses[-1],lengthSc[-1])) - self.loss=losses - self.lengthSc=lengthSc + optimizer.zero_grad() + output = model(xTrain) + loss = -mll(output, yTrain, xTrain) + loss.backward() + optimizer.step() + losses.append(loss.item()) + lengthSc.append(model.covar_module.base_kernel.lengthscale.item()) + if i == 0 or (i+1) % 100 == 0: + print('...... GPR-hyperparameters Optimization, iter %d/%d - loss: %.3f - lengthsc: %.3f' % + (i + 1, nIter, losses[-1], lengthSc[-1])) + self.loss = losses + self.lengthSc = lengthSc # Plot convergence of hyperparameters optimization if gprOpts['convPlot']: - self.optim_conv_plot() + self.optim_conv_plot() - #(6) Posteriors of GPR model with optimized hyperparameters + # (6) Posteriors of GPR model with optimized hyperparameters model.eval() - likelihood.eval() - #(7) Evaluate the posteriors at the test points + likelihood.eval() + # (7) Evaluate the posteriors at the test points with torch.no_grad(): - xTest=torch.from_numpy(xTest) + xTest = torch.from_numpy(xTest) post_f = model(xTest) post_obs = likelihood(post_f, xTest) - self.post_f=post_f - self.post_y=post_obs + self.post_f = post_f + self.post_y = post_obs def gprTorch_pd(self): """ @@ -252,59 +271,62 @@ def gprTorch_pd(self): Observations (X_i,Y_i) are assumed to be independent but their noise variance can be either the same (iid=homoscedastic) or different (heteroscedastic). """ - if self.nResp==1: #single-task (=single-response) - self.gprTorch_pd_singleTask() - else: #multi-task (=multi-response) - raise ValueError('Multitask version of GPR is not available yet!') + if self.nResp == 1: # single-task (=single-response) + self.gprTorch_pd_singleTask() + else: # multi-task (=multi-response) + raise ValueError('Multitask version of GPR is not available yet!') def gprTorch_pd_singleTask(self): """ GPR for p>1 uncertain parameter and single-variate response y """ - xTrain=self.xTrain - yTrain=self.yTrain[:,0] - noiseSdev=self.noiseV - xTest=self.xTest - gprOpts=self.gprOpts - p=self.p - #(0) Assignments - nIter=gprOpts['nIter'] #number of iterations in optimization of hyperparameters - lr_ =gprOpts['lr'] #learning rate for the optimizaer of the hyperparameters - torch.set_printoptions(precision=8) #to avoid losing accuracy in print after converting to torch - #(1) convert numpy arrays to torch tensors - xTrain=torch.from_numpy(xTrain) - yTrain=torch.from_numpy(yTrain) - yLogVar=torch.from_numpy(np.log(noiseSdev**2.)) - - #(2) Construct the GPR for the noise + xTrain = self.xTrain + yTrain = self.yTrain[:, 0] + noiseSdev = self.noiseV + xTest = self.xTest + gprOpts = self.gprOpts + p = self.p + # (0) Assignments + # number of iterations in optimization of hyperparameters + nIter = gprOpts['nIter'] + # learning rate for the optimizaer of the hyperparameters + lr_ = gprOpts['lr'] + # to avoid losing accuracy in print after converting to torch + torch.set_printoptions(precision=8) + # (1) convert numpy arrays to torch tensors + xTrain = torch.from_numpy(xTrain) + yTrain = torch.from_numpy(yTrain) + yLogVar = torch.from_numpy(np.log(noiseSdev**2.)) + + # (2) Construct the GPR for the noise log_noise_model = SingletaskGPModel_mIn( - xTrain, - yLogVar, - GaussianLikelihood(), - ) - #(3) Construct GPR for f(q) + xTrain, + yLogVar, + GaussianLikelihood(), + ) + # (3) Construct GPR for f(q) # (a) Likelihood likelihood = _GaussianLikelihoodBase( - noise_covar=HeteroskedasticNoise(log_noise_model), - ) - # likelihood = GaussianLikelihood(noise=noiseSdev**2.) - ##common Gaussian likelihood with no inference for heteroscedastic noise levels + noise_covar=HeteroskedasticNoise(log_noise_model), + ) + # likelihood = GaussianLikelihood(noise=noiseSdev**2.) + # common Gaussian likelihood with no inference for heteroscedastic noise levels # (b) prior GPR model model = SingletaskGPModel_mIn(xTrain, yTrain, likelihood) - #(3) Optimize the hyperparameters + # (3) Optimize the hyperparameters model.train().double() likelihood.train().double() optimizer = torch.optim.Adam([ - {'params': model.parameters()}, # Includes Likelihood parameters -# {'params': list(model.parameters()) + list(likelihood.parameters())}, - ], - lr=lr_) #lr: learning rate + {'params': model.parameters()}, # Includes Likelihood parameters + # {'params': list(model.parameters()) + list(likelihood.parameters())}, + ], + lr=lr_) # lr: learning rate # "Loss" for GPs - the marginal log likelihood mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model) - losses=[] - lengthSc=[[] for _ in range(p)] + losses = [] + lengthSc = [[] for _ in range(p)] for i in range(nIter): # Zero gradients from previous iteration optimizer.zero_grad() @@ -316,62 +338,66 @@ def gprTorch_pd_singleTask(self): # optimize optimizer.step() # info on optimization - loss_=loss.item() - lengthSc_=[] + loss_ = loss.item() + lengthSc_ = [] for j in range(p): - lengthSc_.append(model.covar_module.base_kernel.lengthscale.squeeze()[j].item()) - #lengthSC_.append(model.covar_module.base_kernel.lengthscale.item()) - ##if all lengthscales are the same (see the definition of self.covar_module, above) - if i==0 or (i+1) % 100 == 0: - print('...... GPR-hyperparameters Optimization, iter %d/%d - loss: %.3f' %((i + 1), nIter, loss_),end=" ") - print('lengthscales='+'%.3f '*p %(tuple(lengthSc_))) + lengthSc_.append( + model.covar_module.base_kernel.lengthscale.squeeze()[j].item()) + # lengthSC_.append(model.covar_module.base_kernel.lengthscale.item()) + # if all lengthscales are the same (see the definition of self.covar_module, above) + if i == 0 or (i+1) % 100 == 0: + print('...... GPR-hyperparameters Optimization, iter %d/%d - loss: %.3f' % + ((i + 1), nIter, loss_), end=" ") + print('lengthscales='+'%.3f '*p % (tuple(lengthSc_))) losses.append(loss_) for j in range(p): lengthSc[j].append(lengthSc_[j]) - self.loss=losses - self.lengthSc=lengthSc - #print('lr=',optimizer.param_groups[0]['lr']) - #print('pars',optimizer.param_groups[0]['params']) + self.loss = losses + self.lengthSc = lengthSc + # print('lr=',optimizer.param_groups[0]['lr']) + # print('pars',optimizer.param_groups[0]['params']) # Plot convergence of hyperparameters optimization if gprOpts['convPlot']: - self.optim_conv_plot() - #(4) Posteriors of GPR model with optimized hyperparameters + self.optim_conv_plot() + # (4) Posteriors of GPR model with optimized hyperparameters model.eval() likelihood.eval() - #(3) Prediction at test inputs + # (3) Prediction at test inputs with torch.no_grad(): - xTest=torch.from_numpy(xTest) + xTest = torch.from_numpy(xTest) post_f = model(xTest) post_obs = likelihood(post_f, xTest) - self.post_f=post_f - self.post_y=post_obs + self.post_f = post_f + self.post_y = post_obs def optim_conv_plot(self): """ Plot convergence of loss and length-scale during the optimization """ - plt.figure(figsize=(12,3)) + plt.figure(figsize=(12, 3)) plt.subplot(121) - plt.plot(self.loss,'-r') - plt.ylabel('Loss',fontsize=16) - plt.xlabel('Iteration',fontsize=16) + plt.plot(self.loss, '-r') + plt.ylabel('Loss', fontsize=16) + plt.xlabel('Iteration', fontsize=16) plt.xticks(fontsize=13) plt.yticks(fontsize=13) plt.subplot(122) - if self.p==1: - plt.plot(self.lengthSc,'-b') - elif self.p>1: - for j in range(self.p): - plt.plot(self.lengthSc[j],'-',label='Lengthscale, param%d'%(j+1)) - plt.ylabel('Lengthscale',fontsize=16) - plt.xlabel('Iteration',fontsize=16) + if self.p == 1: + plt.plot(self.lengthSc, '-b') + elif self.p > 1: + for j in range(self.p): + plt.plot(self.lengthSc[j], '-', + label='Lengthscale, param%d' % (j+1)) + plt.ylabel('Lengthscale', fontsize=16) + plt.xlabel('Iteration', fontsize=16) plt.xticks(fontsize=13) plt.yticks(fontsize=13) - if self.p>1: - plt.legend(loc='best',fontsize=14) + if self.p > 1: + plt.legend(loc='best', fontsize=14) plt.suptitle("Convergence of GPR's hyperparam optimization") plt.show() + class gprPost: """ Post-processing a constructed GPR @@ -399,33 +425,36 @@ class gprPost: `ciU`: pD numpy array of size (nTest1,nTest2,...,nTestp) Upper 95% CI of the GP-posterior """ - def __init__(self,gpPost,nTest,shift=0.0,scale=1.0): - self.gpPost=gpPost - self.nTest=nTest - self.shift=shift - self.scale=scale + + def __init__(self, gpPost, nTest, shift=0.0, scale=1.0): + self.gpPost = gpPost + self.nTest = nTest + self.shift = shift + self.scale = scale def torchPost(self): """ Computes mean, standard-deviation, and CI of the GPR-posterior created by GPyTorch - """ + """ with torch.no_grad(): - post_=self.gpPost - nTest=self.nTest - post_mean_=post_.mean.detach().numpy() - post_mean =post_mean_.reshape(nTest,order='F') #posterior mean - lower_, upper_ = post_.confidence_region() #\pm 2*sdev of posterior mean - lower_=lower_.detach().numpy().reshape(nTest,order='F') - upper_=upper_.detach().numpy().reshape(nTest,order='F') - post_sdev = (post_mean-lower_)/2.0 #sdev of the posterior mean of f(q) - post_sdev *= self.scale #de-standardize sdev prediction - post_mean =post_mean * self.scale + self.shift #de-standardize mean prediction + post_ = self.gpPost + nTest = self.nTest + post_mean_ = post_.mean.detach().numpy() + post_mean = post_mean_.reshape(nTest, order='F') # posterior mean + lower_, upper_ = post_.confidence_region() # \pm 2*sdev of posterior mean + lower_ = lower_.detach().numpy().reshape(nTest, order='F') + upper_ = upper_.detach().numpy().reshape(nTest, order='F') + # sdev of the posterior mean of f(q) + post_sdev = (post_mean-lower_)/2.0 + post_sdev *= self.scale # de-standardize sdev prediction + post_mean = post_mean * self.scale + self.shift # de-standardize mean prediction lower_ = post_mean-1.96*post_sdev upper_ = post_mean+1.96*post_sdev - self.mean=post_mean - self.sdev=post_sdev - self.ciL=lower_ - self.ciU=upper_ + self.mean = post_mean + self.sdev = post_sdev + self.ciL = lower_ + self.ciU = upper_ + class gprPlot: """ @@ -452,66 +481,71 @@ class gprPlot: `torch2d_3dSurf()`: 3D plot of the GPR surface (mean+CI) constructed for a 2D input. """ - def __init__(self,pltOpts={}): - self.pltOpts=pltOpts + + def __init__(self, pltOpts={}): + self.pltOpts = pltOpts self._set_vals() def _set_vals(self): - keys_=self.pltOpts.keys() #provided keys + keys_ = self.pltOpts.keys() # provided keys if 'title' in keys_: - self.title=self.pltOpts['title'] + self.title = self.pltOpts['title'] if 'xlab' in keys_: - self.xlab=self.pltOpts['xlab'] + self.xlab = self.pltOpts['xlab'] if 'ylab' in keys_: - self.ylab=self.pltOpts['ylab'] - #default values for fontsizes - self.titleFS=15 #title fs - self.legFS=15 #legend fs - self.labFS=[17,17] #axes-label fs - self.ticksFS=[18,18] #axes-label fs + self.ylab = self.pltOpts['ylab'] + # default values for fontsizes + self.titleFS = 15 # title fs + self.legFS = 15 # legend fs + self.labFS = [17, 17] # axes-label fs + self.ticksFS = [18, 18] # axes-label fs if 'titleFS' in keys_: - self.titleFS=self.pltOpts['titleFS'] + self.titleFS = self.pltOpts['titleFS'] if 'legFS' in keys_: - self.legFS=self.pltOpts['legFS'] + self.legFS = self.pltOpts['legFS'] if 'labFS' in keys_: - self.labFS=self.pltOpts['labFS'] - if len(self.labFS)!=2: - raise ValueError("Value of 'labFS' should have length 2 (x,y-axes label fontsize).") + self.labFS = self.pltOpts['labFS'] + if len(self.labFS) != 2: + raise ValueError( + "Value of 'labFS' should have length 2 (x,y-axes label fontsize).") if 'ticksFS' in keys_: - self.ticksFS=self.pltOpts['ticksFS'] - if len(self.ticksFS)!=2: - raise ValueError("Value of 'ticksFS' should have length 2 (x,y-axes label fontsize).") - #options for saving the plot - self.figSave=False + self.ticksFS = self.pltOpts['ticksFS'] + if len(self.ticksFS) != 2: + raise ValueError( + "Value of 'ticksFS' should have length 2 (x,y-axes label fontsize).") + # options for saving the plot + self.figSave = False if 'save' in keys_ and self.pltOpts['save']: - list_=['figName','figDir','figSize'] - for L_ in list_: - if L_ not in keys_: - raise KeyError("%s is required in pltOpts since 'save' is True." %L_) - else: - globals()[L_]=self.pltOpts[L_] - self.figName=figName - self.figDir=figDir - self.figSize=figSize - self.figSave=True - + list_ = ['figName', 'figDir', 'figSize'] + for L_ in list_: + if L_ not in keys_: + raise KeyError( + "%s is required in pltOpts since 'save' is True." % L_) + else: + globals()[L_] = self.pltOpts[L_] + self.figName = figName + self.figDir = figDir + self.figSize = figSize + self.figSave = True + def _figSaver(self): """ Saves figures """ fig = plt.gcf() if not os.path.exists(self.figDir): - os.makedirs(self.figDir) - figSave_=figDir+figName - figOut=figDir+figName + os.makedirs(self.figDir) + figSave_ = figDir+figName + figOut = figDir+figName DPI = fig.get_dpi() - fig.set_size_inches(self.figSize[0]/float(DPI),self.figSize[1]/float(DPI)) - plt.savefig(figOut+'.pdf',bbox_inches='tight') - - def torch1d(self,post_f,post_obs,xTrain,yTrain,xTest,fExTest,shift=0.0,scale=1.0): + fig.set_size_inches( + self.figSize[0]/float(DPI), self.figSize[1]/float(DPI)) + plt.savefig(figOut+'.pdf', bbox_inches='tight') + + def torch1d(self, post_f, post_obs, xTrain, yTrain, xTest, fExTest, shift=0.0, scale=1.0): """ Plots the GPR constructed by GPyToch for a 1D input. - + Args: `post_f`: GpyTorch object Posterior density of the model function f(q) @@ -531,45 +565,52 @@ def torch1d(self,post_f,post_obs,xTrain,yTrain,xTest,fExTest,shift=0.0,scale=1.0 Value by which the original training data 'yTrain' were scaled for standardization """ with torch.no_grad(): - lower_f, upper_f = post_f.confidence_region() - lower_obs, upper_obs = post_obs.confidence_region() - post_obs_sample_=post_obs.sample().numpy() - #de-standardize mean and sdev - yTrain_=yTrain*scale+shift - post_f_mean_ = post_f.mean[:].numpy() * scale + shift - post_obs_mean_ = post_obs.mean[:].numpy() * scale + shift - #NOTE: confidence interval = 2* sdev, see - #https://github.com/cornellius-gp/gpytorch/blob/4a1ba02d2367e4e9dd03eb1ccbfa4707da02dd08/gpytorch/distributions/multivariate_normal.py - post_f_sdev_ = ((upper_f.numpy() - lower_f.numpy())/4.0)*scale - post_obs_sdev_ = ((upper_obs.numpy() - lower_obs.numpy())/4.0)*scale - lower_f_=post_f_mean_ - 1.96 * post_f_sdev_ - upper_f_=post_f_mean_ + 1.96 * post_f_sdev_ - lower_obs_=post_obs_mean_ - 1.96 * post_obs_sdev_ - upper_obs_=post_obs_mean_ + 1.96 * post_obs_sdev_ - post_obs_sample_=post_obs_sample_*scale+shift - - plt.figure(figsize=(10,6)) - plt.plot(xTest,fExTest,'--b',label='Exact Output') - plt.plot(xTrain, yTrain_, 'ok',markersize=4,label='Training obs. y') - plt.plot(xTest, post_f_mean_, '-r',lw=2,label='Mean Model') - plt.plot(xTest, post_obs_mean_, ':m',lw=2,label='Mean Posterior Pred') - plt.plot(xTest, post_obs_sample_, '-k',lw=1,label='Sample Posterior Pred') - plt.fill_between(xTest, lower_f_, upper_f_, alpha=0.3,label='CI for f(q)') - plt.fill_between(xTest, lower_obs_, upper_obs_, alpha=0.15, color='r',label='CI for obs. y') - plt.legend(loc='best',fontsize=self.legFS) - if hasattr(self,'title'): - plt.title(self.title,fontsize=self.titleFS) - else: - plt.title('Single-task GP, 1D parameter',fontsize=self.titleFS) - plt.xticks(fontsize=self.ticksFS[0]) - plt.yticks(fontsize=self.ticksFS[1]) - plt.xlabel(r'$\mathbf{q}$',fontsize=self.labFS[0]) - plt.ylabel(r'$y$',fontsize=self.labFS[1]) - if self.figSave: + lower_f, upper_f = post_f.confidence_region() + lower_obs, upper_obs = post_obs.confidence_region() + post_obs_sample_ = post_obs.sample().numpy() + # de-standardize mean and sdev + yTrain_ = yTrain*scale+shift + post_f_mean_ = post_f.mean[:].numpy() * scale + shift + post_obs_mean_ = post_obs.mean[:].numpy() * scale + shift + # NOTE: confidence interval = 2* sdev, see + # https://github.com/cornellius-gp/gpytorch/blob/4a1ba02d2367e4e9dd03eb1ccbfa4707da02dd08/gpytorch/distributions/multivariate_normal.py + post_f_sdev_ = ((upper_f.numpy() - lower_f.numpy())/4.0)*scale + post_obs_sdev_ = ( + (upper_obs.numpy() - lower_obs.numpy())/4.0)*scale + lower_f_ = post_f_mean_ - 1.96 * post_f_sdev_ + upper_f_ = post_f_mean_ + 1.96 * post_f_sdev_ + lower_obs_ = post_obs_mean_ - 1.96 * post_obs_sdev_ + upper_obs_ = post_obs_mean_ + 1.96 * post_obs_sdev_ + post_obs_sample_ = post_obs_sample_*scale+shift + + plt.figure(figsize=(10, 6)) + plt.plot(xTest, fExTest, '--b', label='Exact Output') + plt.plot(xTrain, yTrain_, 'ok', markersize=4, + label='Training obs. y') + plt.plot(xTest, post_f_mean_, '-r', lw=2, label='Mean Model') + plt.plot(xTest, post_obs_mean_, ':m', + lw=2, label='Mean Posterior Pred') + plt.plot(xTest, post_obs_sample_, '-k', + lw=1, label='Sample Posterior Pred') + plt.fill_between(xTest, lower_f_, upper_f_, + alpha=0.3, label='CI for f(q)') + plt.fill_between(xTest, lower_obs_, upper_obs_, + alpha=0.15, color='r', label='CI for obs. y') + plt.legend(loc='best', fontsize=self.legFS) + if hasattr(self, 'title'): + plt.title(self.title, fontsize=self.titleFS) + else: + plt.title('Single-task GP, 1D parameter', + fontsize=self.titleFS) + plt.xticks(fontsize=self.ticksFS[0]) + plt.yticks(fontsize=self.ticksFS[1]) + plt.xlabel(r'$\mathbf{q}$', fontsize=self.labFS[0]) + plt.ylabel(r'$y$', fontsize=self.labFS[1]) + if self.figSave: self._figSaver() - plt.show() + plt.show() - def torch2d_2dcont(self,xTrain,qTest,fTestGrid): + def torch2d_2dcont(self, xTrain, qTest, fTestGrid): """ Planar contour plots of a GPR constructed over a 2D input space. @@ -584,24 +625,25 @@ def torch2d_2dcont(self,xTrain,qTest,fTestGrid): `fTestGrid`: 2D numpy array of shape (nTest_1,nTest_2) Response values at a tensor-product grid constructed from `qTest` """ - fig = plt.figure(figsize=(5,5)) - plt.plot(xTrain[:,0],xTrain[:,1],'or') - CS1=plt.contour(qTest[0],qTest[1],fTestGrid.T,levels=40) - plt.clabel(CS1, inline=True, fontsize=15,colors='k',fmt='%0.2f',rightside_up=True,manual=False) - plt.legend(loc='best',fontsize=self.legFS) + fig = plt.figure(figsize=(5, 5)) + plt.plot(xTrain[:, 0], xTrain[:, 1], 'or') + CS1 = plt.contour(qTest[0], qTest[1], fTestGrid.T, levels=40) + plt.clabel(CS1, inline=True, fontsize=15, colors='k', + fmt='%0.2f', rightside_up=True, manual=False) + plt.legend(loc='best', fontsize=self.legFS) plt.xticks(fontsize=self.ticksFS[0]) plt.yticks(fontsize=self.ticksFS[1]) - if hasattr(self,'xlab'): - plt.xlabel(self.xlab,fontsize=self.labFS[0]) - if hasattr(self,'ylab'): - plt.ylabel(self.ylab,fontsize=self.labFS[1]) - if hasattr(self,'title'): - plt.title(self.title,fontsize=self.titleFS) + if hasattr(self, 'xlab'): + plt.xlabel(self.xlab, fontsize=self.labFS[0]) + if hasattr(self, 'ylab'): + plt.ylabel(self.ylab, fontsize=self.labFS[1]) + if hasattr(self, 'title'): + plt.title(self.title, fontsize=self.titleFS) if self.figSave: - self._figSaver() + self._figSaver() plt.show() - - def torch2d_3dSurf(self,xTrain,yTrain,qTest,post_,shift=0.0,scale=1.0): + + def torch2d_3dSurf(self, xTrain, yTrain, qTest, post_, shift=0.0, scale=1.0): """ 3D plot of the GPR surface (mean+CI) constructed for a 2D input (parameter). @@ -620,35 +662,38 @@ def torch2d_3dSurf(self,xTrain,yTrain,qTest,post_,shift=0.0,scale=1.0): `scale`: Float scalar (optional) Value by which the original training data 'yTrain' were scaled for standardization """ - nTest=[len(qTest[i]) for i in range(len(qTest))] - #Predicted mean and variance at the test grid - fP_=gprPost(post_,nTest) + nTest = [len(qTest[i]) for i in range(len(qTest))] + # Predicted mean and variance at the test grid + fP_ = gprPost(post_, nTest) fP_.torchPost() - #de-standardized the mean and sdev - yTrain=yTrain * scale + shift - post_mean=fP_.mean * scale + shift - post_sdev=fP_.sdev * scale - lower_=post_mean-1.96*post_sdev - upper_=post_mean+1.96*post_sdev - - xTestGrid1,xTestGrid2=np.meshgrid(qTest[0],qTest[1], sparse=False, indexing='ij') - fig = plt.figure(figsize=(10,10)) + # de-standardized the mean and sdev + yTrain = yTrain * scale + shift + post_mean = fP_.mean * scale + shift + post_sdev = fP_.sdev * scale + lower_ = post_mean-1.96*post_sdev + upper_ = post_mean+1.96*post_sdev + + xTestGrid1, xTestGrid2 = np.meshgrid( + qTest[0], qTest[1], sparse=False, indexing='ij') + fig = plt.figure(figsize=(10, 10)) ax = fig.gca(projection='3d') - mean_surf = ax.plot_surface(xTestGrid1, xTestGrid2, post_mean,cmap='jet', - antialiased=True,rstride=1,cstride=1,linewidth=0,alpha=0.4) - - ax.view_init(elev=30., azim=45) #default view - - upper_surf_obs = ax.plot_wireframe(xTestGrid1, xTestGrid2, upper_, linewidth=1,alpha=0.25,color='r') - lower_surf_obs = ax.plot_wireframe(xTestGrid1, xTestGrid2, lower_, linewidth=1,alpha=0.25,color='b') - plt.plot(xTrain[:,0],xTrain[:,1],yTrain,'ok',ms='5') - if hasattr(self,'xlab'): - plt.xlabel(self.xlab,fontsize=self.labFS[0]) - if hasattr(self,'ylab'): - plt.ylabel(self.ylab,fontsize=self.labFS[1]) - if hasattr(self,'title'): - plt.title(self.title,fontsize=self.titleFS) + mean_surf = ax.plot_surface(xTestGrid1, xTestGrid2, post_mean, cmap='jet', + antialiased=True, rstride=1, cstride=1, linewidth=0, alpha=0.4) + + ax.view_init(elev=30., azim=45) # default view + + upper_surf_obs = ax.plot_wireframe( + xTestGrid1, xTestGrid2, upper_, linewidth=1, alpha=0.25, color='r') + lower_surf_obs = ax.plot_wireframe( + xTestGrid1, xTestGrid2, lower_, linewidth=1, alpha=0.25, color='b') + plt.plot(xTrain[:, 0], xTrain[:, 1], yTrain, 'ok', ms='5') + if hasattr(self, 'xlab'): + plt.xlabel(self.xlab, fontsize=self.labFS[0]) + if hasattr(self, 'ylab'): + plt.ylabel(self.ylab, fontsize=self.labFS[1]) + if hasattr(self, 'title'): + plt.title(self.title, fontsize=self.titleFS) if self.figSave: - self._figSaver() + self._figSaver() plt.show() # diff --git a/UQit/lagInt.py b/UQit/lagInt.py index 6cc9a13..7745e64 100644 --- a/UQit/lagInt.py +++ b/UQit/lagInt.py @@ -1,10 +1,10 @@ ####################################################################### -# Lagrange Interpolation is constructed based on tensor-product rule -# using the response values at a given nodal set and then it is used to +# Lagrange Interpolation is constructed based on tensor-product rule +# using the response values at a given nodal set and then it is used to # interpolate at some test points within the parameter space. ####################################################################### # Saleh Rezaeiravesh, salehr@kth.se -#----------------------------------------------------------------------- +# ----------------------------------------------------------------------- # """ Note: @@ -17,194 +17,206 @@ import UQit.reshaper as reshaper # # + + class lagInt: - """ - Lagrange interpolation over a p-D parameter space, where p=1,2,... - The interpolation order is :math:`(n_1-1)*(n_2-1)*...*(n_p-1)`, where, - :math:`n_k, k=1,2,...,p` refers to the number of training nodes in the i-th dimension of the parameter space. - - .. math:: - F(\mathbf{Q})=\sum_{k_1=1}^{n_1} ... \sum_{k_p=1}^{n_p} [fNodes(k_1,k_2,...,k_p) L_{k_1}(Q_1) L_{k_2}(Q_2) ... L_{k_p}(Q_p)] - - where, :math:`L_{k_i}(Q_i)` is the single-variate Lagrange basis in the i-th dimension. - - Args: - `qNodes`: A list of length p - `qNodes=[qNodes_1,qNodes_2,...,qNodes_p]`, - where `qNodes_k` is a 1D numpy array of `n_k` nodes in the k-th parameter dimension. - `fNodes`: A p-D numpy array of shape `(n_1,n_2,...,n_p)` (or a 1D array of size `n=n_1*n_2*...*n_p` if p>1) - Response values at the training nodes. - `qTest`: A list of length p - Containing test samples from the p-D parameter space, i.e. `qTest=[Q_0,Q_1,...,Q_{p-1}]` - where `Q_i` is a 1D numpy array of size `m_i, i=0,1,...,p-1`. - Note that `Q_i` must be in `[min(qNodes_i),max(qNodes_i)]` for `i=1,2,...,p`. - `liDict`: dict (optional) - To set different options for Lagrange interpolation over the p-D space when p>1, with the following key: - * 'testRule': The rule for treating the multi-dimensional test points with the values: - - 'tensorProd': A tensor-product grid is constructed from `Q_i`s in list `qTest`. - Hence, the total number of sample points is `m=m0*m1*m_{p-1}`. - - 'set': All `m_i` for `i=0,1,...p-1` should be equal. As a result, - a set of `m=m_i` test points is considered. - Returns: - `val`: f(qTest), Interpolated values of f(q) at `qTest` - * If p==1: 1D numpy array of size m1 - * If p>1 and 'testRule'=='tensorProd' => `val`: pD numpy array of shape `(m1,m2,...,mp)` - * If p>1 and 'testRule'=='set' => 'val': 1D numpy array of size `m1*m2*...*mp` - """ - def __init__(self,qNodes,fNodes,qTest,liDict=[]): - self.qNodes=qNodes - self.fNodes=fNodes - self.qTest=qTest - self.liDict=liDict - self._info() - self.interp() - - def _info(self): - """ - Checks the consistency in the inputs - """ - self.p=len(self.qNodes) - if self.p>1: - if 'testRule' not in self.liDict.keys(): - raise KeyError("For p>1, 'testRule' ('tensorProd' or 'set') should be set in 'liDict.") - else: - testRule=self.liDict['testRule'] - if testRule!='tensorProd' and testRule!='set': - raise ValueError("Invalid value for 'tensorProd'") - - def interp(self): - """ - Lagrange interpolation, for p=1,2, ... - """ - if self.p==1: - self.interp_1d() - elif self.p>1: - self.interp_pd() - - def basis1d(self,qNodes_,k,Q_): - """ - Constructs single-variate Lagrange bases :math:`L_k(q)` using `n` nodes `qNodes_` for - `k=0,1,...,n-1`. The bases are evaluated at `m` samples `Q_` of a single-variate parameter - - Args: - `qNodes_`: 1D numpy array of size n - The single-variate training nodal set - `Q_`: 1D numpy array of size m - Test samples - `k`: 1D numpy array of int values - Order of the polynomial bases - - Returns: - `prod`: n-by-m numpy array - Values of :math:`L_k(Q)` for `k=0,1,...,n-1` evaluated at `Q_` - """ - n=qNodes_.size - m=Q_.shape[-1] - prod=np.zeros((k.shape[-1],m)) - for k_ in k: - prod_=1.0 - for j in range(n): - if j!=k_: - prod_*=(Q_-qNodes_[j])/(qNodes_[k_]-qNodes_[j]) - prod[k_,:]=prod_ - return prod - - def interp_1d(self): - R""" - Lagrange interpolation of order (n-1) constructed over a 1D parameter space. - Bases are constructed from n nodes `qNodes` and are evaluated at test points `Q` in `[min(qNodes_),max(qNodes_)]`. - - .. math:: - F(Q)=\sum_{k=0}^{n-1} fNodes_{k} L_k(Q) - - where, Lagrange Bases L_k(q) are constructed from the nodal set. - """ - qNodes=self.qNodes[0] - fNodes=self.fNodes - Q=self.qTest[0] - if (np.all(Q>max(qNodes))): - raise ValueError('qTest cannot be larger than max(qNodes).') - if (np.all(Qnp.amax(qNodes[i],axis=0))): - raise ValueError('qTest cannot be larger than max(qNodes) in %d-th dim.'%i) - if (np.all(Q[i]2: - mulInd=[[i for i in range(p-1,-1,-1)]]*(p-1) #list of indices for tensordot - else: - mulInd=p - - for j in range(mTot): #test points - idxTest_=idxTestGrid[j] - if testRule=='tensorProd': - Lk_prod=Lk[0][:,int(idxTest_[0])] - for i in range(1,p): - Lk_prod=np.tensordot(Lk_prod,Lk[i][:,int(idxTest_[i])],0) - elif testRule=='set': - Lk_prod=Lk[0][:,int(idxTest_)] - for i in range(1,p): - Lk_prod=np.tensordot(Lk_prod,Lk[i][:,int(idxTest_)],0) - fInterp[j]=np.tensordot(Lk_prod,fNodes,mulInd) - if testRule=='tensorProd': - fInterp=fInterp.reshape(mList,order='F') - self.val=fInterp + """ + Lagrange interpolation over a p-D parameter space, where p=1,2,... + The interpolation order is :math:`(n_1-1)*(n_2-1)*...*(n_p-1)`, where, + :math:`n_k, k=1,2,...,p` refers to the number of training nodes in the i-th dimension of the parameter space. + + .. math:: + F(\mathbf{Q})=\sum_{k_1=1}^{n_1} ... \sum_{k_p=1}^{n_p} [fNodes(k_1,k_2,...,k_p) L_{k_1}(Q_1) L_{k_2}(Q_2) ... L_{k_p}(Q_p)] + + where, :math:`L_{k_i}(Q_i)` is the single-variate Lagrange basis in the i-th dimension. + + Args: + `qNodes`: A list of length p + `qNodes=[qNodes_1,qNodes_2,...,qNodes_p]`, + where `qNodes_k` is a 1D numpy array of `n_k` nodes in the k-th parameter dimension. + `fNodes`: A p-D numpy array of shape `(n_1,n_2,...,n_p)` (or a 1D array of size `n=n_1*n_2*...*n_p` if p>1) + Response values at the training nodes. + `qTest`: A list of length p + Containing test samples from the p-D parameter space, i.e. `qTest=[Q_0,Q_1,...,Q_{p-1}]` + where `Q_i` is a 1D numpy array of size `m_i, i=0,1,...,p-1`. + Note that `Q_i` must be in `[min(qNodes_i),max(qNodes_i)]` for `i=1,2,...,p`. + `liDict`: dict (optional) + To set different options for Lagrange interpolation over the p-D space when p>1, with the following key: + * 'testRule': The rule for treating the multi-dimensional test points with the values: + - 'tensorProd': A tensor-product grid is constructed from `Q_i`s in list `qTest`. + Hence, the total number of sample points is `m=m0*m1*m_{p-1}`. + - 'set': All `m_i` for `i=0,1,...p-1` should be equal. As a result, + a set of `m=m_i` test points is considered. + Returns: + `val`: f(qTest), Interpolated values of f(q) at `qTest` + * If p==1: 1D numpy array of size m1 + * If p>1 and 'testRule'=='tensorProd' => `val`: pD numpy array of shape `(m1,m2,...,mp)` + * If p>1 and 'testRule'=='set' => 'val': 1D numpy array of size `m1*m2*...*mp` + """ + + def __init__(self, qNodes, fNodes, qTest, liDict=[]): + self.qNodes = qNodes + self.fNodes = fNodes + self.qTest = qTest + self.liDict = liDict + self._info() + self.interp() + + def _info(self): + """ + Checks the consistency in the inputs + """ + self.p = len(self.qNodes) + if self.p > 1: + if 'testRule' not in self.liDict.keys(): + raise KeyError( + "For p>1, 'testRule' ('tensorProd' or 'set') should be set in 'liDict.") + else: + testRule = self.liDict['testRule'] + if testRule != 'tensorProd' and testRule != 'set': + raise ValueError("Invalid value for 'tensorProd'") + + def interp(self): + """ + Lagrange interpolation, for p=1,2, ... + """ + if self.p == 1: + self.interp_1d() + elif self.p > 1: + self.interp_pd() + + def basis1d(self, qNodes_, k, Q_): + """ + Constructs single-variate Lagrange bases :math:`L_k(q)` using `n` nodes `qNodes_` for + `k=0,1,...,n-1`. The bases are evaluated at `m` samples `Q_` of a single-variate parameter + + Args: + `qNodes_`: 1D numpy array of size n + The single-variate training nodal set + `Q_`: 1D numpy array of size m + Test samples + `k`: 1D numpy array of int values + Order of the polynomial bases + + Returns: + `prod`: n-by-m numpy array + Values of :math:`L_k(Q)` for `k=0,1,...,n-1` evaluated at `Q_` + """ + n = qNodes_.size + m = Q_.shape[-1] + prod = np.zeros((k.shape[-1], m)) + for k_ in k: + prod_ = 1.0 + for j in range(n): + if j != k_: + prod_ *= (Q_-qNodes_[j])/(qNodes_[k_]-qNodes_[j]) + prod[k_, :] = prod_ + return prod + + def interp_1d(self): + R""" + Lagrange interpolation of order (n-1) constructed over a 1D parameter space. + Bases are constructed from n nodes `qNodes` and are evaluated at test points `Q` in `[min(qNodes_),max(qNodes_)]`. + + .. math:: + F(Q)=\sum_{k=0}^{n-1} fNodes_{k} L_k(Q) + + where, Lagrange Bases L_k(q) are constructed from the nodal set. + """ + qNodes = self.qNodes[0] + fNodes = self.fNodes + Q = self.qTest[0] + if (np.all(Q > max(qNodes))): + raise ValueError('qTest cannot be larger than max(qNodes).') + if (np.all(Q < min(qNodes))): + raise ValueError('qTest cannot be smaller than min(qNodes).') + if (fNodes.size != qNodes.size): + raise ValueError('qNodes and fNodes should have the same size.') + n = fNodes.size + k = np.arange(n) + Lk = self.basis1d(qNodes, k, Q) + fInterp = np.matmul(fNodes[None, :], Lk).T + self.val = fInterp[:, 0] + + def interp_pd(self): + R""" + Lagrange interpolation of order :math:`(n_1-1)*(n_2-1)*...*(n_p-1)` constructed over a + p-D parameter space. Here, :math:`n_k, k=1,2,...,p` refers to the number of training nodes in the i-th dimension of the parameter space. + + .. math:: + F(\mathbf{Q})=\sum_{k_1=1}^{n_1} ... \sum_{k_p=1}^{n_p} [fNodes(k_1,k_2,...,k_p) L_{k_1}(Q_1) L_{k_2}(Q_2) ... L_{k_p}(Q_p)] + + where, :math:`L_{k_i}(Q_p)` is the single-variate Lagrange basis in the i-th dimension. + """ + Q = np.asarray(self.qTest) + qNodes = self.qNodes + fNodes = self.fNodes + p = self.p + nList = [] # list of the number of nodes in each dimension + mList = [] # list of the number of test points + for i in range(p): + n_ = qNodes[i].shape[0] + nList.append(n_) + mList.append(Q[i].shape[0]) + if fNodes.ndim == 1: + # NOTE: the smaller index changes slowest (fortran like) + fNodes = np.reshape(fNodes, nList, order='F') + # check the arguments + testRule = self.liDict['testRule'] + if (p != len(Q)): + raise ValueError( + 'qNodes and qTest should be of the same dimension.') + for i in range(p): + if (np.all(Q[i] > np.amax(qNodes[i], axis=0))): + raise ValueError( + 'qTest cannot be larger than max(qNodes) in %d-th dim.' % i) + if (np.all(Q[i] < np.amin(qNodes[i], axis=0))): + raise ValueError( + 'qTest cannot be smaller than min(qNodes) in %d-th dim.' % i) + if (fNodes.shape[i] != nList[i]): + raise ValueError( + 'qNodes and fNodes should be of the same size.') + # Construct and evaluate Lagrange interpolation + idxTest = [] # List of indices counting the test points + Lk = [] # List of Lagrange bases + for i in range(p): + idxTest.append(np.arange(mList[i])) + k = np.arange(nList[i]) + Lk_ = self.basis1d(qNodes[i], k, Q[i]) # Basis at the i-th dim + Lk.append(Lk_) + + if testRule == 'tensorProd': + idxTestGrid = reshaper.vecs2grid(idxTest) + elif testRule == 'set': + idxTestGrid = idxTest[0] # same for all dimensions + + mTot = idxTestGrid.shape[0] # total number of test points + fInterp = np.zeros(mTot) + if p > 2: + # list of indices for tensordot + mulInd = [[i for i in range(p-1, -1, -1)]]*(p-1) + else: + mulInd = p + + for j in range(mTot): # test points + idxTest_ = idxTestGrid[j] + if testRule == 'tensorProd': + Lk_prod = Lk[0][:, int(idxTest_[0])] + for i in range(1, p): + Lk_prod = np.tensordot( + Lk_prod, Lk[i][:, int(idxTest_[i])], 0) + elif testRule == 'set': + Lk_prod = Lk[0][:, int(idxTest_)] + for i in range(1, p): + Lk_prod = np.tensordot(Lk_prod, Lk[i][:, int(idxTest_)], 0) + fInterp[j] = np.tensordot(Lk_prod, fNodes, mulInd) + if testRule == 'tensorProd': + fInterp = fInterp.reshape(mList, order='F') + self.val = fInterp # -def lagInt_Quads2Line(fNodes,qNodes,lineDef): + + +def lagInt_Quads2Line(fNodes, qNodes, lineDef): """ Constructing a Lagrange interpolation from tensor-product nodal set in a 2D parameter space and then evaluating it at the test points located on a straight line lying in the same parameter plane. @@ -222,32 +234,34 @@ def lagInt_Quads2Line(fNodes,qNodes,lineDef): `val`: f(qTest), Interpolated values of f(q) at `qTest` 1D numpy array of size `m1` """ - p=len(qNodes) - nQ=[qNodes[0].shape[0],qNodes[1].shape[0]] - if (fNodes.ndim==1): - fNodes=np.reshape(fNodes,nQ,'F') - #limits of the parameters space - qBound=[[min(qNodes[0]),max(qNodes[0])], - [min(qNodes[1]),max(qNodes[1])]] - #(1) Check if the line is relying on the parameters plane - lineStart=lineDef['start'] - lineEnd =lineDef['end'] + p = len(qNodes) + nQ = [qNodes[0].shape[0], qNodes[1].shape[0]] + if (fNodes.ndim == 1): + fNodes = np.reshape(fNodes, nQ, 'F') + # limits of the parameters space + qBound = [[min(qNodes[0]), max(qNodes[0])], + [min(qNodes[1]), max(qNodes[1])]] + # (1) Check if the line is relying on the parameters plane + lineStart = lineDef['start'] + lineEnd = lineDef['end'] for i in range(p): - if (lineStart[i]qBound[i][1]): - print(i,lineStart[i],qBound[i][0],lineStart[i],qBound[i][1]) - raise ValueError('Test line cannot be outside of the training plane. Check lineStart in dim %d' %i) - if (lineEnd[i]qBound[i][1]): - raise ValueError('Test line cannot be outside of the training plane. Check lineEnd in dim %d' %i) - - #(2) Generate interpolation points over the line - noPtsLine=lineDef['noPtsLine'] - q1=np.linspace(lineStart[0],lineEnd[0],noPtsLine) - slope=(lineEnd[1]-lineStart[1])/(lineEnd[0]-lineStart[0]) - q2=slope*(q1-lineEnd[0])+lineEnd[1] - - #(3) Lagerange interpolation from quadratures to the points on the line - qLine=[q1,q2] - fLine=lagInt(fNodes=fNodes,qNodes=qNodes,qTest=qLine,liDict={'testRule':'set'}).val - return qLine,fLine -# + if (lineStart[i] < qBound[i][0] or lineStart[i] > qBound[i][1]): + print(i, lineStart[i], qBound[i][0], lineStart[i], qBound[i][1]) + raise ValueError( + 'Test line cannot be outside of the training plane. Check lineStart in dim %d' % i) + if (lineEnd[i] < qBound[i][0] or lineEnd[i] > qBound[i][1]): + raise ValueError( + 'Test line cannot be outside of the training plane. Check lineEnd in dim %d' % i) + + # (2) Generate interpolation points over the line + noPtsLine = lineDef['noPtsLine'] + q1 = np.linspace(lineStart[0], lineEnd[0], noPtsLine) + slope = (lineEnd[1]-lineStart[1])/(lineEnd[0]-lineStart[0]) + q2 = slope*(q1-lineEnd[0])+lineEnd[1] + # (3) Lagerange interpolation from quadratures to the points on the line + qLine = [q1, q2] + fLine = lagInt(fNodes=fNodes, qNodes=qNodes, qTest=qLine, + liDict={'testRule': 'set'}).val + return qLine, fLine +# diff --git a/UQit/linAlg.py b/UQit/linAlg.py index a1bf042..640cd6a 100644 --- a/UQit/linAlg.py +++ b/UQit/linAlg.py @@ -1,22 +1,24 @@ ########################################### -# Tools for Linear Algebra +# Tools for Linear Algebra ########################################### # Saleh Rezaeiravesh, salehr@kth.se -#------------------------------------------ +# ------------------------------------------ # import numpy as np import cvxpy as cvx # -def myLinearRegress(A,R,L_=1,solver_='OSQP',max_iter_=100000): + + +def myLinearRegress(A, R, L_=1, solver_='OSQP', max_iter_=100000): """ Solves the linear system of equations Af=R in its normalized form A'Af=A'R This solver works for uniquely-, over-, and under-determined linear systems. If system is under-determined, the compressed sensing method with L1 or L2 regularization is used. For this purpose, the library cvxpy is used (https://www.cvxpy.org). - + Args: `A`: numpy array of shape (n,K) - + `R`: numpy array of size n `L_`: int (optional) @@ -34,35 +36,38 @@ def myLinearRegress(A,R,L_=1,solver_='OSQP',max_iter_=100000): `f`: 1D numpy array of size K The solution of the linear system A.f=R """ - def linearSysSolver(A,R,L_=1): + def linearSysSolver(A, R, L_=1): """ Linear Regression """ - n=A.shape[0] #number of data - K=A.shape[1] #number of unknowns - #make the system normal - M=np.dot(A.T,A) - R=np.dot(A.T,R) + n = A.shape[0] # number of data + K = A.shape[1] # number of unknowns + # make the system normal + M = np.dot(A.T, A) + R = np.dot(A.T, R) # if (K>=n): #only under-determined system => use compressed sensing - if (K!=n): #always use compressed sensing (regularization) - f = cvx.Variable(K) - objective = cvx.Minimize(cvx.norm(f, L_)) #L1/L2-regularization - constraints = [M*f == R] - prob = cvx.Problem(objective, constraints) - if solver_ not in cvx.installed_solvers(): - raise ValueError("Chosen solver=%s is not in cvxpy list of solvers." %solver_) - if solver_ == 'OSQP': - object_value = prob.solve(solver=solver_,verbose=True,max_iter=max_iter_) - elif solver_ in ['SCS','ECOS','ECOS_BB']: - object_value = prob.solve(solver=solver_,verbose=True,max_iters=max_iter_) - else: - object_value = prob.solve(solver=solver_,verbose=True) - print('...... Compressed sensing (regularization) is done.') - print(' Min objective value=||fHat||= %g in L%d-sense.'%(object_value,L_)) - fHat=f.value + if (K != n): # always use compressed sensing (regularization) + f = cvx.Variable(K) + objective = cvx.Minimize(cvx.norm(f, L_)) # L1/L2-regularization + constraints = [M*f == R] + prob = cvx.Problem(objective, constraints) + if solver_ not in cvx.installed_solvers(): + raise ValueError( + "Chosen solver=%s is not in cvxpy list of solvers." % solver_) + if solver_ == 'OSQP': + object_value = prob.solve( + solver=solver_, verbose=True, max_iter=max_iter_) + elif solver_ in ['SCS', 'ECOS', 'ECOS_BB']: + object_value = prob.solve( + solver=solver_, verbose=True, max_iters=max_iter_) + else: + object_value = prob.solve(solver=solver_, verbose=True) + print('...... Compressed sensing (regularization) is done.') + print(' Min objective value=||fHat||= %g in L%d-sense.' % + (object_value, L_)) + fHat = f.value else: - fHat=np.linalg.solve(M,R) + fHat = np.linalg.solve(M, R) return fHat - f=linearSysSolver(A,R,1) + f = linearSysSolver(A, R, 1) return f - diff --git a/UQit/nodes.py b/UQit/nodes.py index 912c5c7..4474a09 100644 --- a/UQit/nodes.py +++ b/UQit/nodes.py @@ -2,12 +2,14 @@ # Spectral nodes from different types of polynomials ###################################################### # Saleh Rezaeiravesh, salehr@kth.se -#----------------------------------------------------- +# ----------------------------------------------------- # import numpy as np import math as mt # # + + def Clenshaw_pts(n): """ Generates Clenshaw points over range [-1,1] @@ -20,12 +22,14 @@ def Clenshaw_pts(n): `x_`: 1D numpy array of size `n` Clenshaw points """ - x=np.zeros(n) + x = np.zeros(n) for i in range(n): - x[i]=np.cos(i*mt.pi/(n-1)) - x_=x[::-1] + x[i] = np.cos(i*mt.pi/(n-1)) + x_ = x[::-1] return x_ # + + def ClenshawCurtis_pts(l): """ Generates Clenshaw-Curtis nodes at level l over range [0,1] @@ -38,18 +42,20 @@ def ClenshawCurtis_pts(l): `x`: 1D numpy array of size `n` Contains Clenshaw-Curtis points """ - if l>1: - n=2**(l-1)+1 - x=np.zeros(n) - for i in range(n): - x[i]=(1-np.cos(i*mt.pi/(n-1)))/2 + if l > 1: + n = 2**(l-1)+1 + x = np.zeros(n) + for i in range(n): + x[i] = (1-np.cos(i*mt.pi/(n-1)))/2 else: - n=1 - x=np.zeros(n) - x[0]=0.5 + n = 1 + x = np.zeros(n) + x[0] = 0.5 return x # -def gllPts(n,eps=10**-8,maxIter=1000): + + +def gllPts(n, eps=10**-8, maxIter=1000): """ Generating Gauss-Lobatto-Legendre (GLL) nodes of order n using the Newton-Raphson iteration. @@ -73,26 +79,26 @@ def gllPts(n,eps=10**-8,maxIter=1000): "Spectral Methods in Fluid Dynamics," Section 2.3. Springer-Verlag 1987. https://link.springer.com/book/10.1007/978-3-642-84108-8 """ - V=np.zeros((n,n)) #Legendre Vandermonde Matrix - #Initial guess for the nodes: Clenshaw points - xi=Clenshaw_pts(n) - iter_=0 - err=1000 - xi_old=xi - while iter_eps: - iter_+=1 - #Update the Legendre-Vandermonde matrix - V[:,0]=1. - V[:,1]=xi - for j in range(2,n): - V[:,j]=((2.*j-1)*xi*V[:,j-1] - (j-1)*V[:,j-2])/float(j) - #Newton-Raphson iteration - xi=xi_old-(xi*V[:,n-1]-V[:,n-2])/(n*V[:,n-1]) - err=max(abs(xi-xi_old).flatten()) - xi_old=xi - if (iter_>maxIter and err>eps): - print('gllPts(): max iterations reached without convergence!') - #Weights - w=2./(n*(n-1)*V[:,n-1]**2.) - return xi,w + V = np.zeros((n, n)) # Legendre Vandermonde Matrix + # Initial guess for the nodes: Clenshaw points + xi = Clenshaw_pts(n) + iter_ = 0 + err = 1000 + xi_old = xi + while iter_ < maxIter and err > eps: + iter_ += 1 + # Update the Legendre-Vandermonde matrix + V[:, 0] = 1. + V[:, 1] = xi + for j in range(2, n): + V[:, j] = ((2.*j-1)*xi*V[:, j-1] - (j-1)*V[:, j-2])/float(j) + # Newton-Raphson iteration + xi = xi_old-(xi*V[:, n-1]-V[:, n-2])/(n*V[:, n-1]) + err = max(abs(xi-xi_old).flatten()) + xi_old = xi + if (iter_ > maxIter and err > eps): + print('gllPts(): max iterations reached without convergence!') + # Weights + w = 2./(n*(n-1)*V[:, n-1]**2.) + return xi, w # diff --git a/UQit/pce.py b/UQit/pce.py index acbb19f..7d367f4 100644 --- a/UQit/pce.py +++ b/UQit/pce.py @@ -1,9 +1,9 @@ ############################################# # generalized Polynomial Chaos Expansion ############################################# -#-------------------------------------------- +# -------------------------------------------- # Saleh Rezaeiravesh, salehr@kth.se -#-------------------------------------------- +# -------------------------------------------- # import os import sys @@ -15,791 +15,822 @@ import UQit.linAlg as linAlg # # + + class pce: - R""" - Constructs non-intrusive generalized Polynomial Chaos Expansion (PCE). - The parameter space has dimension p. - We have taken n samples from the p-D parameter space and evaluated the - simulator at each sample. - The samples can be taken using any approach, it is only enough to set the options correctly. - The general aim is to estimate :math:`\hat{f}_k` in - - .. math:: - f(q)=\sum_{k=0}^K \hat{f}_k \psi_k(\xi) - - where - - .. math:: - \psi_k(\xi)=\psi_{k_1}(\xi_1)\psi_{k_2}(\xi_2)\cdots\psi_{k_p}(\xi_p) - - Args: - `fVal`: 1D numpy array of size n - Simulator's response values at n training samples. - `xi`: 2D numpy array of shape (n,p) - Training parameter samples over the mapped space. - NOTE: Always have to be provided unless `'sampleType':'GQ'` . - `nQList`: (optional) List of length p, - `nQList=[nQ1,nQ2,...,nQp]`, where `nQi`: number of samples in i-th direction - * nQList=[] (default) if p==1 or p>1 and `'truncMethod':'TO'` - * Required only if p>1 and `'truncMethod':'TP'` - `pceDict`: dict - Contains settings required for constructing the PCE. The keys are: - * `'p'`: int - Dimension of the parameter - * `'distType'`: List of length p, - The i-th value specifies the distribution type of the i-th parameter (based on the gPCE rule) - * `'sampleType'`: string - Type of parameter samples at which observations are made - - `'GQ'` (Gauss quadrature nodes) - - `' '` (other nodal sets, see `class trainSample` in `sampling.py`) - * `'pceSolveMethod'`: string - Method of solving for the PCE coefficients - - `'Projection'`: Projection method; samples have to be Gauss-quadrature nodes. - - `'Regression'`: Regression method for uniquely-, over-, and under-determined systems. - If under-determined, compressed sensing with L1/L2 regularization is automatically used. - * `'truncMethod'`: string (mandatory only if p>1) - Method of truncating the PCE - - `'TP'`: Tensor-Product method - - `'TO'`: Total-Order method - * `'LMax'`: int (optional) - Maximum order of the PCE in each of the parameter dimensions. - It is mandatory for p>1 and `'TuncMethod'=='TO'` - - `'LMax'` can be used only with `'pceSolveMethod':'Regression'` - - If p==1 and `'LMax'` is not provided, it will be assumed to be equal to n. - - If p>1 and `'LMax'` is not provided, it will be assumed to a default value. - `verbose`: bool (optional) - If True (default), info is printed about the PCE being constructed - - Attributes: - `coefs`: 1D numpy array of size K - Coefficients in the PCE - `fMean`: scalar - PCE estimation for E[f(q)] - `fVar`: scalar - PCE estimation for V[f(q)] - `kSet`: List (size K) of p-D lists, p>1 - Index set :math:`[[k_{1,1},k_{2,1},...k_{p,1}],...,[k_{1,K},k_{2,K},..,k_{p,K}]]` - If p==1, then kSet=[] - """ - def __init__(self,fVal,xi,pceDict,nQList=[],verbose=True): - self.fVal=fVal - self.xi=xi - self.nQList=nQList - self.pceDict=pceDict - self.verbose=verbose - self._info() - self._pceDict_corrector() - self.cnstrct() - - def _info(self): - """ - Checks consistency of the inputs - """ - obligKeyList=['p','distType','sampleType','pceSolveMethod'] #obligatory keys in pceDict - optKeyList=['truncMethod','LMax'] #optional keys in pceDict - self.obligKeyList=obligKeyList - self.optKeyList=optKeyList - self.LMax_def=10 #default value of LMax (if not provided) - if self.fVal.ndim >1: - raise ValueError("fVal should be a 1D numpy array of size n.") - self.availDist=['Unif','Norm'] #Currently available distributions - - def _pceDict_corrector(self): - """ - Checks and corrects `pceDict` to ensure the consistency of the options. - * For 'GQ' samples+'TP' truncation scheme: either 'Projection' or 'Regression' can be used - * For all combination of samples and truncation schemes, 'Projection' can be used to compute PCE coefficients - """ - for key_ in self.obligKeyList: - if key_ not in self.pceDict: - raise KeyError("%s is missing from pceDict." %key_) - self.p=self.pceDict['p'] - self.distType=self.pceDict['distType'] - self.sampleType=self.pceDict['sampleType'] - self.pceSolveMethod=self.pceDict['pceSolveMethod'] - - if len(self.distType)!=self.p: - raise ValueError("Provide %d 'distType' for the random parameter" %self.p) - if len(self.nQList)>0 and len(self.nQList)!=self.p: - raise ValueError("For 'truncMethod':'TP', nQList of the length p=%d is needed."%self.pceDict['p']) - if len(self.xi)>0: #if xi is provided - if self.xi.ndim!=2: - raise ValueError("Wrong dimension: xi must be a 2D numpy array.") - else: - if self.xi.shape[1]!=self.p: - raise ValueError("Second dimension of xi should be equal to p="%self.p) - - if 'truncMethod' not in self.pceDict: - if self.pceDict['p']!=1: - raise KeyError("'truncMethod' is missing from pceDict.") - else: #p==1 - if self.pceDict['pceSolveMethod']=='Projection': - if self.pceDict['sampleType'] !='GQ': - self.pceDict['pceSolveMethod']='Regression' - if self.verbose: - print("... Original 'Projection' method for PCE is replaced by 'Regression'.") - else: - if self.pceDict['p']>1: - if self.pceDict['truncMethod']=='TO': - if 'pceSolveMethod' not in self.pceDict or self.pceDict['pceSolveMethod']!='Regression': - self.pceDict['pceSolveMethod']='Regression' - if self.verbose: - print("... Original method for PCE is replaced by 'Regression'.") - if self.pceDict['truncMethod']=='TP': - if 'sampleType' not in self.pceDict or self.pceDict['sampleType']!='GQ': - self.pceDict['pceSolveMethod']='Regression' - if self.verbose: - print("... Original method for PCE is replaced by 'Regression'.") - if self.pceDict['pceSolveMethod']=='Regression' and self.pceDict['truncMethod']!='TP': - if 'LMax' not in self.pceDict: - print("WARNING in pceDict: 'LMax' should be set when Total-Order method is used.") - print("Here 'LMax' is set to default value %d" %self.LMax_def) - self.pceDict.update({'LMax':self.LMax_def}) - #Values associated to pceDict's obligatory keys - if 'LMax' in self.pceDict.keys(): - self.LMax=self.pceDict['LMax'] - if self.p>1: - self.truncMethod=self.pceDict['truncMethod'] - #Check the validity of the distType - for distType_ in self.distType: - if distType_ not in self.availDist: - raise ValueError("Invalid 'distType'! Availabe list: ",self.availDist) - - @classmethod - def mapToUnit(self,x_,xBound_): - R""" - Linearly maps `x_` in `xBound_` to `xi_` in `[-1,1]` - - Args: - `x_`: Either a scalar or a 1D numpy array - - `xBound_`: A list of length 2 specifying the range of `x_` - - Returns: - `xi_`: Mapped value of `x_` - """ - xi_=(2.*(x_-xBound_[0])/(xBound_[1]-xBound_[0])-1.) - return xi_ - - @classmethod - def mapFromUnit(self,xi_,xBound_): - R""" - Linearly maps `xi_` in `[-1,1]` to `x` in `xBound_` - - Args: - `xi_`: Either a scalar or a 1D numpy array - - `xBound_`: A list of length 2 specifying the range of `x_` - - Returns: - `x_`: Mapped value of `xi_` - """ - xi_ = np.array(xi_, copy=False, ndmin=1) - x_=(0.5*(xi_+1.0)*(xBound_[1]-xBound_[0])+xBound_[0]) - return x_ - - @classmethod - def gqPtsWts(self,n_,type_): - R""" - Gauss quadrature nodes and weights associated to distribution type `type_` - based on the gPCE rule. - - Args: - `n_`: int, - Order of the gPCE polynomial. - `type_`: string, - Distribution of the random variable according to the gPCE rule - - Returns: - `quads`: 1D numpy array of size `n_` - Gauss quadrature nodes - `weights`: 1D numpy array of size `n_` - Gauss quadrature weights - """ - if type_=='Unif': - x_=np.polynomial.legendre.leggauss(n_) - elif type_=='Norm': - x_=np.polynomial.hermite_e.hermegauss(n_) - quads=x_[0] - weights=x_[1] - return quads,weights - - @classmethod - def map_xi2q(self,xi_,xAux_,distType_): - R""" - Maps `xi_` in :math:`\Gamma` to `x_` in :math:`\mathbb{Q}`, where :math:`\Gamma` is the - mapped space corresponding to the standard gPCE. - - Args: - `xi_`: 1D numpy array of size n - Samples taken from the mapped parameter space - `distType_`: string - Distribution type of the parameter - `xAux`: List of length 2 - * xAux==xBound (admissible range) if distType_=='Unif', hence :math:`\Gamma=[-1,1]` - * xAux==[m,sdev] if disType_=='Norm', where x_~N(m,sdev^2) and :math:`\Gamma=[-\infty,\infty]` - - Returns: - `x_`: 1D numpy array of size n - Mapped parameter value in :math:`\mathbb{Q}` - """ - if distType_=='Unif': - x_=mapFromUnit(xi_,xAux_) - elif distType_=='Norm': - x_=xAux_[0]+xAux_[1]*xi_ - return x_ - - @classmethod - def basis(self,n_,xi_,distType_): - R""" - Evaluates gPCE polynomial basis of order `n_` at `xi_` points taken from the mapped - space :math:`\Gamma`. - The standard polynomials are chosen based on the gPCE rules. - - Args: - `n_`: int - Order of the basis - `xi_`: 1D numpy array of size m - Points taken from the mapped space - `distType_`: string - Distribution type of the random parameter (based on the gPCE rule) - - Returns: - `psi`: 1D numpy array of size m - Values of the gPCE basis at `xi_` - """ - if distType_=='Unif': - psi=np.polynomial.legendre.legval(xi_,[0]*n_+[1]) - elif distType_=='Norm': - psi=np.polynomial.hermite_e.hermeval(xi_,[0]*n_+[1]) - return psi - - @classmethod - def basisNorm(self,k_,distType_,nInteg=10000): - """ - Evaluates the L2-norm of the gPCE polynomial basis of order `k_` at `nInteg` points taken from the - mapped space :math:`\Gamma`. - - Args: - `k_`: int - Order of the gPCE basis - `distType_`: string - Distribution type of the random parameter (according to the gPCE rule) - `nInteg`: int (optional) - Number of points to evaluate the L2-norm integral - - Returns: - `xi_`: scalar - L2-norm of the gPCE basis of order `k_` - """ - if distType_=='Unif': - xi_=np.linspace(-1,1,nInteg) - elif distType_=='Norm': - xi_=np.linspace(-5,5,nInteg) - psi_k_=self.basis(k_,xi_,distType_) - psiNorm=np.linalg.norm(psi_k_,2) - return psiNorm - - @classmethod - def density(self,xi_,distType_): - R""" - Evaluates the PDF of the standard gPCE random variables with distribution - type `distType_` at points `xi_` taken from the mapped space :math:`\Gamma`. - - Args: - `xi_`: 1D numpy array of size n - Samples taken from the mapped space - `distType_`: string - Distribution type of the random parameter (according to the gPCE rule) - - Returns: - `pdf_`: 1D numpy array of size n - PDF of the random parameter (according to the gPCE rule) - """ - if distType_=='Unif': - pdf_=0.5*np.ones(xi_.shape[-1]) - elif distType_=='Norm': - pdf_=np.exp(-0.5*xi_**2.)/mt.sqrt(2*mt.pi) - return pdf_ - - def _gqInteg_fac(self,distType_): - """ - Multipliers for the Gauss-quadrature integration rule, when using the - weights provided by numpy. - """ - if distType_=='Unif': - fac_=0.5 - elif distType_=='Norm': - fac_=1./mt.sqrt(2*mt.pi) - return fac_ - - def cnstrct(self): - """ - Constructs a PCE over a p-D parameter space (p=1,2,3,...) - """ - if self.p==1: - self.cnstrct_1d() - elif self.p>1: - self.cnstrct_pd() - - def cnstrct_1d(self): - R""" - Constructs a PCE over a 1D parameter space, i.e. p=1. - """ - if self.sampleType=='GQ' and self.pceSolveMethod=='Projection': - self.cnstrct_GQ_1d() - else: #Regression method - if 'LMax' in self.pceDict.keys(): - self.LMax=self.pceDict['LMax'] - else: - self.LMax=len(self.fVal) - self.pceDict['LMax']=self.LMax - if self.verbose: - print("...... No 'LMax' existed, so 'LMax=n='",self.LMax) - self.cnstrct_nonGQ_1d() - - def cnstrct_GQ_1d(self): - R""" - Constructs a PCE over a 1D parameter space using Projection method with Gauss-quadrature nodes. - - Args: - `fVal`: 1D numpy array of size `n` - Simulator's response values at `n` training samples - `pceDict['distType']`: List of length 1 - =[distType1], where distType1:string specifies the distribution type of the parameter - based on the gPCE rule - """ - nQ=len(self.fVal) - distType_=self.distType[0] - xi,w=self.gqPtsWts(nQ,distType_) - K=nQ - #Find the coefficients in the expansion - fCoef=np.zeros(nQ) - sum2=[] - fac_=self._gqInteg_fac(distType_) - for k in range(K): #k-th coeff in PCE - psi_k=self.basis(k,xi,distType_) - sum1=np.sum(self.fVal[:K]*psi_k[:K]*w[:K]*fac_) - sum2_=np.sum((psi_k[:K])**2.*w[:K]*fac_) - fCoef[k]=sum1/sum2_ - sum2.append(sum2_) - #Estimate the mean and variance of f(q) - fMean=fCoef[0] - fVar=np.sum(fCoef[1:]**2.*sum2[1:]) - self.coefs=fCoef - self.fMean=fMean - self.fVar=fVar - - def cnstrct_nonGQ_1d(self): - """ - Constructs a PCE over a 1D parameter space using 'Regression' method - with arbitrary truncation `K=LMax` to compute the PCE coefficients for arbitrarily - chosen parameter samples. - - Args: + R""" + Constructs non-intrusive generalized Polynomial Chaos Expansion (PCE). + The parameter space has dimension p. + We have taken n samples from the p-D parameter space and evaluated the + simulator at each sample. + The samples can be taken using any approach, it is only enough to set the options correctly. + The general aim is to estimate :math:`\hat{f}_k` in + + .. math:: + f(q)=\sum_{k=0}^K \hat{f}_k \psi_k(\xi) + + where + + .. math:: + \psi_k(\xi)=\psi_{k_1}(\xi_1)\psi_{k_2}(\xi_2)\cdots\psi_{k_p}(\xi_p) + + Args: + `fVal`: 1D numpy array of size n + Simulator's response values at n training samples. + `xi`: 2D numpy array of shape (n,p) + Training parameter samples over the mapped space. + NOTE: Always have to be provided unless `'sampleType':'GQ'` . + `nQList`: (optional) List of length p, + `nQList=[nQ1,nQ2,...,nQp]`, where `nQi`: number of samples in i-th direction + * nQList=[] (default) if p==1 or p>1 and `'truncMethod':'TO'` + * Required only if p>1 and `'truncMethod':'TP'` + `pceDict`: dict + Contains settings required for constructing the PCE. The keys are: + * `'p'`: int + Dimension of the parameter + * `'distType'`: List of length p, + The i-th value specifies the distribution type of the i-th parameter (based on the gPCE rule) + * `'sampleType'`: string + Type of parameter samples at which observations are made + - `'GQ'` (Gauss quadrature nodes) + - `' '` (other nodal sets, see `class trainSample` in `sampling.py`) + * `'pceSolveMethod'`: string + Method of solving for the PCE coefficients + - `'Projection'`: Projection method; samples have to be Gauss-quadrature nodes. + - `'Regression'`: Regression method for uniquely-, over-, and under-determined systems. + If under-determined, compressed sensing with L1/L2 regularization is automatically used. + * `'truncMethod'`: string (mandatory only if p>1) + Method of truncating the PCE + - `'TP'`: Tensor-Product method + - `'TO'`: Total-Order method + * `'LMax'`: int (optional) + Maximum order of the PCE in each of the parameter dimensions. + It is mandatory for p>1 and `'TuncMethod'=='TO'` + - `'LMax'` can be used only with `'pceSolveMethod':'Regression'` + - If p==1 and `'LMax'` is not provided, it will be assumed to be equal to n. + - If p>1 and `'LMax'` is not provided, it will be assumed to a default value. + `verbose`: bool (optional) + If True (default), info is printed about the PCE being constructed + + Attributes: + `coefs`: 1D numpy array of size K + Coefficients in the PCE + `fMean`: scalar + PCE estimation for E[f(q)] + `fVar`: scalar + PCE estimation for V[f(q)] + `kSet`: List (size K) of p-D lists, p>1 + Index set :math:`[[k_{1,1},k_{2,1},...k_{p,1}],...,[k_{1,K},k_{2,K},..,k_{p,K}]]` + If p==1, then kSet=[] + """ + + def __init__(self, fVal, xi, pceDict, nQList=[], verbose=True): + self.fVal = fVal + self.xi = xi + self.nQList = nQList + self.pceDict = pceDict + self.verbose = verbose + self._info() + self._pceDict_corrector() + self.cnstrct() + + def _info(self): + """ + Checks consistency of the inputs + """ + obligKeyList = ['p', 'distType', 'sampleType', + 'pceSolveMethod'] # obligatory keys in pceDict + optKeyList = ['truncMethod', 'LMax'] # optional keys in pceDict + self.obligKeyList = obligKeyList + self.optKeyList = optKeyList + self.LMax_def = 10 # default value of LMax (if not provided) + if self.fVal.ndim > 1: + raise ValueError("fVal should be a 1D numpy array of size n.") + self.availDist = ['Unif', 'Norm'] # Currently available distributions + + def _pceDict_corrector(self): + """ + Checks and corrects `pceDict` to ensure the consistency of the options. + * For 'GQ' samples+'TP' truncation scheme: either 'Projection' or 'Regression' can be used + * For all combination of samples and truncation schemes, 'Projection' can be used to compute PCE coefficients + """ + for key_ in self.obligKeyList: + if key_ not in self.pceDict: + raise KeyError("%s is missing from pceDict." % key_) + self.p = self.pceDict['p'] + self.distType = self.pceDict['distType'] + self.sampleType = self.pceDict['sampleType'] + self.pceSolveMethod = self.pceDict['pceSolveMethod'] + + if len(self.distType) != self.p: + raise ValueError( + "Provide %d 'distType' for the random parameter" % self.p) + if len(self.nQList) > 0 and len(self.nQList) != self.p: + raise ValueError( + "For 'truncMethod':'TP', nQList of the length p=%d is needed." % self.pceDict['p']) + if len(self.xi) > 0: # if xi is provided + if self.xi.ndim != 2: + raise ValueError( + "Wrong dimension: xi must be a 2D numpy array.") + else: + if self.xi.shape[1] != self.p: + raise ValueError( + "Second dimension of xi should be equal to p=" % self.p) + + if 'truncMethod' not in self.pceDict: + if self.pceDict['p'] != 1: + raise KeyError("'truncMethod' is missing from pceDict.") + else: # p==1 + if self.pceDict['pceSolveMethod'] == 'Projection': + if self.pceDict['sampleType'] != 'GQ': + self.pceDict['pceSolveMethod'] = 'Regression' + if self.verbose: + print( + "... Original 'Projection' method for PCE is replaced by 'Regression'.") + else: + if self.pceDict['p'] > 1: + if self.pceDict['truncMethod'] == 'TO': + if 'pceSolveMethod' not in self.pceDict or self.pceDict['pceSolveMethod'] != 'Regression': + self.pceDict['pceSolveMethod'] = 'Regression' + if self.verbose: + print( + "... Original method for PCE is replaced by 'Regression'.") + if self.pceDict['truncMethod'] == 'TP': + if 'sampleType' not in self.pceDict or self.pceDict['sampleType'] != 'GQ': + self.pceDict['pceSolveMethod'] = 'Regression' + if self.verbose: + print( + "... Original method for PCE is replaced by 'Regression'.") + if self.pceDict['pceSolveMethod'] == 'Regression' and self.pceDict['truncMethod'] != 'TP': + if 'LMax' not in self.pceDict: + print( + "WARNING in pceDict: 'LMax' should be set when Total-Order method is used.") + print("Here 'LMax' is set to default value %d" % + self.LMax_def) + self.pceDict.update({'LMax': self.LMax_def}) + # Values associated to pceDict's obligatory keys + if 'LMax' in self.pceDict.keys(): + self.LMax = self.pceDict['LMax'] + if self.p > 1: + self.truncMethod = self.pceDict['truncMethod'] + # Check the validity of the distType + for distType_ in self.distType: + if distType_ not in self.availDist: + raise ValueError( + "Invalid 'distType'! Availabe list: ", self.availDist) + + @classmethod + def mapToUnit(self, x_, xBound_): + R""" + Linearly maps `x_` in `xBound_` to `xi_` in `[-1,1]` + + Args: + `x_`: Either a scalar or a 1D numpy array + + `xBound_`: A list of length 2 specifying the range of `x_` + + Returns: + `xi_`: Mapped value of `x_` + """ + xi_ = (2.*(x_-xBound_[0])/(xBound_[1]-xBound_[0])-1.) + return xi_ + + @classmethod + def mapFromUnit(self, xi_, xBound_): + R""" + Linearly maps `xi_` in `[-1,1]` to `x` in `xBound_` + + Args: + `xi_`: Either a scalar or a 1D numpy array + + `xBound_`: A list of length 2 specifying the range of `x_` + + Returns: + `x_`: Mapped value of `xi_` + """ + xi_ = np.array(xi_, copy=False, ndmin=1) + x_ = (0.5*(xi_+1.0)*(xBound_[1]-xBound_[0])+xBound_[0]) + return x_ + + @classmethod + def gqPtsWts(self, n_, type_): + R""" + Gauss quadrature nodes and weights associated to distribution type `type_` + based on the gPCE rule. + + Args: + `n_`: int, + Order of the gPCE polynomial. + `type_`: string, + Distribution of the random variable according to the gPCE rule + + Returns: + `quads`: 1D numpy array of size `n_` + Gauss quadrature nodes + `weights`: 1D numpy array of size `n_` + Gauss quadrature weights + """ + if type_ == 'Unif': + x_ = np.polynomial.legendre.leggauss(n_) + elif type_ == 'Norm': + x_ = np.polynomial.hermite_e.hermegauss(n_) + quads = x_[0] + weights = x_[1] + return quads, weights + + @classmethod + def map_xi2q(self, xi_, xAux_, distType_): + R""" + Maps `xi_` in :math:`\Gamma` to `x_` in :math:`\mathbb{Q}`, where :math:`\Gamma` is the + mapped space corresponding to the standard gPCE. + + Args: + `xi_`: 1D numpy array of size n + Samples taken from the mapped parameter space + `distType_`: string + Distribution type of the parameter + `xAux`: List of length 2 + * xAux==xBound (admissible range) if distType_=='Unif', hence :math:`\Gamma=[-1,1]` + * xAux==[m,sdev] if disType_=='Norm', where x_~N(m,sdev^2) and :math:`\Gamma=[-\infty,\infty]` + + Returns: + `x_`: 1D numpy array of size n + Mapped parameter value in :math:`\mathbb{Q}` + """ + if distType_ == 'Unif': + x_ = mapFromUnit(xi_, xAux_) + elif distType_ == 'Norm': + x_ = xAux_[0]+xAux_[1]*xi_ + return x_ + + @classmethod + def basis(self, n_, xi_, distType_): + R""" + Evaluates gPCE polynomial basis of order `n_` at `xi_` points taken from the mapped + space :math:`\Gamma`. + The standard polynomials are chosen based on the gPCE rules. + + Args: + `n_`: int + Order of the basis + `xi_`: 1D numpy array of size m + Points taken from the mapped space + `distType_`: string + Distribution type of the random parameter (based on the gPCE rule) + + Returns: + `psi`: 1D numpy array of size m + Values of the gPCE basis at `xi_` + """ + if distType_ == 'Unif': + psi = np.polynomial.legendre.legval(xi_, [0]*n_+[1]) + elif distType_ == 'Norm': + psi = np.polynomial.hermite_e.hermeval(xi_, [0]*n_+[1]) + return psi + + @classmethod + def basisNorm(self, k_, distType_, nInteg=10000): + """ + Evaluates the L2-norm of the gPCE polynomial basis of order `k_` at `nInteg` points taken from the + mapped space :math:`\Gamma`. + + Args: + `k_`: int + Order of the gPCE basis + `distType_`: string + Distribution type of the random parameter (according to the gPCE rule) + `nInteg`: int (optional) + Number of points to evaluate the L2-norm integral + + Returns: + `xi_`: scalar + L2-norm of the gPCE basis of order `k_` + """ + if distType_ == 'Unif': + xi_ = np.linspace(-1, 1, nInteg) + elif distType_ == 'Norm': + xi_ = np.linspace(-5, 5, nInteg) + psi_k_ = self.basis(k_, xi_, distType_) + psiNorm = np.linalg.norm(psi_k_, 2) + return psiNorm + + @classmethod + def density(self, xi_, distType_): + R""" + Evaluates the PDF of the standard gPCE random variables with distribution + type `distType_` at points `xi_` taken from the mapped space :math:`\Gamma`. + + Args: + `xi_`: 1D numpy array of size n + Samples taken from the mapped space + `distType_`: string + Distribution type of the random parameter (according to the gPCE rule) + + Returns: + `pdf_`: 1D numpy array of size n + PDF of the random parameter (according to the gPCE rule) + """ + if distType_ == 'Unif': + pdf_ = 0.5*np.ones(xi_.shape[-1]) + elif distType_ == 'Norm': + pdf_ = np.exp(-0.5*xi_**2.)/mt.sqrt(2*mt.pi) + return pdf_ + + def _gqInteg_fac(self, distType_): + """ + Multipliers for the Gauss-quadrature integration rule, when using the + weights provided by numpy. + """ + if distType_ == 'Unif': + fac_ = 0.5 + elif distType_ == 'Norm': + fac_ = 1./mt.sqrt(2*mt.pi) + return fac_ + + def cnstrct(self): + """ + Constructs a PCE over a p-D parameter space (p=1,2,3,...) + """ + if self.p == 1: + self.cnstrct_1d() + elif self.p > 1: + self.cnstrct_pd() + + def cnstrct_1d(self): + R""" + Constructs a PCE over a 1D parameter space, i.e. p=1. + """ + if self.sampleType == 'GQ' and self.pceSolveMethod == 'Projection': + self.cnstrct_GQ_1d() + else: # Regression method + if 'LMax' in self.pceDict.keys(): + self.LMax = self.pceDict['LMax'] + else: + self.LMax = len(self.fVal) + self.pceDict['LMax'] = self.LMax + if self.verbose: + print("...... No 'LMax' existed, so 'LMax=n='", self.LMax) + self.cnstrct_nonGQ_1d() + + def cnstrct_GQ_1d(self): + R""" + Constructs a PCE over a 1D parameter space using Projection method with Gauss-quadrature nodes. + + Args: `fVal`: 1D numpy array of size `n` Simulator's response values at `n` training samples - `pceDict['distType']`: List of length 1, + `pceDict['distType']`: List of length 1 =[distType1], where distType1:string specifies the distribution type of the parameter based on the gPCE rule - `xi`: 2D numpy array of size (n,1) - Training parameter samples over the mapped space - `pceDict['LMax']`: int - Maximum order of the PCE. - `LMax` is required since `'pceSolveMethod'=='Regression'`. - """ - nQ=len(self.fVal) #number of quadratures (collocation samples) - K=self.LMax #truncation in the PCE - distType_=self.distType[0] - if self.verbose: - print('...... Number of terms in PCE, K= ',K) - nData=len(self.fVal) #number of observations - if self.verbose: - print('...... Number of Data point, n= ',nData) - #(2) Find the coefficients in the expansion:Only Regression method can be used. - A=np.zeros((nData,K)) - sum2=[] - xi_aux,w_aux=self.gqPtsWts(K+1,distType_) #auxiliary GQ rule for computing gamma_k - fac_=self._gqInteg_fac(distType_) - for k in range(K): - psi_aux=self.basis(k,xi_aux,distType_) - for j in range(nData): - A[j,k]=self.basis(k,self.xi[j],distType_) - sum2.append(np.sum((psi_aux[:K+1])**2.*w_aux[:K+1]*fac_)) - #Find the PCE coeffs by Linear Regression - fCoef=linAlg.myLinearRegress(A,self.fVal) - #(3) Find the mean and variance of f(q) estimated by PCE - fMean=fCoef[0] - fVar=np.sum(fCoef[1:]**2.*sum2[1:]) - self.coefs=fCoef - self.fMean=fMean - self.fVar=fVar - - def cnstrct_pd(self): - """ - Construct a PCE over a p-D parameter space, where p>1. - """ - if self.sampleType=='GQ' and self.truncMethod=='TP': - #'GQ'+'TP': use either 'Projection' or 'Regression' - self.cnstrct_GQTP_pd() - else: #'Regression' - self.cnstrct_nonGQTP_pd() - - def cnstrct_GQTP_pd(self): - R""" - Constructs a PCE over a pD parameter space (p>1) using the following settings: - * `'sampType':'GQ'` (Gauss-Quadrature nodes) - * `'truncMethod': 'TP'` (Tensor-product) - * `'pceSolveMethod':'Projection'` or 'Regression' - """ - if self.verbose: - print('... A gPCE for a %d-D parameter space is constructed.'%self.p) - print('...... Samples in each direction are Gauss Quadrature nodes (User should check this!).') - print('...... PCE truncation method: TP') - print('...... Method of computing PCE coefficients: %s' %self.pceSolveMethod) - distType=self.distType - p=self.p - #(1) Quadrature rule - xi=[] - w=[] - fac=[] - K=1 - for i in range(p): - xi_,w_=self.gqPtsWts(self.nQList[i],distType[i]) - xi.append(xi_) - w.append(w_) - K*=self.nQList[i] - fac.append(self._gqInteg_fac(distType[i])) - if self.verbose: - print('...... Number of terms in PCE, K= ',K) - nData=len(self.fVal) #number of observations - if self.verbose: - print('...... Number of Data point, n= ',nData) - if K!=nData: - raise ValueError("K=%d is not equal to nData=%d"%(K,nData)) - #(2) Index set - kSet=[] #index set for the constructed PCE - kGlob=np.arange(K) #Global index - kLoc=kGlob.reshape(self.nQList,order='F') #Local index - for i in range(K): - k_=np.where(kLoc==kGlob[i]) - kSet_=[] - for j in range(p): - kSet_.append(k_[j][0]) - kSet.append(kSet_) - #(3) Find the coefficients in the expansion - #By default, Projection method is used (assuming samples are Gauss-Quadrature points) - fCoef=np.zeros(K) - sum2=[] - fVal_=self.fVal.reshape(self.nQList,order='F').T #global to local index - for k in range(K): - psi_k=[] - for j in range(p): - psi_k.append(self.basis(kSet[k][j],xi[j],distType[j])) - sum1 =np.matmul(fVal_,(psi_k[0]*w[0]))*fac[0] - sum2_=np.sum(psi_k[0]**2*w[0])*fac[0] - for i in range(1,p): - num_=(psi_k[i]*w[i]) - sum1=np.matmul(sum1,num_)*fac[i] - sum2_*=np.sum(psi_k[i]**2*w[i])*fac[i] - fCoef[k]=sum1/sum2_ - sum2.append(sum2_) - #(3b) Compute fCoef via Regression - if self.pceDict['pceSolveMethod']=='Regression': - xiGrid=reshaper.vecs2grid(xi) - A=np.zeros((nData,K)) - for k in range(K): - aij_=self.basis(kSet[k][0],xiGrid[:,0],distType[0]) - for i in range(1,p): - aij_*=self.basis(kSet[k][i],xiGrid[:,i],distType[i]) - A[:,k]=aij_ - fCoef=linAlg.myLinearRegress(A,self.fVal) #This is a uniquely determined system - #(4) Find the mean and variance of f(q) as estimated by PCE - fMean=fCoef[0] - fVar=np.sum(fCoef[1:]**2.*sum2[1:]) - self.coefs=fCoef - self.fMean=fMean - self.fVar=fVar - self.kSet=kSet - - def cnstrct_nonGQTP_pd(self): - R""" - Constructs a PCE over a pD parameter space (p>1), for the following settings: - * `'truncMethod'`: `'TO'` or `'TP'` - * `'pceSolveMethod'`: `'Regression'` (only allowed method) - * This method is used for any combination of `'sampleType'` and `'truncMethod'` - but `'GQ'`+`'TP'` - """ - p=self.p - distType=self.distType - xiGrid=self.xi - if self.verbose: - print('... A gPCE for a %d-D parameter space is constructed.' %p) - print('...... PCE truncation method: %s' %self.truncMethod) - print('...... Method of computing PCE coefficients: %s' %self.pceSolveMethod) - if self.truncMethod=='TO': - LMax=self.LMax #max order of polynomial in each direction - if self.verbose: - print(' with LMax=%d as the max polynomial order in each direction.' %LMax) - if self.pceSolveMethod!='Regression': - raise ValueError("only Regression method can be used for PCE with Total Order truncation method.") - #(1) Preliminaries - #Number of terms in PCE - if self.truncMethod=='TO': - K=int(mt.factorial(LMax+p)/(mt.factorial(LMax)*mt.factorial(p))) - Nmax=(LMax+1)**p - kOrderList=[LMax+1]*p - elif self.truncMethod=='TP': #'TP' needs nQList - K=np.prod(np.asarray(self.nQList)) - Nmax=K - kOrderList=self.nQList - # Quadrature rule only to compute \gamma_k - xiAux=[] - wAux=[] - fac=[] - for i in range(p): - xi_,w_=self.gqPtsWts(self.nQList[i],distType[i]) - xiAux.append(xi_) - wAux.append(w_) - fac.append(self._gqInteg_fac(distType[i])) - # Index set - kSet=[] #index set for the constructed PCE - kGlob=np.arange(Nmax) #Global index - kLoc=kGlob.reshape(kOrderList,order='F') #Local index - for i in range(Nmax): - k_=np.where(kLoc==kGlob[i]) - kSet_=[] - for j in range(p): - kSet_.append(k_[j][0]) - if (self.truncMethod=='TO' and sum(kSet_)<=LMax) or self.truncMethod=='TP': - kSet.append(kSet_) - if self.verbose: - print('...... Number of terms in PCE, K= ',K) - nData=len(self.fVal) #number of observations - if self.verbose: - print('...... Number of Data point, n= ',nData) - #(2) Find the coefficients in the expansion:Only Regression method can be used. - A=np.zeros((nData,K)) - sum2=[] - for k in range(K): - psi_k_=self.basis(kSet[k][0],xiAux[0],distType[0]) - sum2_=np.sum(psi_k_**2*wAux[0])*fac[0] - aij_=self.basis(kSet[k][0],xiGrid[:,0],distType[0]) - for i in range(1,p): - aij_*=self.basis(kSet[k][i],xiGrid[:,i],distType[i]) - psi_k_=self.basis(kSet[k][i],xiAux[i],distType[i]) - sum2_*=np.sum(psi_k_**2*wAux[i])*fac[i] - A[:,k]=aij_ - sum2.append(sum2_) - #Find the PCE coeffs by Linear Regression - fCoef=linAlg.myLinearRegress(A,self.fVal) - #(3) Find the mean and variance of f(q) as estimated by PCE - fMean=fCoef[0] - fVar=0.0 - for k in range(1,K): - fVar+=fCoef[k]*fCoef[k]*sum2[k] - self.coefs=fCoef - self.fMean=fMean - self.fVar=fVar - self.kSet=kSet + """ + nQ = len(self.fVal) + distType_ = self.distType[0] + xi, w = self.gqPtsWts(nQ, distType_) + K = nQ + # Find the coefficients in the expansion + fCoef = np.zeros(nQ) + sum2 = [] + fac_ = self._gqInteg_fac(distType_) + for k in range(K): # k-th coeff in PCE + psi_k = self.basis(k, xi, distType_) + sum1 = np.sum(self.fVal[:K]*psi_k[:K]*w[:K]*fac_) + sum2_ = np.sum((psi_k[:K])**2.*w[:K]*fac_) + fCoef[k] = sum1/sum2_ + sum2.append(sum2_) + # Estimate the mean and variance of f(q) + fMean = fCoef[0] + fVar = np.sum(fCoef[1:]**2.*sum2[1:]) + self.coefs = fCoef + self.fMean = fMean + self.fVar = fVar + + def cnstrct_nonGQ_1d(self): + """ + Constructs a PCE over a 1D parameter space using 'Regression' method + with arbitrary truncation `K=LMax` to compute the PCE coefficients for arbitrarily + chosen parameter samples. + + Args: + `fVal`: 1D numpy array of size `n` + Simulator's response values at `n` training samples + `pceDict['distType']`: List of length 1, + =[distType1], where distType1:string specifies the distribution type of the parameter + based on the gPCE rule + `xi`: 2D numpy array of size (n,1) + Training parameter samples over the mapped space + `pceDict['LMax']`: int + Maximum order of the PCE. + `LMax` is required since `'pceSolveMethod'=='Regression'`. + """ + nQ = len(self.fVal) # number of quadratures (collocation samples) + K = self.LMax # truncation in the PCE + distType_ = self.distType[0] + if self.verbose: + print('...... Number of terms in PCE, K= ', K) + nData = len(self.fVal) # number of observations + if self.verbose: + print('...... Number of Data point, n= ', nData) + # (2) Find the coefficients in the expansion:Only Regression method can be used. + A = np.zeros((nData, K)) + sum2 = [] + # auxiliary GQ rule for computing gamma_k + xi_aux, w_aux = self.gqPtsWts(K+1, distType_) + fac_ = self._gqInteg_fac(distType_) + for k in range(K): + psi_aux = self.basis(k, xi_aux, distType_) + for j in range(nData): + A[j, k] = self.basis(k, self.xi[j], distType_) + sum2.append(np.sum((psi_aux[:K+1])**2.*w_aux[:K+1]*fac_)) + # Find the PCE coeffs by Linear Regression + fCoef = linAlg.myLinearRegress(A, self.fVal) + # (3) Find the mean and variance of f(q) estimated by PCE + fMean = fCoef[0] + fVar = np.sum(fCoef[1:]**2.*sum2[1:]) + self.coefs = fCoef + self.fMean = fMean + self.fVar = fVar + + def cnstrct_pd(self): + """ + Construct a PCE over a p-D parameter space, where p>1. + """ + if self.sampleType == 'GQ' and self.truncMethod == 'TP': + # 'GQ'+'TP': use either 'Projection' or 'Regression' + self.cnstrct_GQTP_pd() + else: # 'Regression' + self.cnstrct_nonGQTP_pd() + + def cnstrct_GQTP_pd(self): + R""" + Constructs a PCE over a pD parameter space (p>1) using the following settings: + * `'sampType':'GQ'` (Gauss-Quadrature nodes) + * `'truncMethod': 'TP'` (Tensor-product) + * `'pceSolveMethod':'Projection'` or 'Regression' + """ + if self.verbose: + print('... A gPCE for a %d-D parameter space is constructed.' % self.p) + print( + '...... Samples in each direction are Gauss Quadrature nodes (User should check this!).') + print('...... PCE truncation method: TP') + print('...... Method of computing PCE coefficients: %s' % + self.pceSolveMethod) + distType = self.distType + p = self.p + # (1) Quadrature rule + xi = [] + w = [] + fac = [] + K = 1 + for i in range(p): + xi_, w_ = self.gqPtsWts(self.nQList[i], distType[i]) + xi.append(xi_) + w.append(w_) + K *= self.nQList[i] + fac.append(self._gqInteg_fac(distType[i])) + if self.verbose: + print('...... Number of terms in PCE, K= ', K) + nData = len(self.fVal) # number of observations + if self.verbose: + print('...... Number of Data point, n= ', nData) + if K != nData: + raise ValueError("K=%d is not equal to nData=%d" % (K, nData)) + # (2) Index set + kSet = [] # index set for the constructed PCE + kGlob = np.arange(K) # Global index + kLoc = kGlob.reshape(self.nQList, order='F') # Local index + for i in range(K): + k_ = np.where(kLoc == kGlob[i]) + kSet_ = [] + for j in range(p): + kSet_.append(k_[j][0]) + kSet.append(kSet_) + # (3) Find the coefficients in the expansion + # By default, Projection method is used (assuming samples are Gauss-Quadrature points) + fCoef = np.zeros(K) + sum2 = [] + # global to local index + fVal_ = self.fVal.reshape(self.nQList, order='F').T + for k in range(K): + psi_k = [] + for j in range(p): + psi_k.append(self.basis(kSet[k][j], xi[j], distType[j])) + sum1 = np.matmul(fVal_, (psi_k[0]*w[0]))*fac[0] + sum2_ = np.sum(psi_k[0]**2*w[0])*fac[0] + for i in range(1, p): + num_ = (psi_k[i]*w[i]) + sum1 = np.matmul(sum1, num_)*fac[i] + sum2_ *= np.sum(psi_k[i]**2*w[i])*fac[i] + fCoef[k] = sum1/sum2_ + sum2.append(sum2_) + # (3b) Compute fCoef via Regression + if self.pceDict['pceSolveMethod'] == 'Regression': + xiGrid = reshaper.vecs2grid(xi) + A = np.zeros((nData, K)) + for k in range(K): + aij_ = self.basis(kSet[k][0], xiGrid[:, 0], distType[0]) + for i in range(1, p): + aij_ *= self.basis(kSet[k][i], xiGrid[:, i], distType[i]) + A[:, k] = aij_ + # This is a uniquely determined system + fCoef = linAlg.myLinearRegress(A, self.fVal) + # (4) Find the mean and variance of f(q) as estimated by PCE + fMean = fCoef[0] + fVar = np.sum(fCoef[1:]**2.*sum2[1:]) + self.coefs = fCoef + self.fMean = fMean + self.fVar = fVar + self.kSet = kSet + + def cnstrct_nonGQTP_pd(self): + R""" + Constructs a PCE over a pD parameter space (p>1), for the following settings: + * `'truncMethod'`: `'TO'` or `'TP'` + * `'pceSolveMethod'`: `'Regression'` (only allowed method) + * This method is used for any combination of `'sampleType'` and `'truncMethod'` + but `'GQ'`+`'TP'` + """ + p = self.p + distType = self.distType + xiGrid = self.xi + if self.verbose: + print('... A gPCE for a %d-D parameter space is constructed.' % p) + print('...... PCE truncation method: %s' % self.truncMethod) + print('...... Method of computing PCE coefficients: %s' % + self.pceSolveMethod) + if self.truncMethod == 'TO': + LMax = self.LMax # max order of polynomial in each direction + if self.verbose: + print( + ' with LMax=%d as the max polynomial order in each direction.' % LMax) + if self.pceSolveMethod != 'Regression': + raise ValueError( + "only Regression method can be used for PCE with Total Order truncation method.") + # (1) Preliminaries + # Number of terms in PCE + if self.truncMethod == 'TO': + K = int(mt.factorial(LMax+p)/(mt.factorial(LMax)*mt.factorial(p))) + Nmax = (LMax+1)**p + kOrderList = [LMax+1]*p + elif self.truncMethod == 'TP': # 'TP' needs nQList + K = np.prod(np.asarray(self.nQList)) + Nmax = K + kOrderList = self.nQList + # Quadrature rule only to compute \gamma_k + xiAux = [] + wAux = [] + fac = [] + for i in range(p): + xi_, w_ = self.gqPtsWts(self.nQList[i], distType[i]) + xiAux.append(xi_) + wAux.append(w_) + fac.append(self._gqInteg_fac(distType[i])) + # Index set + kSet = [] # index set for the constructed PCE + kGlob = np.arange(Nmax) # Global index + kLoc = kGlob.reshape(kOrderList, order='F') # Local index + for i in range(Nmax): + k_ = np.where(kLoc == kGlob[i]) + kSet_ = [] + for j in range(p): + kSet_.append(k_[j][0]) + if (self.truncMethod == 'TO' and sum(kSet_) <= LMax) or self.truncMethod == 'TP': + kSet.append(kSet_) + if self.verbose: + print('...... Number of terms in PCE, K= ', K) + nData = len(self.fVal) # number of observations + if self.verbose: + print('...... Number of Data point, n= ', nData) + # (2) Find the coefficients in the expansion:Only Regression method can be used. + A = np.zeros((nData, K)) + sum2 = [] + for k in range(K): + psi_k_ = self.basis(kSet[k][0], xiAux[0], distType[0]) + sum2_ = np.sum(psi_k_**2*wAux[0])*fac[0] + aij_ = self.basis(kSet[k][0], xiGrid[:, 0], distType[0]) + for i in range(1, p): + aij_ *= self.basis(kSet[k][i], xiGrid[:, i], distType[i]) + psi_k_ = self.basis(kSet[k][i], xiAux[i], distType[i]) + sum2_ *= np.sum(psi_k_**2*wAux[i])*fac[i] + A[:, k] = aij_ + sum2.append(sum2_) + # Find the PCE coeffs by Linear Regression + fCoef = linAlg.myLinearRegress(A, self.fVal) + # (3) Find the mean and variance of f(q) as estimated by PCE + fMean = fCoef[0] + fVar = 0.0 + for k in range(1, K): + fVar += fCoef[k]*fCoef[k]*sum2[k] + self.coefs = fCoef + self.fMean = fMean + self.fVar = fVar + self.kSet = kSet # + + class pceEval: - R""" - Evaluates a constructed PCE at test samples taken from the parameter space. - The parameter space has dimension p. - The number of test samples is m. - - Args: - `coefs`: 1D numpy array of size K - PCE coefficients - `xi`: A list of length p - `xi=[xi_1,xi_2,..,xi_p]`, where `xi_k` is a 1D numpy array containing - `m_k` test samples taken from the mapped space of the k-th parameter. - Always a tensor-product grid of the test samples is constructed over the p-D space, - therefore, `m=m_1*m_2*...*m_p`. - `distType`: List of length p of strings, - The i-th value specifies the distribution type of the i-th parameter (based on the gPCE rule) - `kSet`: (optional, required only if p>1) List of length `K` - The index set produced when constructing the PCE with a specified truncation scheme. - :math:`kSet=[[k_{1,1},k_{2,1},...k_{p,1}],...,[k_{1,K},k_{2,K},..,k_{p,K}]]` with :math:`k_{i,j}` being integer - - Attribute: - `pceVal`: 1D numpy array of size m, - Response values predicted (interpolated) by the PCE at the test samples - """ - def __init__(self,coefs,xi,distType,kSet=[]): - self.coefs=coefs - self.xi=xi - self.distType=distType - self.kSet=kSet - self._info() - self.eval() - - def _info(self): - p=len(self.distType) - self.p=p #param dimension - K=len(self.coefs) - self.K=K #PCE truncation - self.availDist=['Unif','Norm'] #available distributions - #Check the validity of the distType - distType_=self.distType - for distType__ in distType_: - if distType__ not in self.availDist: - raise ValueError("Invalid 'distType'! Availabe list: ",self.availDist) - - def eval(self): - if self.p==1: - self.eval_1d() - elif self.p>1: - self.eval_pd() - - def eval_1d(self): - """ - Evaluates a PCE over a 1D parameter space at a set of test samples xi - taken from the mapped parameter space. - """ - fpce=[] - xi_=self.xi[0] - for i in range(xi_.size): - sum1=0.0 - for k in range(self.K): - sum1+=self.coefs[k]*pce.basis(k,xi_[i],self.distType[0]) - fpce.append(sum1) - self.pceVal=np.asarray(fpce) - - def eval_pd(self): - """ - Evaluates a PCE over a p-D (p>1) parameter space at a set of test samples xi - taken from the mapped parameter space. - """ - p=self.p - K=self.K - for k in range(K): - a_=self.coefs[k]*pce.basis(self.kSet[k][0],self.xi[0],self.distType[0]) - for i in range(1,p): - b_=pce.basis(self.kSet[k][i],self.xi[i],self.distType[i]) - a_=np.matmul(np.expand_dims(a_,axis=a_.ndim),b_[None,:]) - if k>0: - A+=a_ - else: - A=a_ - fpce=A - self.pceVal=fpce + R""" + Evaluates a constructed PCE at test samples taken from the parameter space. + The parameter space has dimension p. + The number of test samples is m. + + Args: + `coefs`: 1D numpy array of size K + PCE coefficients + `xi`: A list of length p + `xi=[xi_1,xi_2,..,xi_p]`, where `xi_k` is a 1D numpy array containing + `m_k` test samples taken from the mapped space of the k-th parameter. + Always a tensor-product grid of the test samples is constructed over the p-D space, + therefore, `m=m_1*m_2*...*m_p`. + `distType`: List of length p of strings, + The i-th value specifies the distribution type of the i-th parameter (based on the gPCE rule) + `kSet`: (optional, required only if p>1) List of length `K` + The index set produced when constructing the PCE with a specified truncation scheme. + :math:`kSet=[[k_{1,1},k_{2,1},...k_{p,1}],...,[k_{1,K},k_{2,K},..,k_{p,K}]]` with :math:`k_{i,j}` being integer + + Attribute: + `pceVal`: 1D numpy array of size m, + Response values predicted (interpolated) by the PCE at the test samples + """ + + def __init__(self, coefs, xi, distType, kSet=[]): + self.coefs = coefs + self.xi = xi + self.distType = distType + self.kSet = kSet + self._info() + self.eval() + + def _info(self): + p = len(self.distType) + self.p = p # param dimension + K = len(self.coefs) + self.K = K # PCE truncation + self.availDist = ['Unif', 'Norm'] # available distributions + # Check the validity of the distType + distType_ = self.distType + for distType__ in distType_: + if distType__ not in self.availDist: + raise ValueError( + "Invalid 'distType'! Availabe list: ", self.availDist) + + def eval(self): + if self.p == 1: + self.eval_1d() + elif self.p > 1: + self.eval_pd() + + def eval_1d(self): + """ + Evaluates a PCE over a 1D parameter space at a set of test samples xi + taken from the mapped parameter space. + """ + fpce = [] + xi_ = self.xi[0] + for i in range(xi_.size): + sum1 = 0.0 + for k in range(self.K): + sum1 += self.coefs[k]*pce.basis(k, xi_[i], self.distType[0]) + fpce.append(sum1) + self.pceVal = np.asarray(fpce) + + def eval_pd(self): + """ + Evaluates a PCE over a p-D (p>1) parameter space at a set of test samples xi + taken from the mapped parameter space. + """ + p = self.p + K = self.K + for k in range(K): + a_ = self.coefs[k]*pce.basis(self.kSet[k] + [0], self.xi[0], self.distType[0]) + for i in range(1, p): + b_ = pce.basis(self.kSet[k][i], self.xi[i], self.distType[i]) + a_ = np.matmul(np.expand_dims(a_, axis=a_.ndim), b_[None, :]) + if k > 0: + A += a_ + else: + A = a_ + fpce = A + self.pceVal = fpce + class convPlot: - R""" - Computes and plots the convergence indicator of a PCE that is defined as, - - .. math:: - \vartheta_k = ||f_k \Psi_k||/|f_0| - - versus :math:`|k|=\sum_{i=1}^p k_i`. - The dimension of the parameter space, p, is arbitrary. But for p>1, kSet have to - be provided. - - Args: - `coefs`: 1D numpy array of length K - The PCE coefficients - `distType`: List of length p of strings, - The i-th value specifies the distribution type of the i-th parameter (based on the gPCE rule) - `kSet`: (optional, required only if p>1) List of length `K` - The index set produced when constructing the PCE with a specified truncation scheme. - :math:`kSet=[[k_{1,1},k_{2,1},...k_{p,1}],...,[k_{1,K},k_{2,K},..,k_{p,K}]]` with :math:`k_{i,j}` being integer. - `convPltOpts`: (optional) dict - Containing the options to save the figure. It includes the following keys: - * 'figDir': - Path to the directory at which the figure is saved (if not exists, is created) - * 'figName': - Name of the figure - - Attributes: - `kMag`: List of K integers - `=|k|`, sum of the PCE uni-directional indices - `pceConvIndic`: 1D numpy array of size K - The PCE convergence indicator - - Methods: - `pceConv()`: - Computes the convergence indicator - - `pceConvPlot()`: - Plots the PCE convergence indicator - """ - def __init__(self,coefs,distType,kSet=[],convPltOpts=[]): - self.coefs=coefs - self.distType=distType - self.kSet=kSet - self.convPltOpts=convPltOpts - self._get_info() - self.pceConv() - self.pceConvPlot() - - def _get_info(self): - if len(self.kSet)==0: - p=1 - else: - p=len(self.distType) - self.p=p #param dimension - K=len(self.coefs) - self.K=K #PCE truncation - self.availDist=['Unif','Norm'] #available distributions - #Check the validity of the distType - if self.p==1: - distType_=[self.distType] - else: - distType_=self.distType - for distType__ in distType_: - if distType__ not in self.availDist: - raise ValueError("Invalid 'distType'! Availabe list: ",self.availDist) - - def pceConv(self): - """ - Computes the convergence indicator - """ - K=self.K - if self.p==1: - distType_=[self.distType] - kSet_=[[i] for i in range(K)] - else: - kSet_=self.kSet - distType_=self.distType - #magnitude of the pce indices - kMag=[sum(kSet_[i]) for i in range(K)] - #compute norm of the PCE bases - termNorm=[] - for ik in range(K): #over PCE terms - PsiNorm=1.0 - for ip in range(self.p): #over parameter dimension - k_=kSet_[ik][ip] - PsiNorm*=pce.basisNorm(k_,distType_[ip]) - termNorm.append(abs(self.coefs[ik])*PsiNorm) - termNorm0=termNorm[0] - self.kMag=kMag - self.pceConvIndic=termNorm/termNorm0 - - def pceConvPlot(self): - """ - Plots the PCE convergence indicator - """ - plt.figure(figsize=(10,5)) - plt.semilogy(self.kMag,self.pceConvIndic,'ob',fillstyle='none') - plt.xlabel(r'$|\mathbf{k}|$',fontsize=18) - plt.ylabel(r'$|\hat{f}_\mathbf{k}|\, ||\Psi_{\mathbf{k}(\mathbf{\xi})}||_2/|\hat{f}_0|$',fontsize=18) - plt.xticks(ticks=self.kMag,fontsize=17) - plt.yticks(fontsize=17) - plt.grid(alpha=0.3) - if self.convPltOpts: - if 'ylim' in self.convPltOpts: - plt.ylim(self.convPltOpts['ylim']) - fig = plt.gcf() - DPI = fig.get_dpi() - fig.set_size_inches(800/float(DPI),400/float(DPI)) - figDir=self.convPltOpts['figDir'] - if not os.path.exists(figDir): - os.makedirs(figDir) - outName=self.convPltOpts['figName'] - plt.savefig(figDir+outName+'.pdf',bbox_inches='tight') - plt.show() - else: - plt.show() -# + R""" + Computes and plots the convergence indicator of a PCE that is defined as, + + .. math:: + \vartheta_k = ||f_k \Psi_k||/|f_0| + + versus :math:`|k|=\sum_{i=1}^p k_i`. + The dimension of the parameter space, p, is arbitrary. But for p>1, kSet have to + be provided. + + Args: + `coefs`: 1D numpy array of length K + The PCE coefficients + `distType`: List of length p of strings, + The i-th value specifies the distribution type of the i-th parameter (based on the gPCE rule) + `kSet`: (optional, required only if p>1) List of length `K` + The index set produced when constructing the PCE with a specified truncation scheme. + :math:`kSet=[[k_{1,1},k_{2,1},...k_{p,1}],...,[k_{1,K},k_{2,K},..,k_{p,K}]]` with :math:`k_{i,j}` being integer. + `convPltOpts`: (optional) dict + Containing the options to save the figure. It includes the following keys: + * 'figDir': + Path to the directory at which the figure is saved (if not exists, is created) + * 'figName': + Name of the figure + + Attributes: + `kMag`: List of K integers + `=|k|`, sum of the PCE uni-directional indices + `pceConvIndic`: 1D numpy array of size K + The PCE convergence indicator + + Methods: + `pceConv()`: + Computes the convergence indicator + + `pceConvPlot()`: + Plots the PCE convergence indicator + """ + + def __init__(self, coefs, distType, kSet=[], convPltOpts=[]): + self.coefs = coefs + self.distType = distType + self.kSet = kSet + self.convPltOpts = convPltOpts + self._get_info() + self.pceConv() + self.pceConvPlot() + + def _get_info(self): + if len(self.kSet) == 0: + p = 1 + else: + p = len(self.distType) + self.p = p # param dimension + K = len(self.coefs) + self.K = K # PCE truncation + self.availDist = ['Unif', 'Norm'] # available distributions + # Check the validity of the distType + if self.p == 1: + distType_ = [self.distType] + else: + distType_ = self.distType + for distType__ in distType_: + if distType__ not in self.availDist: + raise ValueError( + "Invalid 'distType'! Availabe list: ", self.availDist) + + def pceConv(self): + """ + Computes the convergence indicator + """ + K = self.K + if self.p == 1: + distType_ = [self.distType] + kSet_ = [[i] for i in range(K)] + else: + kSet_ = self.kSet + distType_ = self.distType + # magnitude of the pce indices + kMag = [sum(kSet_[i]) for i in range(K)] + # compute norm of the PCE bases + termNorm = [] + for ik in range(K): # over PCE terms + PsiNorm = 1.0 + for ip in range(self.p): # over parameter dimension + k_ = kSet_[ik][ip] + PsiNorm *= pce.basisNorm(k_, distType_[ip]) + termNorm.append(abs(self.coefs[ik])*PsiNorm) + termNorm0 = termNorm[0] + self.kMag = kMag + self.pceConvIndic = termNorm/termNorm0 + + def pceConvPlot(self): + """ + Plots the PCE convergence indicator + """ + plt.figure(figsize=(10, 5)) + plt.semilogy(self.kMag, self.pceConvIndic, 'ob', fillstyle='none') + plt.xlabel(r'$|\mathbf{k}|$', fontsize=18) + plt.ylabel( + r'$|\hat{f}_\mathbf{k}|\, ||\Psi_{\mathbf{k}(\mathbf{\xi})}||_2/|\hat{f}_0|$', fontsize=18) + plt.xticks(ticks=self.kMag, fontsize=17) + plt.yticks(fontsize=17) + plt.grid(alpha=0.3) + if self.convPltOpts: + if 'ylim' in self.convPltOpts: + plt.ylim(self.convPltOpts['ylim']) + fig = plt.gcf() + DPI = fig.get_dpi() + fig.set_size_inches(800/float(DPI), 400/float(DPI)) + figDir = self.convPltOpts['figDir'] + if not os.path.exists(figDir): + os.makedirs(figDir) + outName = self.convPltOpts['figName'] + plt.savefig(figDir+outName+'.pdf', bbox_inches='tight') + plt.show() + else: + plt.show() +# diff --git a/UQit/ppce.py b/UQit/ppce.py index 6e22fcd..c49404c 100644 --- a/UQit/ppce.py +++ b/UQit/ppce.py @@ -1,9 +1,9 @@ ############################################################### # Probabilistic generalized Polynomial Chaos Expansion (PPCE) ############################################################### -#-------------------------------------------------------------- +# -------------------------------------------------------------- # Saleh Rezaeiravesh, salehr@kth.se -#-------------------------------------------------------------- +# -------------------------------------------------------------- # import os import sys @@ -15,6 +15,8 @@ import UQit.sampling as sampling # # + + class ppce: """ Probabilistic Polynomial Chaos Expansion (PPCE) over a p-D parameter space, @@ -49,7 +51,7 @@ class ppce: * 'convPlot_gpr': bool If true, values of the hyper-parameters are plotted vs. iteration during the optimization. - + Attributes: `fMean_list`: 1D numpy array of size `nMC` PCE estimates for the mean of f(q) @@ -64,176 +66,186 @@ class ppce: `ppceDict['nGQtest'][i]` containing the GQ test samples in the i-th direction of the parameter. """ - def __init__(self,qTrain,yTrain,noiseV,ppceDict): - self.qTrain=qTrain - self.yTrain=yTrain - self.noiseV=noiseV - self.ppceDict=ppceDict + + def __init__(self, qTrain, yTrain, noiseV, ppceDict): + self.qTrain = qTrain + self.yTrain = yTrain + self.noiseV = noiseV + self.ppceDict = ppceDict self.info() self.cnstrct() def info(self): - if self.qTrain.ndim==1: - self.p=1 - elif self.qTrain.ndim==2: - self.p=self.qTrain.shape[-1] + if self.qTrain.ndim == 1: + self.p = 1 + elif self.qTrain.ndim == 2: + self.p = self.qTrain.shape[-1] - if self.qTrain.shape[0]==self.yTrain.shape: - raise ValueError("Size of qTrain and yTrain should be the same.") + if self.qTrain.shape[0] == self.yTrain.shape: + raise ValueError("Size of qTrain and yTrain should be the same.") - obligKeys=['nGQtest','qInfo','nMC','nIter_gpr','lr_gpr','convPlot_gpr'] + obligKeys = ['nGQtest', 'qInfo', 'nMC', + 'nIter_gpr', 'lr_gpr', 'convPlot_gpr'] for key_ in obligKeys: if key_ not in self.ppceDict.keys(): - raise KeyError("%s is required in ppceDict." %key_) - + raise KeyError("%s is required in ppceDict." % key_) + def cnstrct(self): - if self.p==1: - self.ppce_cnstrct_1d() - elif self.p>1: - self.ppce_cnstrct_pd() + if self.p == 1: + self.ppce_cnstrct_1d() + elif self.p > 1: + self.ppce_cnstrct_pd() def ppce_cnstrct_1d(self): - """ - Constructing a probabilistic PCE over a 1D parameter space - """ - print('... Probabilistic PCE for 1D input parameter.') - p=self.p - ppceDict=self.ppceDict - qTrain=self.qTrain - yTrain=self.yTrain - noiseSdev=self.noiseV - #(0) Assignments - nGQ=ppceDict['nGQtest'] - qInfo=ppceDict['qInfo'] - nMC=ppceDict['nMC'] - nw_=int(nMC/10) - distType=ppceDict['distType'] - #Make a dict for GPR - gprOpts={'nIter':ppceDict['nIter_gpr'], - 'lr':ppceDict['lr_gpr'], - 'convPlot':ppceDict['convPlot_gpr'] - } - standardizeYTrain_=False - if 'standardizeYTrain_gpr' in ppceDict.keys(): - gprOpts.update({'standardizeYTrain':ppceDict['standardizeYTrain_gpr']}) - standardizeYTrain_=True - - #(1) Generate test points that are Gauss quadratures chosen based on the - # distribution of q (gPCE rule) - sampsGQ=sampling.trainSample(sampleType='GQ',GQdistType=distType,qInfo=qInfo,nSamp=nGQ) - qTest=sampsGQ.q - - #(2) Construct GPR surrogate based on training data - gpr_=gpr_torch.gpr(qTrain[:,None],yTrain[:,None],noiseSdev,qTest[:,None],gprOpts) - post_f=gpr_.post_f - post_obs=gpr_.post_y - shift_=0.0 - scale_=1.0 - if standardizeYTrain_: - shift_=gpr_.shift[0] #0: single-response - scale_=gpr_.scale[0] - - #(3) Use samples of GPR tested at GQ nodes to construct a PCE - fMean_list=[] - fVar_list =[] - pceDict={'p':p,'sampleType':'GQ','pceSolveMethod':'Projection','distType':[distType]} - for j in range(nMC): - # Draw a sample for f(q) from GPR surrogate - f_=post_obs.sample().numpy()*scale_+shift_ - # Construct PCE for the drawn sample - pce_=pce(fVal=f_,xi=[],pceDict=pceDict,verbose=False) #GP+TP - fMean_list.append(pce_.fMean) - fVar_list.append(pce_.fVar) - if ((j+1)%nw_==0): - print("...... ppce repetition for finding samples of the PCE coefficients, iter = %d/%d" - %(j+1,nMC)) - - #(4) Outputs - fMean_list=np.asarray(fMean_list) - fVar_list=np.asarray(fVar_list) - #Optional outputs: only used for gprPlot - optOut={'post_f':post_f,'post_obs':post_obs,'qTest':[qTest]} - self.fMean_samps=fMean_list - self.fVar_samps=fVar_list - self.optOut=optOut + """ + Constructing a probabilistic PCE over a 1D parameter space + """ + print('... Probabilistic PCE for 1D input parameter.') + p = self.p + ppceDict = self.ppceDict + qTrain = self.qTrain + yTrain = self.yTrain + noiseSdev = self.noiseV + # (0) Assignments + nGQ = ppceDict['nGQtest'] + qInfo = ppceDict['qInfo'] + nMC = ppceDict['nMC'] + nw_ = int(nMC/10) + distType = ppceDict['distType'] + # Make a dict for GPR + gprOpts = {'nIter': ppceDict['nIter_gpr'], + 'lr': ppceDict['lr_gpr'], + 'convPlot': ppceDict['convPlot_gpr'] + } + standardizeYTrain_ = False + if 'standardizeYTrain_gpr' in ppceDict.keys(): + gprOpts.update( + {'standardizeYTrain': ppceDict['standardizeYTrain_gpr']}) + standardizeYTrain_ = True + + # (1) Generate test points that are Gauss quadratures chosen based on the + # distribution of q (gPCE rule) + sampsGQ = sampling.trainSample( + sampleType='GQ', GQdistType=distType, qInfo=qInfo, nSamp=nGQ) + qTest = sampsGQ.q + + # (2) Construct GPR surrogate based on training data + gpr_ = gpr_torch.gpr( + qTrain[:, None], yTrain[:, None], noiseSdev, qTest[:, None], gprOpts) + post_f = gpr_.post_f + post_obs = gpr_.post_y + shift_ = 0.0 + scale_ = 1.0 + if standardizeYTrain_: + shift_ = gpr_.shift[0] # 0: single-response + scale_ = gpr_.scale[0] + + # (3) Use samples of GPR tested at GQ nodes to construct a PCE + fMean_list = [] + fVar_list = [] + pceDict = {'p': p, 'sampleType': 'GQ', + 'pceSolveMethod': 'Projection', 'distType': [distType]} + for j in range(nMC): + # Draw a sample for f(q) from GPR surrogate + f_ = post_obs.sample().numpy()*scale_+shift_ + # Construct PCE for the drawn sample + pce_ = pce(fVal=f_, xi=[], pceDict=pceDict, verbose=False) # GP+TP + fMean_list.append(pce_.fMean) + fVar_list.append(pce_.fVar) + if ((j+1) % nw_ == 0): + print("...... ppce repetition for finding samples of the PCE coefficients, iter = %d/%d" + % (j+1, nMC)) + + # (4) Outputs + fMean_list = np.asarray(fMean_list) + fVar_list = np.asarray(fVar_list) + # Optional outputs: only used for gprPlot + optOut = {'post_f': post_f, 'post_obs': post_obs, 'qTest': [qTest]} + self.fMean_samps = fMean_list + self.fVar_samps = fVar_list + self.optOut = optOut # + def ppce_cnstrct_pd(self): - """ - Constructing a probabilistic PCE over a p-D parameter space, p>1 - """ - p=self.p - print('... Probabilistic PCE for %d-D input parameter.' %p) - ppceDict=self.ppceDict - qTrain=self.qTrain - yTrain=self.yTrain - noiseSdev=self.noiseV - #(0) Assignments - nGQ=ppceDict['nGQtest'] - qInfo=ppceDict['qInfo'] - nMC=ppceDict['nMC'] - nw_=int(nMC/10) - distType=ppceDict['distType'] - #Make a dict for gpr (do NOT change) - gprOpts={'nIter':ppceDict['nIter_gpr'], - 'lr':ppceDict['lr_gpr'], - 'convPlot':ppceDict['convPlot_gpr'] - } - standardizeYTrain_=False - if 'standardizeYTrain_gpr' in ppceDict.keys(): - gprOpts.update({'standardizeYTrain':ppceDict['standardizeYTrain_gpr']}) - standardizeYTrain_=True - - #Make a dict for PCE (do NOT change) - #Always use TP truncation with GQ sampling (hence Projection method) - pceDict={'p':p, - 'truncMethod':'TP', - 'sampleType':'GQ', - 'pceSolveMethod':'Projection', - 'distType':distType - } - - #(1) Generate test points that are Gauss quadratures chosen based on - # the distribution of q (gPCE rule) - qTestList=[] - for i in range(p): - sampsGQ=sampling.trainSample(sampleType='GQ',GQdistType=distType[i], - qInfo=qInfo[i],nSamp=nGQ[i]) - qTestList.append(sampsGQ.q) - qTestGrid=reshaper.vecs2grid(qTestList) - - #(2) Construct GPR surrogate based on training data - gpr_=gpr_torch.gpr(qTrain,yTrain[:,None],noiseSdev,qTestGrid,gprOpts) - post_f=gpr_.post_f - post_obs=gpr_.post_y - shift_=0.0 - scale_=1.0 - if standardizeYTrain_: - shift_=gpr_.shift[0] #0: single-response - scale_=gpr_.scale[0] - - #optional: plot constructed response surface - #gpr_torch.gprPlot().torch2d_3dSurf(qTrain,yTrain,qTestList,post_obs,shift=shift_,scale=scale_) - - #(3) Use samples of GPR tested at GQ nodes to construct a PCE - fMean_list=[] - fVar_list =[] - for j in range(nMC): - # Draw a sample for f(q) from GPR surrogate - f_=post_obs.sample().numpy()*scale_+shift_ - # Construct PCE for the drawn sample - pce_=pce(fVal=f_,nQList=nGQ,xi=[],pceDict=pceDict,verbose=False) - fMean_list.append(pce_.fMean) - fVar_list.append(pce_.fVar) - if ((j+1)%nw_==0): - print("...... ppce repetition for finding samples of the PCE coefficients, iter = %d/%d" - %(j+1,nMC)) - - #(4) Outputs - fMean_list=np.asarray(fMean_list) - fVar_list=np.asarray(fVar_list) - #Optional outputs: only used for gprPlot - optOut={'post_f':post_f,'post_obs':post_obs,'qTest':qTestList} - self.optOut=optOut - self.fMean_samps=fMean_list - self.fVar_samps=fVar_list + """ + Constructing a probabilistic PCE over a p-D parameter space, p>1 + """ + p = self.p + print('... Probabilistic PCE for %d-D input parameter.' % p) + ppceDict = self.ppceDict + qTrain = self.qTrain + yTrain = self.yTrain + noiseSdev = self.noiseV + # (0) Assignments + nGQ = ppceDict['nGQtest'] + qInfo = ppceDict['qInfo'] + nMC = ppceDict['nMC'] + nw_ = int(nMC/10) + distType = ppceDict['distType'] + # Make a dict for gpr (do NOT change) + gprOpts = {'nIter': ppceDict['nIter_gpr'], + 'lr': ppceDict['lr_gpr'], + 'convPlot': ppceDict['convPlot_gpr'] + } + standardizeYTrain_ = False + if 'standardizeYTrain_gpr' in ppceDict.keys(): + gprOpts.update( + {'standardizeYTrain': ppceDict['standardizeYTrain_gpr']}) + standardizeYTrain_ = True + + # Make a dict for PCE (do NOT change) + # Always use TP truncation with GQ sampling (hence Projection method) + pceDict = {'p': p, + 'truncMethod': 'TP', + 'sampleType': 'GQ', + 'pceSolveMethod': 'Projection', + 'distType': distType + } + + # (1) Generate test points that are Gauss quadratures chosen based on + # the distribution of q (gPCE rule) + qTestList = [] + for i in range(p): + sampsGQ = sampling.trainSample(sampleType='GQ', GQdistType=distType[i], + qInfo=qInfo[i], nSamp=nGQ[i]) + qTestList.append(sampsGQ.q) + qTestGrid = reshaper.vecs2grid(qTestList) + + # (2) Construct GPR surrogate based on training data + gpr_ = gpr_torch.gpr( + qTrain, yTrain[:, None], noiseSdev, qTestGrid, gprOpts) + post_f = gpr_.post_f + post_obs = gpr_.post_y + shift_ = 0.0 + scale_ = 1.0 + if standardizeYTrain_: + shift_ = gpr_.shift[0] # 0: single-response + scale_ = gpr_.scale[0] + + # optional: plot constructed response surface + # gpr_torch.gprPlot().torch2d_3dSurf(qTrain,yTrain,qTestList,post_obs,shift=shift_,scale=scale_) + + # (3) Use samples of GPR tested at GQ nodes to construct a PCE + fMean_list = [] + fVar_list = [] + for j in range(nMC): + # Draw a sample for f(q) from GPR surrogate + f_ = post_obs.sample().numpy()*scale_+shift_ + # Construct PCE for the drawn sample + pce_ = pce(fVal=f_, nQList=nGQ, xi=[], + pceDict=pceDict, verbose=False) + fMean_list.append(pce_.fMean) + fVar_list.append(pce_.fVar) + if ((j+1) % nw_ == 0): + print("...... ppce repetition for finding samples of the PCE coefficients, iter = %d/%d" + % (j+1, nMC)) + + # (4) Outputs + fMean_list = np.asarray(fMean_list) + fVar_list = np.asarray(fVar_list) + # Optional outputs: only used for gprPlot + optOut = {'post_f': post_f, 'post_obs': post_obs, 'qTest': qTestList} + self.optOut = optOut + self.fMean_samps = fMean_list + self.fVar_samps = fVar_list # diff --git a/UQit/psobol.py b/UQit/psobol.py index 051f8e5..aa39675 100644 --- a/UQit/psobol.py +++ b/UQit/psobol.py @@ -1,16 +1,22 @@ ############################################################### # Probabilistic Sobol Sensitivity Indices ############################################################### -#-------------------------------------------------------------- +# -------------------------------------------------------------- # Saleh Rezaeiravesh, salehr@kth.se -#-------------------------------------------------------------- -# ToDo: For p>=3, choosing large number of test points to evaluate the -# integrals in the Sobol indices can be a bit problematic since the -# samples generated from GPR (GpyTorch) may lead to memory/CPU issues. -# A solution will be evaluating the integrals for Sobol indices by a -# quadrature rule. Since we have a complete freedom, using GLL rule would be good. -#-------------------------------------------------------------- +# -------------------------------------------------------------- +# ToDo: For p>=3, choosing large number of test points to evaluate the +# integrals in the Sobol indices can be a bit problematic since the +# samples generated from GPR (GpyTorch) may lead to memory/CPU issues. +# A solution will be evaluating the integrals for Sobol indices by a +# quadrature rule. Since we have a complete freedom, using GLL rule would be good. +# -------------------------------------------------------------- # +from math import pi +import UQit.sampling as sampling +import UQit.analyticTestFuncs as analyticTestFuncs +import matplotlib.pyplot as plt +import copy +import math as mt import os import sys import numpy as np @@ -21,6 +27,8 @@ ##import UQit.sampling as sampling # # + + class psobol: """ Probabilistic Sobol (psobol) Sensitivity Indices over a p-D parameter space, @@ -55,7 +63,7 @@ class psobol: * 'convPlot_gpr': bool If true, values of the hyper-parameters are plotted vs. iteration during the optimization. - + Attributes: `Si_samps`: list of length p =[S1_samps,S2_samps,...,Sp_samps] where Si_samps is a 1d numpy array of size nMC containing @@ -73,281 +81,299 @@ class psobol: * `post_f`: Posterior density of f(q) * `post_obs`: Posterior density of y """ - def __init__(self,qTrain,yTrain,noiseV,psobolDict): - self.qTrain=qTrain - self.yTrain=yTrain - self.noiseV=noiseV - self.psobolDict=psobolDict + + def __init__(self, qTrain, yTrain, noiseV, psobolDict): + self.qTrain = qTrain + self.yTrain = yTrain + self.noiseV = noiseV + self.psobolDict = psobolDict self.info() self.cnstrct() def info(self): - if self.qTrain.ndim==1: - raise ValueError("qTrain must be a 2D numpy array.") - elif self.qTrain.ndim==2: - self.p=self.qTrain.shape[-1] + if self.qTrain.ndim == 1: + raise ValueError("qTrain must be a 2D numpy array.") + elif self.qTrain.ndim == 2: + self.p = self.qTrain.shape[-1] - if self.qTrain.shape[0]==self.yTrain.shape: - raise ValueError("Size of qTrain and yTrain should be the same.") + if self.qTrain.shape[0] == self.yTrain.shape: + raise ValueError("Size of qTrain and yTrain should be the same.") - obligKeys=['qTest','pdf','nMC','nIter_gpr','lr_gpr','convPlot_gpr'] + obligKeys = ['qTest', 'pdf', 'nMC', + 'nIter_gpr', 'lr_gpr', 'convPlot_gpr'] for key_ in obligKeys: if key_ not in self.psobolDict.keys(): - raise KeyError("%s is required in psobolDict." %key_) + raise KeyError("%s is required in psobolDict." % key_) - nTest=[] + nTest = [] for i in range(self.p): - nTest.append(self.psobolDict['qTest'][i].shape[0]) - self.nTest=nTest + nTest.append(self.psobolDict['qTest'][i].shape[0]) + self.nTest = nTest def cnstrct(self): - self.psobol_cnstrct() + self.psobol_cnstrct() def psobol_cnstrct(self): - """ - Constructing probabilistic Sobol indices over a p-D parameter space, p>1 - """ - p=self.p - print('... Probabilistic Sobol indices for %d-D input parameter.' %p) - psobolDict=self.psobolDict - qTrain=self.qTrain - yTrain=self.yTrain - noiseSdev=self.noiseV - #(0) Assignments - qTest=psobolDict['qTest'] - pdf=psobolDict['pdf'] - nMC=psobolDict['nMC'] - nw_=int(nMC/10) - #Make a dict for gpr (do NOT change) - gprOpts={'nIter':psobolDict['nIter_gpr'], - 'lr':psobolDict['lr_gpr'], - 'convPlot':psobolDict['convPlot_gpr'] - } - standardizeYTrain_=False - if 'standardizeYTrain_gpr' in psobolDict.keys(): - gprOpts.update({'standardizeYTrain':psobolDict['standardizeYTrain_gpr']}) - standardizeYTrain_=True - - #(1) Generate a tensor product grid from qTest. At the grid samples, the gpr is sampled. - qTestGrid=reshaper.vecs2grid(qTest) - - #(2) Construct GPR surrogate based on training data - gpr_=gpr_torch.gpr(qTrain,yTrain[:,None],noiseSdev,qTestGrid,gprOpts) - post_f=gpr_.post_f - post_obs=gpr_.post_y - shift_=0.0 - scale_=1.0 - if standardizeYTrain_: - shift_=gpr_.shift[0] #0: single-response - scale_=gpr_.scale[0] - - #optional: plot constructed response surface only for p==2 - #gpr_torch.gprPlot().torch2d_3dSurf(qTrain,yTrain,qTest,post_obs,shift=shift_,scale=scale_) - - #(3) Compute Sobol indices for samples of GPR generated at qTestGrid - Si_list_=[] - Sij_list_=[] - STi_list_=[] - for j in range(nMC): - # Draw a sample for f(q) from GPR surrogate - f_=post_obs.sample().numpy()*scale_+shift_ - f_=np.reshape(f_,self.nTest,'F') - # Compute the Sobol indices - sobol_=sobol(qTest,f_,pdf) - Si_list_.append(sobol_.Si) - Sij_list_.append(sobol_.Sij) - STi_list_.append(sobol_.STi) - if ((j+1)%nw_==0): - print("...... psobol repetition for finding samples of the Sobol indices, iter = %d/%d" - %(j+1,nMC)) - #reshape lists and arrays - S_=np.zeros(nMC) - ST_=np.zeros(nMC) - Si_list=[] - Sij_list=[] - STi_list=[] - for i in range(p): - for j in range(nMC): - S_[j]=Si_list_[j][i] - ST_[j]=STi_list_[j][i] - Si_list.append(S_.copy()) - STi_list.append(ST_.copy()) - - for i in range(len(Sij_list_[0])): - for j in range(nMC): - S_[j]=Sij_list_[j][i] - Sij_list.append(S_.copy()) - - self.Si_samps=Si_list - self.Sij_samps=Sij_list - self.STi_samps=STi_list - self.SijName=sobol_.SijName - - #(4) Outputs - #Optional outputs: can only be used for gprPlot - optOut={'post_f':post_f,'post_obs':post_obs} - self.optOut=optOut + """ + Constructing probabilistic Sobol indices over a p-D parameter space, p>1 + """ + p = self.p + print('... Probabilistic Sobol indices for %d-D input parameter.' % p) + psobolDict = self.psobolDict + qTrain = self.qTrain + yTrain = self.yTrain + noiseSdev = self.noiseV + # (0) Assignments + qTest = psobolDict['qTest'] + pdf = psobolDict['pdf'] + nMC = psobolDict['nMC'] + nw_ = int(nMC/10) + # Make a dict for gpr (do NOT change) + gprOpts = {'nIter': psobolDict['nIter_gpr'], + 'lr': psobolDict['lr_gpr'], + 'convPlot': psobolDict['convPlot_gpr'] + } + standardizeYTrain_ = False + if 'standardizeYTrain_gpr' in psobolDict.keys(): + gprOpts.update( + {'standardizeYTrain': psobolDict['standardizeYTrain_gpr']}) + standardizeYTrain_ = True + + # (1) Generate a tensor product grid from qTest. At the grid samples, the gpr is sampled. + qTestGrid = reshaper.vecs2grid(qTest) + + # (2) Construct GPR surrogate based on training data + gpr_ = gpr_torch.gpr( + qTrain, yTrain[:, None], noiseSdev, qTestGrid, gprOpts) + post_f = gpr_.post_f + post_obs = gpr_.post_y + shift_ = 0.0 + scale_ = 1.0 + if standardizeYTrain_: + shift_ = gpr_.shift[0] # 0: single-response + scale_ = gpr_.scale[0] + + # optional: plot constructed response surface only for p==2 + # gpr_torch.gprPlot().torch2d_3dSurf(qTrain,yTrain,qTest,post_obs,shift=shift_,scale=scale_) + + # (3) Compute Sobol indices for samples of GPR generated at qTestGrid + Si_list_ = [] + Sij_list_ = [] + STi_list_ = [] + for j in range(nMC): + # Draw a sample for f(q) from GPR surrogate + f_ = post_obs.sample().numpy()*scale_+shift_ + f_ = np.reshape(f_, self.nTest, 'F') + # Compute the Sobol indices + sobol_ = sobol(qTest, f_, pdf) + Si_list_.append(sobol_.Si) + Sij_list_.append(sobol_.Sij) + STi_list_.append(sobol_.STi) + if ((j+1) % nw_ == 0): + print("...... psobol repetition for finding samples of the Sobol indices, iter = %d/%d" + % (j+1, nMC)) + # reshape lists and arrays + S_ = np.zeros(nMC) + ST_ = np.zeros(nMC) + Si_list = [] + Sij_list = [] + STi_list = [] + for i in range(p): + for j in range(nMC): + S_[j] = Si_list_[j][i] + ST_[j] = STi_list_[j][i] + Si_list.append(S_.copy()) + STi_list.append(ST_.copy()) + + for i in range(len(Sij_list_[0])): + for j in range(nMC): + S_[j] = Sij_list_[j][i] + Sij_list.append(S_.copy()) + + self.Si_samps = Si_list + self.Sij_samps = Sij_list + self.STi_samps = STi_list + self.SijName = sobol_.SijName + + # (4) Outputs + # Optional outputs: can only be used for gprPlot + optOut = {'post_f': post_f, 'post_obs': post_obs} + self.optOut = optOut + # # ######### -##MAIN +# MAIN ######### -import math as mt -import copy -import matplotlib.pyplot as plt -import UQit.analyticTestFuncs as analyticTestFuncs -import UQit.sampling as sampling + + def psobol_2d_test(): """ Test for psobol for 2 parameters """ ## - def plot_trainData(n,fSamples,noiseSdev,yTrain): + def plot_trainData(n, fSamples, noiseSdev, yTrain): """ Plot the noisy training data which are used in GPR. """ - plt.figure(figsize=(10,5)) - x_=np.zeros(n) + plt.figure(figsize=(10, 5)) + x_ = np.zeros(n) for i in range(n): - x_[i]=i+1 - for i in range(500): #only for plottig possible realizations - noise_=noiseSdev*np.random.randn(n) - plt.plot(x_,fSamples+noise_,'.',color='steelblue',alpha=0.4,markersize=1) - plt.errorbar(x_,fSamples,yerr=1.96*abs(noiseSdev),ls='none',capsize=5,ecolor='k', - elinewidth=4,label=r'$95\%$ CI in Obs.') - plt.plot(x_,fSamples,'o' ,markersize=6,markerfacecolor='lime', - markeredgecolor='salmon',label='Mean Observation') - plt.plot(x_,yTrain ,'xr' ,markersize=6,label='Sample Observation') - plt.legend(loc='best',fontsize=15) - plt.ylabel('QoI',fontsize=17) - plt.xlabel('Simulation Index',fontsize=17) + x_[i] = i+1 + for i in range(500): # only for plottig possible realizations + noise_ = noiseSdev*np.random.randn(n) + plt.plot(x_, fSamples+noise_, '.', + color='steelblue', alpha=0.4, markersize=1) + plt.errorbar(x_, fSamples, yerr=1.96*abs(noiseSdev), ls='none', capsize=5, ecolor='k', + elinewidth=4, label=r'$95\%$ CI in Obs.') + plt.plot(x_, fSamples, 'o', markersize=6, markerfacecolor='lime', + markeredgecolor='salmon', label='Mean Observation') + plt.plot(x_, yTrain, 'xr', markersize=6, label='Sample Observation') + plt.legend(loc='best', fontsize=15) + plt.ylabel('QoI', fontsize=17) + plt.xlabel('Simulation Index', fontsize=17) plt.xticks(fontsize=15) - plt.yticks(fontsize=15) + plt.yticks(fontsize=15) plt.title('Training data with associated confidence') plt.show() ## - def trainDataGen(p,sampleType,n,qBound,fExName,noiseType): + + def trainDataGen(p, sampleType, n, qBound, fExName, noiseType): """ Generate Training Data """ - # (a) xTrain - if sampleType=='grid': - nSamp=n[0]*n[1] - gridList=[] - for i in range(p): - grid_=np.linspace(qBound[i][0],qBound[i][1],n[i]) - gridList.append(grid_) - xTrain=reshaper.vecs2grid(gridList) - elif sampleType=='random': - nSamp=n #number of random samples - xTrain=sampling.LHS_sampling(n,qBound) - # (b) Observation noise - noiseSdev=noiseGen(nSamp,noiseType,xTrain,fExName) + # (a) xTrain + if sampleType == 'grid': + nSamp = n[0]*n[1] + gridList = [] + for i in range(p): + grid_ = np.linspace(qBound[i][0], qBound[i][1], n[i]) + gridList.append(grid_) + xTrain = reshaper.vecs2grid(gridList) + elif sampleType == 'random': + nSamp = n # number of random samples + xTrain = sampling.LHS_sampling(n, qBound) + # (b) Observation noise + noiseSdev = noiseGen(nSamp, noiseType, xTrain, fExName) # (c) Training response - yTrain=analyticTestFuncs.fEx2D(xTrain[:,0],xTrain[:,1],fExName,'comp').val - yTrain_noiseFree=yTrain - yTrain=yTrain_noiseFree+noiseSdev*np.random.randn(nSamp) - return xTrain,yTrain,noiseSdev,yTrain_noiseFree - ## - def noiseGen(n,noiseType,xTrain,fExName): - """ - Generate a 1D numpy array of standard deviations of the observation noise - """ - if noiseType=='homo': - sd=0.5 # noise standard deviation (Note: non-zero, to avoid instabilities) - sdV=sd*np.ones(n) - elif noiseType=='hetero': - sdV=0.5*(analyticTestFuncs.fEx2D(xTrain[:,0],xTrain[:,1],fExName,'comp').val+0.001) - return sdV + yTrain = analyticTestFuncs.fEx2D( + xTrain[:, 0], xTrain[:, 1], fExName, 'comp').val + yTrain_noiseFree = yTrain + yTrain = yTrain_noiseFree+noiseSdev*np.random.randn(nSamp) + return xTrain, yTrain, noiseSdev, yTrain_noiseFree + ## + + def noiseGen(n, noiseType, xTrain, fExName): + """ + Generate a 1D numpy array of standard deviations of the observation noise + """ + if noiseType == 'homo': + # noise standard deviation (Note: non-zero, to avoid instabilities) + sd = 0.5 + sdV = sd*np.ones(n) + elif noiseType == 'hetero': + sdV = 0.5 * \ + (analyticTestFuncs.fEx2D( + xTrain[:, 0], xTrain[:, 1], fExName, 'comp').val+0.001) + return sdV # #----- SETTINGS - #definition of the parameters - qInfo=[[-3,2],[0.5,0.7]] #information about the parameters - distType=['Unif','Norm'] #type of distribution of the parameters - - #options for generating training samples - fExName='type1' #Type of simulator in analyticTestFuncs.fEx2D - #'type1', 'type2', 'type3', 'Rosenbrock' - sampleType='random' #'random' or 'grid': type of training samples - if sampleType=='grid': - n=[9,9] #number of training samples in each input dimension - elif sampleType=='random': - n=100 #total number of training samples drawn randomly - noiseType='hetero' #noise type: 'homo'=homoscedastic, 'hetero'=heterscedastic - - #options for Sobol indices - nMC=1000 #number of MC samples to compute psobol indices - nTest=[41,40] #number of test points in each parameter dimension to compute integrals in Sobol indices - - #options for GPR - nIter_gpr_=500 #number of iterations in optimization of GPR hyperparameters - lr_gpr_ =0.05 #learning rate in the optimization of the hyperparameters - convPlot_gpr_=True #plot convergence of optimization of GPR hyperparameters - #------------------------------------------------ - qBound=copy.deepcopy(qInfo) #default - #(1) Generate training data - p=len(qInfo) #dimension of the input + # definition of the parameters + qInfo = [[-3, 2], [0.5, 0.7]] # information about the parameters + distType = ['Unif', 'Norm'] # type of distribution of the parameters + + # options for generating training samples + fExName = 'type1' # Type of simulator in analyticTestFuncs.fEx2D + #'type1', 'type2', 'type3', 'Rosenbrock' + sampleType = 'random' # 'random' or 'grid': type of training samples + if sampleType == 'grid': + n = [9, 9] # number of training samples in each input dimension + elif sampleType == 'random': + n = 100 # total number of training samples drawn randomly + noiseType = 'hetero' # noise type: 'homo'=homoscedastic, 'hetero'=heterscedastic + + # options for Sobol indices + nMC = 1000 # number of MC samples to compute psobol indices + # number of test points in each parameter dimension to compute integrals in Sobol indices + nTest = [41, 40] + + # options for GPR + nIter_gpr_ = 500 # number of iterations in optimization of GPR hyperparameters + lr_gpr_ = 0.05 # learning rate in the optimization of the hyperparameters + convPlot_gpr_ = True # plot convergence of optimization of GPR hyperparameters + # ------------------------------------------------ + qBound = copy.deepcopy(qInfo) # default + # (1) Generate training data + p = len(qInfo) # dimension of the input for i in range(p): - if distType[i]=='Norm': - qBound[i][0]=qInfo[i][0]-5.*qInfo[i][1] - qBound[i][1]=qInfo[i][0]+5.*qInfo[i][1] - xTrain,yTrain,noiseSdev,yTrain_noiseFree=trainDataGen(p,sampleType,n,qBound,fExName,noiseType) - nSamp=yTrain.shape[0] - plot_trainData(nSamp,yTrain_noiseFree,noiseSdev,yTrain) - - #(2) Create the test samples and PDF - qTest=[] - pdf=[] + if distType[i] == 'Norm': + qBound[i][0] = qInfo[i][0]-5.*qInfo[i][1] + qBound[i][1] = qInfo[i][0]+5.*qInfo[i][1] + xTrain, yTrain, noiseSdev, yTrain_noiseFree = trainDataGen( + p, sampleType, n, qBound, fExName, noiseType) + nSamp = yTrain.shape[0] + plot_trainData(nSamp, yTrain_noiseFree, noiseSdev, yTrain) + + # (2) Create the test samples and PDF + qTest = [] + pdf = [] for i in range(p): - q_=np.linspace(qBound[i][0],qBound[i][1],nTest[i]) - if distType[i]=='Norm': - pdf.append(np.exp(-(q_-qInfo[i][0])**2/(2*qInfo[i][1]**2))/(qInfo[i][1]*mt.sqrt(2*mt.pi))) - elif distType[i]=='Unif': - pdf.append(np.ones(nTest[i])/(qBound[i][1]-qBound[i][0])) + q_ = np.linspace(qBound[i][0], qBound[i][1], nTest[i]) + if distType[i] == 'Norm': + pdf.append( + np.exp(-(q_-qInfo[i][0])**2/(2*qInfo[i][1]**2))/(qInfo[i][1]*mt.sqrt(2*mt.pi))) + elif distType[i] == 'Unif': + pdf.append(np.ones(nTest[i])/(qBound[i][1]-qBound[i][0])) else: - raise ValueError("distType of the %d-th parameter can be 'Unif' or 'Norm'"%(i+1)) + raise ValueError( + "distType of the %d-th parameter can be 'Unif' or 'Norm'" % (i+1)) qTest.append(q_) - #plot PDFs + # plot PDFs for i in range(p): - plt.plot(qTest[i],pdf[i],label='pdf of q'+str(i+1)) + plt.plot(qTest[i], pdf[i], label='pdf of q'+str(i+1)) plt.legend(loc='best') - plt.show() + plt.show() - #(3) Assemble the psobolOpts dict - psobolDict={'nMC':nMC,'qTest':qTest,'pdf':pdf, - 'nIter_gpr':nIter_gpr_,'lr_gpr':lr_gpr_,'convPlot_gpr':convPlot_gpr_} + # (3) Assemble the psobolOpts dict + psobolDict = {'nMC': nMC, 'qTest': qTest, 'pdf': pdf, + 'nIter_gpr': nIter_gpr_, 'lr_gpr': lr_gpr_, 'convPlot_gpr': convPlot_gpr_} # (4) Construct probabilistic Sobol indices - psobol_=psobol(xTrain,yTrain,noiseSdev,psobolDict) - Si_samps=psobol_.Si_samps - Sij_samps=psobol_.Sij_samps - STi_samps=psobol_.STi_samps - SijName=psobol_.SijName + psobol_ = psobol(xTrain, yTrain, noiseSdev, psobolDict) + Si_samps = psobol_.Si_samps + Sij_samps = psobol_.Sij_samps + STi_samps = psobol_.STi_samps + SijName = psobol_.SijName # Results # (a) Compute mean and sdev of the estimated indices for i in range(p): - print('Main Sobol index S%d: mean=%g, sdev=%g',str(i+1),np.mean(Si_samps[i]),np.std(Si_samps[i])) - print('Total Sobol index ST%d: mean=%g, sdev=%g',str(i+1),np.mean(STi_samps[i]),np.std(STi_samps[i])) - print('Interacting Sobol indices S12: mean=%g, sdev=%g',np.mean(Sij_samps[0]),np.std(Sij_samps[0])) + print('Main Sobol index S%d: mean=%g, sdev=%g', str( + i+1), np.mean(Si_samps[i]), np.std(Si_samps[i])) + print('Total Sobol index ST%d: mean=%g, sdev=%g', str( + i+1), np.mean(STi_samps[i]), np.std(STi_samps[i])) + print('Interacting Sobol indices S12: mean=%g, sdev=%g', + np.mean(Sij_samps[0]), np.std(Sij_samps[0])) # (b) plot histogram and fitted pdf to different Sobol indices - statsUQit.pdfFit_uniVar(Si_samps[0],True,[]) - statsUQit.pdfFit_uniVar(Si_samps[1],True,[]) - statsUQit.pdfFit_uniVar(Sij_samps[0],True,[]) + statsUQit.pdfFit_uniVar(Si_samps[0], True, []) + statsUQit.pdfFit_uniVar(Si_samps[1], True, []) + statsUQit.pdfFit_uniVar(Sij_samps[0], True, []) # (c) plot samples of the Sobol indices - plt.figure(figsize=(10,5)) + plt.figure(figsize=(10, 5)) for i in range(p): - plt.plot(Si_samps[i],'-',label='S'+str(i+1)) - plt.plot(STi_samps[i],'--',label='ST'+str(i+1)) - plt.plot(Si_samps[i],':',label='S'+str(12)) - plt.legend(loc='best',fontsize=14) - plt.xlabel('Sample',fontsize=14) - plt.ylabel('Sobol indices',fontsize=14) + plt.plot(Si_samps[i], '-', label='S'+str(i+1)) + plt.plot(STi_samps[i], '--', label='ST'+str(i+1)) + plt.plot(Si_samps[i], ':', label='S'+str(12)) + plt.legend(loc='best', fontsize=14) + plt.xlabel('Sample', fontsize=14) + plt.ylabel('Sobol indices', fontsize=14) plt.show() + # -from math import pi + + def psobol_Ishigami_test(): """ Test for psobol for 3 uncertain parameters. @@ -355,75 +381,82 @@ def psobol_Ishigami_test(): The resulting psobol indices can be compared to the standard Sobol indices in order to verify the implementation of psobol. """ - #-------------------------- + # -------------------------- #------- SETTINGS - #definition of the parameters - qBound=[[-pi,pi], #admissible range of parameters - [-pi,pi], - [-pi,pi]] - #options for Training data + # definition of the parameters + qBound = [[-pi, pi], # admissible range of parameters + [-pi, pi], + [-pi, pi]] + # options for Training data # n=[100, 70, 80] #number of samples for q1, q2, q3 - n=500 #number of LHS random samples - a=7 #parameters in Ishigami function - b=0.1 - noise_sdev=0.2 #standard-deviation of observation noise - #options for Sobol indices - nMC=500 #number of MC samples to compute psobol indices - nTest=[20,21,22] #number of test points in each parameter dimension to compute integrals in Sobol indices - #options for GPR - nIter_gpr_=100 #number of iterations in optimization of GPR hyperparameters - lr_gpr_ =0.05 #learning rate in the optimization of the hyperparameters - convPlot_gpr_=True #plot convergence of optimization of GPR hyperparameters - #-------------------------- - p=len(qBound) - #(1) Generate training data - qTrain=sampling.LHS_sampling(n,qBound) #LHS random samples - fEx_=analyticTestFuncs.fEx3D(qTrain[:,0],qTrain[:,1],qTrain[:,2],'Ishigami','comp',{'a':a,'b':b}) - yTrain=fEx_.val+noise_sdev*np.random.randn(n) - noiseSdev=noise_sdev*np.ones(n) - - #(2) Create the test samples and associated PDF - qTest=[] - pdf=[] - #qBound=[[-1,1],[-1,1],[-1,1]] - #for i in range(3): + n = 500 # number of LHS random samples + a = 7 # parameters in Ishigami function + b = 0.1 + noise_sdev = 0.2 # standard-deviation of observation noise + # options for Sobol indices + nMC = 500 # number of MC samples to compute psobol indices + # number of test points in each parameter dimension to compute integrals in Sobol indices + nTest = [20, 21, 22] + # options for GPR + nIter_gpr_ = 100 # number of iterations in optimization of GPR hyperparameters + lr_gpr_ = 0.05 # learning rate in the optimization of the hyperparameters + convPlot_gpr_ = True # plot convergence of optimization of GPR hyperparameters + # -------------------------- + p = len(qBound) + # (1) Generate training data + qTrain = sampling.LHS_sampling(n, qBound) # LHS random samples + fEx_ = analyticTestFuncs.fEx3D( + qTrain[:, 0], qTrain[:, 1], qTrain[:, 2], 'Ishigami', 'comp', {'a': a, 'b': b}) + yTrain = fEx_.val+noise_sdev*np.random.randn(n) + noiseSdev = noise_sdev*np.ones(n) + + # (2) Create the test samples and associated PDF + qTest = [] + pdf = [] + # qBound=[[-1,1],[-1,1],[-1,1]] + # for i in range(3): # qTrain[:,i]=(qTrain[:,i]+pi)/(2*pi)*2-1 for i in range(p): - qTest.append(np.linspace(qBound[i][0],qBound[i][1],nTest[i])) + qTest.append(np.linspace(qBound[i][0], qBound[i][1], nTest[i])) pdf.append(np.ones(nTest[i])/(qBound[i][1]-qBound[i][0])) - #(3) Assemble the psobolOpts dict - psobolDict={'nMC':nMC,'qTest':qTest,'pdf':pdf, - 'nIter_gpr':nIter_gpr_,'lr_gpr':lr_gpr_,'convPlot_gpr':convPlot_gpr_} + # (3) Assemble the psobolOpts dict + psobolDict = {'nMC': nMC, 'qTest': qTest, 'pdf': pdf, + 'nIter_gpr': nIter_gpr_, 'lr_gpr': lr_gpr_, 'convPlot_gpr': convPlot_gpr_} # (4) Construct probabilistic Sobol indices - psobol_=psobol(qTrain,yTrain,noiseSdev,psobolDict) - Si_samps=psobol_.Si_samps - Sij_samps=psobol_.Sij_samps - STi_samps=psobol_.STi_samps - SijName=psobol_.SijName + psobol_ = psobol(qTrain, yTrain, noiseSdev, psobolDict) + Si_samps = psobol_.Si_samps + Sij_samps = psobol_.Sij_samps + STi_samps = psobol_.STi_samps + SijName = psobol_.SijName # (a) Compute mean and sdev of the estimated indices for i in range(p): - print('Main Sobol index S%d: mean=%g, sdev=%g',str(i+1),np.mean(Si_samps[i]),np.std(Si_samps[i])) - print('Total Sobol index ST%d: mean=%g, sdev=%g',str(i+1),np.mean(STi_samps[i]),np.std(STi_samps[i])) - print('Interacting Sobol indices S12: mean=%g, sdev=%g',np.mean(Sij_samps[0]),np.std(Sij_samps[0])) - print('Interacting Sobol indices S13: mean=%g, sdev=%g',np.mean(Sij_samps[1]),np.std(Sij_samps[1])) - print('Interacting Sobol indices S23: mean=%g, sdev=%g',np.mean(Sij_samps[2]),np.std(Sij_samps[2])) + print('Main Sobol index S%d: mean=%g, sdev=%g', str( + i+1), np.mean(Si_samps[i]), np.std(Si_samps[i])) + print('Total Sobol index ST%d: mean=%g, sdev=%g', str( + i+1), np.mean(STi_samps[i]), np.std(STi_samps[i])) + print('Interacting Sobol indices S12: mean=%g, sdev=%g', + np.mean(Sij_samps[0]), np.std(Sij_samps[0])) + print('Interacting Sobol indices S13: mean=%g, sdev=%g', + np.mean(Sij_samps[1]), np.std(Sij_samps[1])) + print('Interacting Sobol indices S23: mean=%g, sdev=%g', + np.mean(Sij_samps[2]), np.std(Sij_samps[2])) # (b) plot histogram and fitted pdf to different Sobol indices - statsUQit.pdfFit_uniVar(Si_samps[0],True,[]) - statsUQit.pdfFit_uniVar(Si_samps[1],True,[]) - statsUQit.pdfFit_uniVar(Sij_samps[0],True,[]) + statsUQit.pdfFit_uniVar(Si_samps[0], True, []) + statsUQit.pdfFit_uniVar(Si_samps[1], True, []) + statsUQit.pdfFit_uniVar(Sij_samps[0], True, []) # (c) plot samples of the Sobol indices - plt.figure(figsize=(10,5)) + plt.figure(figsize=(10, 5)) for i in range(p): - plt.plot(Si_samps[i],'-',label='S'+str(i+1)) - plt.plot(STi_samps[i],'--',label='ST'+str(i+1)) - plt.plot(Si_samps[i],':',label='S'+str(12)) - plt.legend(loc='best',fontsize=14) - plt.xlabel('sample',fontsize=14) - plt.ylabel('Sobol indices',fontsize=14) + plt.plot(Si_samps[i], '-', label='S'+str(i+1)) + plt.plot(STi_samps[i], '--', label='ST'+str(i+1)) + plt.plot(Si_samps[i], ':', label='S'+str(12)) + plt.legend(loc='best', fontsize=14) + plt.xlabel('sample', fontsize=14) + plt.ylabel('Sobol indices', fontsize=14) plt.show() -# +# diff --git a/UQit/reshaper.py b/UQit/reshaper.py index 7f31c3a..073258a 100644 --- a/UQit/reshaper.py +++ b/UQit/reshaper.py @@ -2,23 +2,27 @@ # Tools for converting and reshaping arrays and lists ####################################################### # Saleh Rezaeiravesh, salehr@kth.se -#------------------------------------------------------ +# ------------------------------------------------------ # import numpy as np # # + + def lengthVector(x): """ Returns the length of vector x which can be a list or a numpy array """ - if (isinstance(x,np.ndarray)): - nx=x.size - elif (isinstance(x,(list))): - nx=len(x) + if (isinstance(x, np.ndarray)): + nx = x.size + elif (isinstance(x, (list))): + nx = len(x) else: - raise ValueError("Unknown object x") + raise ValueError("Unknown object x") return int(nx) # + + def vecs2grid(x): """ Makes a p-D (p>1) tensor-product grid from a list of length p containg 1D numpy arrays of points in each dimension. @@ -26,20 +30,22 @@ def vecs2grid(x): Args: `x`: A list of length p>1 x=[x1,x2,...,xp] where xi is a 1D numpy array of size ni - + Returns: 'z': A numpy array of shape (n1*n2*...*np,p) """ - p=len(x) - if p<=1: - raise ValueError("Import a list of length p>1.") - z_=np.meshgrid(*x,copy=True,sparse=False,indexing='ij') - n=z_[-1].size - z=np.zeros((n,p)) + p = len(x) + if p <= 1: + raise ValueError("Import a list of length p>1.") + z_ = np.meshgrid(*x, copy=True, sparse=False, indexing='ij') + n = z_[-1].size + z = np.zeros((n, p)) for i in range(p): - z[:,i]=z_[i].reshape((n,1),order='F')[:,0] + z[:, i] = z_[i].reshape((n, 1), order='F')[:, 0] return z # + + def vecsGlue(*x): """ Makes a set by gluing p>1 1D numpy arrays x0,x1,...,xp of the same size (=n) @@ -51,16 +57,17 @@ def vecsGlue(*x): `z`: numpy array of shape (n,p) z[:,i]=xi where i=1,2,...,p """ - p=len(x) - if p<=1: - raise ValueError("More than one numpy array must be imported.") - n=lengthVector(x[0]) - for i in range(1,p): - n_=len(x[i]) - if n_!=n: - raise ValueError("Imported numpy arrays should have the same size.") - z=np.zeros((n,p)) + p = len(x) + if p <= 1: + raise ValueError("More than one numpy array must be imported.") + n = lengthVector(x[0]) + for i in range(1, p): + n_ = len(x[i]) + if n_ != n: + raise ValueError( + "Imported numpy arrays should have the same size.") + z = np.zeros((n, p)) for j in range(p): - for i in range(n): - z[i,j]=x[j][i] + for i in range(n): + z[i, j] = x[j][i] return z diff --git a/UQit/sampling.py b/UQit/sampling.py index 8c73092..c4847b2 100644 --- a/UQit/sampling.py +++ b/UQit/sampling.py @@ -1,9 +1,9 @@ ############################################# # Sampling from parameter space ############################################# -#-------------------------------------------- +# -------------------------------------------- # Saleh Rezaeiravesh, salehr@kth.se -#-------------------------------------------- +# -------------------------------------------- # import sys import os @@ -12,6 +12,8 @@ import UQit.nodes as nodes import UQit.pce as pce # + + class trainSample: R""" Generating training samples from a 1D paramter space using different methods. @@ -49,105 +51,109 @@ class trainSample: Admissible range of `q` `w`: 1D numpy array of size `nSamp` Weights in Gauss-Quadrature rule only if sampleType=='GQ' - + Examples: ts1=trainSample(sampleType='GQ',GQdistType='Unif',qInfo=[2,3],nSamp=10) ts2=trainSample(sampleType='NormRand',qInfo=[2,3],nSamp=10) ts3=trainSample(sampleType='GLL',qInfo=[2,3],nSamp=10) """ - def __init__(self,sampleType='',GQdistType='',qInfo=[],nSamp=0): + + def __init__(self, sampleType='', GQdistType='', qInfo=[], nSamp=0): self.info() - self.sampleType=sampleType - self.GQdistType=GQdistType + self.sampleType = sampleType + self.GQdistType = GQdistType self.check() - self.qInfo=qInfo - self.nSamp=nSamp - self.w=[[]]*self.nSamp + self.qInfo = qInfo + self.nSamp = nSamp + self.w = [[]]*self.nSamp self.genSamples() def info(self): - sampleTypeList=['GQ','GLL','unifSpaced','unifRand','normRand','Clenshaw',\ - 'Clenshaw-Curtis'] - GQdistList=['Unif','Norm'] #list of available distributions for gpce - self.sampleTypeList=sampleTypeList - self.GQdistList=GQdistList + sampleTypeList = ['GQ', 'GLL', 'unifSpaced', 'unifRand', 'normRand', 'Clenshaw', + 'Clenshaw-Curtis'] + # list of available distributions for gpce + GQdistList = ['Unif', 'Norm'] + self.sampleTypeList = sampleTypeList + self.GQdistList = GQdistList def check(self): if self.sampleType not in self.sampleTypeList: - raise KeyError('#ERROR @ parSample: Invalid sampleType! Choose from'\ - ,self.sampleTypeList) - if self.sampleType=='GQ' and self.GQdistType not in self.GQdistList: - raise KeyError('#ERROR @ parSample: Invalid GQdistType! Choose from'\ - ,self.GQdistList) + raise KeyError( + '#ERROR @ parSample: Invalid sampleType! Choose from', self.sampleTypeList) + if self.sampleType == 'GQ' and self.GQdistType not in self.GQdistList: + raise KeyError( + '#ERROR @ parSample: Invalid GQdistType! Choose from', self.GQdistList) - def genSamples(self): - n=self.nSamp - if self.sampleType=='GQ' and self.GQdistType in self.GQdistList: - self.gqPtsWts() - elif self.sampleType=='normRand': - self.xi=np.random.normal(size=n) - self.xiBound=[min(self.xi),max(self.xi)] - self.mean=self.qInfo[0] - self.sdev=self.qInfo[1] - self.xi2q_map() - else: - if self.sampleType=='GLL': - self.xiBound=[-1,1] - xi_,w_=nodes.gllPts(n) - self.xi=xi_ - self.w=w_ - if self.sampleType=='unifSpaced': - xiBound_=[0,1] - self.xiBound=xiBound_ - self.xi=np.linspace(xiBound_[0],xiBound_[1],n) - elif self.sampleType=='unifRand': - self.xiBound=[0,1] - self.xi=np.random.rand(n) - elif self.sampleType=='Clenshaw': - self.xiBound=[-1,1] - self.xi=nodes.Clenshaw_pts(n) - elif self.sampleType=='Clenshaw-Curtis': - self.xiBound=[0,1] - l_=1+int(mt.log(n-1)/mt.log(2)) - self.xi=nodes.ClenshawCurtis_pts(l_) - self.qBound=self.qInfo - self.xi2q_map() + def genSamples(self): + n = self.nSamp + if self.sampleType == 'GQ' and self.GQdistType in self.GQdistList: + self.gqPtsWts() + elif self.sampleType == 'normRand': + self.xi = np.random.normal(size=n) + self.xiBound = [min(self.xi), max(self.xi)] + self.mean = self.qInfo[0] + self.sdev = self.qInfo[1] + self.xi2q_map() + else: + if self.sampleType == 'GLL': + self.xiBound = [-1, 1] + xi_, w_ = nodes.gllPts(n) + self.xi = xi_ + self.w = w_ + if self.sampleType == 'unifSpaced': + xiBound_ = [0, 1] + self.xiBound = xiBound_ + self.xi = np.linspace(xiBound_[0], xiBound_[1], n) + elif self.sampleType == 'unifRand': + self.xiBound = [0, 1] + self.xi = np.random.rand(n) + elif self.sampleType == 'Clenshaw': + self.xiBound = [-1, 1] + self.xi = nodes.Clenshaw_pts(n) + elif self.sampleType == 'Clenshaw-Curtis': + self.xiBound = [0, 1] + l_ = 1+int(mt.log(n-1)/mt.log(2)) + self.xi = nodes.ClenshawCurtis_pts(l_) + self.qBound = self.qInfo + self.xi2q_map() def gqPtsWts(self): """ Gauss-Quadrature nodes and weights according to the gPCE rule """ - n=self.nSamp - if self.GQdistType=='Unif': - x=np.polynomial.legendre.leggauss(n) - self.xi=x[0] - self.w=x[1] - self.xiBound=[-1,1] - self.qBound=self.qInfo - elif self.GQdistType=='Norm': - x=np.polynomial.hermite_e.hermegauss(n) - self.xi=x[0] - self.w=x[1] - self.xiBound=[min(x[0]),max(x[0])] - self.mean=self.qInfo[0] - self.sdev=self.qInfo[1] + n = self.nSamp + if self.GQdistType == 'Unif': + x = np.polynomial.legendre.leggauss(n) + self.xi = x[0] + self.w = x[1] + self.xiBound = [-1, 1] + self.qBound = self.qInfo + elif self.GQdistType == 'Norm': + x = np.polynomial.hermite_e.hermegauss(n) + self.xi = x[0] + self.w = x[1] + self.xiBound = [min(x[0]), max(x[0])] + self.mean = self.qInfo[0] + self.sdev = self.qInfo[1] self.xi2q_map() def xi2q_map(self): """ Linearly map xi in Gamma to q in Q """ - xi_=self.xi - if (self.sampleType=='GQ' and self.GQdistType=='Norm') or \ - self.sampleType=='normRand': - self.q=self.qInfo[0]+self.qInfo[1]*xi_ - self.qBound=[min(self.q),max(self.q)] + xi_ = self.xi + if (self.sampleType == 'GQ' and self.GQdistType == 'Norm') or \ + self.sampleType == 'normRand': + self.q = self.qInfo[0]+self.qInfo[1]*xi_ + self.qBound = [min(self.q), max(self.q)] else: - xiBound_=self.xiBound - qBound_=self.qBound - self.q=(xi_-xiBound_[0])/(xiBound_[1]-xiBound_[0])*\ - (qBound_[1]-qBound_[0])+qBound_[0] + xiBound_ = self.xiBound + qBound_ = self.qBound + self.q = (xi_-xiBound_[0])/(xiBound_[1]-xiBound_[0]) *\ + (qBound_[1]-qBound_[0])+qBound_[0] # + + class testSample: R""" Generating test samples from a 1D paramter space using different methods. @@ -170,7 +176,7 @@ class testSample: Admissible range of `q` `nSamp`: int Number of samples to draw - + Attributes: `xi`: 1D numpy array of size `nSamp` Samples on the mapped space Gamma @@ -189,72 +195,80 @@ class testSample: ts5=testSample(sampleType='unifSpaced',GQdistType='Unif',qBound=[-1,3],nSamp=10) ts6=testSample(sampleType='GLL',qBound=[-1,3],nSamp=10) """ - def __init__(self,sampleType,qBound,nSamp,GQdistType='Unif',qInfo=[]): + + def __init__(self, sampleType, qBound, nSamp, GQdistType='Unif', qInfo=[]): self.info() - self.sampleType=sampleType - self.GQdistType=GQdistType - self.qInfo=qInfo + self.sampleType = sampleType + self.GQdistType = GQdistType + self.qInfo = qInfo self.check() - self.qBound=qBound - self.nSamp=nSamp + self.qBound = qBound + self.nSamp = nSamp self.genTestSamples() def info(self): - sampleTypeList=['GLL','unifSpaced','unifRand','normRand'] - GQdistList=['Unif','Norm'] #list of available distributions for gpce - self.sampleTypeList=sampleTypeList - self.GQdistList=GQdistList + sampleTypeList = ['GLL', 'unifSpaced', 'unifRand', 'normRand'] + # list of available distributions for gpce + GQdistList = ['Unif', 'Norm'] + self.sampleTypeList = sampleTypeList + self.GQdistList = GQdistList def check(self): if self.sampleType not in self.sampleTypeList: - raise KeyError('#ERROR @ testSample: Invalid sampleType! Choose from'\ - ,self.sampleTypeList) - if self.GQdistType=='Norm' and len(self.qInfo)==0: - raise KeyError("#ERROR @ testSample: qInfo is mandatory for GQdistType='Norm'") + raise KeyError( + '#ERROR @ testSample: Invalid sampleType! Choose from', self.sampleTypeList) + if self.GQdistType == 'Norm' and len(self.qInfo) == 0: + raise KeyError( + "#ERROR @ testSample: qInfo is mandatory for GQdistType='Norm'") - def genTestSamples(self): - n=self.nSamp - if self.GQdistType=='Unif': - self.xiBound=[-1,1] + def genTestSamples(self): + n = self.nSamp + if self.GQdistType == 'Unif': + self.xiBound = [-1, 1] + + if self.sampleType == 'unifSpaced': + q_ = np.linspace(self.qBound[0], self.qBound[1], n) + elif self.sampleType == 'GLL': + self.xiBound = [-1, 1] + xi_, w_ = nodes.gllPts(n) + q_ = (xi_ - self.xiBound[0]) / (self.xiBound[1] - self.xiBound[0]) * \ + (self.qBound[1] - self.qBound[0]) + self.qBound[0] + elif self.sampleType == 'unifRand': + if self.GQdistType != 'Unif': + raise ValueError( + "#ERROR @ testSample: sampleType 'unifRand' should be with GQdistType 'Unif' or ''.") + q_ = np.random.rand( + n)*(self.qBound[1]-self.qBound[0])+self.qBound[0] + q_ = np.sort(q_) + elif self.sampleType == 'normRand': + if self.GQdistType != 'Norm': + raise ValueError("#ERROR @ testSample: sampleType 'normRand' should be with\ + GQdistType 'Norm'.") + else: + xi_ = np.random.normal(size=n) + xi_ = np.sort(xi_) + self.xi = xi_ + q_ = self.qInfo[0]+xi_*self.qInfo[1] + self.q = q_ + self.q2xi_map() - if self.sampleType=='unifSpaced': - q_=np.linspace(self.qBound[0],self.qBound[1],n) - elif self.sampleType=='GLL': - self.xiBound=[-1,1] - xi_,w_=nodes.gllPts(n) - q_=xi_*(self.qBound[1]-self.qBound[0])+self.qBound[0] - elif self.sampleType=='unifRand': - if self.GQdistType!='Unif': - raise ValueError("#ERROR @ testSample: sampleType 'unifRand' should be with GQdistType 'Unif' or ''.") - q_=np.random.rand(n)*(self.qBound[1]-self.qBound[0])+self.qBound[0] - q_=np.sort(q_) - elif self.sampleType=='normRand': - if self.GQdistType!='Norm': - raise ValueError("#ERROR @ testSample: sampleType 'normRand' should be with\ - GQdistType 'Norm'.") - else: - xi_=np.random.normal(size=n) - xi_=np.sort(xi_) - self.xi=xi_ - q_=self.qInfo[0]+xi_*self.qInfo[1] - self.q=q_ - self.q2xi_map() - def q2xi_map(self): """ Linearly map q in Q to xi in Gamma """ - q_=self.q - qBound_=self.qBound - if self.GQdistType=='Norm': - xi_=(q_-self.qInfo[0])/self.qInfo[1] - self.xiBound=[min(xi_),max(xi_)] + q_ = self.q + qBound_ = self.qBound + if self.GQdistType == 'Norm': + xi_ = (q_-self.qInfo[0])/self.qInfo[1] + self.xiBound = [min(xi_), max(xi_)] else: - xi_=(self.xiBound[1]-self.xiBound[0])*(q_-qBound_[0])\ - /(qBound_[1]-qBound_[0])+self.xiBound[0] - self.xi=xi_ + xi_ = (self.xiBound[1]-self.xiBound[0])*(q_-qBound_[0])\ + / (qBound_[1]-qBound_[0])+self.xiBound[0] + self.xi = xi_ # -def LHS_sampling(n,xBound): + + +def LHS_sampling(n, xBound): """ LHS (Latin Hypercube) sampler from a p-D random variable distributed uniformly. Credits: https://zmurchok.github.io/2019/03/15/Latin-Hypercube-Sampling.html @@ -270,10 +284,10 @@ def LHS_sampling(n,xBound): `x`: 2D numpy array of size (n,p) Samples taken from the p-D space with ranges `xBound` """ - p=len(xBound) - x = np.random.uniform(size=[n,p]) + p = len(xBound) + x = np.random.uniform(size=[n, p]) for i in range(p): - x_ = (np.argsort(x[:,i])+0.5)/float(n) - x[:,i]=x_*(xBound[i][1]-xBound[i][0])+xBound[i][0] + x_ = (np.argsort(x[:, i])+0.5)/float(n) + x[:, i] = x_*(xBound[i][1]-xBound[i][0])+xBound[i][0] return x # diff --git a/UQit/sobol.py b/UQit/sobol.py index 4ca7d30..3b958c2 100644 --- a/UQit/sobol.py +++ b/UQit/sobol.py @@ -1,9 +1,9 @@ ########################################## # Sobol Sensitivity Indices ########################################## -#----------------------------------------- +# ----------------------------------------- # Saleh Rezaeiravesh, salehr@kth.se -#----------------------------------------- +# ----------------------------------------- """ * Parameters are assumed to be independent from each other * Parameters can have any arbitrary distribution. The discrete PDF should be imported to sobol(). @@ -17,6 +17,8 @@ import numpy as np from scipy.integrate import simps # + + class sobol: """ Computes Sobol sensitivity indices for a p-D parameter (p>1). @@ -24,7 +26,7 @@ class sobol: Assumptions: * Parameters are independent from each other * Up to second-order interactions are considered. - + Args: `q`: A list of length p q=[q1,q2,...,qp] where qi: 1D numpy array of size ni containing uniformly-spaced parameter samples @@ -32,7 +34,7 @@ class sobol: Response values at samples `q`. The tensor product with ordering='F' (Fortran-like) is presumed. `pdf`: List of length p of 1D numpy arrays The i-th array in the list contains the values of the PDF of q_i, where i=1,2,...,p - + Methods: compute(): Computes the Sobol indices @@ -47,10 +49,11 @@ class sobol: `STi`: A 1D numpy array of size p Total Sobol indices """ - def __init__(self,q,f,pdf): - self.q=q - self.f=f - self.pdf=pdf + + def __init__(self, q, f, pdf): + self.q = q + self.f = f + self.pdf = pdf self._info() self.comp() @@ -58,201 +61,204 @@ def _info(self): """ Settings and Checking """ - self.p=len(self.q) - if self.p<2: - raise ValueError("For Sobol indices at Least p==2 is required.") - self.n=[self.q[i].shape[0] for i in range(self.p)] + self.p = len(self.q) + if self.p < 2: + raise ValueError("For Sobol indices at Least p==2 is required.") + self.n = [self.q[i].shape[0] for i in range(self.p)] - def _permuteIndex(self,i): + def _permuteIndex(self, i): """ Make a permutation list for index `i` over index set {0,1,...,p-1} Example: for p=4, i=0: I=[1,2,3] i=2: I=[3,0,1] """ - indexList=np.arange(self.p) - I1=np.where(indexListi) - I=np.append(I2,I1) - return I + indexList = np.arange(self.p) + I1 = np.where(indexList < i) + I2 = np.where(indexList > i) + I = np.append(I2, I1) + return I def _dualInteractIndex(self): """ Makes a list of pairs of indices required for dual interactions * The indices in each pair are always ascending and non-repeated. """ - pairIndex=[] #Pairs of indices whose interactions to be found - indexList=np.arange(self.p) - while indexList.size>1: - iHead=indexList[0] - indexList=np.delete(indexList,0) + pairIndex = [] # Pairs of indices whose interactions to be found + indexList = np.arange(self.p) + while indexList.size > 1: + iHead = indexList[0] + indexList = np.delete(indexList, 0) for i in range(indexList.size): - pairIndex.append([iHead,indexList[i]]) - compIndex=[] #Complement to the pairs (The dropped indices on which integration is done) - indexList=np.arange(self.p) + pairIndex.append([iHead, indexList[i]]) + # Complement to the pairs (The dropped indices on which integration is done) + compIndex = [] + indexList = np.arange(self.p) for pair_ in pairIndex: - comp_=[] + comp_ = [] for ind_ in indexList: if ind_ not in pair_: - comp_.append(ind_) - compIndex.append(comp_) - return pairIndex,compIndex + comp_.append(ind_) + compIndex.append(comp_) + return pairIndex, compIndex + + def doubleInteg(self, g, x1, x2): + """ + Numerical double integration :math:`\int_{x2} \int_{x1} g(x1,x2) dx1 dx2` - def doubleInteg(self,g,x1,x2): - """ - Numerical double integration :math:`\int_{x2} \int_{x1} g(x1,x2) dx1 dx2` - - Args: - `g`: numpy array of shape (n1,n2) - Values of integrand over the grid of `x1`-`x2` - `x1`,`x2`: 1D numpy arrays of sizes n1 and n2 - Integrated variables + Args: + `g`: numpy array of shape (n1,n2) + Values of integrand over the grid of `x1`-`x2` + `x1`,`x2`: 1D numpy arrays of sizes n1 and n2 + Integrated variables - Returns: - `T`: scalar - Integral value - """ - T_=simps(g,x2,axis=1) - T=simps(T_,x1) - return T + Returns: + `T`: scalar + Integral value + """ + T_ = simps(g, x2, axis=1) + T = simps(T_, x1) + return T - def dualInteractTerm(self,fij_,fi_,fj_): + def dualInteractTerm(self, fij_, fi_, fj_): """ Computes 2nd-order interaction terms in the HDMR decomposition based on :math:`f_{ij}(q_i,q_j)=\int f(q)dq_{\sim{ij}}-f_i(q_i)-f_j(q_j)-f_0` """ - ni=fi_.size - nj=fj_.size - fij=np.zeros((ni,nj)) - f0=self.f0 + ni = fi_.size + nj = fj_.size + fij = np.zeros((ni, nj)) + f0 = self.f0 for i2 in range(nj): for i1 in range(ni): - fij[i1,i2]=fij_[i1,i2]-fi_[i1]-fj_[i2]-f0 + fij[i1, i2] = fij_[i1, i2]-fi_[i1]-fj_[i2]-f0 return fij def hdmr_0_1(self): """ Zeroth- and first-order HDMR decomposition """ - p=self.p - q=self.q - pdf=self.pdf + p = self.p + q = self.q + pdf = self.pdf - fi=[] + fi = [] for i in range(p): - parIndexList=np.arange(p) #Index of remaining parameters to integrate over - #permut list I=[] - I_permute=self._permuteIndex(i) #Permutation index list - f_=self.f - while I_permute.size>0: - iInteg_=I_permute[0] - iAxis_=np.where(parIndexList==iInteg_)[0][0] - #reshape the pdf for multiplying it with f_ - n_=[1]*f_.ndim - n_[iAxis_]=pdf[iInteg_].shape[0] - pdf_=pdf[iInteg_].reshape(n_) + # Index of remaining parameters to integrate over + parIndexList = np.arange(p) + # permut list I=[] + I_permute = self._permuteIndex(i) # Permutation index list + f_ = self.f + while I_permute.size > 0: + iInteg_ = I_permute[0] + iAxis_ = np.where(parIndexList == iInteg_)[0][0] + # reshape the pdf for multiplying it with f_ + n_ = [1]*f_.ndim + n_[iAxis_] = pdf[iInteg_].shape[0] + pdf_ = pdf[iInteg_].reshape(n_) - f_=f_*pdf_ - f_=simps(f_,q[iInteg_],axis=iAxis_) - parIndexList=np.delete(parIndexList,iAxis_) - I_permute=np.delete(I_permute,0) - fi.append(f_) - #compute the mean - f0=simps(fi[0]*pdf[0],q[0],axis=0) - fi=[fi[i]-f0 for i in range(p)] - self.f0=f0 - self.fi=fi + f_ = f_*pdf_ + f_ = simps(f_, q[iInteg_], axis=iAxis_) + parIndexList = np.delete(parIndexList, iAxis_) + I_permute = np.delete(I_permute, 0) + fi.append(f_) + # compute the mean + f0 = simps(fi[0]*pdf[0], q[0], axis=0) + fi = [fi[i]-f0 for i in range(p)] + self.f0 = f0 + self.fi = fi def hdmr_2(self): """ 2nd-order HDMR decomposition """ - pairIndex,compIndex=self._dualInteractIndex() - self.dualInteractIndex=pairIndex - pairNum=len(pairIndex) #number of dual interaction pairs - compLen=len(compIndex[0]) #number of indices in each complement set - q=self.q - pdf=self.pdf - fij=[] + pairIndex, compIndex = self._dualInteractIndex() + self.dualInteractIndex = pairIndex + pairNum = len(pairIndex) # number of dual interaction pairs + compLen = len(compIndex[0]) # number of indices in each complement set + q = self.q + pdf = self.pdf + fij = [] for k in range(pairNum): - indexList=np.arange(self.p) - f_=self.f + indexList = np.arange(self.p) + f_ = self.f for l in range(compLen): - iInteg_=compIndex[k][l] - iAxis_=np.where(indexList==iInteg_)[0][0] - if (self.p>2): - #reshape the pdf for multiplying it with f_ - n_=[1]*f_.ndim - n_[iAxis_]=pdf[iInteg_].shape[0] - pdf_=pdf[iInteg_].reshape(n_) - f_=f_*pdf_ - f_=simps(f_,q[iInteg_],axis=iAxis_) + iInteg_ = compIndex[k][l] + iAxis_ = np.where(indexList == iInteg_)[0][0] + if (self.p > 2): + # reshape the pdf for multiplying it with f_ + n_ = [1]*f_.ndim + n_[iAxis_] = pdf[iInteg_].shape[0] + pdf_ = pdf[iInteg_].reshape(n_) + f_ = f_*pdf_ + f_ = simps(f_, q[iInteg_], axis=iAxis_) else: - f_=self.f - indexList=np.delete(indexList,iAxis_) - i=pairIndex[k][0] - j=pairIndex[k][1] - fi_=self.fi[i] - fj_=self.fi[j] - fij_=self.dualInteractTerm(f_,fi_,fj_) + f_ = self.f + indexList = np.delete(indexList, iAxis_) + i = pairIndex[k][0] + j = pairIndex[k][1] + fi_ = self.fi[i] + fj_ = self.fi[j] + fij_ = self.dualInteractTerm(f_, fi_, fj_) fij.append(fij_) - self.fij=fij + self.fij = fij - def hdmr(self): + def hdmr(self): """ HDMR decomposition """ - #Mean and 1st-order Sobol decomposition terms + # Mean and 1st-order Sobol decomposition terms self.hdmr_0_1() - #2nd-order Interactions + # 2nd-order Interactions self.hdmr_2() def partVariance(self): """ Partial variances and Sobol indices """ - q=self.q - pdf=self.pdf - fi=self.fi - fij=self.fij - #1st-order terms, Di - Di=[] + q = self.q + pdf = self.pdf + fi = self.fi + fij = self.fij + # 1st-order terms, Di + Di = [] for i in range(self.p): - Di_=simps(fi[i]**2.*pdf[i],q[i]) + Di_ = simps(fi[i]**2.*pdf[i], q[i]) Di.append(Di_) - D=sum(Di) - #2nd-order terms, Dij - Dij=[] - SijName=[] - #Main Sobol indices + D = sum(Di) + # 2nd-order terms, Dij + Dij = [] + SijName = [] + # Main Sobol indices for k in range(len(fij)): - i=self.dualInteractIndex[k][0] - j=self.dualInteractIndex[k][1] - pdf_=pdf[i][:,None]*pdf[j] - Dij.append(self.doubleInteg(np.multiply(fij[k]**2.,pdf_),q[i],q[j])) + i = self.dualInteractIndex[k][0] + j = self.dualInteractIndex[k][1] + pdf_ = pdf[i][:, None]*pdf[j] + Dij.append(self.doubleInteg( + np.multiply(fij[k]**2., pdf_), q[i], q[j])) SijName.append('S'+str(i+1)+str(j+1)) - D+=sum(Dij) - self.Si=Di/D - self.Sij=Dij/D - self.SijName=SijName - #Total Sobol indices - Sij_sum=[] + D += sum(Dij) + self.Si = Di/D + self.Sij = Dij/D + self.SijName = SijName + # Total Sobol indices + Sij_sum = [] for l in range(self.p): - sum_=0 + sum_ = 0 for k in range(len(fij)): - i=self.dualInteractIndex[k][0] - j=self.dualInteractIndex[k][1] - if i==l or j==l: - sum_+=self.Sij[k] - Sij_sum.append(sum_) - self.STi=self.Si+Sij_sum + i = self.dualInteractIndex[k][0] + j = self.dualInteractIndex[k][1] + if i == l or j == l: + sum_ += self.Sij[k] + Sij_sum.append(sum_) + self.STi = self.Si+Sij_sum def comp(self): """ Computes Sobol indices based on HDMR """ - #HDMR (Sobol) decomposition of the model function + # HDMR (Sobol) decomposition of the model function self.hdmr() - #Partial Variances & Sobol indices + # Partial Variances & Sobol indices self.partVariance() # diff --git a/UQit/stats.py b/UQit/stats.py index 7947e8d..ac32304 100644 --- a/UQit/stats.py +++ b/UQit/stats.py @@ -11,11 +11,13 @@ import UQit.write as writeUQ # # -def pdfFit_uniVar(f,doPlot,pwOpts): + + +def pdfFit_uniVar(f, doPlot, pwOpts): """ Fits a PDF to samples `f` and plots both histogram and the fitted PDF. As an option, the plots and data can be saved on disk. - + Args: `f`: 1D numpy array of size n Samples @@ -28,60 +30,62 @@ def pdfFit_uniVar(f,doPlot,pwOpts): * 'header': string, the header of the dumped file * 'iLoc': int, After converting to string will be added to the `figName` """ - if f.ndim>1: - print('Note: input f to pdfFit_uniVar(f) must be a 1D numpy array. We reshape f!') - nTot=1 - for i in range(f.ndim): - nTot*=f.shape[i] - f=np.reshape(f,nTot) - #fit kde + if f.ndim > 1: + print('Note: input f to pdfFit_uniVar(f) must be a 1D numpy array. We reshape f!') + nTot = 1 + for i in range(f.ndim): + nTot *= f.shape[i] + f = np.reshape(f, nTot) + # fit kde kde = sm.nonparametric.KDEUnivariate(f) kde.fit() - #plot + # plot if doPlot: - plt.figure(figsize=(10,4)); - ax=plt.gca(); - plt.plot(kde.support,kde.density,'-r',lw=2) - binsNum='auto' #70 - BIN=plt.hist(f,bins=binsNum,density=True,color='steelblue',alpha=0.4,edgecolor='b') - plt.xticks(fontsize=15) - plt.yticks(fontsize=15) - plt.grid(alpha=0.4) - if pwOpts: #if not empty - #Dump the data of the PDf - if pwOpts['figDir']: - figDir=pwOpts['figDir'] - wrtDir=figDir+'/dumpData/' - if pwOpts['figName']: - outName=pwOpts['figName'] - if pwOpts['iLoc']: - outName+='_'+str(pwOpts['iLoc']) - if not os.path.exists(figDir): - os.makedirs(figDir) - if not os.path.exists(wrtDir): - os.makedirs(wrtDir) - F1=open(wrtDir+outName+'.dat','w') - if pwOpts['header']: - F1.write('# '+pwOpts['header']+'\n') - F1.write('# kde.support \t\t kde.density \n') - F1.write('# '+writeUQ.printRepeated('-',50)+'\n') - for i in range(len(kde.support)): - F1.write('%g \t %g \n' %(kde.support[i],kde.density[i])) - F1.write('# '+writeUQ.printRepeated('-',50)+'\n') - F1.write('# bin.support \t\t bin.density \n') - F1.write('# '+writeUQ.printRepeated('-',50)+'\n') - for i in range(len(BIN[0])): - F1.write('%g \t %g \n' %(BIN[1][i],BIN[0][i])) - #Save the figure of the PDF - fig = plt.gcf() - DPI = fig.get_dpi() - fig.set_size_inches(800/float(DPI),400/float(DPI)) - plt.savefig(figDir+outName+'.pdf',bbox_inches='tight') - else: - plt.show() + plt.figure(figsize=(10, 4)) + ax = plt.gca() + plt.plot(kde.support, kde.density, '-r', lw=2) + binsNum = 'auto' # 70 + BIN = plt.hist(f, bins=binsNum, density=True, + color='steelblue', alpha=0.4, edgecolor='b') + plt.xticks(fontsize=15) + plt.yticks(fontsize=15) + plt.grid(alpha=0.4) + if pwOpts: # if not empty + # Dump the data of the PDf + if pwOpts['figDir']: + figDir = pwOpts['figDir'] + wrtDir = figDir+'/dumpData/' + if pwOpts['figName']: + outName = pwOpts['figName'] + if pwOpts['iLoc']: + outName += '_'+str(pwOpts['iLoc']) + if not os.path.exists(figDir): + os.makedirs(figDir) + if not os.path.exists(wrtDir): + os.makedirs(wrtDir) + F1 = open(wrtDir+outName+'.dat', 'w') + if pwOpts['header']: + F1.write('# '+pwOpts['header']+'\n') + F1.write('# kde.support \t\t kde.density \n') + F1.write('# '+writeUQ.printRepeated('-', 50)+'\n') + for i in range(len(kde.support)): + F1.write('%g \t %g \n' % (kde.support[i], kde.density[i])) + F1.write('# '+writeUQ.printRepeated('-', 50)+'\n') + F1.write('# bin.support \t\t bin.density \n') + F1.write('# '+writeUQ.printRepeated('-', 50)+'\n') + for i in range(len(BIN[0])): + F1.write('%g \t %g \n' % (BIN[1][i], BIN[0][i])) + # Save the figure of the PDF + fig = plt.gcf() + DPI = fig.get_dpi() + fig.set_size_inches(800/float(DPI), 400/float(DPI)) + plt.savefig(figDir+outName+'.pdf', bbox_inches='tight') + else: + plt.show() return kde -def pdfPredict_uniVar(f,fTest,doPlot): + +def pdfPredict_uniVar(f, fTest, doPlot): """ Evaluates the continuous PDF fitted to `f` at `fTest`. @@ -89,16 +93,16 @@ def pdfPredict_uniVar(f,fTest,doPlot): `f`: 1D numpy array `fTest`: List of length m - + Returns: `pdfPred`: 1D numpy array of size m """ - #Fit the PDF to f - kde=pdfFit_uniVar(f,doPlot,{}) - #Evaluate the PDF at f0 - pdfPred=[] + # Fit the PDF to f + kde = pdfFit_uniVar(f, doPlot, {}) + # Evaluate the PDF at f0 + pdfPred = [] for i in range(len(fTest)): pdfPred.append(kde.evaluate(fTest[i])[0]) - pdfPred=np.asarray(pdfPred) + pdfPred = np.asarray(pdfPred) return pdfPred -# +# diff --git a/UQit/surr2surr.py b/UQit/surr2surr.py index e700a90..58169fe 100644 --- a/UQit/surr2surr.py +++ b/UQit/surr2surr.py @@ -1,9 +1,9 @@ ############################################################### # Interpolation from a surrogate to another surrogate ############################################################### -#-------------------------------------------------------------- +# -------------------------------------------------------------- # Saleh Rezaeiravesh, salehr@kth.se -#-------------------------------------------------------------- +# -------------------------------------------------------------- # import os import sys @@ -12,7 +12,9 @@ import UQit.reshaper as reshaper # # -def lagIntAtGQs(fValM1,qM1,spaceM1,nM2,spaceM2,distType): + + +def lagIntAtGQs(fValM1, qM1, spaceM1, nM2, spaceM2, distType): """ Given response values `fValM1` at `nM1` arbitrary samples over the p-D `spaceM1`, the values at `nM2` Gauss quadrature points over `spaceM2` are computed using Lagrange interpolation. @@ -51,35 +53,37 @@ def lagIntAtGQs(fValM1,qM1,spaceM1,nM2,spaceM2,distType): `fValM2`: 1D numpy array of size (nM1_1*nM1_2*...*nM1_p). Interpolated response values at `xiM2` """ - #(1) Check the inputs - ndim=len(spaceM1) - if (ndim!=len(spaceM2) or ndim!=len(qM1)): - raise ValueError('SpaceM1 and SpaceM2 should have the same dimension, p.') + # (1) Check the inputs + ndim = len(spaceM1) + if (ndim != len(spaceM2) or ndim != len(qM1)): + raise ValueError( + 'SpaceM1 and SpaceM2 should have the same dimension, p.') for idim in range(ndim): - d1=spaceM1[idim][1]-spaceM1[idim][0] - d2=spaceM2[idim][1]-spaceM2[idim][0] - if (d2>d1): - print("Wrong parameter range in direction ",ldim+1) - raise ValueError("||spaceM2|| should be smaller than ||spaceM1||.") - #(2) Construct the Gauss-quadrature stochastic samples for model2 - qM2=[] - xiM2=[] + d1 = spaceM1[idim][1]-spaceM1[idim][0] + d2 = spaceM2[idim][1]-spaceM2[idim][0] + if (d2 > d1): + print("Wrong parameter range in direction ", ldim+1) + raise ValueError("||spaceM2|| should be smaller than ||spaceM1||.") + # (2) Construct the Gauss-quadrature stochastic samples for model2 + qM2 = [] + xiM2 = [] for i in range(ndim): - xi_,w=pce.gqPtsWts(nM2[i],distType[i]) - qM2_=pce.mapFromUnit(xi_,spaceM2[i]) + xi_, w = pce.gqPtsWts(nM2[i], distType[i]) + qM2_ = pce.mapFromUnit(xi_, spaceM2[i]) qM2.append(qM2_) xiM2.append(xi_) - if (ndim==1): - qM2=qM2[0] - xiM2=xiM2[0] - elif (ndim>1): - xiM2=reshaper.vecs2grid(xiM2) - #(3) Use lagrange interpolation to find values at q2, given fVal1 at q1 - if ndim==1: - fVal2Interp=lagInt(fNodes=fValM1,qNodes=[qM1[0]],qTest=[qM2]).val - elif (ndim>1): - fVal2Interp_=lagInt(fNodes=fValM1,qNodes=qM1,qTest=qM2,liDict={'testRule':'tensorProd'}).val - nM2_=fVal2Interp_.size - fVal2Interp=fVal2Interp_.reshape(nM2_,order='F') - return qM2,xiM2,fVal2Interp + if (ndim == 1): + qM2 = qM2[0] + xiM2 = xiM2[0] + elif (ndim > 1): + xiM2 = reshaper.vecs2grid(xiM2) + # (3) Use lagrange interpolation to find values at q2, given fVal1 at q1 + if ndim == 1: + fVal2Interp = lagInt(fNodes=fValM1, qNodes=[qM1[0]], qTest=[qM2]).val + elif (ndim > 1): + fVal2Interp_ = lagInt(fNodes=fValM1, qNodes=qM1, qTest=qM2, liDict={ + 'testRule': 'tensorProd'}).val + nM2_ = fVal2Interp_.size + fVal2Interp = fVal2Interp_.reshape(nM2_, order='F') + return qM2, xiM2, fVal2Interp # diff --git a/UQit/write.py b/UQit/write.py index 11a4b84..004cddf 100644 --- a/UQit/write.py +++ b/UQit/write.py @@ -2,7 +2,7 @@ # Tools for printing or writing data in file ##################################################### # Saleh Rezaeiravesh, salehr@kth.se -#---------------------------------------------------- +# ---------------------------------------------------- # def printRepeated(string_to_expand, length): """