Changeset 5166 for branches/uq/puq/puq_analyze.py

Timestamp:

Mar 23, 2015, 10:48:24 AM (9 years ago)

Author:

mmh

Message:

snapshot of uq work

File:

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

branches/uq/puq/puq_analyze.py

@@ -62 +62 @@
     xy.text = pts
 
+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)
+    # 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)
+
+    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
+    curve = xml.SubElement(dout, 'curve', {'id': 'curve_pdf-%s-%s' % (vname, percent)})
+    about = xml.SubElement(curve, 'about')
+    if percent == 0:
+        xml.SubElement(about, 'label').text = "mean"
+    else:
+        xml.SubElement(about, 'label').text = "middle %s%%" % percent
+    xml.SubElement(about, 'group').text = '%s (Probability)' % desc
+    comp = xml.SubElement(curve, 'component')
+    xy = xml.SubElement(comp, 'xy')
+    pts = ""
+    for x, y in zip(xarr, yp):
+        pts += "%s %s " % (x, y)
+    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])
+    xy.text = pts
+
 
 def plot_pdf(v, pdf, desc):
@@ -91 +165 @@
 pcurves = {}
 xvals = {}
-
-regexp = re.compile('([ \da-z]+)\[([ \d]+)\]')
+acurves = {}
+
+reg1 = re.compile('([ \da-zA-Z]+)\[([ \d]+)\]')
+reg2 = re.compile('([ \da-zA-Z]+)\.([xy])\[([ \d]+)\]')
 
 uqtype = h5.attrs['UQtype']
@@ -102 +178 @@
     # For curves built from pdfs, just put them in a dict for now
     if 'x' in odata.attrs:
-        matches = regexp.findall(v)
+        matches = reg1.findall(v)
         vname, vindex = matches[0]
         vindex = int(vindex)
@@ -110 +186 @@
         xvals[vname][vindex] = odata.attrs['x']
         pcurves[vname][vindex] = pdf
+    elif reg2.findall(v) != []:
+        matches = reg2.findall(v)
+        vname, xy, vindex = matches[0]
+        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
     else:
         desc = h5['/output/data/%s' % v].attrs['description']
@@ -121 +208 @@
     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)
+
 with open('run_uq.xml', 'w') as f:
     f.write("<?xml version=\"1.0\"?>\n")