# 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 | |
---|---|

1 | import openturns as ot |

2 | from math import sqrt |

3 | import itertools as it |

4 | import numpy as np |

5 | from scipy import linalg |

6 | |

7 | |

8 | def 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 | |

14 | def reprn(x): |

15 | |

16 | if x>=0: |

17 | return "%6.2e" %x |

18 | else: |

19 | return "%7.2e" %x |

20 | |

21 | def 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 | |

31 | def 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 | |

44 | def 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 | |

76 | def 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 | |

172 | def 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 |