Madagascar development blog

Open Data/Open Source: Seismic Unix scripts to process a 2D land line by Karl Schleicher is an example of a reproducible paper created with Seismic Unix. It is also a unique example of generating a full processing sequence with a publicly-available dataset, from field data to the final image. Work is under way for generating more examples and for translating the processing flow from SU (suproj) to Madagascar (proj).

Scientists' Elusive Goal: Reproducing Study Results , an article published on the first page of the Wall Street Journal, describes the crisis of scientific reproducibility in bio-medical research.




Many of the issues described in the article sound familiar.

...Reproducibility is the foundation of all modern research, the standard by which scientific claims are evaluated. In the U.S. alone, biomedical research is a $100-billion-year enterprise. So when published medical findings can't be validated by others, there are major consequences...

...There is also a more insidious and pervasive problem: a preference for positive results...

...Some studies can't be redone for a more prosaic reason: the authors won't make all their raw data available to rival scientists...

Geophysical research does not affect human lives directly but its quality can suffer from non-reproducibility in very much the same way.

tkvpconvert provides a Tkinter GUI (graphical user interface) for vpconvert. It is a silly exercise, which simply demonstrates an ability to generate specialized GUIs for Madagascar scripts and programs.

sfsmooth, one of the most useful Madagascar programs, implements smoothing by triangle filtering.

The idea of triangle filtering is explained in Jon Claerbout's books Processing versus Inversion and Image Estimation by Example. You can find the explanation in the chapter Smoothing with box and triangle.

The triangle filter of radius is a correlation of two box filters. In the Z-transform notation,
where and k is the smoothing radius.

Triangle filtering is a very efficient operation, because it requires only N multiplications for the input of size N.

Triangle filtering serves as a good approximation to more expensive Gaussian filtering. Repeated triangle filtering rapidly approaches Gaussian filtering (a consequence of the central limit theorem). The following example from jsg/shape/smoo demonstrates how repeated application of triangle smoothing turns a triangle filter into a Gaussian filter


The following example from gee/ajt/triangle demonstrated how sfsmooth handles boundary conditions: the energy is reflected from the side in order to preserve the zero-frequency (DC) component of the input signal. Repeated smoothing or smoothing with a very large radius turns every signal into a constant.


In multiple dimensions, triangle smoothing is applied sequentially in different directions, possible with different radii (controlled by rect#= ) parameter. The following example from gee/ajt/mound shows shapes created by 2-D triangle smoothing after one pass and two passes (repeat=2).


Triangle smoothing is theoretically a self-adjoint operator. Numerically, the adjoint and forward operators are different. The action is controlled by the adj= parameter.

For smoothing with a box instead of a triangle, use sfboxsmooth.

For a different approach to smoothing, try sfgaussmooth, which implements smoothing by recursive Gaussian filtering.

Previous programs of the month:

• sfnoise
• sfgraph
• sfclip
• sfagc
• sfenvelope
• sfcontour

A short user survey was conducted after the 2011 Madagascar School in Beijing.



The results are overwhelmingly positive. 100% of those who responded to the survey stated that they would be interested in attending Madagascar events in the future, and 100% would recommend it to their colleagues. As for the location of a future event, 95% suggested Beijing again.

Many of those surveyed liked the excellent organization of the school. They disliked the fact that the school was too short and did not have enough space to accommodate all interested students. Organizers of future schools should take all survey suggestions into account.

A new paper is added to the collection of reproducible documents:

Time-frequency analysis of seismic data using local attributes

sfcontour is a program for generating contour plots of 2-D data. It shares many of its parameters with other 2-D plotting programs, such as sfgrey and sfgraph. You can access these parameters by running sfdoc stdplot.

Let us discuss some of the parameters specific to sfcontour.

The number of contours, their origin and increment are given by nc=, c0=, dc=. If these parameters are not specified, they are determined from the data values. The contour lines in the following picture from jsg/agath/radon were created with sfcontour c0=10 dc=10 nc=10



Alternatively, the contours can be specified in a list with c= or in a file with cfile=.

By default, sfcontour plots only contours for the positive values of the input. If you want both positive and negative values to be represented, use allpos=n.

A 3-D version of sfcontour is sfcontour3.

The contour lines in the following picture from jsg/flat/comaz were created with sfcontour3 cfile=

A new paper is added to the collection of reproducible documents:

Seismic data interpolation beyond aliasing using regularized nonstationary autoregression

A new paper is added to the collection of reproducible documents:

Fourier finite-difference wave propagation

Complex trace attributes were introduced into geophysics by the paper

Taner, M. T., F. Koehler, and R. E. Sheriff, 1979, Complex seismic trace analysis: Geophysics, 44, 1041-1063.

If is the input seismic trace, then the analytical trace is defined as the complex-valued signal
where is the Hilbert transform of :

The signal envelope is the positive signal

A phase-rotated seismic signal is
where is the phase of rotation.

By default, sfenvelope computes the signal envelope. It can also produce a phase-rotated signal if given hilb=y and phase=. If phase=90 (the default value), the phase-rotated signal will be simply the Hilbert transform of the input.

The following figure from book/rsf/rsf/sfenvelope illustrate an application of sfenvelope:


Computing the discrete Hilbert transform is not a trivial task. In the Fourier domain, the continuous Hilbert transform is given by where sgn is the sign function. The discontinuity of the sign function in the frequency domain at is related to the slow decay of the filter impulse response in the time domain. The discontinuity at the Nyquist frequency creates additional oscillations. Different practical implementations shorten the filter impulse response by effectively smoothing the Fourier-domain discontinuities.

The Madagascar implementation of the discrete Hilbert transform follows the algorithm described in

Pei, S.-C., and P.-H. Wang, 2001, Closed-form design of maximally flat FIR Hilbert transformers, differentiators, and fractional delayers by power series expansion: IEEE Trans. on Circuits and Systems, v. 48, No. 4, 389-398.

The accuracy/cost trade-off is controlled by two parameters: order= and ref=. The following figures frombook/rsf/rsf/sfenvelope illustrate the effect of the order= parameter:



The Seismic Unix implementation (suhilb program) applies a Hamming window in the time domain. For some reason, it has the filter polarity reversed:



A multidimensional analog of the Hilbert transform is the Riesz transform. It is implemented in the sfriesz program.

New users often ask for a brief comprehensible introduction to Madagascar.

Jeff Gowdin contributes A brief introduction to Madagascar, a tutorial book designed to provide such an introduction.

The purpose of this document is to teach new users how to use the powerful tools in Madagascar to: process data, produce documents and build your own Madagascar programs.

You can access the tutorial in the following forms:

• Wiki: http://www.ahay.org/wiki/Tutorial
• PDF: http://www.ahay.org/RSF/book/tutorial/book.pdf
• HTML chapters: http://www.ahay.org/RSF/book/tutorial/
• LaTeX Source: $RSFSRC/book/tutorial

You can endorse or discuss Science Code Manifesto published this week at http://sciencecodemanifesto.org/

Software is a cornerstone of science. Without software, twenty-first century science would be impossible. Without better software, science cannot progress.

But the culture and institutions of science have not yet adjusted to this reality. We need to reform them to address this challenge, by adopting these five principles:

Code

All source code written specifically to process data for a published paper must be available to the reviewers and readers of the paper.

Copyright

The copyright ownership and license of any released source code must be clearly stated.

Citation

Researchers who use or adapt science source code in their research must credit the code’s creators in resulting publications.

Credit

Software contributions must be included in systems of scientific assessment, credit, and recognition.

Curation

Source code must remain available, linked to related materials, for the useful lifetime of the publication.

Nick Barnes, the author of the Manifesto, explains its creation as follows:

I wrote it for the Climate Code Foundation, initially as a response and contribution to the Royal Society’s policy study on “Science as a Public Enterprise”. It is partly inspired by the Panton Principles, a bold statement of ideals in scientific data sharing. It refines the ideas I laid out in an opinion piece for Nature in 2010.

However, I did not originate these ideas. They are simply extensions of the core principle of science: publication. Publication is what distinguishes science from alchemy, and is what has propelled science – and human society – so far and so fast in the last 300 years. The Manifesto is the natural application of this principle to the relatively new, and increasingly important, area of science software.

My own ideals, influenced by the Free and Open Source Software movement, go beyond those stated in the Manifesto: I believe that Open Source publication of all science software will be one outcome of the current revolution in scientific methods, a revolution in which I hope this Manifesto will play a part.

sfagc implements Automatic Gain Control, a popular technique for normalizing signal amplitudes.

The algorithm is simple: divide the signal by its smoothed absolute value. The smoothing is controlled by rect#= and repeat= parameters, similar to the ones used by sfsmooth.

The following example from rsf/rsf/sfagc illustrates the application of sfagc in comparison with the application of sfpow, which applies a gain based on a power of time. The gains are applied to a shot gather from Alaska from the collection of shot gathers by Yilmaz and Cumro. A similar example appears on page 236 in Jon Claerbout's Imaging the Earth's Interior.



For a more accurate algorithm, try sfshapeagc, which computes the gain function using shaping regularization.

During the last four years, there have been nearly 170,000 visits to the Madagascar website (including 126 visits from the island of Madagascar). The top ten countries, as counted by Google Analytics, are: USA, China, Canada, UK, Brazil, Germany, Italy, Saudi Arabia, France, and India.



If the top ten are normalized by population (to compute visits per capita), they become: Canada, Saudi Arabia, USA, UK, Italy, France, Germany, Brazil, China, and India.

sfclip is very simple yet very useful program. It "clips" the data to the specified maximum by the absolute value.

Here is a simple test. First, let us make some data.

bash$ sfmath n1=10 output="sin(x1)" > data.rsf

bash$ < data.rsf sfdisfil

   0:             0       0.8415       0.9093       0.1411      -0.7568

   5:       -0.9589      -0.2794        0.657       0.9894       0.4121

Now clip it to 0.5 by maximum absolute value.

bash$ < data.rsf sfclip clip=0.5 > clip.rsf

bash$ < clip.rsf sfdisfil

   0:           0.5          0.5          0.5       0.1411         -0.5

   5:          -0.5      -0.2794          0.5          0.5       0.4121 

What if you need to clip the data not by the maximum value but to a specified range? Use sfclip2.

bash$ < data.rsf sfclip2 lower=0 upper=0.9 > clip2.rsf

bash$ < clip2.rsf sfdisfil

   0:             0       0.8415          0.9       0.1411            0

   5:             0            0        0.657          0.9       0.4121

sfclip should handle correctly infinite values, for example those resulting from division by zero.

bash$ sfmath n1=10 output=1/x1 > data.rsf

bash$ < data.rsf sfdisfil

   0:           inf            1          0.5       0.3333         0.25

   5:           0.2       0.1667       0.1429        0.125       0.1111

bash$ < data.rsf sfclip clip=0.3 > clip.rsf

bash$ < clip.rsf sfdisfil

   0:           0.3          0.3          0.3          0.3         0.25

   5:           0.2       0.1667       0.1429        0.125       0.1111

A prototype of sfclip is used as an example in Guide to madagascar API. The actual program is a little different.