Blog: Amélioration du chaos - Postprocessing - « Pretty print » des résultats: exemple_chaos_printchaos.py

File exemple_chaos_printchaos.py, 4.8 KB (added by souchaud@…, 8 years ago)
Line 
1import openturns as ot
2
3import printchaos as pc
4
5# Definir le nombre de variables aleatoires
6NbVA = 5
7
8# Liste pour sauvegarder le plan d'experience en attendant getPEX
9Plan_exp = []
10
11
12class myfunction(ot.OpenTURNSPythonFunction):
13    def __init__(self):
14        ot.OpenTURNSPythonFunction.__init__(self, NbVA, 1)
15    def f(self, X):
16        I=X[1]*X[2]**3/12.
17        v=X[4]*X[0]**3/(48.*X[3]*I)
18        #Sauvegarde du plan d'experience en attendant getPEX
19        Plan_exp.append([X[0],X[1],X[2],X[3],X[4],v])
20        return [v]
21
22####
23# Modele physique = NumericalMathFunction
24####
25PhysModel = ot.NumericalMathFunction(myfunction())
26
27
28# Definition de listes de variables et d'une liste de leur nom respectif
29L = []
30des = []
31L.append(ot.LogNormal(5.,0.02,0.,ot.LogNormal.MU_SIGMAOVERMU))
32des.append("L")
33L.append(ot.LogNormal(0.2,0.05,0.,ot.LogNormal.MU_SIGMAOVERMU))
34des.append("b")
35L.append(ot.LogNormal(0.4,0.05,0.,ot.LogNormal.MU_SIGMAOVERMU))
36des.append("h")
37L.append(ot.LogNormal(30000.,0.12,0.,ot.LogNormal.MU_SIGMAOVERMU))
38des.append("E")
39L.append(ot.LogNormal(0.1,0.2,0.,ot.LogNormal.MU_SIGMAOVERMU))
40des.append("F")
41
42# Automatisation de la creation de la DistributionCollection
43CollectionD=ot.DistributionCollection(NbVA)
44for i in range(NbVA):
45    L[i].setName(des[i])
46    CollectionD[i]=ot.Distribution(L[i])
47
48# Copule
49IC = ot.IndependentCopula(NbVA)
50
51####
52# Modele probabiliste = ComposedDistribution
53####
54ProbaModel = ot.ComposedDistribution(CollectionD,ot.Copula(IC))
55
56
57# Ordre de developpement du chaos
58ordre = 3
59
60# Strategie :
61# --> 0 FixedStrategy
62# --> 1 SequentialStrategy
63# --> 2 CleaningStrategy
64Strategie = 0
65
66# Type de plan d'experience & nombre de simulations
67# --> 0 Monte Carlo
68# --> 1 LHS
69PEX = 1
70N1 = 100
71
72
73
74# Base des polynomes
75PolynomialFamilyColl = ot.PolynomialFamilyCollection(NbVA)
76for i in range(NbVA):
77# Exponential, Gumbel, Normal, LogNormal, Weibull ==> Hermite
78    if ((ProbaModel.getMarginal(i).getClassName() == "Exponential") or \
79        (ProbaModel.getMarginal(i).getClassName() == "Gumbel") or \
80        (ProbaModel.getMarginal(i).getClassName() == "Normal") or \
81        (ProbaModel.getMarginal(i).getClassName() == "LogNormal") or \
82        (ProbaModel.getMarginal(i).getClassName() == "Weibull")):
83        PolynomialFamilyColl[i] = ot.OrthogonalUniVariatePolynomialFamily(ot.HermiteFactory())
84# Uniform ==> Legendre
85    if (ProbaModel.getMarginal(i).getClassName() == "Uniform"):
86        PolynomialFamilyColl[i] = ot.OrthogonalUniVariatePolynomialFamily(ot.LegendreFactory())
87# Beta(b+1,a+b+2,binf,bsup) ==> Jacobi(a,b)
88    if (ProbaModel.getMarginal(i).getClassName() == "Beta"):
89        #b = Modelproba.getMarginal(i).getR()-1.0
90        #a = Modelproba.getMarginal(i).getT()-b-2.0
91        b = ProbaModel.getMarginal(i).getParametersCollection()[0][0]-1.0
92        a = ProbaModel.getMarginal(i).getParametersCollection()[0][1]-b-2.0
93        PolynomialFamilyColl[i] = ot.OrthogonalUniVariatePolynomialFamily(ot.JacobiFactory(a,b))
94# Gamma(k+1,lambda,gamma) ==> Laguerre(k)
95    if (ProbaModel.getMarginal(i).getClassName() == "Gamma"):
96        k = ProbaModel.getMarginal(i).getK()-1.0
97        PolynomialFamilyColl[i] = ot.OrthogonalUniVariatePolynomialFamily(ot.LaguerreFactory(k))
98
99multivariateBasis = ot.OrthogonalProductPolynomialFactory(PolynomialFamilyColl, ot.EnumerateFunction(NbVA))
100
101
102
103# Nombre de polynomes de degre <= ordre
104if (Strategie == 0):
105    p = ot.EnumerateFunction(NbVA).getStrateCumulatedCardinal(ordre)
106    truncatureBasisStrategy = ot.FixedStrategy(ot.OrthogonalBasis(multivariateBasis), p)
107if (Strategie == 1):
108    maximumCardinalBasis = N1/2
109    truncatureBasisStrategy = ot.SequentialStrategy(ot.OrthogonalBasis(multivariateBasis), maximumCardinalBasis)
110if (Strategie == 2):
111    maximumConsideredTerms = N1/2
112    mostSignificant = 10
113    significanceFactor = 1.0e-4
114    truncatureBasisStrategy = ot.CleaningStrategy(ot.OrthogonalBasis(multivariateBasis), maximumConsideredTerms, mostSignificant, significanceFactor, True)
115
116
117
118if (PEX == 0):
119    EvaluationCoeffStrategy = ot.LeastSquaresStrategy(ot.MonteCarloExperiment(N1))
120if (PEX == 1):
121    EvaluationCoeffStrategy = ot.LeastSquaresStrategy(ot.LHSExperiment(N1))
122
123
124
125polynomialChaosAlgorithm = ot.FunctionalChaosAlgorithm(PhysModel, ot.Distribution(ProbaModel),
126                                                 ot.AdaptiveStrategy(truncatureBasisStrategy),
127                                                 ot.ProjectionStrategy(EvaluationCoeffStrategy))
128polynomialChaosAlgorithm.run()
129polynomialChaosResult = polynomialChaosAlgorithm.getResult()
130
131####
132# Results
133####
134metaModel = polynomialChaosResult.getMetaModel()
135
136# pretty print
137pc.print_moments(polynomialChaosResult)
138pc.print_coefficients(polynomialChaosResult, PM = ProbaModel)
139pc.print_plot_sensitivities(polynomialChaosResult, PM = ProbaModel)
140# recuperation de l'echantillon de sortie dans la fonction (utilisable 0.13.2)
141pc.print_plot_regression(Plan_exp, polynomialChaosResult)
142