[cig-commits] r20265 - in short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns: . reverse strikeslip
brad at geodynamics.org
brad at geodynamics.org
Fri Jun 1 09:05:22 PDT 2012
Author: brad
Date: 2012-06-01 09:05:21 -0700 (Fri, 01 Jun 2012)
New Revision: 20265
Added:
short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/README
short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/invert_slip.py
short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/reverse/plot_invresults.py
short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/plot_invresults.py
Removed:
short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/reverse/make-plot.sh
short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/reverse/plot_results.py
short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/reverse/run-inv.sh
short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/reverse/simple_invert.py
short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/make-plot.sh
short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/plot_results.py
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
Modified:
short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/Makefile.am
short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/reverse/Makefile.am
short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/reverse/README
short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/Makefile.am
short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/README
Log:
Factored out common stuff (README, invert_slip.py) and small cleanup of python scripts.
Modified: short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/Makefile.am
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/Makefile.am 2012-06-01 15:21:56 UTC (rev 20264)
+++ short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/Makefile.am 2012-06-01 16:05:21 UTC (rev 20265)
@@ -16,6 +16,11 @@
# ----------------------------------------------------------------------
#
+dist_noinst_DATA = \
+ README \
+ invert_slip.py
+
+
SUBDIRS = \
strikeslip \
reverse
Added: short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/README
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/README (rev 0)
+++ short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/README 2012-06-01 16:05:21 UTC (rev 20265)
@@ -0,0 +1,117 @@
+This directory contains two examples that demonstrate computing static
+Green's functions via PyLith and using them in a simple inversion for
+fault slip. The directory strikeslip contains a 2-D strike-slip fault
+example and the directory reverse contains a 2-D reverse fault
+example. The two examples are run using exactly the same commands.
+
+The main features of these examples are:
+
+ * Generating a finite-element mesh using CUBIT
+ + Variable mesh resolution
+ * Spatially variable slip
+ * Computing Green's functions using PyLith
+
+Each example consists of two simulations. In the first simulation we
+compute the displacements due to a variable slip on the fault fault
+(idealized in a 2-D model with plane-strain); this serves as the
+"observed" displacements. In the second simulation we compute Green's
+functions for the fault. We use these Green's functions in the second
+simulation and the observed displacements from the first simulation to
+invert for the slip (in the first simulation).
+
+The imposed slip in the first simulation includes a region of uniform
+slip with a linear taper at each end. We use uniform linear elastic
+properties with a plane-strain formulation.
+
+The parameters for the bulk constitutive models are defined in
+ mat_elastic.spatialdb
+
+The simulation will output the displacements on the ground surface as
+well as a selection of points (defined in output_points.txt).
+
+For each of the simulations, we recommend examining the displacements,
+stress field, and fault slip.
+
+
+Mesh generation (optional)
+
+ NOTE: The result of this step will overwrite the included file
+ tri3_mesh.exo. You may want to copy/rename this file so that
+ you have a backup copy in case you have difficulty running
+ CUBIT.
+
+ Start CUBIT and play the journal file "mesh_tri3.jou". We highly
+ recommend that you study the contents of the journal files to
+ understand the mesh generation process.
+
+
+Step 1. Forward simulation
+
+ This simulation mimics an earthquake generating an observed
+ displacement field.
+
+ The parameters for the earthquake slip are defined in
+ eqslip.spatialdb
+
+ Run the simulation via the following command:
+ pylith eqsim.cfg
+
+
+Step 2. Generate Green's functions
+
+ This simulation generates Green's functions for a subset of the fault.
+
+ The parameters for the amplitude (and sign) of the slip impulses are
+ defined in
+ impulse_amplitude.spatialdb
+
+ Run the simulation via the following command:
+ pylith --problem=pylith.problems.GreensFns
+
+
+Step 3. Invert for coseismic slip
+
+ Now that we have the forward simulation (data) and Green's functions, we
+ can perform an inversion for the fault slip. We have written a very
+ simple Python script that performs a miminim moment solution using a
+ range of penalty parameters (contained in penalty_params.txt). The inversion
+ code is in invert_slip.py. You can run the inversion by typing:
+
+ ../invert_slip.py --impulses=output/greensfns-fault.h5 --responses=output/greensfns-points.h5 --data=output/eqsim-points.h5 --penalty=penalty_params.txt --output=output/slip_inverted.txt
+
+ This will generate a set of predicted slip values, contained in
+ output/slip_inverted.txt.
+
+
+Step 4. Plot the true and predicted slip values (**requires matplotlib**).
+
+ The python script, plot_invresults.py, compares the inversion with the
+ forward simulation. If you have matplotlib installed, you can run
+ this script by typing:
+
+ plot_invresults.py --solution=output/eqsim-fault.h5 --predicted=output/slip_inverted.txt
+
+ This will display a matplotlib window with the true solution in
+ black and the predicted solution corresponding to different values
+ of the penalty parameter.
+
+
+Suggestions variations
+
+ The list below includes some suggested modifications to the problem
+ that will allow you to become more familiar with PyLith while
+ examining some interesting physics.
+
+ * Create a spatial variation in material properties. For example you
+ might create a velocity contrast across the fault.
+
+ * Adjust the slip distribution.
+
+ * Adjust the location and number of the observation points.
+
+ * Increase/decrease the resolution of the mesh.
+
+ * Try a different penalty function for the inversion. One simple option
+ would be a 1D Laplacian approximation (sets of values of [-1, 2, -1])
+ centered along the diagonal.
+
Copied: short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/invert_slip.py (from rev 20263, short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/reverse/simple_invert.py)
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/invert_slip.py (rev 0)
+++ short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/invert_slip.py 2012-06-01 16:05:21 UTC (rev 20265)
@@ -0,0 +1,159 @@
+#!/usr/bin/env python
+"""
+This is an extremely simple example showing how to set up an inversion
+using PyLith-generated Green's functions. In this simple example,
+there are no data uncertainties, and for a penalty function we use a
+simple minimum moment.
+"""
+
+# The code requires the numpy and h5py packages.
+import numpy
+import h5py
+
+
+# ----------------------------------------------------------------------
+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", driver="sec2")
+ 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", driver="sec2")
+ 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", driver="sec2")
+ dataC = data['geometry/vertices'][:]
+ dataV = data['vertex_fields/displacement'][:]
+ data.close()
+
+ return (dataC, dataV)
+
+
+# ======================================================================
+# The main part of the code is below.
+# Get command-line arguments.
+from optparse import OptionParser
+parser = OptionParser()
+parser.add_option("-i", "--impulses", action="store", type="string",
+ dest="impulse_file", help="HDF5 file with fault GF info")
+parser.add_option("-r", "--responses", action="store", type="string",
+ dest="response_file", help="HDF5 file with GF responses")
+parser.add_option("-d", "--data", action="store", type="string",
+ dest="data_file", help="HDF5 file with data")
+parser.add_option("-p", "--penalty", action="store", type="string",
+ dest="penalty_file", help="text file with penalty parameters")
+parser.add_option("-o", "--output", 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 Y_Coord"
+
+# Loop over number of inversions.
+for inversion in range(numInv):
+ penalty = penalties[inversion]
+ head += " Penalty=%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, fmt="%14.6e")
+f.close()
Modified: short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/reverse/Makefile.am
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/reverse/Makefile.am 2012-06-01 15:21:56 UTC (rev 20264)
+++ short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/reverse/Makefile.am 2012-06-01 16:05:21 UTC (rev 20265)
@@ -29,7 +29,9 @@
eqslip.spatialdb \
impulse_amplitude.spatialdb \
mat_elastic.spatialdb \
- output_points.txt
+ output_points.txt \
+ penalty_params.txt \
+ plot_invresults.py
SUBDIRS = \
Modified: short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/reverse/README
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/reverse/README 2012-06-01 15:21:56 UTC (rev 20264)
+++ short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/reverse/README 2012-06-01 16:05:21 UTC (rev 20265)
@@ -1,114 +1,2 @@
-The example in this directory demonstrates the use of computing static
-Green's functions using PyLith and using them in a simple inversion
-for slip on a reverse fault.
-
-The main features of this example are:
-
- * Generating a finite-element mesh using CUBIT
- + Variable mesh resolution
- * Spatially variable slip
- * Computing Green's functions using PyLith
-
-The examples consists of two simulations. In the first simulation we
-compute the displacements due to a variable slip on a reverse
-fault (idealized in a 2-D model with plane-strain); this serves as the
-"observed" displacements. In the second simulation we compute Green's
-functions for the fault. We use these Green's functions in the second
-simulation and the observed displacements from the first simulation to
-invert for the slip (in the first simulation).
-
-The imposed slip in the first simulation includes a region of uniform
-slip with a linear taper at each end. We use uniform linear elastic
-properties with a plane-strain formulation.
-
-The parameters for the bulk constitutive models are defined in
- mat_elastic.spatialdb
-
-The simulation will output the displacements on the ground surface as
-well as a selection of points (defined in output_points.txt).
-
-For each of the simulations, we recommend examining the displacements,
-stress field, and fault slip.
-
-
-Mesh generation (optional)
-
- NOTE: The result of this step will overwrite the included file
- tri3_mesh.exo. You may want to copy/rename this file so that
- you have a backup copy in case you have difficulty running
- CUBIT.
-
- Start CUBIT and play the journal file "mesh_tri3.jou". We highly
- recommend that you study the contents of the journal files to
- understand the mesh generation process.
-
-
-Step 1. Forward simulation
-
- This simulation mimics an earthquake generating an observed
- displacement field.
-
- The parameters for the earthquake slip are defined in
- eqslip.spatialdb
-
- Run the simulation via the following command:
- pylith eqsim.cfg.
-
-
-Step 2. Generate Green's functions
-
- This simulation generates Green's functions for a subset of the fault.
-
- The parameters for the amplitude (and sign) of the slip impulses are
- defined in
- impulse_amplitude.spatialdb
-
- Run the simulation via the following command:
- pylith --problem=pylith.problems.GreensFns
-
-
-Step 3. Invert for coseismic slip
-
- Now that we have the forward simulation (data) and Green's functions, we
- can perform an inversion for the fault slip. We have written a very
- simple Python script that performs a miminim moment solution using a
- range of penalty parameters (contained in penalty_params.txt). The Python
- code is in simple_invert.py, and a script to run it is contained in
- run-inv.sh. You can run the inversion by simply typing:
- ./run-inv.sh
-
- This will generate a set of predicted slip values, contained in
- strikeslip_inversion.txt.
-
-
-Step 4. Plot the true and predicted slip values.
-
- To see how the predicted solution compares with the applied slip, we
- use a simple Python script (plot_results.py) that makes use of
- matplotlib. You can run this script by simply typing:
- ./make-plot.sh
-
- This will bring up a matplotlib plot window with the true solution (shown
- in black) and the predicted solution corresponding to different values
- of the penalty parameter.
-
-
-Suggestions variations
-
- The list below includes some suggested modifications to the problem
- that will allow you to become more familiar with PyLith while
- examining some interesting physics.
-
- * Create a spatial variation in material properties. For example you
- might create a velocity contrast across the fault.
-
- * Adjust the slip distribution.
-
- * Adjust the location and number of the observation points.
-
- * Increase/decrease the resolution of the mesh.
-
- * Try a different penalty function for the inversion. One simple option
- would be a 1D Laplacian approximation (sets of values of [-1, 2, -1])
- centered along the diagonal.
-
+See the README one level up (examples/2d/greensfns/README) for a
+description and detailed instructions.
Deleted: short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/reverse/make-plot.sh
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/reverse/make-plot.sh 2012-06-01 15:21:56 UTC (rev 20264)
+++ short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/reverse/make-plot.sh 2012-06-01 16:05:21 UTC (rev 20265)
@@ -1,2 +0,0 @@
-#!/bin/bash
-./plot_results.py --solution_file=output/eqsim-fault.h5 --predicted_file=reverse_inversion.txt
Copied: short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/reverse/plot_invresults.py (from rev 20263, short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/reverse/plot_results.py)
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/reverse/plot_invresults.py (rev 0)
+++ short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/reverse/plot_invresults.py 2012-06-01 16:05:21 UTC (rev 20265)
@@ -0,0 +1,83 @@
+#!/usr/bin/env python
+"""
+This script generates a set of line plots comparing the predicted solution to
+the true solution.
+"""
+
+# The code requires the numpy, h5py, and matplotlib packages.
+import numpy
+import h5py
+import matplotlib.pyplot as pyplot
+import math
+
+# Define line colors and other parameters.
+solnColor = "black"
+predColors = ["cyan", "red", "blue", "green", "orange", "purple", "yellow"]
+
+# ----------------------------------------------------------------------
+def getSolution():
+ """
+ Function to get applied slip and fault coordinates.
+ """
+
+ # Open solution file and get slip and coordinates.
+ solution = h5py.File(solutionFile, "r")
+ solC = solution['geometry/vertices'][:]
+ solV = solution['vertex_fields/slip'][:,:,0].flatten()
+
+ # Sort by y-coordinate.
+ solInds = numpy.argsort(solC[:,1])
+ solCSort = solC[solInds,:]
+ solVSort = solV[solInds]
+ # Reverse sign, since this is right-lateral slip.
+ solVSort *= -1.0
+ solution.close()
+
+ return (solCSort, solVSort)
+
+
+# ======================================================================
+# The main part of the code is below.
+# Get command-line arguments.
+from optparse import OptionParser
+parser = OptionParser()
+parser.add_option("-s", "--solution", action="store", type="string",
+ dest="solution_file",
+ help="HDF5 file with applied fault slip")
+parser.add_option("-p", "--predicted", action="store", type="string",
+ dest="predicted_file", help="file with predicted solutions")
+
+(options, args) = parser.parse_args()
+
+if not options.solution_file:
+ parser.error("Solution input file must be specified.")
+
+if not options.predicted_file:
+ parser.error("Predicted input file must be specified.")
+
+solutionFile = options.solution_file
+predictedFile = options.predicted_file
+
+# Get applied fault slip.
+(solCoords, solVals) = getSolution()
+
+# Get predicted fault slip.
+predicted = numpy.loadtxt(predictedFile)
+predCoords = predicted[:,0:2]
+predSlip = predicted[:, 2:]
+
+# Generate figure, starting with true solution.
+pyplot.plot(solVals, -solCoords[:,1]/1.0e+3, linewidth=2, color=solnColor)
+pyplot.xlabel("Reverse slip (m)")
+pyplot.ylabel("Depth (km)")
+pyplot.ylim( (22, 0) )
+
+# Loop over predicted results, using a different line color for each.
+numInv = predSlip.shape[1]
+
+for inversion in range(numInv):
+ pyplot.plot(predSlip[:,inversion], -predCoords[:,1]/1.0e+3,
+ color=predColors[inversion])
+
+
+pyplot.show()
Deleted: short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/reverse/plot_results.py
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/reverse/plot_results.py 2012-06-01 15:21:56 UTC (rev 20264)
+++ short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/reverse/plot_results.py 2012-06-01 16:05:21 UTC (rev 20265)
@@ -1,85 +0,0 @@
-#!/usr/bin/env python
-"""
-This script generates a set of line plots comparing the predicted solution to
-the true solution.
-"""
-
-# The code requires the numpy, h5py, and matplotlib packages.
-import numpy
-import h5py
-import matplotlib.pyplot as pyplot
-# import pdb
-
-# Define line colors and other parameters.
-solnColor = "black"
-predColors = ["cyan", "red", "blue", "green", "orange", "purple", "yellow"]
-style = "lightbg"
-fontsize = 8
-
-
-def getSolution():
- """
- Function to get applied slip and fault coordinates.
- """
-
- # Open solution file and get slip and coordinates.
- solution = h5py.File(solutionFile, "r")
- solC = solution['geometry/vertices'][:]
- solV = solution['vertex_fields/slip'][:,:,0].flatten()
-
- # Sort by y-coordinate.
- solInds = numpy.argsort(solC[:,1])
- solCSort = solC[solInds,:]
- solVSort = solV[solInds]
- # Reverse sign, since this is right-lateral slip.
- solVSort *= -1.0
- solution.close()
-
- return (solCSort, solVSort)
-
-
-# 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", "--solution_file", action="store", type="string",
- dest="solution_file",
- help="HDF5 file with applied fault slip")
-parser.add_option("-r", "--predicted_file", action="store", type="string",
- dest="predicted_file", help="file with predicted solutions")
-
-(options, args) = parser.parse_args()
-
-if (not options.solution_file):
- parser.error("solution input file must be specified")
-
-if (not options.predicted_file):
- parser.error("predicted input file must be specified")
-
-solutionFile = options.solution_file
-predictedFile = options.predicted_file
-
-# Get applied fault slip.
-(solCoords, solVals) = getSolution()
-
-# Get predicted fault slip.
-predicted = numpy.loadtxt(predictedFile, skiprows=1, dtype=numpy.float64)
-predCoords = predicted[:,0:2]
-predSlip = predicted[:, 2:]
-
-# Generate figure, starting with true solution.
-pyplot.plot(solCoords[:,1], solVals, color=solnColor)
-pyplot.xlabel("Y-distance (m)")
-pyplot.ylabel("Reverse slip (m)")
-
-# Loop over predicted results, using a different line color for each.
-numInv = predSlip.shape[1]
-
-for inversion in range(numInv):
- pyplot.plot(predCoords[:,1], predSlip[:,inversion],
- color=predColors[inversion])
-
-
-pyplot.show()
Deleted: short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/reverse/run-inv.sh
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/reverse/run-inv.sh 2012-06-01 15:21:56 UTC (rev 20264)
+++ short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/reverse/run-inv.sh 2012-06-01 16:05:21 UTC (rev 20265)
@@ -1,2 +0,0 @@
-#!/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=reverse_inversion.txt
Deleted: short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/reverse/simple_invert.py
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/reverse/simple_invert.py 2012-06-01 15:21:56 UTC (rev 20264)
+++ short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/reverse/simple_invert.py 2012-06-01 16:05:21 UTC (rev 20265)
@@ -1,157 +0,0 @@
-#!/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()
Modified: short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/Makefile.am
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/Makefile.am 2012-06-01 15:21:56 UTC (rev 20264)
+++ short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/Makefile.am 2012-06-01 16:05:21 UTC (rev 20265)
@@ -29,7 +29,10 @@
eqslip.spatialdb \
impulse_amplitude.spatialdb \
mat_elastic.spatialdb \
- output_points.txt
+ output_points.txt \
+ invert_slip.py \
+ penalty_params.txt \
+ plot_invresults.py
SUBDIRS = \
Modified: short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/README
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/README 2012-06-01 15:21:56 UTC (rev 20264)
+++ short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/README 2012-06-01 16:05:21 UTC (rev 20265)
@@ -1,113 +1,2 @@
-The example in this directory demonstrates the use of computing static
-Green's functions using PyLith and using them in a simple inversion
-for slip on a strike-slip fault.
-
-The main features of this example are:
-
- * Generating a finite-element mesh using CUBIT
- + Variable mesh resolution
- * Spatially variable slip
- * Computing Green's functions using PyLith
-
-The examples consists of two simulations. In the first simulation we
-compute the displacements due to a variable slip on a strike-slip
-fault (idealized in a 2-D model with plane-strain); this serves as the
-"observed" displacements. In the second simulation we compute Green's
-functions for the fault. We use these Green's functions in the second
-simulation and the observed displacements from the first simulation to
-invert for the slip (in the first simulation).
-
-The imposed slip in the first simulation includes a region of uniform
-slip with a linear taper at each end. We use uniform linear elastic
-properties with a plane-strain formulation.
-
-The parameters for the bulk constitutive models are defined in
- mat_elastic.spatialdb
-
-The simulation will output the displacements on the ground surface as
-well as a selection of points (defined in output_points.txt).
-
-For each of the simulations, we recommend examining the displacements,
-stress field, and fault slip.
-
-
-Mesh generation (optional)
-
- NOTE: The result of this step will overwrite the included file
- tri3_mesh.exo. You may want to copy/rename this file so that
- you have a backup copy in case you have difficulty running
- CUBIT.
-
- Start CUBIT and play the journal file "mesh_tri3.jou". We highly
- recommend that you study the contents of the journal files to
- understand the mesh generation process.
-
-
-Step 1. Forward simulation
-
- This simulation mimics an earthquake generating an observed
- displacement field.
-
- The parameters for the earthquake slip are defined in
- eqslip.spatialdb
-
- Run the simulation via the following command:
- pylith eqsim.cfg.
-
-
-Step 2. Generate Green's functions
-
- This simulation generates Green's functions for a subset of the fault.
-
- The parameters for the amplitude (and sign) of the slip impulses are
- defined in
- impulse_amplitude.spatialdb
-
- Run the simulation via the following command:
- pylith --problem=pylith.problems.GreensFns
-
-
-Step 3. Invert for coseismic slip
-
- Now that we have the forward simulation (data) and Green's functions, we
- can perform an inversion for the fault slip. We have written a very
- simple Python script that performs a miminim moment solution using a
- range of penalty parameters (contained in penalty_params.txt). The Python
- code is in simple_invert.py, and a script to run it is contained in
- run-inv.sh. You can run the inversion by simply typing:
- ./run-inv.sh
-
- This will generate a set of predicted slip values, contained in
- strikeslip_inversion.txt.
-
-
-Step 4. Plot the true and predicted slip values.
-
- To see how the predicted solution compares with the applied slip, we
- use a simple Python script (plot_results.py) that makes use of
- matplotlib. You can run this script by simply typing:
- ./make-plot.sh
-
- This will bring up a matplotlib plot window with the true solution (shown
- in black) and the predicted solution corresponding to different values
- of the penalty parameter.
-
-
-Suggestions variations
-
- The list below includes some suggested modifications to the problem
- that will allow you to become more familiar with PyLith while
- examining some interesting physics.
-
- * Create a spatial variation in material properties. For example you
- might create a velocity contrast across the fault.
-
- * Adjust the slip distribution.
-
- * Adjust the location and number of the observation points.
-
- * Increase/decrease the resolution of the mesh.
-
- * Try a different penalty function for the inversion. One simple option
- would be a 1D Laplacian approximation (sets of values of [-1, 2, -1])
- centered along the diagonal.
+See the README one level up (examples/2d/greensfns/README) for a
+description and detailed instructions.
Deleted: short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/make-plot.sh
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/make-plot.sh 2012-06-01 15:21:56 UTC (rev 20264)
+++ short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/make-plot.sh 2012-06-01 16:05:21 UTC (rev 20265)
@@ -1,2 +0,0 @@
-#!/bin/bash
-./plot_results.py --solution_file=output/eqsim-fault.h5 --predicted_file=strikeslip_inversion.txt
Copied: short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/plot_invresults.py (from rev 20263, short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/plot_results.py)
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/plot_invresults.py (rev 0)
+++ short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/plot_invresults.py 2012-06-01 16:05:21 UTC (rev 20265)
@@ -0,0 +1,81 @@
+#!/usr/bin/env python
+"""
+This script generates a plot comparing the predicted solutions to
+the true solution.
+"""
+
+# The code requires the numpy, h5py, and matplotlib packages.
+import numpy
+import h5py
+import matplotlib.pyplot as pyplot
+
+# Define line colors and other parameters.
+solnColor = "black"
+predColors = ["cyan", "red", "blue", "green", "orange", "purple", "yellow"]
+
+
+# ----------------------------------------------------------------------
+def getSolution():
+ """
+ Function to get applied slip and fault coordinates.
+ """
+
+ # Open solution file and get slip and coordinates.
+ solution = h5py.File(solutionFile, "r")
+ solC = solution['geometry/vertices'][:]
+ solV = solution['vertex_fields/slip'][:,:,0].flatten()
+
+ # Sort by y-coordinate.
+ solInds = numpy.argsort(solC[:,1])
+ solCSort = solC[solInds,:]
+ solVSort = solV[solInds]
+ # Reverse sign, since this is right-lateral slip.
+ solVSort *= -1.0
+ solution.close()
+
+ return (solCSort, solVSort)
+
+# ======================================================================
+# The main part of the code is below.
+# Get command-line arguments.
+from optparse import OptionParser
+parser = OptionParser()
+parser.add_option("-s", "--solution", action="store", type="string",
+ dest="solution_file",
+ help="HDF5 file with applied fault slip")
+parser.add_option("-p", "--predicted", action="store", type="string",
+ dest="predicted_file", help="file with predicted solutions")
+
+(options, args) = parser.parse_args()
+
+if not options.solution_file:
+ parser.error("Solution input file must be specified.")
+
+if not options.predicted_file:
+ parser.error("Predicted input file must be specified.")
+
+solutionFile = options.solution_file
+predictedFile = options.predicted_file
+
+# Get applied fault slip.
+(solCoords, solVals) = getSolution()
+
+# Get predicted fault slip.
+predicted = numpy.loadtxt(predictedFile)
+predCoords = predicted[:,0:2]
+predSlip = predicted[:, 2:]
+
+# Generate figure, starting with true solution.
+pyplot.plot(solCoords[:,1]/1.0e+3, solVals, linewidth=2, color=solnColor)
+pyplot.xlabel("Distance Along Strike (km)")
+pyplot.ylabel("Right-Lateral Slip (m)")
+
+# Loop over predicted results, using a different line color for each.
+numInv = predSlip.shape[1]
+
+for inversion in range(numInv):
+ pyplot.plot(predCoords[:,1]/1.0e+3, predSlip[:,inversion],
+ color=predColors[inversion])
+
+
+pyplot.show()
Deleted: short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/plot_results.py
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/plot_results.py 2012-06-01 15:21:56 UTC (rev 20264)
+++ short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/plot_results.py 2012-06-01 16:05:21 UTC (rev 20265)
@@ -1,85 +0,0 @@
-#!/usr/bin/env python
-"""
-This script generates a set of line plots comparing the predicted solution to
-the true solution.
-"""
-
-# The code requires the numpy, h5py, and matplotlib packages.
-import numpy
-import h5py
-import matplotlib.pyplot as pyplot
-# import pdb
-
-# Define line colors and other parameters.
-solnColor = "black"
-predColors = ["cyan", "red", "blue", "green", "orange", "purple", "yellow"]
-style = "lightbg"
-fontsize = 8
-
-
-def getSolution():
- """
- Function to get applied slip and fault coordinates.
- """
-
- # Open solution file and get slip and coordinates.
- solution = h5py.File(solutionFile, "r")
- solC = solution['geometry/vertices'][:]
- solV = solution['vertex_fields/slip'][:,:,0].flatten()
-
- # Sort by y-coordinate.
- solInds = numpy.argsort(solC[:,1])
- solCSort = solC[solInds,:]
- solVSort = solV[solInds]
- # Reverse sign, since this is right-lateral slip.
- solVSort *= -1.0
- solution.close()
-
- return (solCSort, solVSort)
-
-
-# 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", "--solution_file", action="store", type="string",
- dest="solution_file",
- help="HDF5 file with applied fault slip")
-parser.add_option("-r", "--predicted_file", action="store", type="string",
- dest="predicted_file", help="file with predicted solutions")
-
-(options, args) = parser.parse_args()
-
-if (not options.solution_file):
- parser.error("solution input file must be specified")
-
-if (not options.predicted_file):
- parser.error("predicted input file must be specified")
-
-solutionFile = options.solution_file
-predictedFile = options.predicted_file
-
-# Get applied fault slip.
-(solCoords, solVals) = getSolution()
-
-# Get predicted fault slip.
-predicted = numpy.loadtxt(predictedFile, skiprows=1, dtype=numpy.float64)
-predCoords = predicted[:,0:2]
-predSlip = predicted[:, 2:]
-
-# Generate figure, starting with true solution.
-pyplot.plot(solCoords[:,1], solVals, color=solnColor)
-pyplot.xlabel("Y-distance (m)")
-pyplot.ylabel("Right-lateral slip (m)")
-
-# Loop over predicted results, using a different line color for each.
-numInv = predSlip.shape[1]
-
-for inversion in range(numInv):
- pyplot.plot(predCoords[:,1], predSlip[:,inversion],
- color=predColors[inversion])
-
-
-pyplot.show()
Deleted: 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 2012-06-01 15:21:56 UTC (rev 20264)
+++ short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/run-inv.sh 2012-06-01 16:05:21 UTC (rev 20265)
@@ -1,2 +0,0 @@
-#!/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
Deleted: 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 2012-06-01 15:21:56 UTC (rev 20264)
+++ short/3D/PyLith/branches/v1.7-trunk/examples/2d/greensfns/strikeslip/simple_invert.py 2012-06-01 16:05:21 UTC (rev 20265)
@@ -1,157 +0,0 @@
-#!/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()
More information about the CIG-COMMITS
mailing list