source: branches/uq/puq/puq_analyze.py @ 5706

Last change on this file since 5706 was 5706, checked in by mmh, 9 years ago

remove response surface code
  • Property svn:executable set to *
File size: 11.5 KB
Line 
1#!/usr/bin/env python
2"""
3Once submit has finished with the jobs, this function is called to have PUQ
4process the results.
5"""
6from __future__ import print_function
7import sys
8import os
9import numpy as np
10import h5py
11import re
12import puq
13from puq.jpickle import unpickle
14import xml.etree.ElementTree as xml
15import Rappture
16import StringIO
17
18# Redirect stdout and stderr to files for debugging.
19# Append to the files created in get_params.py
20sys.stdout = open("uq_debug.out", 'a')
21sys.stderr = open("uq_debug.err", 'a')
22
23
24# Restore the state of a PUQ session from a HDF5 file.
25def load_from_hdf5(name):
26    h5 = h5py.File(name, 'r+')
27    sw = unpickle(h5['private/sweep'].value)
28    sw.fname = os.path.splitext(name)[0]
29    h5.close()
30
31    sw.psweep._sweep = sw
32
33    if hasattr(sw.psweep, 'reinit'):
34        sw.psweep.reinit()
35    return sw
36
37
38# Plots probability curves
39def plot_pdf_curve(xvals, vname, percent, desc):
40    print('plot_pdf_curve %s %s %s' % (vname, percent, desc))
41    # compute upper and lower percentiles
42    pm = (100 - percent)/200.0
43    pp = 1 - pm
44
45    # collect data into an array
46    xarr = np.empty(len(xvals[vname]))
47    yp = np.empty(len(xvals[vname]))
48    ym = np.empty(len(xvals[vname]))
49    for vindex in sorted(xvals[vname].keys()):
50        xarr[vindex] = xvals[vname][vindex]
51        yp[vindex] = pcurves[vname][vindex].ppf(pp)
52        ym[vindex] = pcurves[vname][vindex].ppf(pm)
53
54    curve = xml.SubElement(dout, 'curve', {'id': 'curve_pdf-%s-%s' % (vname, percent)})
55    about = xml.SubElement(curve, 'about')
56    if percent == 0:
57        xml.SubElement(about, 'label').text = "mean"
58    else:
59        xml.SubElement(about, 'label').text = "middle %s%%" % percent
60    xml.SubElement(about, 'group').text = '%s (Probability)' % desc
61    comp = xml.SubElement(curve, 'component')
62    xy = xml.SubElement(comp, 'xy')
63    pts = ""
64    for x, y in zip(xarr, yp):
65        pts += "%s %s " % (x, y)
66    if percent == 0:
67        pts += '\n'
68    else:
69        for x, y in reversed(zip(xarr, ym)):
70            pts += "%s %s " % (x, y)
71        pts += "%s %s\n" % (xarr[0], yp[0])
72    xy.text = pts
73
74
75def dist(x1, y1, x2, y2, x3, y3): # x3,y3 is the point
76    diffx = x2-x1
77    diffy = y2-y1
78    sqdiff = float(diffx**2 + diffy**2)
79
80    u = ((x3 - x1) * diffx + (y3 - y1) * diffy) / sqdiff
81
82    if u > 1:
83        u = 1
84    elif u < 0:
85        u = 0
86
87    x = x1 + u * diffx
88    y = y1 + u * diffy
89
90    dx = x - x3
91    dy = y - y3
92
93    return np.sqrt(dx*dx + dy*dy)
94
95
96# Plots advanced probability curves
97def plot_pdf_acurve(hf, acurves, vname, percent, desc):
98    print('plot_pdf_acurve %s %s %s' % (vname, percent, desc))
99    # compute upper and lower percentiles
100    pm = (100 - percent)/200.0
101    pp = 1 - pm
102    vlen = len(acurves[vname])
103
104    print('vlen=', vlen)
105    # collect data into an array
106    xp = np.empty(vlen)
107    xm = np.empty(vlen)
108    yp = np.empty(vlen)
109    ym = np.empty(vlen)
110
111    for vindex in sorted(acurves[vname].keys()):
112        print('vindex=', vindex)
113        x1 = acurves[vname][vindex][0].ppf(pp)
114        x2 = acurves[vname][vindex][0].ppf(pm)
115        y1 = acurves[vname][vindex][1].ppf(pp)
116        y2 = acurves[vname][vindex][1].ppf(pm)
117        xd = hf['/output/data/%s.x[%d]' % (vname, int(vindex))].value
118        yd = hf['/output/data/%s.y[%d]' % (vname, int(vindex))].value
119        #p1, p2 = get_closest2()
120        if int(vindex) == 2:
121            print (x1, x2, y1, y2)
122            print ('xd=', repr(xd))
123            print ('yd=', repr(yd))
124
125        #xp[vindex] = x1
126        #xm[vindex] = x2
127        #yp[vindex] = y1
128        #ym[vindex] = y2
129    return
130    curve = xml.SubElement(dout, 'curve', {'id': 'curve_pdf-%s-%s' % (vname, percent)})
131    about = xml.SubElement(curve, 'about')
132    if percent == 0:
133        xml.SubElement(about, 'label').text = "mean"
134    else:
135        xml.SubElement(about, 'label').text = "middle %s%%" % percent
136    xml.SubElement(about, 'group').text = '%s (Probability)' % desc
137    comp = xml.SubElement(curve, 'component')
138    xy = xml.SubElement(comp, 'xy')
139    pts = ""
140    for x, y in zip(xarr, yp):
141        pts += "%s %s " % (x, y)
142    if percent == 0:
143        pts += '\n'
144    else:
145        for x, y in reversed(zip(xarr, ym)):
146            pts += "%s %s " % (x, y)
147        pts += "%s %s\n" % (xarr[0], yp[0])
148    xy.text = pts
149
150
151def plot_pdf(v, pdf, desc):
152    id = '%s (PDF)' % desc
153    elem = xml.SubElement(dout, 'curve', {'id': id})
154    about = xml.SubElement(elem, 'about')
155    xml.SubElement(about, 'label').text = id
156    yaxis = xml.SubElement(elem, 'yaxis')
157    xml.SubElement(yaxis, 'label').text = 'Probability'
158
159    component = xml.SubElement(elem, 'component')
160    xy = xml.SubElement(component, 'xy')
161    pts = "%s 0\n" % pdf.x[0]
162    for x, y in zip(pdf.x, pdf.y):
163        pts += "%s %s\n" % (x, y)
164    pts += "%s 0\n" % pdf.x[-1]
165    xy.text = pts
166
167# Plotting response distance
168# plot_resp1 distance
169# Plotting response maxheight
170# plot_resp1 maxheight
171# desc= The distance the projectile travels.
172# label= Distance
173
174# v=-0.610951707827365*angle**2 + 54.985653704463*angle - 217.610577363519
175
176# SURROGATE MODEL ERROR: 0.442%
177# desc= The maximum height the projectile reaches.
178# label= Maximum Height
179
180# v=3.42068544561648e-12*angle**2 + 8.66308066792632*angle - 134.909576819562
181
182# SURROGATE MODEL ERROR: 0.586%
183# [[  45.          254.92905324]
184#  [  30.19262971  128.95166845]
185#  [  59.80737029  380.90643804]
186#  [  34.52960806  163.81660645]
187#  [  55.47039194  346.04150004]]
188
189
190def write_responses(io, h5):
191    uqtype = h5.attrs['UQtype']
192    for v in h5[uqtype]:
193        if '[' in v:
194            # It is a curve. Ignore.
195            continue
196        try:
197            desc = h5['%s/%s' % (uqtype, v)].attrs['description']
198        except:
199            desc = ''
200        try:
201            label = h5['%s/%s' % (uqtype, v)].attrs['label']
202        except:
203            label = ''
204
205        rsp = h5['/%s/%s/response' % (uqtype, v)].value
206
207        rout = io['output.response(%s)' % label]
208        rout['value'] = rsp
209        rout['about.description'] = desc
210        rout['about.label'] = label + " (Response)"
211
212        rs = unpickle(rsp)
213        rout['variables'] = ' '.join([str(p.name) for p in rs.params])
214        labels = ' '.join([repr(str(p.label)) for p in rs.params])
215        rout['labels'] = labels.replace("'", '"')
216
217        if type(rs) == puq.response.ResponseFunc:
218            rout['equation'] = rs.eqn
219            rout['rmse'] = "{:6.3g}".format(rs.rmse()[1])
220
221        rout['data'] = rs.data
222
223
224def write_params(h5, out):
225    params = map(str, h5['/input/params'].keys())
226    print('#' * 80, file=out)
227    print('INPUT PARAMETERS', file=out)
228
229    for pname in params:
230        print('-' * 80, file=out)
231        p = puq.unpickle(h5['/input/params/' + pname].value)
232        cname = p.__class__.__name__[:-9]
233        pdf_str = '%s [%s - %s] mean=%s dev=%s mode=%s' % (cname, p.pdf.range[0], p.pdf.range[1], p.pdf.mean, p.pdf.dev, p.pdf.mode)
234
235        print("Name:", p.name, file=out)
236        try:
237            print("Label:", p.label, file=out)
238        except:
239            pass
240        print("Desc:", p.description, file=out)
241        print('Value:', pdf_str, file=out)
242    print('#' * 80, file=out)
243    print(file=out)
244
245
246def write_summary(io, h5):
247    outstr = StringIO.StringIO()
248    write_params(h5, outstr)
249    uqtype = h5.attrs['UQtype']
250    for v in h5[uqtype]:
251        if '[' in v:
252            # It is a curve. Ignore.
253            continue
254        desc = h5['%s/%s' % (uqtype, v)].attrs['description']
255        print("QoI: %s (%s)" % (v, desc), file=outstr)
256        rs = unpickle(h5['/%s/%s/response' % (uqtype, v)].value)
257        if type(rs) == puq.response.ResponseFunc:
258            print("\nv=%s\n" % rs.eqn, file=outstr)
259            print("SURROGATE MODEL ERROR:{:6.3g}%".format(rs.rmse()[1]), file=outstr)
260        sens = puq.unpickle(h5['/%s/%s/sensitivity' % (uqtype, v)].value)
261        max_name_len = max(map(len, [p[0] for p in sens]))
262        print("\nSENSITIVITY:", file=outstr)
263        print("Var%s     u*          dev" % (' '*(max_name_len)), file=outstr)
264        print('-'*(28+max_name_len), file=outstr)
265        for item in sens:
266            pad = ' '*(max_name_len - len(item[0]))
267            print("{}{}  {:10.4g}  {:10.4g}".format(pad, item[0],
268                item[1]['ustar'], item[1]['std']), file=outstr)
269        print('-'*(28+max_name_len), file=outstr)
270        print(file=outstr)
271    iostr = io['output.string(UQ Summary)']
272    iostr['about.label'] = 'UQ Summary'
273    iostr['current'] = outstr.getvalue()
274    outstr.close()
275
276sw = load_from_hdf5(sys.argv[1])
277sw.analyze()
278
279h5 = h5py.File(sys.argv[1], 'r+')
280dtree = xml.parse('run_uq.xml')
281droot = dtree.getroot()
282dout = droot.find('output')
283
284# curves built from pdfs
285pcurves = {}
286xvals = {}
287acurves = {}
288
289reg1 = re.compile('([ \da-zA-Z]+)\[([ \d]+)\]')
290reg2 = re.compile('([ \da-zA-Z]+)\.([xy])\[([ \d]+)\]')
291
292uqtype = h5.attrs['UQtype']
293for v in h5[uqtype]:
294    rsp = h5['/%s/%s/response' % (uqtype, v)].value
295    print("rsp=", rsp)
296    rs = unpickle(rsp)
297    pdf = rs.pdf(fit=False)
298    odata = h5['/output/data/%s' % v]
299
300    # For curves built from pdfs, just put them in a dict for now
301    if 'x' in odata.attrs:
302        matches = reg1.findall(v)
303        vname, vindex = matches[0]
304        vindex = int(vindex)
305        if vname not in pcurves:
306            pcurves[vname] = {}
307            xvals[vname] = {}
308        xvals[vname][vindex] = odata.attrs['x']
309        pcurves[vname][vindex] = pdf
310    elif reg2.findall(v) != []:
311        matches = reg2.findall(v)
312        vname, xy, vindex = matches[0]
313        if vname not in acurves:
314            acurves[vname] = {}
315        if vindex not in acurves[vname]:
316            acurves[vname][vindex] = {}
317        if xy == 'x':
318            acurves[vname][vindex][0] = pdf
319        else:
320            acurves[vname][vindex][1] = pdf
321    else:
322        desc = h5['/output/data/%s' % v].attrs['label']
323        plot_pdf(v, pdf, desc)
324
325# now do pdf curves
326for vname in xvals:
327    desc = h5['/output/data/%s[0]' % vname].attrs['description']
328    plot_pdf_curve(xvals, vname, 0, desc)
329    plot_pdf_curve(xvals, vname, 50, desc)
330    plot_pdf_curve(xvals, vname, 95, desc)
331
332# for vname in acurves:
333#    desc = h5['/output/data/%s.x[0]' % vname].attrs['description']
334#    plot_pdf_acurve(h5, acurves, vname, 0, desc)
335#    plot_pdf_acurve(h5, acurves, vname, 50, desc)
336#    plot_pdf_acurve(h5, acurves, vname, 95, desc)
337
338
339# If more than one variable, display sensitivity.
340# Curves have indexed variables, so skip them.
341if len(h5['/input/params']) > 1 and ['[' in x for x in h5[uqtype]].count(False):
342    for v in h5[uqtype]:
343        if '[' in v:
344            continue
345        desc = h5['/output/data/%s' % v].attrs['label']
346        sens = unpickle(h5['/%s/%s/sensitivity' % (uqtype, v)].value)
347        elem = xml.SubElement(dout, 'histogram', {'id': 'sens-%s' % v})
348        about = xml.SubElement(elem, 'about')
349        xml.SubElement(about, 'label').text = '%s (Sensitivity)' % desc
350        xml.SubElement(elem, 'type').text = 'scatter'
351        xaxis = xml.SubElement(elem, 'xaxis')
352        xml.SubElement(xaxis, 'label').text = 'Parameters'
353        yaxis = xml.SubElement(elem, 'yaxis')
354        xml.SubElement(yaxis, 'label').text = 'Sensitivity'
355        comp = xml.SubElement(elem, 'component')
356        xy = xml.SubElement(comp, 'xy')
357        pts = ''
358        for name in sens:
359            n = name[0]
360            try:
361                n = h5['/input/params/%s' % n].attrs['label']
362            except:
363                pass
364            pts += "\"%s\" %s\n" % (n, name[1]['ustar'])
365        xy.text = pts
366
367with open('run_uq.xml', 'w') as f:
368    f.write("<?xml version=\"1.0\"?>\n")
369    dtree.write(f)
370
371io = Rappture.PyXml('run_uq.xml')
372write_responses(io, h5)
373write_summary(io, h5)
374io.close()
375
376h5.close()
Note: See TracBrowser for help on using the repository browser.