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

Timestamp:

May 17, 2015, 6:20:41 PM (9 years ago)

Author:

mmh

Message:

add uq summary output

File:

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

branches/uq/puq/puq_analyze.py

@@ -4 +4 @@
 process the results.
 """
+from __future__ import print_function
 import sys
 import os
@@ -9 +10 @@
 import h5py
 import re
+import puq
 from puq.jpickle import unpickle
 import xml.etree.ElementTree as xml
 from itertools import product, combinations
+import Rappture
+import StringIO
 
 # Redirect stdout and stderr to files for debugging.
@@ -46 +50 @@
 
 def plot_resp1(io, resp, name, h5):
-    print 'plot_resp1', name
+    print('plot_resp1', name)
     varlist = subs_names(resp.vars, h5)
     try:
-        tname = h5['/output/data/%s' % name].attrs['description']
+        tname = h5['/output/data/%s' % name].attrs['label']
     except:
         tname = name
@@ -56 +60 @@
 
     title = '%s (Response)' % (tname)
-    curve = 'output.curve(response-%s)' % name
-    io.put(curve + '.about.label', title)
-    io.put(curve + '.about.group', title)
-    io.put(curve + '.xaxis.label', varlist[0])
-    io.put(curve + '.yaxis.label', tname)
+
+    curve = io['output.curve(response-%s)' % name]
+    curve['about.label'] = title
+    curve['about.group'] = title
+    curve['xaxis.label'] = varlist[0]
+    curve['yaxis.label'] = tname
     x = np.linspace(*resp.vars[0][1], num=50)
     y = np.array(resp.eval(x))
-    for a, b in zip(x, y):
-        io.put(curve + '.component.xy', "%s %s\n" % (a, b), append=1)
+    curve['component.xy'] = (x, y)
 
     # scatter plot sampled data on response surface
-    curve = 'output.curve(response-%s-scatter)' % name
-    io.put(curve + '.about.label', 'Data Points')
-    io.put(curve + '.about.group', title)
-    io.put(curve + '.about.type', 'scatter')
-    io.put(curve + '.xaxis.label', varlist[0])
-    io.put(curve + '.yaxis.label', tname)
-    for a, b in zip(xdata, ydata):
-        io.put(curve + '.component.xy', "%s %s\n" % (a, b), append=1)
+    curve = io['output.curve(response-%s-scatter)' % name]
+    curve['about.label'] = 'Data Points'
+    curve['about.group'] = title
+    curve['about.type'] = 'scatter'
+    curve['xaxis.label'] = varlist[0]
+    curve['yaxis.label'] = tname
+    curve['component.xy'] = (xdata, ydata)
 
 
 def plot_resp2(io, resp, name, varl, h5):
+    print("plot_resp2", name)
     numpoints = 50
-    m2d = 'output.mesh('+varl[0][0]+varl[1][0]+')'
-    io.put(m2d+'.about.label', '2D Mesh')
-    io.put(m2d+'.dim', '2')
-    io.put(m2d+'.hide', 'yes')
-    io.put(m2d+'.grid.xaxis.numpoints', '%s' % numpoints)
-    io.put(m2d+'.grid.yaxis.numpoints', '%s' % numpoints)
-    io.put(m2d+'.grid.xaxis.min', '%g' % (varl[0][1][0]))
-    io.put(m2d+'.grid.xaxis.max', '%g' % (varl[0][1][1]))
-    io.put(m2d+'.grid.yaxis.min', '%g' % (varl[1][1][0]))
-    io.put(m2d+'.grid.yaxis.max', '%g' % (varl[1][1][1]))
+    m2d = io['output.mesh(%s%s)' % (varl[0][0], varl[1][0])]
+
+    m2d['about.label'] = '2D Mesh'
+    m2d['dim'] = 2
+    m2d['hide'] = 'yes'
+    m2d['grid.xaxis.numpoints'] = numpoints
+    m2d['grid.yaxis.numpoints'] = numpoints
+    m2d['grid.xaxis.min'] = varl[0][1][0]
+    m2d['grid.xaxis.max'] = varl[0][1][1]
+    m2d['grid.yaxis.min'] = varl[1][1][0]
+    m2d['grid.yaxis.max'] = varl[1][1][1]
 
     varlist = subs_names(varl, h5)
-    fname = 'output.field(response-'+name+varlist[0]+varlist[1]+')'
-    io.put(fname + '.about.xaxis.label', varlist[0])
-    io.put(fname + '.about.yaxis.label', varlist[1])
+    fname = io['output.field(response-'+name+varlist[0]+varlist[1]+')']
+    fname['about.xaxis.label'] = varlist[0]
+    fname['about.yaxis.label'] = varlist[1]
     try:
-        tname = h5['/output/data/%s' % name].attrs['description']
+        tname = h5['/output/data/%s' % name].attrs['label']
     except:
         tname = name
     title = tname + ' (Response %s)' % ' vs '.join(varlist)
-    io.put(fname + '.about.label', title)
-    io.put(fname + '.component.mesh', m2d)
-    io.put(fname + '.about.view', 'heightmap')
+    fname['about.label'] = title
+    fname['component.mesh'] = m2d.name
+    fname['about.view'] = 'heightmap'
     x = np.linspace(*varl[0][1], num=numpoints)
     y = np.linspace(*varl[1][1], num=numpoints)
@@ -114 +119 @@
             allpts[:, i] = np.mean(v[1])
     pts = np.array(resp.evala(allpts))
-    for z in pts:
-        io.put(fname + '.component.values', '%g\n' % z, append=1)
+    fname['component.values'] = pts
 
 
 def plot_resp3(io, resp, name, varl, h5):
     numpoints = 40
-    m3d = 'output.mesh(m3d)'
-    io.put(m3d+'.about.label', '3D Mesh')
-    io.put(m3d+'.dim', '3')
-    io.put(m3d+'.hide', 'yes')
-    io.put(m3d+'.grid.xaxis.numpoints', '%s' % numpoints)
-    io.put(m3d+'.grid.yaxis.numpoints', '%s' % numpoints)
-    io.put(m3d+'.grid.zaxis.numpoints', '%s' % numpoints)
+    m3d = io['output.mesh(m3d)']
+    m3d['about.label'] = '3D Mesh'
+    m3d['dim'] = 3
+    m3d['hide'] = 'yes'
+    m3d['grid.xaxis.numpoints'] = numpoints
+    m3d['grid.yaxis.numpoints'] = numpoints
+    m3d['grid.zaxis.numpoints'] = numpoints
 
     # for now, scale to [0,1]
-    io.put(m3d+'.grid.xaxis.min', '0')
-    io.put(m3d+'.grid.xaxis.max', '1')
-    io.put(m3d+'.grid.yaxis.min', '0')
-    io.put(m3d+'.grid.yaxis.max', '1')
-    io.put(m3d+'.grid.zaxis.min', '0')
-    io.put(m3d+'.grid.zaxis.max', '1')
+    m3d['grid.xaxis.min'] = 0
+    m3d['grid.xaxis.max'] = 1
+    m3d['grid.yaxis.min'] = 0
+    m3d['grid.yaxis.max'] = 1
+    m3d['grid.zaxis.min'] = 0
+    m3d['grid.zaxis.max'] = 1
 
     varlist = subs_names(varl, h5)
     try:
-        tname = h5['/output/data/%s' % name].attrs['description']
+        tname = h5['/output/data/%s' % name].attrs['label']
     except:
         tname = name
-    fname = 'output.field(response-'+name+varlist[0]+varlist[1]+varlist[2]+')'
-    io.put(fname + '.about.xaxis.label',
-                   '{%s [%.3g - %.3g]}' %
-                   (varlist[0], varl[0][1][0], varl[0][1][1]))
-    io.put(fname + '.about.yaxis.label',
-                   '{%s [%.3g - %.3g]}' %
-                   (varlist[1], varl[1][1][0], varl[1][1][1]))
-    io.put(fname + '.about.zaxis.label',
-                   '{%s [%.3g - %.3g]}' %
-                   (varlist[2], varl[2][1][0], varl[2][1][1]))
+    fname = io['output.field(response-'+name+varlist[0]+varlist[1]+varlist[2]+')']
+    fname['about.xaxis.label'] = '{%s [%.3g - %.3g]}' % (varlist[0], varl[0][1][0], varl[0][1][1])
+    fname['about.yaxis.label'] = '{%s [%.3g - %.3g]}' % (varlist[1], varl[1][1][0], varl[1][1][1])
+    fname['about.zaxis.label'] = '{%s [%.3g - %.3g]}' % (varlist[2], varl[2][1][0], varl[2][1][1])
     title = tname + ' (Response %s)' % ' vs '.join(varlist)
-    io.put(fname + '.about.label', title)
-    io.put(fname + '.component.mesh', m3d)
-    io.put(fname + '.about.description', '3D Field Description')
-    io.put(fname + '.about.view', 'vtkvolume')
+    fname['about.label'] = title
+    fname['component.mesh'] = m3d.name
+    fname['about.description'] = '3D Field Description'
+    fname['about.view'] = 'vtkvolume'
     x = np.linspace(*varl[0][1], num=numpoints)
     y = np.linspace(*varl[1][1], num=numpoints)
@@ -173 +171 @@
             allpts[:, i] = np.mean(v[1])
     pts = np.array(resp.evala(allpts))
-    for z in pts:
-        io.put(fname + '.component.values', '%g\n' % z, append=1)
+    fname['component.values'] = pts
 
 
@@ -193 +190 @@
 # Plots probability curves
 def plot_pdf_curve(xvals, vname, percent, desc):
-    print 'plot_pdf_curve %s %s %s' % (vname, percent, desc)
+    print('plot_pdf_curve %s %s %s' % (vname, percent, desc))
     # compute upper and lower percentiles
     pm = (100 - percent)/200.0
@@ -251 +248 @@
 # Plots advanced probability curves
 def plot_pdf_acurve(hf, acurves, vname, percent, desc):
-    print 'plot_pdf_acurve %s %s %s' % (vname, percent, desc)
+    print('plot_pdf_acurve %s %s %s' % (vname, percent, desc))
     # compute upper and lower percentiles
     pm = (100 - percent)/200.0
@@ -257 +254 @@
     vlen = len(acurves[vname])
 
-    print 'vlen=', vlen
+    print('vlen=', vlen)
     # collect data into an array
     xp = np.empty(vlen)
@@ -265 +262 @@
 
     for vindex in sorted(acurves[vname].keys()):
-        print 'vindex=', vindex
+        print('vindex=', vindex)
         x1 = acurves[vname][vindex][0].ppf(pp)
         x2 = acurves[vname][vindex][0].ppf(pm)
@@ -274 +271 @@
         #p1, p2 = get_closest2()
         if int(vindex) == 2:
-            print x1, x2, y1, y2
-            print 'xd=', repr(xd)
-            print 'yd=', repr(yd)
+            print (x1, x2, y1, y2)
+            print ('xd=', repr(xd))
+            print ('yd=', repr(yd))
 
         #xp[vindex] = x1
@@ -320 +317 @@
     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]]
+
+def write_responses(io, h5):
+    uqtype = h5.attrs['UQtype']
+    for v in h5[uqtype]:
+        if '[' in v:
+            # It is a curve. Ignore.
+            continue
+        try:
+            desc = h5['%s/%s' % (uqtype, v)].attrs['description']
+        except:
+            desc = ''
+        try:
+            label = h5['%s/%s' % (uqtype, v)].attrs['label']
+        except:
+            label = ''
+
+        rsp = h5['/%s/%s/response' % (uqtype, v)].value
+        rout = io['output.response(%s)' % label]
+        rout['value'] = rsp
+        rout['about.description'] = desc
+        rout['about.label'] = label
+
+        rs = unpickle(rsp)
+
+        if type(rs) == puq.response.ResponseFunc:
+            rout['equation'] = rs.eqn
+            rout['rmse'] = "{:6.3g}".format(rs.rmse()[1])
+
+        rout['data'] = rs.data
+
+
+def write_params(h5, out):
+    params = map(str, h5['/input/params'].keys())
+    print('#' * 80, file=out)
+    print('INPUT PARAMETERS', file=out)
+
+    for pname in params:
+        print('-' * 80, file=out)
+        p = puq.unpickle(h5['/input/params/' + pname].value)
+        cname = p.__class__.__name__[:-9]
+        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)
+
+        print("Name:", p.name, file=out)
+        try:
+            print("Label:", p.label, file=out)
+        except:
+            pass
+        print("Desc:", p.description, file=out)
+        print('Value:', pdf_str, file=out)
+    print('#' * 80, file=out)
+    print(file=out)
+
+
+def write_summary(io, h5):
+    outstr = StringIO.StringIO()
+    write_params(h5, outstr)
+    uqtype = h5.attrs['UQtype']
+    for v in h5[uqtype]:
+        if '[' in v:
+            # It is a curve. Ignore.
+            continue
+        desc = h5['%s/%s' % (uqtype, v)].attrs['description']
+        print("QoI: %s (%s)" % (v, desc), file=outstr)
+        rs = unpickle(h5['/%s/%s/response' % (uqtype, v)].value)
+        if type(rs) == puq.response.ResponseFunc:
+            print("\nv=%s\n" % rs.eqn, file=outstr)
+            print("SURROGATE MODEL ERROR:{:6.3g}%".format(rs.rmse()[1]), file=outstr)
+        sens = puq.unpickle(h5['/%s/%s/sensitivity' % (uqtype, v)].value)
+        max_name_len = max(map(len, [p[0] for p in sens]))
+        print("\nSENSITIVITY:", file=outstr)
+        print("Var%s     u*          dev" % (' '*(max_name_len)), file=outstr)
+        print('-'*(28+max_name_len), file=outstr)
+        for item in sens:
+            pad = ' '*(max_name_len - len(item[0]))
+            print("{}{}  {:10.4g}  {:10.4g}".format(pad, item[0],
+                item[1]['ustar'], item[1]['std']), file=outstr)
+        print('-'*(28+max_name_len), file=outstr)
+        print(file=outstr)
+    iostr = io['output.string(UQ Summary)']
+    iostr['about.label'] = 'UQ Summary'
+    iostr['current'] = outstr.getvalue()
+    outstr.close()
 
 sw = load_from_hdf5(sys.argv[1])
@@ -365 +465 @@
             acurves[vname][vindex][1] = pdf
     else:
-        desc = h5['/output/data/%s' % v].attrs['description']
+        desc = h5['/output/data/%s' % v].attrs['label']
         plot_pdf(v, pdf, desc)
 
@@ -388 +488 @@
         if '[' in v:
             continue
-        desc = h5['/output/data/%s' % v].attrs['description']
+        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})
@@ -415 +515 @@
 
 
-import Rappture
-io = Rappture.library('run_uq.xml')
-
+io = Rappture.PyXml('run_uq.xml')
 for v in h5[uqtype]:
     if '[' in v:
         continue
     rs = unpickle(h5['/%s/%s/response' % (uqtype, v)].value)
-    print 'Plotting response', v
+    print ('Plotting response', v)
     plot_resp(io, rs, v, h5)
 
-Rappture.result(io)
-
+# write_responses(io, h5)
+write_summary(io, h5)
+
+io.close()
 h5.close()