Changeset 5148


Ignore:
Timestamp:
Mar 17, 2015, 6:18:04 PM (5 years ago)
Author:
mmh
Message:

uq snap

Location:
branches/uq/puq
Files:
1 deleted
2 edited

Legend:

Unmodified
Added
Removed
  • branches/uq/puq/get_params.py

    r5029 r5148  
    11#!/usr/bin/env python
    2 from puq import *
    3 import csv, sys, os
     2import sys
     3from puq import NormalParameter, UniformParameter, Sweep, RapptureHost, Smolyak
     4import csv, os
    45import numpy as np
    56from puq.jpickle import unpickle
     7
     8sys.stdout = open("gp.out", 'w')
     9sys.stderr = open("gp.err", 'w')
     10
    611
    712dfile, pid, varlist, uq_type, args = sys.argv[1:]
     
    914hname = "puq_%s.hdf5" % pid
    1015
    11 
    12 sys.stdout = open("gp.log", 'w')
    13 sys.stderr = open("gp.err", 'w')
    14 
    1516print sys.argv[1:]
    1617varlist = unpickle(varlist)
    1718print "varlist=", varlist
    18 print type(varlist)
     19
    1920
    2021v = {}
  • branches/uq/puq/puq_analyze.py

    r5103 r5148  
    11#!/usr/bin/env python
    22import sys
     3import os
     4import numpy as np
     5import h5py
     6import re
     7from puq.jpickle import unpickle
     8import xml.etree.ElementTree as xml
     9
    310sys.stdout = open("analyze.out", 'w')
    411sys.stderr = open("analyze.err", 'w')
    512
    6 import puq
    7 import csv
    8 import numpy as np
    9 import os, h5py
    10 from puq.jpickle import unpickle
    11 import xml.etree.ElementTree as xml
    1213
    13 
     14# Restore the state of a PUQ session from a HDF5 file.
    1415def load_from_hdf5(name):
    1516    h5 = h5py.File(name, 'r+')
     
    2425    return sw
    2526
    26 sw = load_from_hdf5(sys.argv[1])
    27 sw.analyze()
    2827
    29 h5 = h5py.File(sys.argv[1], 'r+')
     28# Plots probability curves
     29def plot_pdf_curve(xvals, vname, percent, desc):
     30    print 'plot_pdf_curve %s %s %s' % (vname, percent, desc)
     31    # compute upper and lower percentiles
     32    pm = (100 - percent)/200.0
     33    pp = 1 - pm
     34
     35    # collect data into an array
     36    xarr = np.empty(len(xvals[vname]))
     37    yp = np.empty(len(xvals[vname]))
     38    ym = np.empty(len(xvals[vname]))
     39    for vindex in sorted(xvals[vname].keys()):
     40        xarr[vindex] = xvals[vname][vindex]
     41        yp[vindex] = pcurves[vname][vindex].ppf(pp)
     42        ym[vindex] = pcurves[vname][vindex].ppf(pm)
     43
     44    curve = xml.SubElement(dout, 'curve', {'id': 'curve_pdf-%s-%s' % (vname, percent)})
     45    about = xml.SubElement(curve, 'about')
     46    if percent == 0:
     47        xml.SubElement(about, 'label').text = "mean"
     48    else:
     49        xml.SubElement(about, 'label').text = "middle %s%%" % percent
     50    xml.SubElement(about, 'group').text = '%s (Probability)' % desc
     51    comp = xml.SubElement(curve, 'component')
     52    xy = xml.SubElement(comp, 'xy')
     53    pts = ""
     54    for x, y in zip(xarr, yp):
     55        pts += "%s %s " % (x, y)
     56    if percent == 0:
     57        pts += '\n'
     58    else:
     59        for x, y in reversed(zip(xarr, ym)):
     60            pts += "%s %s " % (x, y)
     61        pts += "%s %s\n" % (xarr[0], yp[0])
     62    xy.text = pts
    3063
    3164
    32 dtree = xml.parse('run_uq.xml')
    33 droot = dtree.getroot()
    34 dout = droot.find('output')
    35 
    36 uqtype = h5.attrs['UQtype']
    37 for v in h5[uqtype]:
    38     print 'PDF', v
    39     rs = unpickle(h5['/%s/%s/response' % (uqtype, v)].value)
    40     pdf = rs.pdf(fit=False)
    41 
    42     id = '%s PDF' % v
     65def plot_pdf(v, pdf, desc):
     66    id = '%s (PDF)' % desc
    4367    elem = xml.SubElement(dout, 'curve', {'id': id})
    4468    about = xml.SubElement(elem, 'about')
    45     label = xml.SubElement(about, 'label')
    46     label.text = id
    47     desc = xml.SubElement(about, 'description')
    48     desc.text = "PDF for %s" % id
     69    xml.SubElement(about, 'label').text = id
    4970    yaxis = xml.SubElement(elem, 'yaxis')
    50     label = xml.SubElement(yaxis, 'label')
    51     label.text = 'Probability'
     71    xml.SubElement(yaxis, 'label').text = 'Probability'
    5272
    5373    component = xml.SubElement(elem, 'component')
     
    5979    xy.text = pts
    6080
     81
     82sw = load_from_hdf5(sys.argv[1])
     83sw.analyze()
     84
     85h5 = h5py.File(sys.argv[1], 'r+')
     86dtree = xml.parse('run_uq.xml')
     87droot = dtree.getroot()
     88dout = droot.find('output')
     89
     90# curves built from pdfs
     91pcurves = {}
     92xvals = {}
     93
     94regexp = re.compile('([ \da-z]+)\[([ \d]+)\]')
     95
     96uqtype = h5.attrs['UQtype']
     97for v in h5[uqtype]:
     98    rs = unpickle(h5['/%s/%s/response' % (uqtype, v)].value)
     99    pdf = rs.pdf(fit=False)
     100    odata = h5['/output/data/%s' % v]
     101
     102    # For curves built from pdfs, just put them in a dict for now
     103    if 'x' in odata.attrs:
     104        matches = regexp.findall(v)
     105        vname, vindex = matches[0]
     106        vindex = int(vindex)
     107        if vname not in pcurves:
     108            pcurves[vname] = {}
     109            xvals[vname] = {}
     110        xvals[vname][vindex] = odata.attrs['x']
     111        pcurves[vname][vindex] = pdf
     112    else:
     113        desc = h5['/output/data/%s' % v].attrs['description']
     114        plot_pdf(v, pdf, desc)
     115
     116# now do pdf curves
     117for vname in xvals:
     118    desc = h5['/output/data/%s[0]' % vname].attrs['description']
     119    plot_pdf_curve(xvals, vname, 0, desc)
     120    plot_pdf_curve(xvals, vname, 50, desc)
     121    plot_pdf_curve(xvals, vname, 95, desc)
     122
    61123with open('run_uq.xml', 'w') as f:
    62124    f.write("<?xml version=\"1.0\"?>\n")
    63125    dtree.write(f)
    64126
    65 
     127with open('run_uq2.xml', 'w') as f:
     128    f.write("<?xml version=\"1.0\"?>\n")
     129    dtree.write(f)
    66130h5.close()
Note: See TracChangeset for help on using the changeset viewer.