Changeset 5746

Timestamp:

Jul 14, 2015 10:11:27 PM (9 years ago)

Author:

mmh

Message:

support parametric probability curves

File:

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

branches/uq/puq/puq_analyze.py

@@ -37 +37 @@
 
 # Plots probability curves
-def plot_pdf_curve(xvals, vname, percent, desc):
-    print('plot_pdf_curve %s %s %s' % (vname, percent, desc))
+def plot_pdf_curve(xvals, vname, percent):
+    print('plot_pdf_curve %s %s' % (vname, percent))
     # compute upper and lower percentiles
     pm = (100 - percent)/200.0
     pp = 1 - pm
+
+    label = None
 
     # collect data into an array
@@ -48 +50 @@
     ym = np.empty(len(xvals[vname]))
     for vindex in sorted(xvals[vname].keys()):
+        if label is None:
+            label = h5['/output/data/%s[%d]' % (vname, vindex)].attrs['label']
         xarr[vindex] = xvals[vname][vindex]
         yp[vindex] = pcurves[vname][vindex].ppf(pp)
@@ -58 +62 @@
     else:
         xml.SubElement(about, 'label').text = "middle %s%%" % percent
-    xml.SubElement(about, 'group').text = '%s:UQ:6' % desc
+    xml.SubElement(about, 'group').text = '%s:UQ:6' % label
     comp = xml.SubElement(curve, 'component')
     xy = xml.SubElement(comp, 'xy')
@@ -73 +77 @@
 
 
-def dist(x1, y1, x2, y2, x3, y3): # x3,y3 is the point
-    diffx = x2-x1
-    diffy = y2-y1
-    sqdiff = float(diffx**2 + diffy**2)
-
-    u = ((x3 - x1) * diffx + (y3 - y1) * diffy) / sqdiff
-
-    if u > 1:
-        u = 1
-    elif u < 0:
-        u = 0
-
-    x = x1 + u * diffx
-    y = y1 + u * diffy
-
-    dx = x - x3
-    dy = y - y3
-
-    return np.sqrt(dx*dx + dy*dy)
-
-
-# Plots advanced probability curves
-def plot_pdf_acurve(hf, acurves, vname, percent, desc):
-    print('plot_pdf_acurve %s %s %s' % (vname, percent, desc))
+def plot_pdf_acurve(hf, acurves, vname, percent):
+    """
+    This function plots the probability curves for parametric
+    PDFs.
+    """
+    print('plot_pdf_acurve %s %s' % (vname, percent))
     # compute upper and lower percentiles
     pm = (100 - percent)/200.0
     pp = 1 - pm
-    vlen = len(acurves[vname])
-
-    print('vlen=', vlen)
-    # collect data into an array
-    xp = np.empty(vlen)
-    xm = np.empty(vlen)
-    yp = np.empty(vlen)
-    ym = np.empty(vlen)
-
+
+    label = None
+    maxvals = []
+    minvals = []
     for vindex in sorted(acurves[vname].keys()):
-        print('vindex=', vindex)
-        x1 = acurves[vname][vindex][0].ppf(pp)
-        x2 = acurves[vname][vindex][0].ppf(pm)
-        y1 = acurves[vname][vindex][1].ppf(pp)
-        y2 = acurves[vname][vindex][1].ppf(pm)
-        xd = hf['/output/data/%s.x[%d]' % (vname, int(vindex))].value
-        yd = hf['/output/data/%s.y[%d]' % (vname, int(vindex))].value
-        #p1, p2 = get_closest2()
-        if int(vindex) == 2:
-            print (x1, x2, y1, y2)
-            print ('xd=', repr(xd))
-            print ('yd=', repr(yd))
-
-        #xp[vindex] = x1
-        #xm[vindex] = x2
-        #yp[vindex] = y1
-        #ym[vindex] = y2
-    return
+        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')
@@ -134 +105 @@
     else:
         xml.SubElement(about, 'label').text = "middle %s%%" % percent
-    xml.SubElement(about, 'group').text = '%s:UQ:6' % desc
+    xml.SubElement(about, 'group').text = '%s:UQ:6' % label
     comp = xml.SubElement(curve, 'component')
     xy = xml.SubElement(comp, 'xy')
-    pts = ""
-    for x, y in zip(xarr, yp):
-        pts += "%s %s " % (x, y)
+
+    pts = ''
+    for v in maxvals:
+        pts += "%s %s " % (v[0], v[1])
+
     if percent == 0:
         pts += '\n'
     else:
-        for x, y in reversed(zip(xarr, ym)):
-            pts += "%s %s " % (x, y)
-        pts += "%s %s\n" % (xarr[0], yp[0])
+        for v in reversed(minvals):
+            pts += "%s %s " % (v[0], v[1])
     xy.text = pts
 
@@ -274 +246 @@
     outstr.close()
 
+
 sw = load_from_hdf5(sys.argv[1])
 sw.analyze()
@@ -288 +261 @@
 
 reg1 = re.compile('([ \da-zA-Z]+)\[([ \d]+)\]')
-reg2 = re.compile('([ \da-zA-Z]+)\.([xy])\[([ \d]+)\]')
 
 uqtype = h5.attrs['UQtype']
 for v in h5[uqtype]:
     rsp = h5['/%s/%s/response' % (uqtype, v)].value
-    print("rsp=", rsp)
     rs = unpickle(rsp)
     pdf = rs.pdf(fit=False)
@@ -308 +279 @@
         xvals[vname][vindex] = odata.attrs['x']
         pcurves[vname][vindex] = pdf
-    elif reg2.findall(v) != []:
-        matches = reg2.findall(v)
-        vname, xy, vindex = matches[0]
+    elif 'curve' in odata.attrs:
+        matches = reg1.findall(v)
+        vname, vindex = matches[0]
+        print('ACURVE: %s - %s' % (vname, vindex))
         if vname not in acurves:
             acurves[vname] = {}
-        if vindex not in acurves[vname]:
-            acurves[vname][vindex] = {}
-        if xy == 'x':
-            acurves[vname][vindex][0] = pdf
-        else:
-            acurves[vname][vindex][1] = pdf
+        acurves[vname][int(vindex)] = pdf
     else:
         desc = h5['/output/data/%s' % v].attrs['label']
         plot_pdf(v, pdf, desc)
 
-# now do pdf curves
+# now do probability curves
 for vname in xvals:
-    desc = h5['/output/data/%s[0]' % vname].attrs['description']
-    plot_pdf_curve(xvals, vname, 0, desc)
-    plot_pdf_curve(xvals, vname, 50, desc)
-    plot_pdf_curve(xvals, vname, 95, desc)
-
-# for vname in acurves:
-#    desc = h5['/output/data/%s.x[0]' % vname].attrs['description']
-#    plot_pdf_acurve(h5, acurves, vname, 0, desc)
-#    plot_pdf_acurve(h5, acurves, vname, 50, desc)
-#    plot_pdf_acurve(h5, acurves, vname, 95, desc)
+    plot_pdf_curve(xvals, vname, 0)
+    plot_pdf_curve(xvals, vname, 50)
+    plot_pdf_curve(xvals, vname, 95)
+
+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)