Changeset 5829

Timestamp:

Aug 21, 2015, 12:21:27 PM (9 years ago)

Author:

mmh

Message:

improved probability curve plotting, bug fixes, cleanup

File:

• branches/uq/puq/puq_analyze.py (modified) (9 diffs)

branches/uq/puq/puq_analyze.py

@@ -12 +12 @@
 import puq
 from puq.jpickle import unpickle
-import xml.etree.ElementTree as xml
 import Rappture
 import StringIO
+from scipy.spatial import ConvexHull
+# geometry library
+from shapely.geometry import Polygon
+from shapely.ops import unary_union
 
 # Redirect stdout and stderr to files for debugging.
@@ -37 +40 @@
 
 # Plots probability curves
-def plot_pdf_curve(xvals, vname, percent):
+def plot_pdf_curve(io, h5, xvals, vname, percent):
     print('plot_pdf_curve %s %s' % (vname, percent))
     # compute upper and lower percentiles
@@ -56 +59 @@
         ym[vindex] = pcurves[vname][vindex].ppf(pm)
 
-    curve = xml.SubElement(dout, 'curve', {'id': 'curve_pdf-%s-%s' % (vname, percent)})
-    about = xml.SubElement(curve, 'about')
+    curve = io['output.curve(curve_pdf-%s-%s)' % (vname, percent)]
     if percent == 0:
-        xml.SubElement(about, 'label').text = "mean"
+        curve['about.label'] = "mean"
     else:
-        xml.SubElement(about, 'label').text = "middle %s%%" % percent
-    xml.SubElement(about, 'group').text = '%s:UQ:6' % label
-    comp = xml.SubElement(curve, 'component')
-    xy = xml.SubElement(comp, 'xy')
+        curve['about.label'] = "middle %s%%" % percent
+    curve['about.group'] = label
+    curve['about.uqtype'] = 'Probability'
+
     pts = ""
     for x, y in zip(xarr, yp):
@@ -74 +76 @@
             pts += "%s %s " % (x, y)
         pts += "%s %s\n" % (xarr[0], yp[0])
-    xy.text = pts
-
-
-def plot_pdf_acurve(hf, acurves, vname, percent):
+
+    curve['component.xy'] = pts
+
+
+def add_pts(f1, percent):
+    # compute upper and lower percentiles
+    pm = (100 - percent) / 200.0
+    pp = 1 - pm
+    prob = np.linspace(pm, pp, 31)
+    x, y = f1.eval(prob)
+    return np.array(zip(x, y))
+
+
+def plot_pdf_acurve(io, h5, acurves, vname, percent):
     """
     This function plots the probability curves for parametric
@@ -83 +95 @@
     """
     print('plot_pdf_acurve %s %s' % (vname, percent))
-    # compute upper and lower percentiles
-    pm = (100 - percent)/200.0
-    pp = 1 - pm
 
     label = None
-    maxvals = []
-    minvals = []
+    prev_pts = None  # last set of points
+
+    poly = []
     for vindex in sorted(acurves[vname].keys()):
         if label is None:
             label = h5['/output/data/%s[%d]' % (vname, vindex)].attrs['label']
-        f1 = unpickle(hf['/output/data/%s[%d]' % (vname, vindex)].attrs['curve'])
-        pdf = acurves[vname][vindex]
-        a, b = f1.eval([pdf.ppf(pm), pdf.ppf(pp)])
-        maxvals.append((a[1], b[1]))
-        minvals.append((a[0], b[0]))
-
-    curve = xml.SubElement(dout, 'curve', {'id': 'curve_pdf-%s-%s' % (vname, percent)})
-    about = xml.SubElement(curve, 'about')
+        f1 = unpickle(h5['/output/data/%s[%d]' % (vname, vindex)].attrs['curve'])
+        bpts = add_pts(f1, percent)
+
+        # first data set? Just remember it.
+        if prev_pts is None:
+            prev_pts = bpts
+            continue
+
+        pts = np.array((prev_pts, bpts)).ravel().reshape(-1, 2)
+        hull = ConvexHull(pts, qhull_options='Pp')
+        p1 = Polygon([hull.points[v] for v in hull.vertices])
+        poly.append(p1)
+        prev_pts = bpts
+
+    u = unary_union(poly)
+
+    curve = io['output.curve(curve_pdf-%s-%s)' % (vname, percent)]
     if percent == 0:
-        xml.SubElement(about, 'label').text = "mean"
+        curve['about.label'] = "mean"
     else:
-        xml.SubElement(about, 'label').text = "middle %s%%" % percent
-    xml.SubElement(about, 'group').text = '%s:UQ:6' % label
-    comp = xml.SubElement(curve, 'component')
-    xy = xml.SubElement(comp, 'xy')
-
-    pts = ''
-    for v in maxvals:
-        pts += "%s %s " % (v[0], v[1])
-
-    if percent == 0:
-        pts += '\n'
-    else:
-        for v in reversed(minvals):
-            pts += "%s %s " % (v[0], v[1])
-    xy.text = pts
-
-
-def plot_pdf(v, pdf, desc):
-    id = '%s:UQ:1' % desc
-    elem = xml.SubElement(dout, 'curve', {'id': id})
-    about = xml.SubElement(elem, 'about')
-    xml.SubElement(about, 'label').text = id
-    yaxis = xml.SubElement(elem, 'yaxis')
-    xml.SubElement(yaxis, 'label').text = 'Probability'
-
-    component = xml.SubElement(elem, 'component')
-    xy = xml.SubElement(component, 'xy')
+        curve['about.label'] = "middle %s%%" % percent
+    curve['about.group'] = label
+    curve['about.uqtype'] = 'Probability'
+    curve['component.xy'] = np.array(u.exterior.xy)
+
+
+def plot_pdf(io, v, pdf, desc):
+    p = io['output.curve(%s)' % desc]
+    p['about.label'] = desc
+    p['about.uqtype'] = "PDF"
+    p['yaxis.label'] = 'Probability'
+
     pts = "%s 0\n" % pdf.x[0]
     for x, y in zip(pdf.x, pdf.y):
         pts += "%s %s\n" % (x, y)
     pts += "%s 0\n" % pdf.x[-1]
-    xy.text = pts
-
-# Plotting response distance
-# plot_resp1 distance
-# Plotting response maxheight
-# plot_resp1 maxheight
-# desc= The distance the projectile travels.
-# label= Distance
-
-# v=-0.610951707827365*angle**2 + 54.985653704463*angle - 217.610577363519
-
-# SURROGATE MODEL ERROR: 0.442%
-# desc= The maximum height the projectile reaches.
-# label= Maximum Height
-
-# v=3.42068544561648e-12*angle**2 + 8.66308066792632*angle - 134.909576819562
-
-# SURROGATE MODEL ERROR: 0.586%
-# [[  45.          254.92905324]
-#  [  30.19262971  128.95166845]
-#  [  59.80737029  380.90643804]
-#  [  34.52960806  163.81660645]
-#  [  55.47039194  346.04150004]]
+    p['component.xy'] = pts
 
 
@@ -180 +162 @@
         rout['value'] = rsp
         rout['about.description'] = desc
-        rout['about.label'] = label + ":UQ:5"
+        rout['about.label'] = label
+        rout['about.uqtype'] = 'Response'
 
         rs = unpickle(rsp)
@@ -247 +230 @@
 
 
+def write_sensitivity(io, h5):
+    # If more than one variable, display sensitivity.
+    # Curves have indexed variables, so skip them.
+    if len(h5['/input/params']) > 1 and ['[' in x for x in h5[uqtype]].count(False):
+        for v in h5[uqtype]:
+            if '[' in v:
+                # curve. skip it.
+                continue
+            desc = h5['/output/data/%s' % v].attrs['label']
+            sens = unpickle(h5['/%s/%s/sensitivity' % (uqtype, v)].value)
+
+            hist = io['output.histogram(sens-%s)' % v]
+            hist['about.label'] = desc
+            hist['about.uqtype'] = 'Sensitivity'
+            hist['about.type'] = 'scatter'
+            hist['xaxis.label'] = 'Parameters'
+            hist['yaxis.label'] = 'Sensitivity'
+            pts = ''
+            for name in sens:
+                n = name[0]
+                try:
+                    n = h5['/input/params/%s' % n].attrs['label']
+                except:
+                    pass
+                pts += "\"%s\" %s\n" % (n, name[1]['ustar'])
+            hist['component.xy'] = pts
+
 sw = load_from_hdf5(sys.argv[1])
 sw.analyze()
 
 h5 = h5py.File(sys.argv[1], 'r+')
-dtree = xml.parse('run_uq.xml')
-droot = dtree.getroot()
-dout = droot.find('output')
+io = Rappture.PyXml('run_uq.xml')
 
 # curves built from pdfs
@@ -260 +268 @@
 acurves = {}
 
-reg1 = re.compile('([ \da-zA-Z]+)\[([ \d]+)\]')
+reg1 = re.compile('([ \da-zA-Z_]+)\[([ \d]+)\]')
 
 uqtype = h5.attrs['UQtype']
@@ -288 +296 @@
     else:
         desc = h5['/output/data/%s' % v].attrs['label']
-        plot_pdf(v, pdf, desc)
+        plot_pdf(io, v, pdf, desc)
 
 # now do probability curves
 for vname in xvals:
-    plot_pdf_curve(xvals, vname, 0)
-    plot_pdf_curve(xvals, vname, 50)
-    plot_pdf_curve(xvals, vname, 95)
+    plot_pdf_curve(io, h5, xvals, vname, 95)
+    plot_pdf_curve(io, h5, xvals, vname, 50)
 
 for vname in acurves:
-    plot_pdf_acurve(h5, acurves, vname, 0)
-    plot_pdf_acurve(h5, acurves, vname, 50)
-    plot_pdf_acurve(h5, acurves, vname, 95)
-
-
-# If more than one variable, display sensitivity.
-# Curves have indexed variables, so skip them.
-if len(h5['/input/params']) > 1 and ['[' in x for x in h5[uqtype]].count(False):
-    for v in h5[uqtype]:
-        if '[' in v:
-            continue
-        desc = h5['/output/data/%s' % v].attrs['label']
-        sens = unpickle(h5['/%s/%s/sensitivity' % (uqtype, v)].value)
-        elem = xml.SubElement(dout, 'histogram', {'id': 'sens-%s' % v})
-        about = xml.SubElement(elem, 'about')
-        xml.SubElement(about, 'label').text = '%s:UQ:2' % desc
-        xml.SubElement(elem, 'type').text = 'scatter'
-        xaxis = xml.SubElement(elem, 'xaxis')
-        xml.SubElement(xaxis, 'label').text = 'Parameters'
-        yaxis = xml.SubElement(elem, 'yaxis')
-        xml.SubElement(yaxis, 'label').text = 'Sensitivity'
-        comp = xml.SubElement(elem, 'component')
-        xy = xml.SubElement(comp, 'xy')
-        pts = ''
-        for name in sens:
-            n = name[0]
-            try:
-                n = h5['/input/params/%s' % n].attrs['label']
-            except:
-                pass
-            pts += "\"%s\" %s\n" % (n, name[1]['ustar'])
-        xy.text = pts
-
-with open('run_uq.xml', 'w') as f:
-    f.write("<?xml version=\"1.0\"?>\n")
-    dtree.write(f)
-
-io = Rappture.PyXml('run_uq.xml')
+    try:
+        plot_pdf_acurve(io, h5, acurves, vname, 95)
+    except:
+        pass
+    try:
+        plot_pdf_acurve(io, h5, acurves, vname, 50)
+    except:
+        pass
+
+write_sensitivity(io, h5)
 write_responses(io, h5)
 write_summary(io, h5)