Appendix: Working with the examples in this report

This report, together with its associated files, constitutes a reproducible research document. Makefiles tie together the various components - text, code source, data, and postscript figures. The principal make rules are the SEP standards: build, view, clean, burn. In this section, I will assume familiarity with the Stanford Exploration Project reproducible research concept, which guides the structure of this document.

This code behind this document is an application of the Hilbert Class Library. So the first thing you need to do is to make HCL available. If HCL is already installed on your system, you do this by adding a line to your .cshrc file or appropriate component file. Otherwise you must install HCL first. The easiest way to do this is to download it from the TRIP web page:

http://www.trip.caam.rice.edu

and follow the installation instructions.

The code also depends on the SU/SEGY vector class package sVector, which therefore must also be installed. It will be part of the next HCL release.

NB: At Rice/TRIP and Stanford/SEP, no installation is necessary: the packages are already installed. Simply add the following lines to your shell environment files:

• at TRIP, add to your .cshrc:

  setenv HCLROOT /import/masc39c/symes/hclr0.9
  setenv KBDAROOT /import/masc39c/symes/kbda
  setenv QPROOT <path to the root directory of this package>
  

• at SEP, add to your Setup/cshrc.generic:

  setenv HCLROOT /jon/symes/hclr0.9
  setenv KBDAROOT /jon/symes/kbda
  setenv QPROOT <path to the root directory of this package>
  

The ``root directory of this package'' referenced in these instructions is the directory you create by downloading the tar file containing this report. In so doing, you create a directory tree with root named qcqm. This is the ``root'' in question. All pathnames in the following discussion are relative to this root.

HCL includes a set of make rules which evolved from the SEP rule set as it was about two years ago. I am sure that SEP's rules have also evolved, and differently. The makefiles for this report are all output of the HCL makefile autowriter, maw, and use HCL's rules. If you get as far as modifying makefiles by hand, bear the possible HCL/SEP incompatibility in mind.

You can rebuild the entire package simply by entering make build in the root directory. You will construct all of the executables and final results (.ps files in the Fig) directories, including this postscript version of this report. You can make individual results in the usual way.

Note that all figures in this report except the first two will vary from build to build, as you choose a new pseudorandom seed each time you execute the commands.

The command for deconvolution is sfilter/decon.x. You execute it by following it with a parameter file name: decon.x par - it reads all of its parameters from the file. Regrettably the ``getpar'' device used in curren HCL programs is not flexible enough to permit specification of parameters on the command line. Probably I should just steal SEP's getpar!

You can use the executable sfilter/decon.x with other data by altering its input parameters, and so explore the capabilities of the algorithm using data other than that supplied with this report. Parameter control requires you to edit the parameter files manually. HCL parameter files are keyword=value lists; the values can be integers, floating point numbers with any size of mantissa, and strings. Parameters to be read only by one part of the program (typically a class constructor) get an identifying string prepended, with a double colon. Thus the parameter Tol for the conjugate gradient algorithm becomes CG::Tol. That is, the parameter file can specify many variables named Tol, so ong as they have been equipped with qualifiers which allow each program unit to choose a unique value.

Files should be in either SU (stripped SEGY) or SEGY formats. sVector currently supports only native binary floating point representation, so if you port data from Linux to SGI etc. you will have to byte-swap it.

Here is the parameter file structure for the deconvolution example:

Sigma=0.5
Lambda=0.0001
Wavelet="rick15.su"
DataTimeSeries="fnd.su"
CG::Tol=1.e-4
CG::MaxItn=50
CG::DispFlag=3
QCQM::Tol=0.01
QCQM::MaxItn=10
QCQM::DispFlag=1

The parameters are

• Sigma: target noise level
• Lambda: initial estimate of $\lambda$. It may be worth thinking about sensible defaults or crude estimates for this, and the scalar model actually suggests such an estimate. For the moment, set by hand.
• Wavelet: name of wavelet data file,in either SU or SEGY format.
• DataTimeSeries: filename for data time series
• CG::Tol: convergence tolerance for HCL conjugate code - iteration terminates if normal residual, i.e. in this case $\lambda A^T(Ax-d) +R^TRx$
• CG::MaxItn: maximum conjugate gradient iterations
• CG::DispFlag: controls verbosity of CG output, as explained in HCL reference manual (through TRIP web page: point your browser at

  http://www.trip.caam.rice.edu
  

  and follow the links to the HCL reference manual, page on the conjugate gradient algorithms). Level 3 is max output, including a summary of the progress at each iteration. Note that the residual reported is the normal residual for this application!
• QCQM::Tol: the Moré-Hebden algorithm considers itself converged when the relative error in the constraint is within this amount
• QCQM::MaxItn: maximum Moré-Hebden iterations
• QCQM::DispFlag: controls verbosity of Moré-Hebden code: 0 = silent, any other value prints diagnostic information about run

Amongst the diagnostics printed out at the end of a run when QCQM::DispFlag is set you will find ``Lagrange cosine''. This is the cosine of the angle between the constraint gradient and the objective gradient. It diagnoses the success of the constrained optimization: if it is very close to 1, then the two gradients are parallel and the first order necessary condition has been satisfied. This occurs in the deconvolution examples in all cases except that depicted in Figure 8, in which nothing converges and you can't fit the data. Apparently failure to get the Lagrange cosine close to 1 in a reasonable number of CG and Moré-Hebden iterations implies that the noise level has been set too small and you are trying to match data components associated with very small singular values.