[cig-commits] r17924 - short/3D/PyLith/trunk/examples/greensfns/hex8

[willic3 at geodynamics.org]

Author: willic3
Date: 2011-02-20 17:18:37 -0800 (Sun, 20 Feb 2011)
New Revision: 17924

Modified:
   short/3D/PyLith/trunk/examples/greensfns/hex8/designdata.py
   short/3D/PyLith/trunk/examples/greensfns/hex8/getcoords.py
   short/3D/PyLith/trunk/examples/greensfns/hex8/lininvert.py
Log:
Fixed some bugs and added more output info.



Modified: short/3D/PyLith/trunk/examples/greensfns/hex8/designdata.py
===================================================================
--- short/3D/PyLith/trunk/examples/greensfns/hex8/designdata.py	2011-02-21 01:17:43 UTC (rev 17923)
+++ short/3D/PyLith/trunk/examples/greensfns/hex8/designdata.py	2011-02-21 01:18:37 UTC (rev 17924)
@@ -2,11 +2,17 @@
 #
 # ----------------------------------------------------------------------
 #
-#                           Brad T. Aagaard
-#                        U.S. Geological Survey
+# Brad T. Aagaard, U.S. Geological Survey
+# Charles A. Williams, GNS Science
+# Matthew G. Knepley, University of Chicago
 #
-# <LicenseText>
+# This code was developed as part of the Computational Infrastructure
+# for Geodynamics (http://geodynamics.org).
 #
+# Copyright (c) 2010 University of California, Davis
+#
+# See COPYING for license information.
+#
 # ----------------------------------------------------------------------
 #
 

Modified: short/3D/PyLith/trunk/examples/greensfns/hex8/getcoords.py
===================================================================
--- short/3D/PyLith/trunk/examples/greensfns/hex8/getcoords.py	2011-02-21 01:17:43 UTC (rev 17923)
+++ short/3D/PyLith/trunk/examples/greensfns/hex8/getcoords.py	2011-02-21 01:18:37 UTC (rev 17924)
@@ -105,7 +105,8 @@
     lines = f.readlines()
     for line in lines:
       for entry in line.split(', '):
-        self.idList.append(int(entry.rstrip(',\n')))
+        idEntry = int(entry.rstrip(',\n')) - 1
+        self.idList.append(idEntry)
 
     self.numIds = len(self.idList)
     f.close()
@@ -118,23 +119,10 @@
     Function to read a list of coordinates and output the coordinates if
     they correspond to one of the requested ID's.
     """
-    f = open(self.coordFile, 'r')
-    o = open(self.outputFile, 'w')
-    lines = f.readlines()
-    for line in lines:
-      vertexId = int(line.split(', ')[0].rstrip(',\n'))
-      for requestedId in self.idList:
-        if vertexId == requestedId:
-          x = line.split(', ')[1].rstrip(',')
-          y = line.split(', ')[2].rstrip(',')
-          z = line.split(', ')[3].rstrip(',\n')
-          outLine = x + '   ' + y + '   ' + z + '\n'
-          o.write(outLine)
-          break
+    meshCoords = numpy.loadtxt(self.coordFile, dtype=numpy.float64)
+    outputCoords = meshCoords[self.idList,:]
+    numpy.savetxt(self.outputFile, outputCoords)
 
-    f.close()
-    o.close()
-
     return
   
 

Modified: short/3D/PyLith/trunk/examples/greensfns/hex8/lininvert.py
===================================================================
--- short/3D/PyLith/trunk/examples/greensfns/hex8/lininvert.py	2011-02-21 01:17:43 UTC (rev 17923)
+++ short/3D/PyLith/trunk/examples/greensfns/hex8/lininvert.py	2011-02-21 01:18:37 UTC (rev 17924)
@@ -160,9 +160,9 @@
 
     self.numObs = self.dataVals.shape[0]
     self.numParams = self.aprioriVals.shape[0]
-    # print "Number of observations: %i" % self.numObs
-    # print "Number of parameters: %i" % self.numParams
-    # sys.stdout.flush()
+    print "Number of observations: %i" % self.numObs
+    print "Number of parameters: %i" % self.numParams
+    sys.stdout.flush()
 
     return
 
@@ -195,21 +195,29 @@
     self.fData = numpy.diag(fDataVec)
     self.predicted = numpy.dot(self.dataDesign, self.solution)
     dataMisfit = self.dataVals - self.predicted
-    weightMisfit = dataMisfit * fDataVec
-    self.chiSquare = numpy.sum(weightMisfit * weightMisfit)
+    dataWeightMisfit = dataMisfit * fDataVec
+    self.chiSquare = numpy.sum(dataWeightMisfit * dataWeightMisfit)
+    dataMisfitNorm = numpy.linalg.norm(dataMisfit)
+    dataWeightMisfitNorm = numpy.linalg.norm(dataWeightMisfit)
     
     # Write out inversion info.
     string1 = "\nNumber of observations:  %i\n" % self.numObs
     string2 = "Number of parameters:  %i\n" % self.numParams
-    string3 = "Chi-square value:  %e\n" % self.chiSquare
+    string3 = "Data Chi-square value:  %e\n" % self.chiSquare
+    string4 = "Data residual norm:  %e\n" % dataMisfitNorm
+    string5 = "Weighted data residual norm:  %e\n" % dataWeightMisfitNorm
     print string1
     print string2
     print string3
+    print string4
+    print string5
     sys.stdout.flush()
     i = open(self.infoOutputFile, 'w')
     i.write(string1)
     i.write(string2)
     i.write(string3)
+    i.write(string4)
+    i.write(string5)
     i.close()
                                
     # Write out data info.
@@ -222,7 +230,7 @@
     
     dataOut = numpy.transpose(numpy.vstack((self.dataVals, self.predicted,
                                             fDataVecInv, dataMisfit,
-                                            weightMisfit)))
+                                            dataWeightMisfit)))
     numpy.savetxt(self.dataOutputFile, dataOut)
 
     return