source: branches/1.7/puq/analyze.py @ 6226

Last change on this file since 6226 was 5981, checked in by gah, 7 years ago

add simsets, fix puq execution to not use use

  • Property svn:executable set to *
File size: 9.6 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 Rappture
15import StringIO
16from scipy.spatial import ConvexHull
17# geometry library
18from shapely.geometry import Polygon
19from shapely.ops import unary_union
20
21# Redirect stdout and stderr to files for debugging.
22# Append to the files created in get_params.py
23sys.stdout = open("uq_debug.out", 'a')
24sys.stderr = open("uq_debug.err", 'a')
25
26
27# Restore the state of a PUQ session from a HDF5 file.
28def load_from_hdf5(name):
29    h5 = h5py.File(name, 'r+')
30    sw = unpickle(h5['private/sweep'].value)
31    sw.fname = os.path.splitext(name)[0]
32    h5.close()
33
34    sw.psweep._sweep = sw
35
36    if hasattr(sw.psweep, 'reinit'):
37        sw.psweep.reinit()
38    return sw
39
40
41# Plots probability curves
42def plot_pdf_curve(io, h5, xvals, vname, percent):
43    print('plot_pdf_curve %s %s' % (vname, percent))
44    # compute upper and lower percentiles
45    pm = (100 - percent)/200.0
46    pp = 1 - pm
47
48    label = None
49
50    # collect data into an array
51    xarr = np.empty(len(xvals[vname]))
52    yp = np.empty(len(xvals[vname]))
53    ym = np.empty(len(xvals[vname]))
54
55    for vindex in sorted(xvals[vname].keys()):
56        if label is None:
57            label = h5['/output/data/%s[%d]' % (vname, vindex)].attrs['label']
58        xarr[vindex] = xvals[vname][vindex]
59        yp[vindex] = pcurves[vname][vindex].ppf(pp)
60        ym[vindex] = pcurves[vname][vindex].ppf(pm)
61
62    curve = io['output.curve(curve_pdf-%s-%s)' % (vname, percent)]
63    if percent == 0:
64        curve['about.label'] = "mean"
65    else:
66        curve['about.label'] = "middle %s%%" % percent
67    curve['about.group'] = label
68    curve['about.uqtype'] = 'Probability'
69
70    pts = ""
71    for x, y in zip(xarr, yp):
72        pts += "%s %s " % (x, y)
73    if percent == 0:
74        pts += '\n'
75    else:
76        for x, y in reversed(zip(xarr, ym)):
77            pts += "%s %s " % (x, y)
78        pts += "%s %s\n" % (xarr[0], yp[0])
79
80    curve['component.xy'] = pts
81
82
83def add_pts(f1, percent):
84    # compute upper and lower percentiles
85    pm = (100 - percent) / 200.0
86    pp = 1 - pm
87    prob = np.linspace(pm, pp, 31)
88    x, y = f1.eval(prob)
89    return np.array(zip(x, y))
90
91
92def plot_pdf_acurve(io, h5, acurves, vname, percent):
93    """
94    This function plots the probability curves for parametric
95    PDFs.
96    """
97    print('plot_pdf_acurve %s %s' % (vname, percent))
98
99    label = None
100    prev_pts = None  # last set of points
101
102    poly = []
103    for vindex in sorted(acurves[vname].keys()):
104        if label is None:
105            label = h5['/output/data/%s[%d]' % (vname, vindex)].attrs['label']
106        f1 = unpickle(h5['/output/data/%s[%d]' % (vname, vindex)].attrs['curve'])
107        bpts = add_pts(f1, percent)
108
109        # first data set? Just remember it.
110        if prev_pts is None:
111            prev_pts = bpts
112            continue
113
114        pts = np.array((prev_pts, bpts)).ravel().reshape(-1, 2)
115        hull = ConvexHull(pts, qhull_options='Pp')
116        p1 = Polygon([hull.points[v] for v in hull.vertices])
117        poly.append(p1)
118        prev_pts = bpts
119
120    u = unary_union(poly)
121
122    curve = io['output.curve(curve_pdf-%s-%s)' % (vname, percent)]
123    if percent == 0:
124        curve['about.label'] = "mean"
125    else:
126        curve['about.label'] = "middle %s%%" % percent
127    curve['about.group'] = label
128    curve['about.uqtype'] = 'Probability'
129    curve['component.xy'] = np.array(u.exterior.xy)
130
131
132def plot_pdf(io, v, pdf, desc):
133    print("plot_pdf %s desc=%s" % (v, desc))
134    p = io['output.curve(pdf-%s)' % v]
135    p['about.label'] = desc
136    p['about.uqtype'] = "PDF"
137    p['yaxis.label'] = 'Probability'
138
139    pts = "%s 0\n" % pdf.x[0]
140    for x, y in zip(pdf.x, pdf.y):
141        pts += "%s %s\n" % (x, y)
142    pts += "%s 0\n" % pdf.x[-1]
143    p['component.xy'] = pts
144
145
146def write_responses(io, h5):
147    uqtype = h5.attrs['UQtype']
148    for v in h5[uqtype]:
149        print("write_responses", v)
150        if '[' in v:
151            # It is a curve. Ignore.
152            continue
153        try:
154            desc = h5['%s/%s' % (uqtype, v)].attrs['description']
155        except:
156            desc = ''
157        try:
158            label = h5['%s/%s' % (uqtype, v)].attrs['label']
159        except:
160            label = ''
161
162        rsp = h5['/%s/%s/response' % (uqtype, v)].value
163        rout = io['output.response(%s)' % v]
164        rout['value'] = rsp
165        rout['about.description'] = desc
166        rout['about.label'] = label
167        rout['about.uqtype'] = 'Response'
168
169        rs = unpickle(rsp)
170        rout['variables'] = ' '.join([str(p.name) for p in rs.params])
171        labels = ' '.join([repr(str(p.label)) for p in rs.params])
172        rout['labels'] = labels.replace("'", '"')
173
174        if type(rs) == puq.response.ResponseFunc:
175            rout['equation'] = rs.eqn
176            rout['rmse'] = "{:6.3g}".format(rs.rmse()[1])
177
178        rout['data'] = rs.data
179
180
181def write_params(h5, out):
182    params = map(str, h5['/input/params'].keys())
183    print('#' * 80, file=out)
184    print('INPUT PARAMETERS', file=out)
185
186    for pname in params:
187        print('-' * 80, file=out)
188        p = puq.unpickle(h5['/input/params/' + pname].value)
189        cname = p.__class__.__name__[:-9]
190        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)
191
192        print("Name:", p.name, file=out)
193        try:
194            print("Label:", p.label, file=out)
195        except:
196            pass
197        print("Desc:", p.description, file=out)
198        print('Value:', pdf_str, file=out)
199    print('#' * 80, file=out)
200    print(file=out)
201
202
203def write_summary(io, h5):
204    outstr = StringIO.StringIO()
205    write_params(h5, outstr)
206    uqtype = h5.attrs['UQtype']
207    for v in h5[uqtype]:
208        if '[' in v:
209            # It is a curve. Ignore.
210            continue
211        desc = h5['%s/%s' % (uqtype, v)].attrs['description']
212        print("QoI: %s (%s)" % (v, desc), file=outstr)
213        rs = unpickle(h5['/%s/%s/response' % (uqtype, v)].value)
214        if type(rs) == puq.response.ResponseFunc:
215            print("\nv=%s\n" % rs.eqn, file=outstr)
216            print("SURROGATE MODEL ERROR:{:6.3g}%".format(rs.rmse()[1]), file=outstr)
217        sens = puq.unpickle(h5['/%s/%s/sensitivity' % (uqtype, v)].value)
218        max_name_len = max(map(len, [p[0] for p in sens]))
219        print("\nSENSITIVITY:", file=outstr)
220        print("Var%s     u*          dev" % (' '*(max_name_len)), file=outstr)
221        print('-'*(28+max_name_len), file=outstr)
222        for item in sens:
223            pad = ' '*(max_name_len - len(item[0]))
224            print("{}{}  {:10.4g}  {:10.4g}".format(pad, item[0],
225                item[1]['ustar'], item[1]['std']), file=outstr)
226        print('-'*(28+max_name_len), file=outstr)
227        print(file=outstr)
228    iostr = io['output.string(UQ Summary)']
229    iostr['about.label'] = 'UQ Summary'
230    iostr['current'] = outstr.getvalue()
231    outstr.close()
232
233
234def write_sensitivity(io, h5):
235    # If more than one variable, display sensitivity.
236    # Curves have indexed variables, so skip them.
237    if len(h5['/input/params']) > 1 and ['[' in x for x in h5[uqtype]].count(False):
238        for v in h5[uqtype]:
239            if '[' in v:
240                # curve. skip it.
241                continue
242            desc = h5['/output/data/%s' % v].attrs['label']
243            sens = unpickle(h5['/%s/%s/sensitivity' % (uqtype, v)].value)
244
245            hist = io['output.histogram(sens-%s)' % v]
246            hist['about.label'] = desc
247            hist['about.uqtype'] = 'Sensitivity'
248            hist['about.type'] = 'scatter'
249            hist['xaxis.label'] = 'Parameters'
250            hist['yaxis.label'] = 'Sensitivity'
251            pts = ''
252            for name in sens:
253                n = name[0]
254                try:
255                    n = h5['/input/params/%s' % n].attrs['label']
256                except:
257                    pass
258                pts += "\"%s\" %s\n" % (n, name[1]['ustar'])
259            hist['component.xy'] = pts
260
261sw = load_from_hdf5(sys.argv[1])
262sw.analyze()
263
264h5 = h5py.File(sys.argv[1], 'r+')
265io = Rappture.PyXml('run_uq.xml')
266
267# curves built from pdfs
268pcurves = {}
269xvals = {}
270acurves = {}
271
272reg1 = re.compile('([ \da-zA-Z_]+)\[([ \d]+)\]')
273
274uqtype = h5.attrs['UQtype']
275for v in h5[uqtype]:
276    print('v=', v)
277    rsp = h5['/%s/%s/response' % (uqtype, v)].value
278    rs = unpickle(rsp)
279    pdf = rs.pdf(fit=False)
280    odata = h5['/output/data/%s' % v]
281
282    # For curves built from pdfs, just put them in a dict for now
283    if 'x' in odata.attrs:
284        matches = reg1.findall(v)
285        vname, vindex = matches[0]
286        print('CURVE: vname=%s   vindex=%s' % (vname, vindex))
287        vindex = int(vindex)
288        if vname not in pcurves:
289            pcurves[vname] = {}
290            xvals[vname] = {}
291        xvals[vname][vindex] = odata.attrs['x']
292        pcurves[vname][vindex] = pdf
293    elif 'curve' in odata.attrs:
294        matches = reg1.findall(v)
295        vname, vindex = matches[0]
296        print('ACURVE: %s - %s' % (vname, vindex))
297        if vname not in acurves:
298            acurves[vname] = {}
299        acurves[vname][int(vindex)] = pdf
300    else:
301        desc = h5['/output/data/%s' % v].attrs['label']
302        plot_pdf(io, v, pdf, desc)
303
304# now do probability curves
305for vname in xvals:
306    plot_pdf_curve(io, h5, xvals, vname, 95)
307    plot_pdf_curve(io, h5, xvals, vname, 50)
308
309for vname in acurves:
310    try:
311        plot_pdf_acurve(io, h5, acurves, vname, 95)
312    except:
313        pass
314    try:
315        plot_pdf_acurve(io, h5, acurves, vname, 50)
316    except:
317        pass
318
319write_sensitivity(io, h5)
320write_responses(io, h5)
321write_summary(io, h5)
322
323io.close()
324h5.close()
Note: See TracBrowser for help on using the repository browser.