Posts by author M. Berveiller EDF

Convergence des moments statistiques de la fonction d'Ishigami en fonction de la taille du plan d'expériences avec le chaos polynomial

Testé avec la version 0.12.3 d'OpenTURNS

# Chargement des librairies
from openturns import *
from math import *
from openturns_viewer import ViewImage
 
# Initialisation de la graîne pour les simulations de Monte-Carlo
RandomGenerator().SetSeed(0)
 
# Définition de la fonction factorielle
def factorielle(x):
	if x==0:
		return 1
	else:
		return x*factorielle(x-1)
 
# Définition des constantes intervenant dans le fonction d'Ishigami
A = 7.0
B = 0.1 
 
#Modele physique : fonction d'Ishigami
class myfunction(OpenTURNSPythonFunction):
	def __init__(self):
		OpenTURNSPythonFunction.__init__(self, 3, 1)
	def f(self, X):
		H = sin(X[0]) + A*(sin(X[1])**2) + B*(X[2]**4)*sin(X[0])
		Y=[H]
		return Y
 
function=myfunction()
g = NumericalMathFunction(function)
 
# Nombre de variables aléatoires d'entrées
dim = 3
 
#Modele probababiliste
mean=NumericalPoint(3)
mean[0]=0.0
mean[1]=0.0
mean[2]=0.0
 
L1=Uniform(-pi,pi)
L2=Uniform(-pi,pi)
L3=Uniform(-pi,pi)
 
Collection=DistributionCollection(dim)
Collection[0]=Distribution(L1)
Collection[1]=Distribution(L2)
Collection[2]=Distribution(L3)
 
Copule=IndependentCopula(dim)
Modelproba=ComposedDistribution(Collection,Copula(Copule))
 
# Definition des polynômes unidimensionels du chaos polynomial
legendre = OrthogonalUniVariatePolynomialFamily(LegendreFactory())
polyColl = PolynomialFamilyCollection(dim)
polyColl[0]=legendre
polyColl[1]=legendre
polyColl[2]=legendre
 
# Affichage d'un polynôme de la base
enumFunc = EnumerateFunction(dim)
MultiIndice_1 = enumFunc(0)
print "Polynome 1 de la base : ", MultiIndice_1
 
productPoly = OrthogonalProductPolynomialFactory(polyColl,enumFunc)
basis = OrthogonalBasis(productPoly)
 
mesure = basis.getMeasure()
################################################################	
# Calcul du nombre de termes dans le dvpt sur CP
# On fixe le degré du chaos polynomial
p=5
P=factorielle(dim+p)/(factorielle(dim)*factorielle(p))
print "Taille du chaos polynomial pour p=3 : ", P
 
 
################################################################	
nb_calculs=range(0,3)
 
moyenne=[0]
variance=[0]
taille_plan=[0]
taille_base_finale=[0]
for i in nb_calculs:
	sampleSize = P+10*i
	taille_plan.append(sampleSize)
	print "Taille du plan d experiences : ", sampleSize
	regressionStrategy = RegressionStrategy(MonteCarloExperiment(sampleSize))
	#regressionStrategy = RegressionStrategy(LHSExperiment(sampleSize))		
 
	#Fixed Strategy
	basisSize = P
	#print "Taille de la base", basisSize
	fixedStrategy= FixedStrategy(basis,basisSize)
	#polynomialChaosAlgorithm = FunctionalChaosAlgorithm(g,Distribution(Modelproba),AdaptiveStrategy(fixedStrategy),ProjectionStrategy(regressionStrategy))
 
	#SequentialStrategy:	
	sequentialStrategy=SequentialStrategy(basis,basisSize)	
	#polynomialChaosAlgorithm = FunctionalChaosAlgorithm(g,Distribution(Modelproba),AdaptiveStrategy(sequentialStrategy),ProjectionStrategy(regressionStrategy))
 
	#CleaningStrategy:		
 
	maxIndex = 3*P #Numero du dernier polynome Psi que l on peut utiliser
	basisSize_cleaning = P #ici le basisSise_cleaning est la taille max de la base toleree
	significanceFactor = 1.0e-3	
	cleaningStrategy = CleaningStrategy(basis,maxIndex,basisSize_cleaning , significanceFactor,True)	
	polynomialChaosAlgorithm = FunctionalChaosAlgorithm(g,Distribution(Modelproba),AdaptiveStrategy(cleaningStrategy),ProjectionStrategy(regressionStrategy))
 
 
	polynomialChaosAlgorithm.run()
 
	#result	
	polynomialChaosResult = polynomialChaosAlgorithm.getResult()	
 
	#Get the polynomial chaos coefficients	
	coefficients = polynomialChaosResult.getCoefficients()	
 
	#Get the metamodel as a NumericalMathFunction	
	metaModel = polynomialChaosResult.getMetaModel()	
 
	#Get the basis		
	basisCollection = polynomialChaosResult.getBasis()	
	print "Taille de la base finale : ", polynomialChaosResult.getBasis().getSize()
 
	taille_base_finale.append(polynomialChaosResult.getBasis().getSize())
	#Get the distribution of variables U	
	mesure = polynomialChaosResult.getMeasure()	
 
	#Get the composedmodel which is the model of the reduced variables U	
	composedModel = polynomialChaosResult.getComposedModel	
 
 
	# Calcul des moments a partir des coefficients
	print "Moyenne a partir des coefficients : " , polynomialChaosResult.getCoefficients()[0]
	print "Erreur sur la moyenne : " , (polynomialChaosResult.getCoefficients()[0]-3.5)/3.5
	moyenne.append(polynomialChaosResult.getCoefficients()[0])
	# Calcul de la variance a partir des coefficients
	variance_coeffs=polynomialChaosResult.getCoefficients().norm2()-	polynomialChaosResult.getCoefficients()[0]**2
	print  "variance a partir des coefficients : " , variance_coeffs
	VarianceTheorique = (A**2)/8.0 + (B)*(pi**4)/5.0 + (B**2)*(pi**8)/18.0 +1.0/2.0
	print  "Erreur sur la variance : " , (variance_coeffs-VarianceTheorique)/VarianceTheorique,"\n"
	#print  "ecart-type a partir des coefficients : " , sqrt(variance_coeffs)
	variance.append(variance_coeffs)
 
fic_cv=open('cv_Ishigami_cleaning_5','w')
for i in nb_calculs:
	fic_cv.write(str(taille_plan[i+1])+' '+str(taille_base_finale[i+1])+' '+str(moyenne[i+1])+' '+str(variance[i+1])+'\n')
fic_cv.close()
 
 
print moyenne
print variance
# Variance theorique
VarianceTheorique = (A**2)/8.0 + (B)*(pi**4)/5.0 + (B**2)*(pi**8)/18.0 +1.0/2.0
print "variance theorique : ", VarianceTheorique
 
fic_coeff=open('coefficients_Ishigami_cleaning_5','w')
nb_coeffs =  polynomialChaosResult.getCoefficients().getDimension()
affichage_coefficients=range(0,nb_coeffs)
for i in affichage_coefficients:
	#print enumFunc(i), polynomialChaosResult.getCoefficients()[i]
	fic_coeff.write(str(enumFunc(i))+' '+str(polynomialChaosResult.getCoefficients()[i])+'\n')
fic_coeff.close()