[cig-commits] r20261 - short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip
willic3 at geodynamics.org
willic3 at geodynamics.org
Thu May 31 18:15:09 PDT 2012
Author: willic3
Date: 2012-05-31 18:15:09 -0700 (Thu, 31 May 2012)
New Revision: 20261
Added:
short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/penalty_params.txt
short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/run-inv.sh
short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/simple_invert.py
Log:
Added simple inversion code to make use of GF.
Still need to update README and include code to plot results.
Added: short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/penalty_params.txt
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/penalty_params.txt (rev 0)
+++ short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/penalty_params.txt 2012-06-01 01:15:09 UTC (rev 20261)
@@ -0,0 +1,7 @@
+0.00001
+0.0001
+0.001
+0.01
+0.1
+1.0
+10.0
Added: short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/run-inv.sh
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/run-inv.sh (rev 0)
+++ short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/run-inv.sh 2012-06-01 01:15:09 UTC (rev 20261)
@@ -0,0 +1,2 @@
+#!/bin/bash
+./simple_invert.py --impulse_file=output/greensfns-fault.h5 --response_file=output/greensfns-points.h5 --data_file=output/eqsim-points.h5 --penalty_file=penalty_params.txt --output_file=strikeslip_inversion.txt
Property changes on: short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/run-inv.sh
___________________________________________________________________
Name: svn:executable
+ *
Added: short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/simple_invert.py
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/simple_invert.py (rev 0)
+++ short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/simple_invert.py 2012-06-01 01:15:09 UTC (rev 20261)
@@ -0,0 +1,157 @@
+#!/usr/bin/env python
+"""
+This is an extremely simple example showing how to set up and inversion using
+PyLith-generated Green's functions. In this simple example, there are no data
+uncertainties, and as a penalty function we use a simple minimum moment.
+"""
+
+# The code requires the numpy and h5py packages.
+import numpy
+import h5py
+# import pdb
+
+
+def getImpResp():
+ """
+ Function to get impulse and response coordinates and values. Both are sorted
+ along the y-dimension of the fault.
+ """
+
+ # Open impulse file and determine which fault vertices were used.
+ impulses = h5py.File(impulseFile, "r")
+ impC = impulses['geometry/vertices'][:]
+ impV = impulses['vertex_fields/slip'][:,:,0]
+ impInds = numpy.nonzero(impV != 0.0)
+ impCUsed = impC[impInds[1]]
+ impVUsed = impV[impInds[0], impInds[1]]
+
+ # Sort by y-coordinate.
+ impInds = numpy.argsort(impCUsed[:,1])
+ impCSort = impCUsed[impInds,:]
+ impVSort = impVUsed[impInds]
+
+ # Close impulse file and get responses.
+ impulses.close()
+ responses = h5py.File(responseFile, "r")
+ respC = responses['geometry/vertices'][:]
+ respV = responses['vertex_fields/displacement'][:]
+ respVIsort = respV[impInds,:,:]
+ responses.close()
+
+ return (impCSort, impVSort, respC, respVIsort)
+
+
+def getData():
+ """
+ Function to get data values and coordinates.
+ """
+
+ # Open data file. Since the response and data coordinates are in the same
+ # order, we don't have to worry about sorting.
+ data = h5py.File(dataFile, "r")
+ dataC = data['geometry/vertices'][:]
+ dataV = data['vertex_fields/displacement'][:]
+ data.close()
+
+ return (dataC, dataV)
+
+# pdb.set_trace()
+
+# The main part of the code is below.
+# Get command-line arguments.
+from optparse import OptionParser
+parser = OptionParser()
+parser.add_option("-i", "--impulse_file", action="store", type="string",
+ dest="impulse_file", help="HDF5 file with fault GF info")
+parser.add_option("-r", "--response_file", action="store", type="string",
+ dest="response_file", help="HDF5 file with GF responses")
+parser.add_option("-d", "--data_file", action="store", type="string",
+ dest="data_file", help="HDF5 file with data")
+parser.add_option("-p", "--penalty_file", action="store", type="string",
+ dest="penalty_file", help="text file with penalty parameters")
+parser.add_option("-o", "--output_file", action="store", type="string",
+ dest="output_file", help="text file with estimated slip")
+
+(options, args) = parser.parse_args()
+
+if (not options.impulse_file):
+ parser.error("impulse input file must be specified")
+
+if (not options.response_file):
+ parser.error("response input file must be specified")
+
+if (not options.data_file):
+ parser.error("data input file must be specified")
+
+if (not options.penalty_file):
+ parser.error("penalty input file must be specified")
+
+if (not options.output_file):
+ parser.error("output file must be specified")
+
+impulseFile = options.impulse_file
+responseFile = options.response_file
+dataFile = options.data_file
+penaltyFile = options.penalty_file
+outputFile = options.output_file
+
+# Get GF info.
+(impCoords, impVals, respCoords, respVals) = getImpResp()
+
+# Get observed displacements and observation locations.
+(dataCoords, dataVals) = getData()
+
+# Get penalty parameters.
+penalties = numpy.loadtxt(penaltyFile, dtype=numpy.float64)
+
+# Determine matrix sizes and set up A-matrix.
+numParams = impVals.shape[0]
+numObs = 2 * dataVals.shape[1]
+aMat = respVals.reshape((numParams, numObs)).transpose()
+
+# Create diagonal matrix to use as the penalty.
+parDiag = numpy.eye(numParams, dtype=numpy.float64)
+
+# Data vector is a flattened version of the dataVals, plus the a priori
+# values of the parameters (assumed to be zero).
+dataVec = numpy.concatenate((dataVals.flatten(),
+ numpy.zeros(numParams, dtype=numpy.float64)))
+
+# Determine number of inversions and create empty array to hold results.
+numInv = penalties.shape[0]
+invResults = numpy.zeros((numParams, 2 + numInv))
+invResults[:,0:2] = impCoords
+head = "X_Coord\tY_Coord"
+
+# Loop over number of inversions.
+for inversion in range(numInv):
+ penalty = penalties[inversion]
+ head += "\tPenalty=%g" % penalty
+
+ # Scale diagonal by penalty parameter, and stack A-matrix with penalty matrix.
+ penMat = penalty * parDiag
+ designMat = numpy.vstack((aMat, penMat))
+ designMatTrans = designMat.transpose()
+
+ # Form generalized inverse matrix.
+ genInv = numpy.dot(numpy.linalg.inv(numpy.dot(designMatTrans, designMat)),
+ designMatTrans)
+
+ # Solution is matrix-vector product of generalized inverse with data vector.
+ solution = numpy.dot(genInv, dataVec)
+ invResults[:,2 + inversion] = solution
+
+ # Compute predicted results and residual.
+ predicted = numpy.dot(aMat, solution)
+ residual = dataVals.flatten() - predicted
+ residualNorm = numpy.linalg.norm(residual)
+ print "Penalty parameter: %g" % penalty
+ print " Residual norm: %g" % residualNorm
+
+head += "\n"
+
+# Output results.
+f = open(outputFile, "w")
+f.write(head)
+numpy.savetxt(f, invResults, delimiter="\t")
+f.close()
Property changes on: short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/simple_invert.py
___________________________________________________________________
Name: svn:executable
+ *
More information about the CIG-COMMITS
mailing list