| 1 | #!/usr/bin/env python |
|---|
| 2 | import sys |
|---|
| 3 | import os |
|---|
| 4 | import numpy as np |
|---|
| 5 | import h5py |
|---|
| 6 | import re |
|---|
| 7 | from puq.jpickle import unpickle |
|---|
| 8 | import xml.etree.ElementTree as xml |
|---|
| 9 | |
|---|
| 10 | sys.stdout = open("analyze.out", 'w') |
|---|
| 11 | sys.stderr = open("analyze.err", 'w') |
|---|
| 12 | |
|---|
| 13 | |
|---|
| 14 | # Restore the state of a PUQ session from a HDF5 file. |
|---|
| 15 | def load_from_hdf5(name): |
|---|
| 16 | h5 = h5py.File(name, 'r+') |
|---|
| 17 | sw = unpickle(h5['private/sweep'].value) |
|---|
| 18 | sw.fname = os.path.splitext(name)[0] |
|---|
| 19 | h5.close() |
|---|
| 20 | |
|---|
| 21 | sw.psweep._sweep = sw |
|---|
| 22 | |
|---|
| 23 | if hasattr(sw.psweep, 'reinit'): |
|---|
| 24 | sw.psweep.reinit() |
|---|
| 25 | return sw |
|---|
| 26 | |
|---|
| 27 | |
|---|
| 28 | # Plots probability curves |
|---|
| 29 | def plot_pdf_curve(xvals, vname, percent, desc): |
|---|
| 30 | print 'plot_pdf_curve %s %s %s' % (vname, percent, desc) |
|---|
| 31 | # compute upper and lower percentiles |
|---|
| 32 | pm = (100 - percent)/200.0 |
|---|
| 33 | pp = 1 - pm |
|---|
| 34 | |
|---|
| 35 | # collect data into an array |
|---|
| 36 | xarr = np.empty(len(xvals[vname])) |
|---|
| 37 | yp = np.empty(len(xvals[vname])) |
|---|
| 38 | ym = np.empty(len(xvals[vname])) |
|---|
| 39 | for vindex in sorted(xvals[vname].keys()): |
|---|
| 40 | xarr[vindex] = xvals[vname][vindex] |
|---|
| 41 | yp[vindex] = pcurves[vname][vindex].ppf(pp) |
|---|
| 42 | ym[vindex] = pcurves[vname][vindex].ppf(pm) |
|---|
| 43 | |
|---|
| 44 | curve = xml.SubElement(dout, 'curve', {'id': 'curve_pdf-%s-%s' % (vname, percent)}) |
|---|
| 45 | about = xml.SubElement(curve, 'about') |
|---|
| 46 | if percent == 0: |
|---|
| 47 | xml.SubElement(about, 'label').text = "mean" |
|---|
| 48 | else: |
|---|
| 49 | xml.SubElement(about, 'label').text = "middle %s%%" % percent |
|---|
| 50 | xml.SubElement(about, 'group').text = '%s (Probability)' % desc |
|---|
| 51 | comp = xml.SubElement(curve, 'component') |
|---|
| 52 | xy = xml.SubElement(comp, 'xy') |
|---|
| 53 | pts = "" |
|---|
| 54 | for x, y in zip(xarr, yp): |
|---|
| 55 | pts += "%s %s " % (x, y) |
|---|
| 56 | if percent == 0: |
|---|
| 57 | pts += '\n' |
|---|
| 58 | else: |
|---|
| 59 | for x, y in reversed(zip(xarr, ym)): |
|---|
| 60 | pts += "%s %s " % (x, y) |
|---|
| 61 | pts += "%s %s\n" % (xarr[0], yp[0]) |
|---|
| 62 | xy.text = pts |
|---|
| 63 | |
|---|
| 64 | |
|---|
| 65 | def plot_pdf(v, pdf, desc): |
|---|
| 66 | id = '%s (PDF)' % desc |
|---|
| 67 | elem = xml.SubElement(dout, 'curve', {'id': id}) |
|---|
| 68 | about = xml.SubElement(elem, 'about') |
|---|
| 69 | xml.SubElement(about, 'label').text = id |
|---|
| 70 | yaxis = xml.SubElement(elem, 'yaxis') |
|---|
| 71 | xml.SubElement(yaxis, 'label').text = 'Probability' |
|---|
| 72 | |
|---|
| 73 | component = xml.SubElement(elem, 'component') |
|---|
| 74 | xy = xml.SubElement(component, 'xy') |
|---|
| 75 | pts = "%s 0\n" % pdf.x[0] |
|---|
| 76 | for x, y in zip(pdf.x, pdf.y): |
|---|
| 77 | pts += "%s %s\n" % (x, y) |
|---|
| 78 | pts += "%s 0\n" % pdf.x[-1] |
|---|
| 79 | xy.text = pts |
|---|
| 80 | |
|---|
| 81 | |
|---|
| 82 | sw = load_from_hdf5(sys.argv[1]) |
|---|
| 83 | sw.analyze() |
|---|
| 84 | |
|---|
| 85 | h5 = h5py.File(sys.argv[1], 'r+') |
|---|
| 86 | dtree = xml.parse('run_uq.xml') |
|---|
| 87 | droot = dtree.getroot() |
|---|
| 88 | dout = droot.find('output') |
|---|
| 89 | |
|---|
| 90 | # curves built from pdfs |
|---|
| 91 | pcurves = {} |
|---|
| 92 | xvals = {} |
|---|
| 93 | |
|---|
| 94 | regexp = re.compile('([ \da-z]+)\[([ \d]+)\]') |
|---|
| 95 | |
|---|
| 96 | uqtype = h5.attrs['UQtype'] |
|---|
| 97 | for v in h5[uqtype]: |
|---|
| 98 | rs = unpickle(h5['/%s/%s/response' % (uqtype, v)].value) |
|---|
| 99 | pdf = rs.pdf(fit=False) |
|---|
| 100 | odata = h5['/output/data/%s' % v] |
|---|
| 101 | |
|---|
| 102 | # For curves built from pdfs, just put them in a dict for now |
|---|
| 103 | if 'x' in odata.attrs: |
|---|
| 104 | matches = regexp.findall(v) |
|---|
| 105 | vname, vindex = matches[0] |
|---|
| 106 | vindex = int(vindex) |
|---|
| 107 | if vname not in pcurves: |
|---|
| 108 | pcurves[vname] = {} |
|---|
| 109 | xvals[vname] = {} |
|---|
| 110 | xvals[vname][vindex] = odata.attrs['x'] |
|---|
| 111 | pcurves[vname][vindex] = pdf |
|---|
| 112 | else: |
|---|
| 113 | desc = h5['/output/data/%s' % v].attrs['description'] |
|---|
| 114 | plot_pdf(v, pdf, desc) |
|---|
| 115 | |
|---|
| 116 | # now do pdf curves |
|---|
| 117 | for vname in xvals: |
|---|
| 118 | desc = h5['/output/data/%s[0]' % vname].attrs['description'] |
|---|
| 119 | plot_pdf_curve(xvals, vname, 0, desc) |
|---|
| 120 | plot_pdf_curve(xvals, vname, 50, desc) |
|---|
| 121 | plot_pdf_curve(xvals, vname, 95, desc) |
|---|
| 122 | |
|---|
| 123 | with open('run_uq.xml', 'w') as f: |
|---|
| 124 | f.write("<?xml version=\"1.0\"?>\n") |
|---|
| 125 | dtree.write(f) |
|---|
| 126 | |
|---|
| 127 | with open('run_uq2.xml', 'w') as f: |
|---|
| 128 | f.write("<?xml version=\"1.0\"?>\n") |
|---|
| 129 | dtree.write(f) |
|---|
| 130 | h5.close() |
|---|
