Changeset 5169 for branches/uq

Timestamp:

Mar 26, 2015, 6:48:40 AM (9 years ago)

Author:

mmh

Message:

fix multiple simulations

Location:

branches/uq

Files:

• gui/scripts/analyzer.tcl (modified) (5 diffs)
• gui/scripts/resultset.tcl (modified) (8 diffs)
• lang/tcl/scripts/library.tcl (modified) (11 diffs)
• lang/tcl/scripts/task.tcl (modified) (4 diffs)
• puq/get_params.py (modified) (1 diff)
• puq/puq_analyze.py (modified) (6 diffs)

branches/uq/gui/scripts/analyzer.tcl

@@ -107 +107 @@
 itcl::body Rappture::Analyzer::constructor {tool args} {
     set _tool $tool
-    puts "Analyzer init"
+
     # use this to store all simulation results
     set _resultset [Rappture::ResultSet ::\#auto]
@@ -420 +420 @@
             return
         }
-    }
-
-    lappend args -uq_type [$uq type]
-    lappend args -uq_args [$uq args]
-    puts "args=$args"
+        lappend args -uq_type [$uq type]
+        lappend args -uq_args [$uq args]
+    }
+
+    puts "simulate args=$args"
 
     # simulation is needed -- go to simulation page
@@ -501 +501 @@
 # ----------------------------------------------------------------------
 itcl::body Rappture::Analyzer::reset {{when -eventually}} {
+    puts "RESET $when"
     if {$when == "-eventually"} {
         after cancel [list catch [itcl::code $this reset -now]]
@@ -506 +507 @@
         return
     }
+
+    puts "RESET contains [$_resultset contains [$_tool xml object]]"
+    puts [$_tool xml object]
 
     # check to see if simulation is really needed
@@ -1123 +1127 @@
     }
 
-
     #
     # Break up the remaining lines according to =RAPPTURE-ERROR=> messages.

branches/uq/gui/scripts/resultset.tcl

@@ -1 @@
-# -*- mode: tcl; indent-tabs-mode: nil -*- 
+# -*- mode: tcl; indent-tabs-mode: nil -*-
 # ----------------------------------------------------------------------
 #  COMPONENT: ResultSet - set of XML objects for simulated results
@@ -40 +40 @@
     private variable _notify         ;# info used for notify command
 }
-                                                                                
+
 # ----------------------------------------------------------------------
 # CONSTRUCTOR
@@ -189 +189 @@
 # ----------------------------------------------------------------------
 itcl::body Rappture::ResultSet::diff {option args} {
+    puts "ResultSet::diff $option $args"
+    flush stdout
     switch -- $option {
         names {
@@ -198 +200 @@
             }
             set col [lindex $args 0]
-
+            puts "col=$col"
             set which "all"
             if {[llength $args] > 1} {
@@ -306 +308 @@
 # ----------------------------------------------------------------------
 itcl::body Rappture::ResultSet::contains {xmlobj} {
+    puts "CONTAINS [$_results size]"
     # no results? then this must be new
     if {[$_results size] == 0} {
@@ -316 +319 @@
     # is a column to represent the quantity with the difference.
     #
+    puts "CONTAINS A"
     set xmlobj0 [$_results get -format xmlobj end]
+    puts "CONTAINS B $xmlobj0 , $xmlobj"
+    puts "DIFF IS [$xmlobj0 diff $xmlobj]"
+    flush stdout
     foreach {op vpath oldval newval} [$xmlobj0 diff $xmlobj] {
+        puts "$op $vpath $oldval $newval"
+        flush stdout
         if {[$xmlobj get $vpath.about.diffs] == "ignore"} {
             continue
@@ -328 +337 @@
         if {[$_results column names $vpath] == ""} {
             # no column for this quantity yet
+            puts "NO COLUMN"
             return 0
         }
     }
-
+    puts "CONTAINS F"
     #
     # If we got this far, then look through existing results for
@@ -357 +367 @@
         }
     }
-
+    puts "CONTAINS Z"
     # must be some differences
     return 0

branches/uq/lang/tcl/scripts/library.tcl

@@ -70 +70 @@
         }
         set fname [file join $Rappture::installdir lib library.xml]
-
         set fid [::open $fname r]
         set info [read $fid]
@@ -215 +214 @@
     public method remove {{path ""}}
     public method xml {{path ""}}
-    public method uq_get_vars {{template 0}}
+    public method uq_get_vars {{tfile ""}}
 
     public method diff {libobj}
@@ -691 +690 @@
 # ----------------------------------------------------------------------
 itcl::body Rappture::LibraryObj::diff {libobj} {
+
+    puts "Library DIFF $libobj"
     set rlist ""
 
@@ -697 +698 @@
     set otherv [Rappture::entities $libobj input]
 
+    puts "thisv=$thisv\n"
+    puts "otherv=$otherv\n"
+
     # scan through values for this object, and compare against other one
     foreach path $thisv {
+      puts "path=$path"
         set i [lsearch -exact $otherv $path]
         if {$i < 0} {
@@ -704 +709 @@
             lappend rlist - $path $raw ""
         } else {
+          puts "tvalue=[value $this $path]"
+          puts "ovalue=[value $libobj $path]"
             foreach {traw tnorm} [value $this $path] break
             foreach {oraw onorm} [value $libobj $path] break
+            puts "tnorm=$tnorm onorm=$onorm traw=$traw oraw=$oraw"
             if {![string equal $tnorm $onorm]} {
                 lappend rlist c $path $traw $oraw
@@ -718 +726 @@
         lappend rlist + $path "" $oraw
     }
+    puts "LDIFF returning $rlist"
     return $rlist
 }
@@ -740 +749 @@
 # ----------------------------------------------------------------------
 itcl::body Rappture::LibraryObj::value {libobj path} {
+  puts "VALUE $path [$libobj element -as type $path]"
     switch -- [$libobj element -as type $path] {
         structure {
@@ -771 +781 @@
                 # then normalize to default units
                 set units [$libobj get $path.units]
+                puts "val=$val units=$units"
                 if {"" != $units} {
-                    set val [Rappture::Units::convert $val \
+                    set val [Rappture::Units::mconvert $val \
                         -context $units -to $units -units off]
                 }
@@ -1040 +1051 @@
 }
 
-itcl::body Rappture::LibraryObj::uq_get_vars {{template 0}} {
+itcl::body Rappture::LibraryObj::uq_get_vars {{tfile ""}} {
     set varlist ""
-    set n [$_node selectNodes /run/input//number]
+
+    if {$tfile == ""} {
+        set node $_node
+    } else {
+        set fid [::open $tfile r]
+        set doc [dom parse [read $fid]]
+        set node [$doc documentElement]
+        close $fid
+    }
+
+    set n [$node selectNodes /run/input//number]
     foreach _n $n {
         set x [$_n selectNodes current/text()]
@@ -1054 +1075 @@
                 -context $units -to $units -units off]
             }
-            if {$template != 0} {
+            if {$tfile != ""} {
               $x nodeValue @@[$_n getAttribute id]
             }
@@ -1060 +1081 @@
         }
     }
+
+    if {$tfile != ""} {
+        set fid [open $tfile w]
+        puts $fid "<?xml version=\"1.0\"?>"
+        puts $fid [$node asXML]
+        close $fid
+        $doc delete
+    }
     return $varlist
 }

branches/uq/lang/tcl/scripts/task.tcl

@@ -173 +173 @@
     # if there are any args, use them to override parameters
     set _outputcb ""
+    set uq_type ""
     foreach {path val} $args {
         if {$path == "-output"} {
@@ -211 +212 @@
     } result]
 
-    # This will turn the driver xml into a template
-    # and return a list of the UQ variables and their PDFs.
-    set uq_varlist [$_xmlobj uq_get_vars 1]
-
-    if {$uq_varlist != ""} {
-        # write out the template file for submit
+    if {$uq_type != ""} {
+        # Copy xml into a new file
         set tfile "template[pid].xml"
-        set status [catch {
-            set fid [open $tfile w]
-            puts $fid "<?xml version=\"1.0\"?>"
-            puts $fid [$_xmlobj xml]
-            close $fid
-        } result]
+        set fid [open $tfile w]
+        puts $fid "<?xml version=\"1.0\"?>"
+        puts $fid [$_xmlobj xml]
+        close $fid
+
+        # Return a list of the UQ variables and their PDFs.
+        # Also turns $tfile into a template file.
+        set uq_varlist [$_xmlobj uq_get_vars $tfile]
     }
 
@@ -239 +238 @@
         }
 
-        if {$uq_varlist == ""} {
+        if {$uq_type == ""} {
             regsub -all @driver $cmd $file cmd
 
@@ -404 +403 @@
     if {$status == 0} {
         set result [string trim $job(output)]
-        puts "result=$result"
-        if {$uq_varlist != ""} {
+
+        if {$uq_type != ""} {
             # UQ. Collect data from all jobs and put it in one xml run file.
-            file delete -force -- run_uq.xml
-            exec puq_analyze.py puq_[pid].hdf5
-            append result "\n" "=RAPPTURE-RUN=>run_uq.xml"
+            #file delete -force -- run_uq.xml
+            set res [exec puq_analyze.py puq_[pid].hdf5]
+            append result "\n" $res
         }
         if {[regexp {=RAPPTURE-RUN=>([^\n]+)} $result match file]} {

branches/uq/puq/get_params.py

@@ -55 +55 @@
 # This just saves PUQ state into HDF5 file, so later we can have PUQ analyze
 # the results and it knows about the input parameters, UQ method, etc.
-sw = Sweep(uq, RapptureHost('puq', dfile), None)
+host = RapptureHost('puq', dfile)
+sw = Sweep(uq, host, None)
+host.psweep = sw.psweep
 sw.run(hname, overwrite=True)
 

branches/uq/puq/puq_analyze.py

@@ -11 +11 @@
 from puq.jpickle import unpickle
 import xml.etree.ElementTree as xml
+from itertools import product, combinations
 
 # Redirect stdout and stderr to files for debugging.
@@ -16 +17 @@
 sys.stdout = open("uq_debug.out", 'a')
 sys.stderr = open("uq_debug.err", 'a')
+
 
 # Restore the state of a PUQ session from a HDF5 file.
@@ -29 +31 @@
         sw.psweep.reinit()
     return sw
+
+
+# do some name substitution for better labels
+def subs_names(varl):
+    varlist = []
+    for v in varl:
+        varlist.append(str(v[0]))
+    return varlist
+
+
+def plot_resp1(io, resp, name):
+    print 'plot_resp1', name
+    varlist = subs_names(resp.vars)
+    curve = 'output.curve(response-%s)' % name
+    io.put(curve + '.about.label', '%s Response for %s' % (name, varlist[0]))
+    io.put(curve + '.xaxis.label', varlist[0])
+    io.put(curve + '.yaxis.label', name)
+    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)
+
+
+def plot_resp2(io, resp, name, varl):
+    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]))
+
+    varlist = subs_names(varl)
+    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])
+    title = name + ' Response (%s)' % ' vs '.join(varlist)
+    io.put(fname + '.about.label', title)
+    io.put(fname + '.component.mesh', m2d)
+    io.put(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))
+    for z in pts:
+        io.put(fname + '.component.values', '%g\n' % z, append=1)
+
+def plot_resp3(io, resp, name, varl):
+    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)
+
+    # 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')
+
+    varlist = subs_names(varl)
+    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]))
+
+    title = name + ' 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')
+    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))
+    for z in pts:
+        io.put(fname + '.component.values', '%g\n' % z, append=1)
+
+
+def plot_resp(io, resp, name):
+    if resp is not None:
+        if len(resp.vars) == 1:
+            plot_resp1(io, resp, name)
+            return
+        if len(resp.vars) >= 3:
+            for v1, v2, v3 in combinations(resp.vars, 3):
+                plot_resp3(io, resp, name, [v1,v2,v3])
+
+        # plot all combinations of 2 variables
+        for v1, v2 in combinations(resp.vars, 2):
+            plot_resp2(io, resp, name, [v1,v2])
 
 
@@ -68 +197 @@
 
 
+
 def dist(x1,y1, x2,y2, x3,y3): # x3,y3 is the point
     diffx = x2-x1
@@ -87 +217 @@
 
     return np.sqrt(dx*dx + dy*dy)
+
 
 # Plots advanced probability curves
@@ -220 +351 @@
 #    plot_pdf_acurve(h5, acurves, vname, 95, desc)
 
+
+# If more than one variable, display sensitivity.
+# Curves have indexed variables, so skip them.
+if ['[' in x for x in h5[uqtype]].count(False):
+    for v in h5[uqtype]:
+        if '[' in v:
+            continue
+        lab = h5['/output/data/%s' % v].attrs['description']
+        sens = unpickle(h5['/%s/%s/sensitivity' % (uqtype, v)].value)
+        elem = xml.SubElement(dout, 'histogram', {'id': 'sens-%s' % v})
+        about = xml.SubElement(elem, 'about')
+        xml.SubElement(about, 'label').text = '%s Sensitivity' % lab
+        xml.SubElement(elem, 'type').text = 'scatter'
+        xaxis = xml.SubElement(elem, 'xaxis')
+        xml.SubElement(xaxis, 'label').text = 'Parameters'
+        yaxis = xml.SubElement(elem, 'yaxis')
+        xml.SubElement(yaxis, 'label').text = 'Sensitivity'
+        comp = xml.SubElement(elem, 'component')
+        xy = xml.SubElement(comp, 'xy')
+        pts = ''
+        for name in sens:
+            n = name[0]
+            try:
+                n = h5['/input/params/%s' % n].attrs['label']
+            except:
+                pass
+            pts += "%s %s\n" % (n, name[1]['ustar'])
+        xy.text = pts
+
 with open('run_uq.xml', 'w') as f:
     f.write("<?xml version=\"1.0\"?>\n")
     dtree.write(f)
 
-with open('run_uq2.xml', 'w') as f:
-    f.write("<?xml version=\"1.0\"?>\n")
-    dtree.write(f)
+
+import Rappture
+io = Rappture.library('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)
+
+Rappture.result(io)
+
 h5.close()