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

Timestamp:

Jun 16, 2015 9:18:59 PM (9 years ago)

Author:

mmh

Message:

remove response surface code

File:

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

branches/uq/puq/puq_analyze.py

@@ -13 +13 @@
 from puq.jpickle import unpickle
 import xml.etree.ElementTree as xml
-from itertools import product, combinations
 import Rappture
 import StringIO
@@ -35 +34 @@
         sw.psweep.reinit()
     return sw
-
-
-# variable names to labels
-def subs_names(varl, h5):
-    varlist = []
-    for v in varl:
-        try:
-            lab = h5['/input/params/%s' % v[0]].attrs['label']
-        except:
-            lab = str(v[0])
-        varlist.append(lab)
-    return varlist
-
-
-def plot_resp1(io, resp, name, h5):
-    print('plot_resp1', name)
-    varlist = subs_names(resp.vars, h5)
-    try:
-        tname = h5['/output/data/%s' % name].attrs['label']
-    except:
-        tname = name
-    ydata = h5['/output/data/%s' % name].value
-    xdata = h5['/input/param_array'].value.T[0]
-
-    title = '%s (Response)' % (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))
-    curve['component.xy'] = (x, y)
-
-    # scatter plot sampled data on response surface
-    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 = 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 = 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['label']
-    except:
-        tname = name
-    title = tname + ' (Response %s)' % ' vs '.join(varlist)
-    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)
-    pts = np.array([(b, a) for a, b in product(y, x)])
-    allpts = np.empty((numpoints**2, len(resp.vars)))
-    for i, v in enumerate(resp.vars):
-        if v[0] == varl[0][0]:
-            allpts[:, i] = pts[:, 0]
-        elif v[0] == varl[1][0]:
-            allpts[:, i] = pts[:, 1]
-        else:
-            allpts[:, i] = np.mean(v[1])
-    pts = np.array(resp.evala(allpts))
-    fname['component.values'] = pts
-
-
-def plot_resp3(io, resp, name, varl, h5):
-    numpoints = 40
-    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]
-    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['label']
-    except:
-        tname = name
-    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)
-    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)
-    z = np.linspace(*varl[2][1], num=numpoints)
-
-    # generate points in the right order, with x changing fastest
-    pts = np.array([(c, b, a) for a, b, c in product(z, y, x)])
-    allpts = np.empty((numpoints**3, len(resp.vars)))
-    for i, v in enumerate(resp.vars):
-        if v[0] == varl[0][0]:
-            allpts[:, i] = pts[:, 0]
-        elif v[0] == varl[1][0]:
-            allpts[:, i] = pts[:, 1]
-        elif v[0] == varl[2][0]:
-            allpts[:, i] = pts[:, 2]
-        else:
-            allpts[:, i] = np.mean(v[1])
-    pts = np.array(resp.evala(allpts))
-    fname['component.values'] = pts
-
-
-def plot_resp(io, resp, name, h5):
-    if resp is not None:
-        if len(resp.vars) == 1:
-            plot_resp1(io, resp, name, h5)
-            return
-        if len(resp.vars) >= 3:
-            for v1, v2, v3 in combinations(resp.vars, 3):
-                plot_resp3(io, resp, name, [v1, v2, v3], h5)
-
-        # plot all combinations of 2 variables
-        for v1, v2 in combinations(resp.vars, 2):
-            plot_resp2(io, resp, name, [v1, v2], h5)
 
 
@@ -338 +186 @@
 #  [  34.52960806  163.81660645]
 #  [  55.47039194  346.04150004]]
+
 
 def write_responses(io, h5):
@@ -355 +204 @@
 
         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
+        rout['about.label'] = label + " (Response)"
 
         rs = unpickle(rsp)
+        rout['variables'] = ' '.join([str(p.name) for p in rs.params])
+        labels = ' '.join([repr(str(p.label)) for p in rs.params])
+        rout['labels'] = labels.replace("'", '"')
 
         if type(rs) == puq.response.ResponseFunc:
@@ -439 +292 @@
 uqtype = h5.attrs['UQtype']
 for v in h5[uqtype]:
-    rs = unpickle(h5['/%s/%s/response' % (uqtype, v)].value)
+    rsp = h5['/%s/%s/response' % (uqtype, v)].value
+    print("rsp=", rsp)
+    rs = unpickle(rsp)
     pdf = rs.pdf(fit=False)
     odata = h5['/output/data/%s' % v]
@@ -514 +369 @@
     dtree.write(f)
 
-
 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)
-    plot_resp(io, rs, v, h5)
-
-# write_responses(io, h5)
+write_responses(io, h5)
 write_summary(io, h5)
-
 io.close()
+
 h5.close()