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

File printchaos.py, 8.3 KB (added by souchaud@…, 8 years ago)
Line 
1import openturns as ot
2from math import sqrt
3import itertools as it
4import numpy as np
5from scipy import linalg
6
7
8def powerset(iterable):
9    s = list(iterable)
10    return it.chain.from_iterable(it.combinations(s, r) for r in range(len(s)+1))
11
12
13
14def reprn(x):
15
16    if x>=0:
17        return "%6.2e" %x
18    else:
19        return "%7.2e" %x
20
21def reppc(x):
22    x = 100.0*x
23    if x>1e-5:
24        return "%1.4f %%" %x
25    elif x<-1e-5:
26        return "%2.4f %%" %x
27    else:
28        return "0 %"
29
30
31def print_moments(res):
32
33    result = ot.FunctionalChaosRandomVector(res)
34    mean = result.getMean()[0]
35    variance = result.getCovariance()[0, 0]
36    stdev = sqrt(variance)
37
38    print "\nMOMENTS:"
39    print "Mean", reprn(mean)
40    print "Standard deviation", reprn(stdev)
41
42
43
44def print_coefficients(res, PM = None):
45
46    Pol = []
47    strpoly = ""
48    NbVA = res.getMetaModel().getInputDimension()
49    if PM != None:
50        for i in range(NbVA):
51            if ((PM.getMarginal(i).getClassName()=="Normal")or(PM.getMarginal(i).getClassName()=="LogNormal")or(PM.getMarginal(i).getClassName()=="Exponential")or(PM.getMarginal(i).getClassName()=="Weibull")or(PM.getMarginal(i).getClassName()=="Gumbel")or(PM.getMarginal(i).getClassName()=="Logistic")or(PM.getMarginal(i).getClassName()=="Rayleigh")):
52                Pol.append('He')
53            if (PM.getMarginal(i).getClassName()=="Uniform"):
54                Pol.append('Le')
55            if (PM.getMarginal(i).getClassName()=="Beta"):
56                Pol.append('Ja')
57            if (PM.getMarginal(i).getClassName()=="Gamma"):
58                Pol.append('La')
59
60        for i in range(len(Pol)):
61            strpoly = strpoly + Pol[i] + " "
62
63    Coefficients = res.getCoefficients()
64    Indices = res.getIndices()
65
66    print "\nPOLYNOMIALS: " + strpoly
67    for i in range(Coefficients.getDimension()):
68        j = Indices[i]
69        if Coefficients[i]>=0:
70            print "Polynomial:", ot.EnumerateFunction(NbVA)(j), "- Coefficient value ", reprn(Coefficients[i])
71        else:
72            print "Polynomial:", ot.EnumerateFunction(NbVA)(j), "- Coefficient value", reprn(Coefficients[i])
73
74
75
76def print_plot_sensitivities(res, PM = None):
77
78    print "\nSENSITIVITIES (Sobol' indices):"
79    FCRV = ot.FunctionalChaosRandomVector(res)
80    NbVA = FCRV.getFunction().getInputDimension()
81    TotalSI = 0
82    SobolTotalIndex = []
83
84    li = []
85    for i in range(NbVA):
86        li.append(i)
87    code = list(powerset(li))
88
89    for i in range(len(code)):
90        if i != 0:
91
92            if len(code[i]) !=len(code[i-1]) and i!=1:
93                print ""
94
95            indices = ot.Indices(len(code[i]))
96            for j in range(len(code[i])):
97                indices[j] = code[i][j]
98
99            if len(code[i]) == 1:
100                if PM != None:
101                    print "Order:",len(code[i]),"- variables",PM.getMarginal(indices[j]).getName(),indices,':',reppc(FCRV.getSobolIndex(indices))
102                else:
103                    print "Order:",len(code[i]),"- variables",str(indices),':',reppc(FCRV.getSobolIndex(indices))
104            else:
105                print "Order:",len(code[i]),"- variables",str(indices),':',reppc(FCRV.getSobolIndex(indices))
106
107    print ""
108    if PM != None:
109        for i in range(NbVA):
110            SobolTotalIndex.append(FCRV.getSobolTotalIndex(i))
111            TotalSI = TotalSI + SobolTotalIndex[i]
112            print "Total Sobol' indice for the variable", PM.getMarginal(i).getName(), ":", reppc(SobolTotalIndex[i])
113    else:
114        for i in range(NbVA):
115            SobolTotalIndex.append(FCRV.getSobolTotalIndex(i))
116            TotalSI = TotalSI + SobolTotalIndex[i]
117            print "Total Sobol' indice for the variable", i ,":", reppc(SobolTotalIndex[i])
118
119    sample_ = []
120    for i in range(NbVA):
121        sample_.append(SobolTotalIndex[i]/TotalSI)
122
123    sample = ot.NumericalSample(NbVA,1)
124    for i in range(NbVA):
125        sample[i][0] = sample_[i]
126
127    toutes_les_couleurs=["green","red","blue","yellow","darkblue","orange",
128    "lightgreen","darkcyan","cyan","magenta","darkgreen","violet","brown",
129    "darkred","pink","ivory","gold","darkgrey","grey","white","aliceblue",
130    "antiquewhite","aquamarine","azure","beige","bisque","black","blanchedalmond",
131    "blueviolet","burlywood","cadetblue","chartreuse","chocolate","coral",
132    "cornflowerblue","cornsilk","darkgoldenrod","darkgray","darkkhaki",
133    "darkmagenta","darkolivegreen","darkorange","darkorchid","darksalmon",
134    "darkseagreen","darkslateblue","darkslategray","darkslategrey","darkturquoise",
135    "darkviolet","deeppink","deepskyblue","dimgray","dimgrey","dodgerblue",
136    "firebrick","floralwhite","forestgreen","gainsboro","ghostwhite","goldenrod",
137    "gray","greenyellow","grey0","honeydew","hotpink","indianred","khaki",
138    "lavender","lavenderblush","lawngreen","lemonchiffon","lightblue",
139    "lightcoral","lightcyan","lightgoldenrod","lightgoldenrodyellow","lightgray",
140    "lightgrey","lightpink","lightsalmon","lightseagreen","lightskyblue","lightslateblue",
141    "lightslategray","lightslategrey","lightsteelblue","lightyellow","limegreen",
142    "linen","maroon","mediumaquamarine","mediumblue","mediumorchid","mediumpurple",
143    "mediumseagreen","mediumslateblue","mediumspringgreen","mediumturquoise",
144    "mediumvioletred","midnightblue","mintcream","mistyrose","moccasin",
145    "navajowhite","navy","navyblue","oldlace","olivedrab","orangered",
146    "orchid","palegoldenrod","palegreen","paleturquoise","palevioletred",
147    "papayawhip","peachpuff","peru","plum","powderblue","purple","rosybrown",
148    "royalblue","saddlebrown","salmon","sandybrown","seagreen","seashell",
149    "sienna","skyblue","slateblue","slategray","slategrey","snow",
150    "springgreen","steelblue","tan","thistle","tomato","turquoise",
151    "violetred","wheat","whitesmoke","yellowgreen"]
152
153    ColorColl = ot.Description(NbVA)
154    for i in range(NbVA):
155        ColorColl[i]=toutes_les_couleurs[i]
156    colors = ot.Description(NbVA)
157
158    legend = ot.Description(NbVA)
159    for i in range(NbVA):
160        if PM != None:
161            legend[i] = PM.getMarginal(i).getName()+": %.1f %%" % (sample[i][0]*100)
162        else:
163            legend[i] = "Marginal "+str(i)+": %2.0f" % (sample[i][0]*100)
164
165    SobolPie=ot.Pie(sample,legend,ot.NumericalPoint(2,0.0),1,ColorColl)
166    myGraph = ot.Graph("Total Sobol Indexes","","",False)
167    myGraph.addDrawable (SobolPie)
168    myGraph.draw("Sobol", 800, 600, ot.GraphImplementation.PNG)
169
170
171
172def print_plot_regression(Pex, res):
173    MM = res.getMetaModel()
174    inputX = ot.NumericalSample(len(Pex),MM.getInputDimension())
175    output_Y = ot.NumericalSample(len(Pex),1)
176    sample2D = ot.NumericalSample(len(Pex),2)
177
178    for i in range(len(Pex)):
179        for j in range(MM.getInputDimension()):
180            inputX[i][j] = Pex[i][j]
181    output_MM_Y = MM(inputX)
182
183    for k in range(len(Pex)):
184        output_Y[k][0] = Pex[k][MM.getInputDimension()]
185
186    for l in range(len(Pex)):
187        sample2D[l][0] = output_Y[l][0]
188        sample2D[l][1] = output_MM_Y[l][0]
189
190    mySampleDrawable = ot.Cloud(sample2D , "blue" , "square" , "Model_vs_MetaModel")
191    graph = ot.Graph("PEX: Model vs. MetaModel" , "Model output" , "MetaModel output" , True)
192    graph.addDrawable(mySampleDrawable)
193    graph.draw("Model_vs_MetaModel", 800, 600, ot.GraphImplementation.PNG)
194
195
196# R2 evaluation
197    meanYi = output_Y.computeMean()[0]
198    stdevYi = output_Y.computeStandardDeviationPerComponent()[0]
199    varYi = output_Y.computeCovariance()
200
201    SStot = 0
202    SSerr = 0
203    for i in range(len(Pex)):
204        SStot = SStot + (output_Y[i][0] - meanYi)**2
205        SSerr = SSerr + (output_Y[i][0] - output_MM_Y[i][0])**2
206    R2 = 1.0 - SSerr/SStot
207
208
209
210# Q2 evaluation
211    TI = res.getTransformation()
212    polynoms = res.getReducedBasis()
213    coeff = res.getCoefficients()
214
215    F = np.matrix(np.zeros((inputX.getSize(),coeff.getSize())))
216    for i in range(F.shape[0]):
217        for j in range(F.shape[1]):
218            F[i,j] = polynoms[j](TI(inputX[i]))[0]
219
220    y_ = np.array(np.zeros(output_Y.getSize()))
221    for i in range(len(y_)):
222        y_[i] = output_Y[i][0]
223
224    coeffs = np.array(np.zeros(coeff.getSize()))
225    for i in range(len(coeffs)):
226        coeffs[i] = coeff[i]
227
228    H = F * linalg.solve((F.T * F), F.T)
229    h = np.array(np.diag(H))
230    delta = (y_ - np.array(coeffs*F.T).ravel()) / (1 - h)
231    errLOO = np.sum(delta**2)/len(y_)
232    Q2 = 1 - errLOO / varYi[0,0]
233
234    print "\nREGRESSION:"
235    print "R^2 coefficient = %1.6f" %R2
236    print "2nd order polynomial Q^2 = %1.6f" %Q2
237