Editing Guide to madagascar programs (section)

=Main programs=
The source files for these programs can be found under [https://github.com/ahay/src/blob/master/system/main system/main] in the Madagascar distribution. The "main" programs perform general-purpose operations on RSF hypercubes regardless of the data dimensionality or physical dimensions.

==sfadd==
{| class="wikitable" align="center" cellspacing="0" border="1"
! colspan="4" style="background:#ffdead;" | Add, multiply, or divide  RSF datasets.
|-
! colspan="4" | sfadd > out.rsf scale= add= sqrt= abs= log= exp= mode= [< file0.rsf] file1.rsf file2.rsf ...
|-
|  colspan="4" | The various operations, if selected, occur in the following order:<br><br>(1) Take absolute value, abs=<br>(2) Add a scalar, add=<br>(3) Take the natural logarithm, log=<br>(4) Take the square root, sqrt=<br>(5) Multiply by a scalar, scale=<br>(6) Compute the base-e exponential, exp=<br>(7) Add, multiply, or divide the data sets, mode=<br><br>sfadd operates on integer, float, or complex data, but all the input<br>and output files must be of the same data type.<br><br>An alternative to sfadd is sfmath, which is more versatile, but may be<br>less efficient.
|-
| ''bools  '' || '''abs=''' ||   || 	If true take absolute value  [nin]
|-
| ''floats '' || '''add=''' ||   || 	Scalar values to add to each dataset  [nin]
|-
| ''bools  '' || '''exp=''' ||   || 	If true compute exponential  [nin]
|-
| ''bools  '' || '''log=''' ||   || 	If true take logarithm  [nin]
|-
| ''string '' || '''mode=''' ||   || 	'a' means add (default), <br>       'p' or 'm' means multiply, <br>       'd' means divide
|-
| ''floats '' || '''scale=''' ||   || 	Scalar values to multiply each dataset with  [nin]
|-
| ''bools  '' || '''sqrt=''' ||   || 	If true take square root  [nin]
|}

<tt>sfadd</tt> is useful for combining (adding, dividing, or
multiplying) several datasets. What if you want to subtract two
datasets? Easy. Use the <tt>scale</tt> parameter as follows:
<pre>
bash&#36; sfadd data1.rsf data2.rsf scale=1,-1 > diff.rsf
</pre>
or
<pre>
bash&#36; sfadd < data1.rsf data2.rsf scale=1,-1 > diff.rsf
</pre>
The same task can be accomplished with the more general <tt>sfmath</tt> program:
<pre>
bash&#36; sfmath one=data1.rsf two=data2.rsf output='one-two' > diff.rsf
</pre>
or
<pre>
bash&#36; sfmath < data1.rsf two=data2.rsf output='input-two' > diff.rsf
</pre>
In both cases, the size and shape of <tt>data1.rsf</tt> and
<tt>data2.rsf</tt> hypercubes should be the same, and a warning message is printed out if the axis sampling parameters (such as
<tt>o1</tt> or <tt>d1</tt>) in these files are different.
====Implementation: [https://github.com/ahay/src/blob/master/system/main/add.c system/main/add.c]====
The first input file is either in the list or in the standard input.
<syntaxhighlight lang="c">
    /* find number of input files */
    if (isatty(fileno(stdin))) { 
        /* no input file in stdin */
	nin=0;
    } else {
        filename[0] = "in";
	nin=1;
    }
</syntaxhighlight>

Collect input files in the <tt>in</tt> array from all command-line
parameters that don't contain an "<tt>=</tt>" sign. The total number
of input files is <tt>nin</tt>.
<syntaxhighlight lang="c">
    for (i=1; i< argc; i++) { /* collect inputs */
	if (NULL != strchr(argv[i],'=')) continue; /* not a file */
	filename[nin] = argv[i];
	nin++;
    }
    if (0==nin) sf_error ("no input");
    /* nin = no of input files*/
</syntaxhighlight>

A helper function <tt>check_compat</tt> checks the compatibility of
input files.
<syntaxhighlight lang="c">
static void 
check_compat (sf_datatype type /* data type */, 
	      size_t      nin  /* number of files */, 
	      sf_file*    in   /* input files [nin] */, 
	      int         dim  /* file dimensionality */, 
	      const int*  n    /* dimensions [dim] */)
/* Check that the input files are compatible. 
   Issue error for type mismatch or size mismatch.
   Issue warning for grid parameters mismatch. */
{
    int ni, id;
    size_t i;
    float d, di, o, oi;
    char key[3];
    const float tol=1.e-5; /* tolerance for comparison */
    
    for (i=1; i < nin; i++) {
	if (sf_gettype(in[i]) != type) 
	    sf_error ("type mismatch: need %d",type);
	for (id=1; id <= dim; id++) {
	    (void) snprintf(key,3,"n%d",id);
	    if (!sf_histint(in[i],key,&ni) || ni != n[id-1])
		sf_error("%s mismatch: need %d",key,n[id-1]);
	    (void) snprintf(key,3,"d%d",id);
	    if (sf_histfloat(in[0],key,&d)) {
		if (!sf_histfloat(in[i],key,&di) || 
		    (fabsf(di-d) > tol*fabsf(d)))
		    sf_warning("%s mismatch: need %g",key,d);
	    } else {
		d = 1.;
	    }
	    (void) snprintf(key,3,"o%d",id);
	    if (sf_histfloat(in[0],key,&o) && 
		(!sf_histfloat(in[i],key,&oi) || 
		 (fabsf(oi-o) > tol*fabsf(d))))
		sf_warning("%s mismatch: need %g",key,o);
	}
    }
}
</syntaxhighlight>

Finally, we enter the main loop, where the input data are getting read
buffer by buffer and combined in the total product depending on the
data type.
<syntaxhighlight lang="c">
    for (nbuf /= sf_esize(in[0]); nsiz > 0; nsiz -= nbuf) {
	if (nbuf > nsiz) nbuf=nsiz;

	for (j=0; j < nin; j++) {
	    collect = (bool) (j != 0);
	    switch(type) {
		case SF_FLOAT:
		    sf_floatread((float*) bufi,
				 nbuf,
				 in[j]);	    
		    add_float(collect, 
			      nbuf,
			      (float*) buf,
			      (const float*) bufi, 
			      cmode, 
			      scale[j], 
			      add[j], 
			      abs_flag[j], 
			      log_flag[j], 
			      sqrt_flag[j], 
			      exp_flag[j]);
		    break;
</syntaxhighlight>

The data combination program for floating point numbers is
<tt>add_float</tt>.  
<syntaxhighlight lang="c">
static void add_float (bool   collect,    /* if collect */
		       size_t nbuf,       /* buffer size */
		       float* buf,        /* output [nbuf] */
		       const float* bufi, /* input  [nbuf] */  
		       char   cmode,      /* operation */
		       float  scale,      /* scale factor */
		       float  add,        /* add factor */
		       bool   abs_flag,   /* if abs */
		       bool   log_flag,   /* if log */
		       bool   sqrt_flag,  /* if sqrt */
		       bool   exp_flag    /* if exp */)
/* Add floating point numbers */
{
    size_t j;
    float f;

    for (j=0; j < nbuf; j++) {
	f = bufi[j];
	if (abs_flag)    f = fabsf(f);
	f += add;
	if (log_flag)    f = logf(f);
	if (sqrt_flag)   f = sqrtf(f);
	if (1. != scale) f *= scale;
	if (exp_flag)    f = expf(f);
	if (collect) {
	    switch (cmode) {
		case 'p': /* product */
		case 'm': /* multiply */
		    buf[j] *= f;
		    break;
		case 'd': /* delete */
		    if (f != 0.) buf[j] /= f;
		    break;
		default:  /* add */
		    buf[j] += f;
		    break;
	    }
	} else {
	    buf[j] = f;
	}
    }
}
</syntaxhighlight>

==sfattr==
{| class="wikitable" align="center" cellspacing="0" border="1"
! colspan="4" style="background:#ffdead;" | Display dataset attributes.
|-
! colspan="4" | sfattr < in.rsf lval=2 want=
|-
|  colspan="4" | <br>Sample output from "sfspike n1=100 | sfbandpass fhi=60 | sfattr"<br>*******************************************<br>rms =      0.992354<br>mean =      0.987576<br>2-norm =       9.92354<br>variance =    0.00955481<br>std dev =     0.0977487<br>max =       1.12735 at 97<br>min =      0.151392 at 100<br>nonzero samples = 100<br>total samples = 100<br>*******************************************<br><br>rms                = sqrt[ sum(data^2) / n ]<br>mean               = sum(data) / n<br>norm               = sum(abs(data)^lval)^(1/lval)<br>variance           = [ sum(data^2) - n*mean^2 ] / [ n-1 ]<br>standard deviation = sqrt [ variance ]
|-
| ''int    '' || '''lval=2''' ||   || 	norm option, lval is a non-negative integer, computes the vector lval-norm
|-
| ''string '' || '''want=''' ||   || 	'all'(default), 'rms', 'mean', 'norm', 'var', <br>       'std', 'max', 'min', 'nonzero', 'samples', 'short' 
:want=   'rms' displays the root mean square
:want=   'norm' displays the square norm, otherwise specified by lval.
:want=   'var' displays the variance
:want=   'std' displays the standard deviation
:want=   'nonzero' displays number of nonzero samples
:want=   'samples' displays total number of samples
:want=   'short' displays a short one-line version
|}


<tt>sfattr</tt> is a useful diagnostic program. It reports certain
statistical values for an RSF dataset: RMS (root-mean-square)
amplitude, mean value, vector norm value, variance, standard deviation,
maximum and minimum values, number of nonzero samples, and the total
number of samples.
If we denote data values as <math>d_i</math> for <math>i=0,1,2,\ldots,n</math>, then the RMS
value is <math>\sqrt{\frac{1}{n}\,\sum\limits_{i=0}^n d_i^2}</math>, the mean
value is <math>\frac{1}{n}\,\sum\limits_{i=0}^n d_i</math>, the <math>L_2</math>-norm value
is <math>\sqrt{\sum\limits_{i=0}^n d_i^2}</math>, the variance is
<math>\frac{1}{n-1}\,\left[\sum\limits_{i=0}^n d_i^2 - \frac{1}{n}\left(\sum\limits_{i=0}^n d_i\right)^2\right]</math>, and the standard
deviation is the square root of the variance. Using <tt>sfattr</tt> is a quick way to see the distribution of data values and check it for anomalies.

The output can be parsed using utilities such as <tt>awk</tt>, to extract only a numeric value for feeding it as a parameter value into a command line interface. Notice the backticks in the example below:
<syntaxhighlight lang="bash">
sfgrey <vel.rsf allpos=y bias=`sfattr <vel.rsf want=min | awk '{print $4}'` | sfpen
</syntaxhighlight>

====Implementation: [https://github.com/ahay/src/blob/master/system/main/attr.c system/main/attr.c]====

Computations start by finding the input data (<tt>in</tt>) size
(<tt>nsiz</tt>) and dimensions (<tt>dim</tt>).
<syntaxhighlight lang="c">
    dim = (size_t) sf_largefiledims (in,n);
    for (nsiz=1, i=0; i < dim; i++) {
	nsiz *= n[i];
    }
</syntaxhighlight>


In the main loop, we read the input data buffer by buffer.
<syntaxhighlight lang="c">
    for (nleft=nsiz; nleft > 0; nleft -= nbuf) {
	nbuf = (bufsiz < nleft)? bufsiz: nleft;
	switch (type) {
	    case SF_FLOAT: 
		sf_floatread((float*) buf,nbuf,in);
		break;
	    case SF_INT:
		sf_intread((int*) buf,nbuf,in);
		break;
            case SF_SHORT:
                sf_shortread((short*) buf,nbuf,in);
		break;
	    case SF_COMPLEX:
		sf_complexread((sf_complex*) buf,nbuf,in);
		break;
	    case SF_UCHAR:
		sf_ucharread((unsigned char*) buf,nbuf,in);
		break;
	    case SF_CHAR:
	    default:
		sf_charread(buf,nbuf,in);
		break;
	}
</syntaxhighlight>


The data attributes are accumulated in corresponding double-precision variables. 
<syntaxhighlight lang="c">
	    fsum += f;
	    fsqr += (double) f*f;
</syntaxhighlight>


Finally, the attributes are reduced and printed out.
<syntaxhighlight lang="c">
    fmean = fsum/nsiz;
    if (lval==2)      fnorm = sqrt(fsqr);
    else if (lval==0) fnorm = nsiz-nzero;
    else              fnorm = pow(flval,1./lval);
    frms = sqrt(fsqr/nsiz);
    if (nsiz > 1) fvar = fabs(fsqr-nsiz*fmean*fmean)/(nsiz-1);
    else          fvar = 0.0;
    fstd = sqrt(fvar);
</syntaxhighlight>

<syntaxhighlight lang="c">
    if(NULL==want || 0==strcmp(want,"rms"))
	printf("rms = %13.6g \n",(float) frms);
    if(NULL==want || 0==strcmp(want,"mean"))
	printf("mean = %13.6g \n",(float) fmean);
    if(NULL==want || 0==strcmp(want,"norm"))
	printf("%d-norm value = %13.6g \n",lval,(float) fnorm);
    if(NULL==want || 0==strcmp(want,"var"))
	printf("variance = %13.6g \n",(float) fvar);
    if(NULL==want || 0==strcmp(want,"std"))
	printf("standard deviation = %13.6g \n",(float) fstd);
</syntaxhighlight>

==sfcat==
{| class="wikitable" align="center" cellspacing="0" border="1"
! colspan="4" style="background:#ffdead;" | Concatenate datasets. 
|-
! colspan="4" | sfcat > out.rsf order= space= axis=3 nspace=(int) (ni/(20*nin) + 1) o= d= [<file0.rsf] file1.rsf file2.rsf ... 
|-
|  colspan="4" | sfmerge inserts additional space between merged data.
|-
| ''int    '' || '''axis=3''' ||   || 	Axis being merged
|-
| ''float  '' || '''d=''' ||   || 	axis sampling
|-
| ''int    '' || '''nspace=(int) (ni/(20*nin) + 1)''' ||   || 	if space=y, number of traces to insert
|-
| ''float  '' || '''o=''' ||   || 	axis origin
|-
| ''ints   '' || '''order=''' ||   || 	concatenation order  [nin]
|-
| ''bool   '' || '''space=''' ||  [y/n] || 	Insert additional space.
:y is default for sfmerge, n is default for sfcat
|}

<tt>sfcat</tt> and <tt>sfmerge</tt> concatenate two or more files
together along a particular axis. It is the same program, only
<tt>sfcat</tt> has the default <tt>space=n</tt> and <tt>sfmerge</tt>
has the default <tt>space=y</tt>.
Example of <tt>sfcat</tt>:
<pre>
bash$ sfspike n1=2 n2=3 > one.rsf
bash$ sfin one.rsf
one.rsf:
    in="/tmp/one.rsf@"
    esize=4 type=float form=native
    n1=2           d1=0.004       o1=0          label1="Time" unit1="s"
    n2=3           d2=0.1         o2=0          label2="Distance" unit2="km"
        6 elements 24 bytes
bash$ sfcat one.rsf one.rsf axis=1 > two.rsf
bash$ sfin two.rsf
two.rsf:
    in="/tmp/two.rsf@"
    esize=4 type=float form=native
    n1=4           d1=0.004       o1=0          label1="Time" unit1="s"
    n2=3           d2=0.1         o2=0          label2="Distance" unit2="km"
        12 elements 48 bytes
</pre>
Example of <tt>sfmerge</tt>:
<pre>
bash$ sfmerge one.rsf one.rsf axis=2 > two.rsf
bash$ sfin two.rsf
two.rsf:
    in="/tmp/two.rsf@"
    esize=4 type=float form=native
    n1=2           d1=0.004       o1=0          label1="Time" unit1="s"
    n2=7           d2=0.1         o2=0          label2="Distance" unit2="km"
        14 elements 56 bytes
</pre>
In this case, an extra empty trace is inserted between the two merged files.
The axes that are not being merged are checked for consistency:
<pre>
bash$ sfcat one.rsf two.rsf > three.rsf
sfcat: n2 mismatch: need 3
</pre>

====Implementation: [https://github.com/ahay/src/blob/master/system/main/cat.c system/main/cat.c]====
The first input file is either in the list or in the standard input.
<syntaxhighlight lang="c">
    in = (sf_file*) sf_alloc ((size_t) argc,sizeof(sf_file));

    if (!sf_stdin()) { /* no input file in stdin */
	nin=0;
    } else {
	filename[0] = "in";
	nin=1;
    }
</syntaxhighlight>

Everything on the command line that does not contain a "=" sign is
treated as a file name, and the corresponding file object is added to the list.
<syntaxhighlight lang="c">
    for (i=1; i< argc; i++) { /* collect inputs */
	if (NULL != strchr(argv[i],'=')) 
	    continue; /* not a file */
	filename[nin] = argv[i];
	nin++;
    }
    if (0==nin) sf_error ("no input");
</syntaxhighlight>

As explained above, if the <tt>space=</tt> parameter is not set, it is
inferred from the program name: <tt>sfmerge</tt> corresponds to
<tt>space=y</tt> and <tt>sfcat</tt> corresponds to <tt>space=n</tt>.
<syntaxhighlight lang="c">
    if (!sf_getbool("space",&space)) {
	/* Insert additional space.
	   y is default for sfmerge, n is default for sfcat */
	prog = sf_getprog();
	if (NULL != strstr (prog, "merge")) {
	    space = true;
	} else if (NULL != strstr (prog, "cat")) {
	    space = false;
	} else {
	    sf_warning("%s is neither merge nor cat,"
		       " assume merge",prog);
	    space = true;
	}
    }
</syntaxhighlight>

Find the axis for the merging (from the command line <tt>axis=</tt>
argument) and figure out two sizes: <tt>n1</tt> for everything after
the axis and <tt>n2</tt> for everything before the axis.
<syntaxhighlight lang="c">
    n1=1;
    n2=1;
    for (i=1; i <= dim; i++) {
	if      (i < axis) n1 *= n[i-1];
	else if (i > axis) n2 *= n[i-1];
    }
</syntaxhighlight>

In the output, the selected axis will get extended.
<syntaxhighlight lang="c">
    /* figure out the length of extended axis */
    ni = 0;
    for (j=0; j < nin; j++) {
	ni += naxis[j];
    }

    if (space) {
	if (!sf_getint("nspace",&nspace)) 
	    nspace = (int) (ni/(20*nin) + 1);
	/* if space=y, number of traces to insert */ 
	ni += nspace*(nin-1);
    } 

    (void) snprintf(key,3,"n%d",axis);
    sf_putint(out,key,(int) ni);
</syntaxhighlight>

The rest is simple: loop through the datasets reading and writing the
data in buffer-size chunks and adding extra empty chunks if
<tt>space=y</tt>.
<syntaxhighlight lang="c">
    for (i2=0; i2 < n2; i2++) {
	for (j=0; j < nin; j++) {
	    for (ni = n1*naxis[j]*esize; ni > 0; ni -= nbuf) {
		nbuf = (BUFSIZ < ni)? BUFSIZ: ni;
		sf_charread (buf,nbuf,in[j]);
		sf_charwrite (buf,nbuf,out);
	    }
	    if (!space || j == nin-1) continue;
	    /* Add spaces */
	    memset(buf,0,BUFSIZ);
	    for (ni = n1*nspace*esize; ni > 0; ni -= nbuf) {
		nbuf = (BUFSIZ < ni)? BUFSIZ: ni;
		sf_charwrite (buf,nbuf,out);
	    }
	}
    }
</syntaxhighlight>

==sfcdottest==
{| class="wikitable" align="center" cellspacing="0" border="1"
! colspan="4" style="background:#ffdead;" | Generic dot-product test for complex linear operators with adjoints 
|-
! colspan="4" | sfcdottest mod=mod.rsf dat=dat.rsf > pip.rsf
|-
| ''file   '' || '''dat=''' ||   || 	auxiliary input file name
|-
| ''file   '' || '''mod=''' ||   || 	auxiliary input file name
|}

A simple demonstration of the program can be made by taking advantage that the complex-to-complex FFT is a linear operator:

<syntaxhighlight lang="bash">
sfspike n1=100 | sfrtoc > spike.rsf
< spike.rsf sffft axis=1 pad=1 > spike2.rsf
sfcdottest sffft mod=spike.rsf dat=spike2.rsf axis=1 pad=1
</syntaxhighlight>

The output should show values identical down to the last decimal:

<pre>
sfcdottest:  L[m]*d=(3.73955,-1.86955)
sfcdottest: L'[d]*m=(3.73955,-1.86955)
</pre>

==sfcmplx==
{| class="wikitable" align="center" cellspacing="0" border="1"
! colspan="4" style="background:#ffdead;" | Create a complex dataset from its real and imaginary parts.
|-
! colspan="4" | sfcmplx < real.rsf > cmplx.rsf real.rsf imag.rsf
|-
|  colspan="4" | There has to be only two input files specified and no additional parameters.
|}


<tt>sfcmplx</tt> simply creates a complex dataset from its real and
imaginary parts. The reverse operation can be accomplished with
<tt>sfreal</tt> and <tt>sfimag</tt>.
Example of <tt>sfcmplx</tt>:
<pre>
bash$ sfspike n1=2 n2=3 > one.rsf
bash$ sfin one.rsf
one.rsf:
    in="/tmp/one.rsf@"
    esize=4 type=float form=native
    n1=2           d1=0.004       o1=0          label1="Time" unit1="s"
    n2=3           d2=0.1         o2=0          label2="Distance" unit2="km"
        6 elements 24 bytes
bash$ sfcmplx one.rsf one.rsf > cmplx.rsf
bash$ sfin cmplx.rsf
cmplx.rsf:
    in="/tmp/cmplx.rsf@"
    esize=8 type=complex form=native
    n1=2           d1=0.004       o1=0          label1="Time" unit1="s"
    n2=3           d2=0.1         o2=0          label2="Distance" unit2="km"
        6 elements 48 bytes
</pre>

====Implementation: [https://github.com/ahay/src/blob/master/system/main/cmplx.c system/main/cmplx.c]====
The program flow is simple. First, get the names of the input files.
<syntaxhighlight lang="c">
    /* the first two non-parameters are real and imag files */
    for (i=1; i< argc; i++) { 
	if (NULL == strchr(argv[i],'=')) {
	    if (NULL == real) {
		real = sf_input (argv[i]);
	    } else {
		imag = sf_input (argv[i]);
		break;
	    }
	}
    }
    if (NULL == imag) {
	if (NULL == real) sf_error ("not enough input");
	/* if only one input, real is in stdin */
	imag = real;
	real = sf_input("in");
    }
</syntaxhighlight>

The main part of the program reads the real and imaginary parts buffer by buffer and assembles and writes out the complex input. 
<syntaxhighlight lang="c">
    for (nleft= (size_t) (rsize*resize); nleft > 0; nleft -= nbuf) {
	nbuf = (BUFSIZ < nleft)? BUFSIZ: nleft;
	sf_charread(rbuf,nbuf,real);
	sf_charread(ibuf,nbuf,imag);
	for (i=0; i < nbuf; i += resize) {
	    memcpy(cbuf+2*i,       rbuf+i,(size_t) resize);
	    memcpy(cbuf+2*i+resize,ibuf+i,(size_t) resize);
	}
	sf_charwrite(cbuf,2*nbuf,cmplx);
    }
</syntaxhighlight>

==sfconjgrad==
{| class="wikitable" align="center" cellspacing="0" border="1"
! colspan="4" style="background:#ffdead;" | Generic conjugate-gradient solver for linear inversion 
|-
! colspan="4" | sfconjgrad < dat.rsf mod=mod.rsf mwt=mwt.rsf known=known.rsf x0=x0.rsf > to.rsf < from.rsf > out.rsf niter=1
|-
| ''string '' || '''known=''' ||   || 	auxiliary input file name
|-
| ''file   '' || '''mod=''' ||   || 	auxiliary input file name
|-
| ''string '' || '''mwt=''' ||   || 	auxiliary input file name
|-
| ''int    '' || '''niter=1''' ||   || 	number of iterations
|-
| ''string '' || '''x0=''' ||   || 	auxiliary input file name
|}

<tt>sfconjgrad</tt> is a generic program for least-squares linear
inversion with the [http://en.wikipedia.org/wiki/Conjugate_gradient_method conjugate-gradient method]. Suppose you have an
executable program <tt><prog></tt> that takes an RSF file from the
standard input and produces an RSF file in the standard output. It may
take any number of additional parameters but one of them must be
<tt>adj=</tt> that sets the forward (<tt>adj=0</tt>) or adjoint
(<tt>adj=1</tt>) operations.  The program <tt><prog></tt> is typically
an RSF program but it could be anything (a script, a multiprocessor
MPI program, etc.) as long as it implements a linear operator
<math>\mathbf{L}</math> and its adjoint. There are no restrictions on the data
size or shape. You can easily test the adjointness with
<tt>sfdottest</tt>. The <tt>sfconjgrad</tt> program searches for a
vector <math>\mathbf{m}</math> that minimizes the least-square misfit 
<math>\|\mathbf{d - L\,m}\|^2</math> for the given input data vector <math>\mathbf{d}</math>. 

The pseudocode for <tt>sfconjgrad</tt> is given at the end of the [https://ahay.org/RSF/book/gee/lsq/paper.pdf "Model fitting with least squares" chapter] of ''Imaging Estimation by Example'' by Jon Claerbout, with the earliest form published in [http://sepwww.stanford.edu/data/media/public/oldreports/sep48/48_25.pdf "Conjugate Gradient Tutorial"] (SEP-48, 1986, same author). A simple toy implementation with a small matrix shows that this algorithm produces the same steps as the algorithm described in equations 45-49 of [http://www.cs.cmu.edu/~quake-papers/painless-conjugate-gradient.pdf "An introduction to the Conjugate Gradient Method Without the Agonizing Pain"] by J.R. Shewchuk, 1994, when the equation <math>A^T A x = A^T b</math> (in Shewchuk's notation) is solved. Multiplying with the transpose ensures a correct solution even when matrix A is square but not symmetric at all. The program [https://ahay.org/RSF/sfcconjgrad.html sfcconjgrad] implements this algorithm for the case when inputs are complex.

Here is an example. The <tt>sfhelicon</tt> program implements Claerbout's multidimensional helical filtering
(Claerbout, 1998<ref>Claerbout, J.,  1998, Multidimensional recursive filters via a helix:  Geophysics, '''63''', 1532--1541.</ref>). It requires a filter to be specified in addition to the input and output vectors. We create a helical  2-D filter using the Unix <tt>echo</tt> command.
<pre>
bash$ echo 1 19 20 n1=3 n=20,20 data_format=ascii_int in=lag.rsf > lag.rsf
bash$ echo 1 1 1 a0=-3 n1=3 data_format=ascii_float in=flt.rsf > flt.rsf
</pre>
Next, we create an example 2-D model and data vector with <tt>sfspike</tt>.
<pre>
bash$ sfspike n1=50 n2=50 > vec.rsf
</pre>
The <tt>sfdottest</tt> program can perform the dot product test to
check that the adjoint mode works correctly.
<pre>
bash$ sfdottest sfhelicon filt=flt.rsf lag=lag.rsf mod=vec.rsf dat=vec.rsf
sfdottest:  L[m]*d=5.28394
sfdottest: L'[d]*m=5.28394
</pre>
Your numbers may be different because <tt>sfdottest</tt> generates new
random input on each run.
Next, let us make some random data with <tt>sfnoise</tt>.
<pre>
bash$ sfnoise seed=2005 rep=y < vec.rsf > dat.rsf
</pre>
and try to invert the filtering operation using <tt>sfconjgrad</tt>:
<pre>
bash$ sfconjgrad sfhelicon filt=flt.rsf lag=lag.rsf mod=vec.rsf < dat.rsf > mod.rsf niter=10
sfconjgrad: iter 1 of 10
sfconjgrad: grad=3253.65
sfconjgrad: iter 2 of 10
sfconjgrad: grad=289.421
sfconjgrad: iter 3 of 10
sfconjgrad: grad=92.3481
sfconjgrad: iter 4 of 10
sfconjgrad: grad=36.9417
sfconjgrad: iter 5 of 10
sfconjgrad: grad=18.7228
sfconjgrad: iter 6 of 10
sfconjgrad: grad=11.1794
sfconjgrad: iter 7 of 10
sfconjgrad: grad=7.26941
sfconjgrad: iter 8 of 10
sfconjgrad: grad=5.15945
sfconjgrad: iter 9 of 10
sfconjgrad: grad=4.23055
sfconjgrad: iter 10 of 10
sfconjgrad: grad=3.57495
</pre>
The output shows that, in 10 iterations, the norm of the gradient vector decreases by almost 1000. 
We can check the residual misfit before
<pre>
bash$ < dat.rsf sfattr want=norm
norm value = 49.7801
</pre>
and after
<pre>
bash$ sfhelicon filt=flt.rsf lag=lag.rsf < mod.rsf | sfadd scale=1,-1 dat.rsf | sfattr want=norm
norm value = 5.73563
</pre>
The misfit decreased by an order of magnitude in 10 iterations. Running the program for more iterations can improve the result.

An equivalent implementation for complex-valued inputs is [https://ahay.org/RSF/sfcconjgrad.html sfcconjgrad]. A lightweight Python implementation can be found in [https://github.com/ahay/src/blob/master/user/fomels/conjgrad.py $PYTHONPATH/rsf/conjgrad.py].

==sfcp==
{| class="wikitable" align="center" cellspacing="0" border="1"
! colspan="4" style="background:#ffdead;" | Copy or move a dataset.
|-
! colspan="4" | sfcp < in.rsf > out.rsf in.rsf out.rsf
|-
|  colspan="4" | sfcp - copy, sfmv - move.<br>Mimics standard Unix commands.
|}


The <tt>sfcp</tt> and <tt>sfmv</tt> command imitate the Unix
<tt>cp</tt> and <tt>mv</tt> commands and serve for copying and moving
RSF files. Example:
<pre>
bash$ sfspike n1=2 n2=3 > one.rsf
bash$ sfin one.rsf
one.rsf:
    in="/tmp/one.rsf@"
    esize=4 type=float form=native
    n1=2           d1=0.004       o1=0          label1="Time" unit1="s"
    n2=3           d2=0.1         o2=0          label2="Distance" unit2="km"
        6 elements 24 bytes
bash$ sfcp one.rsf two.rsf
bash$ sfin two.rsf
two.rsf:
    in="/tmp/two.rsf@"
    esize=4 type=float form=native
    n1=2           d1=0.004       o1=0          label1="Time" unit1="s"
    n2=3           d2=0.1         o2=0          label2="Distance" unit2="km"
        6 elements 24 bytes
</pre>

====Implementation: [https://github.com/ahay/src/blob/master/system/main/cp.c system/main/cp.c]====
First, we look for the two first command-line arguments that don't
have the "=" character in them and consider them as the names of the input and the output files.  
<syntaxhighlight lang="c">
    /* the first two non-parameters are in and out files */
    for (i=1; i< argc; i++) { 
	if (NULL == strchr(argv[i],'=')) {
	    if (NULL == in) {
		infile = argv[i];
		in = sf_input (infile);
	    } else {
		out = sf_output (argv[i]);
		break;
	    }
	}
    }
    if (NULL == in || NULL == out)
	sf_error ("not enough input");
</syntaxhighlight>

Next, we use library functions <font color="#cd4b19">sf_cp</font> and <font color="#cd4b19">sf_rm</font> to do the actual work.
<syntaxhighlight lang="c">
    sf_cp(in,out);
    if (NULL != strstr (prog,"mv")) 
	sf_rm(infile,false,false,false);
</syntaxhighlight>

==sfcut==
{| class="wikitable" align="center" cellspacing="0" border="1"
! colspan="4" style="background:#ffdead;" | Zero a portion of the dataset.
|-
! colspan="4" | sfcut < in.rsf > out.rsf verb=n j#=(1,...) d#=(d1,d2,...) f#=(0,...) min#=(o1,o2,,...) n#=(0,...) max#=(o1+(n1-1)*d1,o2+(n1-1)*d2,,...)
|-
|  colspan="4" | <br>Reverse of window. 
|-
| ''float  '' || '''d#=(d1,d2,...)''' ||   || 	sampling in #-th dimension
|-
| ''largeint'' || '''f#=(0,...)''' ||   || 	window start in #-th dimension
|-
| ''int    '' || '''j#=(1,...)''' ||   || 	jump in #-th dimension
|-
| ''float  '' || '''max#=(o1+(n1-1)*d1,o2+(n1-1)*d2,,...)''' ||   || 	maximum in #-th dimension
|-
| ''float  '' || '''min#=(o1,o2,,...)''' ||   || 	minimum in #-th dimension
|-
| ''int    '' || '''n#=(0,...)''' ||   || 	window size in #-th dimension
|-
| ''bool   '' || '''verb=n''' ||  [y/n] || 	Verbosity flag
|}

The <tt>sfcut</tt> command is related to <tt>sfwindow</tt> and has the same
set of arguments, only instead of extracting the selected window, it fills it
with zeroes. The size of the input data is preserved. 
Examples:
<pre>
bash$ sfspike n1=5 n2=5 > in.rsf
bash$ < in.rsf sfdisfil
   0:             1            1            1            1            1
   5:             1            1            1            1            1
  10:             1            1            1            1            1
  15:             1            1            1            1            1
  20:             1            1            1            1            1
bash$ < in.rsf sfcut n1=2 f1=1 n2=3 f2=2 | sfdisfil
   0:             1            1            1            1            1
   5:             1            1            1            1            1
  10:             1            0            0            1            1
  15:             1            0            0            1            1
  20:             1            0            0            1            1
bash$ < in.rsf sfcut j1=2 | sfdisfil
   0:             0            1            0            1            0
   5:             0            1            0            1            0
  10:             0            1            0            1            0
  15:             0            1            0            1            0
  20:             0            1            0            1            0
</pre>

==sfdd==
{| class="wikitable" align="center" cellspacing="0" border="1"
! colspan="4" style="background:#ffdead;" | Convert between different formats. 
|-
! colspan="4" | sfdd < in.rsf > out.rsf trunc=n line=8 ibm=n form= type= format=
|-
| ''string '' || '''form=''' ||   || 	ascii, native, xdr
|-
| ''string '' || '''format=''' ||   || 	Element format (for conversion to ASCII)
|-
| ''bool   '' || '''ibm=n''' ||  [y/n] || 	Special case - assume integers actually represent IBM floats
|-
| ''int    '' || '''line=8''' ||   || 	Number of numbers per line (for conversion to ASCII)
|-
| ''bool   '' || '''trunc=n''' ||  [y/n] || 	Truncate or round to nearest when converting from float to int/short
|-
| ''string '' || '''type=''' ||   || 	int, float, complex, short, long
|}


The <tt>sfdd</tt> program is used to change either the form (<tt>ascii</tt>,
<tt>xdr</tt>, <tt>native</tt>) or the type (<tt>complex</tt>, <tt>float</tt>,
<tt>int</tt>, <tt>char</tt>) of the input dataset. 
In the example below, we create a plain text (ASCII) file with numbers and
then use <tt>sfdd</tt> to generate an RSF file in <tt>xdr</tt> form with
<tt>complex</tt> numbers. 
<pre>
bash$ cat test.txt
1 2 3 4 5 6
bash$ echo n1=6 data_format=ascii_int in=test.txt > test.rsf
bash$ sfin test.rsf
test.rsf:
    in="test.txt"
    esize=0 type=int form=ascii
    n1=6           d1=?           o1=?
        6 elements
bash$ sfdd < test.rsf form=xdr type=complex > test2.rsf
bash$ sfin test2.rsf
test2.rsf:
    in="/tmp/test2.rsf@"
    esize=8 type=complex form=xdr
    n1=3           d1=?           o1=?
        3 elements 24 bytes
bash$ sfdisfil < test2.rsf
   0:          1,         2i         3,         4i         5,         6i
</pre>
To learn more about the RSF data format, consult the [[Guide to RSF file format| guide to RSF format]].

==sfdisfil==
{| class="wikitable" align="center" cellspacing="0" border="1"
! colspan="4" style="background:#ffdead;" | Print out data values.
|-
! colspan="4" | sfdisfil < in.rsf number=y col=0 format= header= trailer=
|-
|  colspan="4" | <br>Alternatively, use sfdd and convert to ASCII form.
|-
| ''int    '' || '''col=0''' ||   || 	Number of columns.
:The default depends on the data type:
:10 for int and char,
:5 for float,
:3 for complex
|-
| ''string '' || '''format=''' ||   || 	Format for numbers (printf-style).
:The default depends on the data type:<br>       "%4d " for int and char,<br>       "%13.4g" for float,<br>       "%10.4g,%10.4gi" for complex
|-
| ''string '' || '''header=''' ||   || 	Optional header string to output before data
|-
| ''bool   '' || '''number=y''' ||  [y/n] || 	If number the elements
|-
| ''string '' || '''trailer=''' ||   || 	Optional trailer string to output after data
|}


The <tt>sfdisfil</tt> program simply dumps the data contents to the standard
output in a text form. It is used mostly for debugging purposes to examine RSF files quickly. Here is an example:
<pre>
bash$ sfmath o1=0 d1=2 n1=12 output=x1 > test.rsf
bash$ < test.rsf sfdisfil
   0:             0            2            4            6            8
   5:            10           12           14           16           18
  10:            20           22
</pre>
The output format is easily configurable.
<pre>
bash$ < test.rsf sfdisfil col=6 number=n format="%5.1f"
  0.0  2.0  4.0  6.0  8.0 10.0
 12.0 14.0 16.0 18.0 20.0 22.0
</pre>
Along with <tt>sfdd</tt>, <tt>sfdisfil</tt> provides a simple way to convert
RSF data to an ASCII form.

==sfdottest==
{| class="wikitable" align="center" cellspacing="0" border="1"
! colspan="4" style="background:#ffdead;" | Generic dot-product test for linear operators with adjoints 
|-
! colspan="4" | sfdottest mod=mod.rsf dat=dat.rsf > pip.rsf
|-
| ''file   '' || '''dat=''' ||   || 	auxiliary input file name
|-
| ''file   '' || '''mod=''' ||   || 	auxiliary input file name
|}

<tt>sfdottest</tt> is a generic dot-product test program for testing
linear operators. Suppose there is an executable program
<tt><prog></tt> that takes an RSF file from the standard input and
produces an RSF file in the standard output. It may take any number of
additional parameters but one of them must be <tt>adj=</tt> that sets
the forward (<tt>adj=0</tt>) or adjoint (<tt>adj=1</tt>) operations.
The program <tt><prog></tt> is typically an RSF program but it could
be anything (a script, a multiprocessor MPI program, etc.) as long as
it implements a linear operator <math>\mathbf{L}</math> and its adjoint
<math>\mathbf{L}^T</math>. The <tt>sfdottest</tt> program is testing the equality
<center><math>
d^T\,L\,m = m^T\,L^T\,d
</math></center>
by using random vectors <math>\mathbf{m}</math> and <math>\mathbf{d}</math>. You can invoke it with
<pre>
bash$ sfdottest <prog> [optional aruments] mod=mod.rsf dat=dat.rsf
</pre>
where <tt>mod.rsf</tt> and <tt>dat.rsf</tt> are RSF files that
represent vectors from the model and data spaces. Pay attention to the 
dimensions and size of these vectors! If the program does not respond
for a very long time, the dimension and size
of the vectors may be inconsistent with the requirement of the program to be
tested. <tt>sfdottest</tt>
does not create any temporary files and has no restrictive
limitations on the size of the vectors.
Here is an example. We first create a vector with 100 elements using
<tt>sfspike</tt> and then run <tt>sfdottest</tt> to test the
<tt>sfcausint</tt> program. <tt>sfcausint</tt> implements a linear operator of causal integration and its adjoint, the anti-causal
integration.
<pre>
bash$ sfspike n1=100 > vec.rsf
bash$ sfdottest sfcausint mod=vec.rsf dat=vec.rsf
sfdottest:  L[m]*d=1410.2
sfdottest: L'[d]*m=1410.2
bash$ sfdottest sfcausint mod=vec.rsf dat=vec.rsf
sfdottest:  L[m]*d=1165.87
sfdottest: L'[d]*m=1165.87
</pre>
The numbers are different on subsequent runs because of changing
seed in the random number generator.
Here is a somewhat more complicated example. The <tt>sfhelicon</tt>
program implements Claerbout's multidimensional helical filtering
(Claerbout, 1998<ref>Claerbout, J.,  1998, Multidimensional recursive filters via a helix:  Geophysics, '''63''', 1532--1541.</ref>). It requires a filter to be specified in
addition to the input and output vectors. We create a helical 
2-D filter using the Unix <tt>echo</tt> command.
<pre>
bash$ echo 1 19 20 n1=3 n=20,20 data_format=ascii_int in=lag.rsf > lag.rsf
bash$ echo 1 1 1 a0=-3 n1=3 data_format=ascii_float in=flt.rsf > flt.rsf
</pre>
Next, we create an example 2-D model and data vector with <tt>sfspike</tt>.
<pre>
bash$ sfspike n1=50 n2=50 > vec.rsf
</pre>
Now the <tt>sfdottest</tt> program can perform the dot product test.
<pre>
bash$ sfdottest sfhelicon filt=flt.rsf lag=lag.rsf \
> mod=vec.rsf dat=vec.rsf
sfdottest:  L[m]*d=8.97375
sfdottest: L'[d]*m=8.97375
</pre>
Here is the same program tested in the inverse filtering mode:
<pre>
bash$ sfdottest sfhelicon filt=flt.rsf lag=lag.rsf \
> mod=vec.rsf dat=vec.rsf div=y
sfdottest:  L[m]*d=15.0222
sfdottest: L'[d]*m=15.0222
</pre>

==sfget==
{| class="wikitable" align="center" cellspacing="0" border="1"
! colspan="4" style="background:#ffdead;" | Output parameters from the header.
|-
! colspan="4" | sfget parform=y all=n par1 par2 ...
|-
| ''bool   '' || '''all=n''' ||  [y/n] || 	If output all values.
|-
| ''bool   '' || '''parform=y''' ||  [y/n] || 	If y, print out parameter=value. If n, print out value.
|}

The <tt>sfget</tt> program extracts a parameter value from an RSF file. It is
useful mostly for scripting. Here is, for example, a quick calculation of the
maximum value on the first axis in an RSF dataset (the output of
<tt>sfspike</tt>) using the standard Unix <tt>bc</tt> calculator.
<pre>
bash$ ( sfspike n1=100 | sfget n1 d1 o1; echo "o1+(n1-1)*d1" ) | bc
.396
</pre>
See also <tt>sfput</tt>.

====Implementation: [https://github.com/ahay/src/blob/master/system/main/get.c system/main/get.c]====
The implementation is trivial. Loop through all command-line parameters that contain 
the "=" character.
<syntaxhighlight lang="c">
    for (i = 1; i < argc; i++) {
	key = argv[i];
	if (NULL != strchr(key,'=')) continue;
</syntaxhighlight>

Get the parameter value (as a string) and output it as either
<tt>key=value</tt> or <tt>value</tt>, depending on the
<tt>parform</tt> parameter.  
<syntaxhighlight lang="c">
	string = sf_histstring(in,key);
	if (NULL == string) {
	    sf_warning("No key %s",key);
	} else {
	    if (parform) printf ("%s=",key);
	    printf("%s\n",string);
	}
</syntaxhighlight>

==sfheadercut==
{| class="wikitable" align="center" cellspacing="0" border="1"
! colspan="4" style="background:#ffdead;" | Zero a portion of a dataset based on a header mask.
|-
! colspan="4" | sfheadercut mask=head.rsf < in.rsf > out.rsf
|-
|  colspan="4" | <br>The input data is a collection of traces n1xn2,<br>mask is an integer array of size n2.
|-
| ''file   '' || '''mask=''' ||   || 	auxiliary input file name
|}

<tt>sfheadercut</tt> is close to <tt>sfheaderwindow</tt> but instead
of windowing the dataset, it fills the traces specified by the header
mask with zeroes. The size of the input data is preserved.
Here is an example of using <tt>sfheaderwindow</tt> for 
zeroing every other trace in the input file. First, let us create
an input file with ten traces:
<pre>
bash$ sfmath n1=5 n2=10 output=x2+1 > input.rsf
bash$ < input.rsf sfdisfil
   0:             1            1            1            1            1
   5:             2            2            2            2            2
  10:             3            3            3            3            3
  15:             4            4            4            4            4
  20:             5            5            5            5            5
  25:             6            6            6            6            6
  30:             7            7            7            7            7
  35:             8            8            8            8            8
  40:             9            9            9            9            9
  45:            10           10           10           10           10
</pre>
Next, we can create a mask with alternating ones and zeros using
<tt>sfinterleave</tt>.
<pre>
bash$ sfspike n1=5 mag=1 | sfdd type=int > ones.rsf
bash$ sfspike n1=5 mag=0 | sfdd type=int > zeros.rsf
bash$ sfinterleave axis=1 ones.rsf zeros.rsf > mask.rsf
bash$ sfdisfil < mask.rsf
   0:    1    0    1    0    1    0    1    0    1    0
</pre>
Finally, <tt>sfheadercut</tt> zeros the input traces.
<pre>
bash$ sfheadercut < input.rsf mask=mask.rsf > output.rsf
bash$ sfdisfil < output.rsf 
   0:             1            1            1            1            1
   5:             0            0            0            0            0
  10:             3            3            3            3            3
  15:             0            0            0            0            0
  20:             5            5            5            5            5
  25:             0            0            0            0            0
  30:             7            7            7            7            7
  35:             0            0            0            0            0
  40:             9            9            9            9            9
  45:             0            0            0            0            0
</pre>

==sfheadersort==
{| class="wikitable" align="center" cellspacing="0" border="1"
! colspan="4" style="background:#ffdead;" | Sort a dataset according to a header key. 
|-
! colspan="4" | sfheadersort < in.rsf > out.rsf head=
|-
| ''string '' || '''head=''' ||   || 	header file
|}


<tt>sfheadersort</tt> is used to sort traces in the input file
according to trace header information. 
Here is an example of using
<tt>sfheadersort</tt> for randomly shuffling traces in the input
file. First, let us create an input file with seven traces:
<pre>
bash$ sfmath n1=5 n2=7 output=x2+1 > input.rsf
bash$ < input.rsf sfdisfil
   0:             1            1            1            1            1
   5:             2            2            2            2            2
  10:             3            3            3            3            3
  15:             4            4            4            4            4
  20:             5            5            5            5            5
  25:             6            6            6            6            6
  30:             7            7            7            7            7 
</pre>
Next, we can create a random file with seven header values using
<tt>sfnoise</tt>.
<pre>
bash$ sfspike n1=7 | sfnoise rep=y type=n > random.rsf
bash$ < random.rsf sfdisfil
   0:       0.05256      -0.2879       0.1487       0.4097       0.1548
   5:        0.4501       0.2836
</pre>
If you reproduce this example, your numbers will most likely be different,
because, in the absence of <tt>seed=</tt> parameter, <tt>sfnoise</tt>
uses a random seed value to generate pseudo-random numbers. Finally, we
apply <tt>sfheadersort</tt> to shuffle the input traces.
<pre>
bash$ < input.rsf sfheadersort head=random.rsf > output.rsf
bash$ < output.rsf sfdisfil
   0:             2            2            2            2            2
   5:             1            1            1            1            1
  10:             3            3            3            3            3
  15:             5            5            5            5            5
  20:             7            7            7            7            7
  25:             4            4            4            4            4
  30:             6            6            6            6            6
</pre>
As expected, the order of traces in the output file corresponds to the
order of values in the header. Thanks to the separation between
headers and data, the operation of <tt>sfheadersort</tt> is optimally
efficient. It first sorts the headers and only then accesses the data,
reading each data trace only once.

==sfheaderwindow==
{| class="wikitable" align="center" cellspacing="0" border="1"
! colspan="4" style="background:#ffdead;" | Window a dataset based on a header mask.
|-
! colspan="4" | sfheaderwindow mask=head.rsf < in.rsf > out.rsf inv=n
|-
|  colspan="4" | <br>The input data is a collection of traces n1xn2,<br>mask is an integer array os size n2, windowed is n1xm2,<br>where m2 is the number of nonzero elements in mask.
|-
| ''bool   '' || '''inv=n''' ||  [y/n] || 	inversion flag
|-
| ''file   '' || '''mask=''' ||   || 	auxiliary input file name
|}


<tt>sfheaderwindow</tt> is used to window traces in the input file
according to trace header information. 
Here is an example of using <tt>sfheaderwindow</tt> for randomly
selecting part of the traces in the input file. First, let us create
an input file with ten traces:
<pre>
bash$ sfmath n1=5 n2=10 output=x2+1 > input.rsf
bash$ < input.rsf sfdisfil
   0:             1            1            1            1            1
   5:             2            2            2            2            2
  10:             3            3            3            3            3
  15:             4            4            4            4            4
  20:             5            5            5            5            5
  25:             6            6            6            6            6
  30:             7            7            7            7            7
  35:             8            8            8            8            8
  40:             9            9            9            9            9
  45:            10           10           10           10           10
</pre>
Next, we can create a random file with ten header values using
<tt>sfnoise</tt>.
<pre>
bash$ sfspike n1=10 | sfnoise rep=y type=n > random.rsf
bash$ < random.rsf sfdisfil
   0:     -0.005768      0.02258     -0.04331      -0.4129      -0.3909
   5:      -0.03582       0.4595      -0.3326        0.498      -0.3517
</pre>
If you reproduce this example, your numbers will most likely be different,
because, in the absence of <tt>seed=</tt> parameter, <tt>sfnoise</tt>
uses a random seed value to generate pseudo-random numbers. Finally,
we apply <tt>sfheaderwindow</tt> to window the input traces selecting
only those for which the header is greater than zero.
<pre>
bash$ < random.rsf sfmask min=0 > mask.rsf
bash$ < mask.rsf sfdisfil
   0:    0    1    0    0    0    0    1    0    1    0
bash$ < input.rsf sfheaderwindow mask=mask.rsf > output.rsf
bash$ < output.rsf sfdisfil
   0:             2            2            2            2            2
   5:             7            7            7            7            7
  10:             9            9            9            9            9
</pre>
In this case, only three traces are selected for the output. Thanks to
the separation between headers and data, the operation of
<tt>sfheaderwindow</tt> is optimally efficient.

==sfin==
{| class="wikitable" align="center" cellspacing="0" border="1"
! colspan="4" style="background:#ffdead;" | Display basic information about RSF files.
|-
! colspan="4" | sfin info=y check=2. trail=y [<file0.rsf] file1.rsf file2.rsf ...
|-
|  colspan="4" | n1,n2,... are data dimensions<br>o1,o2,... are axis origins<br>d1,d2,... are axis sampling intervals<br>label1,label2,... are axis labels<br>unit1,unit2,... are axis units
|-
| ''float  '' || '''check=2.''' ||   || 	Portion of the data (in Mb) to check for zero values.
|-
| ''bool   '' || '''info=y''' ||  [y/n] || 	If n, only display the name of the data file.
|-
| ''bool   '' || '''trail=y''' ||  [y/n] || 	If n, skip trailing dimensions of  one
|}


<tt>sfin</tt> is one of the most useful programs for operating with
RSF files. It produces quick information on the file hypercube
dimensions and checks the consistency of the associated data file.
Here is an example. Let us create an RSF file and examine it with <tt>sfin</tt>.
<pre>
bash$ sfspike n1=100 n2=20 > spike.rsf
bash$ sfin spike.rsf
spike.rsf:
    in="/tmp/spike.rsf@"
    esize=4 type=float form=native
    n1=100         d1=0.004       o1=0          label1="Time" unit1="s"
    n2=20          d2=0.1         o2=0          label2="Distance" unit2="km"
        2000 elements 8000 bytes
</pre>
<tt>sfin</tt> reports the following information:
  
*location of the data file (<tt>/tmp/spike.rsf\@</tt>) 
*element size (4 bytes) 
*element type (floating point) 
*element form (native) 
*hypercube dimensions (100 by 20) 
*axes scale (0.004 and 0.1) 
*axes origin (0 and 0) 
*axes labels 
*axes units 
*total number of elements 
*total number of bytes in the data file 
Suppose that the file got corrupted by a buggy program and reports
incorrect dimensions. The <tt>sfin</tt> program should be able to
catch the discrepancy.
<pre>
bash$ echo n2=100 >> spike.rsf
bash$ sfin spike.rsf > /dev/null
sfin:           Actually 8000 bytes, 20% of expected.
</pre>
<tt>sfin</tt> also checks the first records in the file for zeros. 
<pre>
bash$ sfspike n1=100 n2=100 k2=99 > spike2.rsf
bash$ sfin spike2.rsf >/dev/null
sfin: The first 32768 bytes are all zeros
</pre>
The number of bytes to check is adjustable
<pre>
bash$ sfin spike2.rsf check=0.01 >/dev/null
sfin: The first 16384 bytes are all zeros
</pre>
You can also output only the location of the data file. This is
sometimes handy in scripts.
<pre>
bash$ sfin spike.rsf spike2.rsf info=n
/tmp/spike.rsf@ /tmp/spike2.rsf@
</pre>
An alternative is to use <tt>sfget</tt>, as follows:
<pre>
bash$ sfget parform=n in < spike.rsf
/tmp/spike.rsf@
</pre>
To actually eliminate annoying trailing dimensions of length one (not just stop displaying them with <tt>trail=n</tt>), you may use <tt>sed</tt>. Example for eliminating axis 6:
<syntaxhighlight lang="bash">
sed -i 's/n6=1//g' file.rsf
</syntaxhighlight>

==sfinterleave==
{| class="wikitable" align="center" cellspacing="0" border="1"
! colspan="4" style="background:#ffdead;" | Combine several datasets by interleaving.
|-
! colspan="4" | sfinterleave > out.rsf axis=3 [< file0.rsf] file1.rsf file2.rsf ...
|-
| ''int    '' || '''axis=3''' ||   || 	Axis for interleaving
|}


<tt>sfinterleave</tt> combines two or more datasets by interleaving them on one
of the axes. Here is a quick example:
<pre>
bash$ sfspike n1=5 n2=5 > one.rsf
bash$ sfdisfil < one.rsf
   0:             1            1            1            1            1
   5:             1            1            1            1            1
  10:             1            1            1            1            1
  15:             1            1            1            1            1
  20:             1            1            1            1            1
bash$ sfscale < one.rsf dscale=2 > two.rsf
bash$ sfdisfil < two.rsf
   0:             2            2            2            2            2
   5:             2            2            2            2            2
  10:             2            2            2            2            2
  15:             2            2            2            2            2
  20:             2            2            2            2            2
bash$ sfinterleave one.rsf two.rsf axis=1 | sfdisfil
   0:             1            2            1            2            1
   5:             2            1            2            1            2
  10:             1            2            1            2            1
  15:             2            1            2            1            2
  20:             1            2            1            2            1
  25:             2            1            2            1            2
  30:             1            2            1            2            1
  35:             2            1            2            1            2
  40:             1            2            1            2            1
  45:             2            1            2            1            2
bash$ sfinterleave < one.rsf two.rsf axis=2 | sfdisfil
   0:             1            1            1            1            1
   5:             2            2            2            2            2
  10:             1            1            1            1            1
  15:             2            2            2            2            2
  20:             1            1            1            1            1
  25:             2            2            2            2            2
  30:             1            1            1            1            1
  35:             2            2            2            2            2
  40:             1            1            1            1            1
  45:             2            2            2            2            2
</pre>

==sfmask==
{| class="wikitable" align="center" cellspacing="0" border="1"
! colspan="4" style="background:#ffdead;" | Create a mask.
|-
! colspan="4" | sfmask < in.rsf > out.rsf min= max=
|-
|  colspan="4" | <br>Mask is an integer data with ones and zeros. <br>Ones correspond to input values between min and max.<br><br>The output can be used with sfheaderwindow.
|-
| ''float  '' || '''max=''' ||   || 	maximum header value
|-
| ''float  '' || '''min=''' ||   || 	minimum header value
|}


<tt>sfmask</tt> creates an integer output of ones and zeros comparing
the values of the input data to specified <tt>min=</tt> and
<tt>max=</tt> parameters. It is useful for <tt>sfheaderwindow</tt> and
in many other applications. Here is a quick example:
<pre>
bash$ sfmath n1=10 output="sin(x1)" > sin.rsf
bash$ < sin.rsf sfdisfil
   0:             0       0.8415       0.9093       0.1411      -0.7568
   5:       -0.9589      -0.2794        0.657       0.9894       0.4121
bash$ < sin.rsf sfmask min=-0.5 max=0.5 | sfdisfil
   0:    1    0    0    1    0    0    1    0    0    1
</pre>

==sfmath==
{| class="wikitable" align="center" cellspacing="0" border="1"
! colspan="4" style="background:#ffdead;" | Mathematical operations on data files.
|-
! colspan="4" | sfmath > out.rsf nostdin=n n#= d#=(1,1,...) o#=(0,0,...) label#= unit#= type= label= unit= output=
|-
|  colspan="4" | <br>Known functions: <br>cos,  sin,  tan,  acos,  asin,  atan, <br>cosh, sinh, tanh, acosh, asinh, atanh,<br>exp,  log,  sqrt, abs,<br>erf,  erfc, sign (for float data),<br>arg,  conj, real, imag (for complex data).<br><br>sfmath will work on float or complex data, but all the input and output<br>files must be of the same data type.<br><br>An alternative to sfmath is sfadd, which may be more efficient, but is<br>less versatile.<br><br>Examples:<br><br>sfmath x=file1.rsf y=file2.rsf power=file3.rsf output='sin((x+2*y)^power)' > out.rsf<br>sfmath < file1.rsf tau=file2.rsf output='exp(tau*input)' > out.rsf<br>sfmath n1=100 type=complex output="exp(I*x1)" > out.rsf<br><br>Arguments which are not treated as variables in mathematical expressions:<br>datapath=, type=, out=<br><br>See also: sfheadermath.
|-
| ''float  '' || '''d#=(1,1,...)''' ||   || 	sampling on #-th axis
|-
| ''string '' || '''label=''' ||   || 	data label
|-
| ''string '' || '''label#=''' ||   || 	label on #-th axis
|-
| ''largeint'' || '''n#=''' ||   || 	size of #-th axis
|-
| ''bool   '' || '''nostdin=n''' ||  [y/n] || 	y - ignore stdin
|-
| ''float  '' || '''o#=(0,0,...)''' ||   || 	origin on #-th axis
|-
| ''string '' || '''output=''' ||   || 	Mathematical description of the output
|-
| ''string '' || '''type=''' ||   || 	output data type [float,complex]
|-
| ''string '' || '''unit=''' ||   || 	data unit
|-
| ''string '' || '''unit#=''' ||   || 	unit on #-th axis
|}

<tt>sfmath</tt> is a versatile program for mathematical operations
with RSF files. It can operate with several input files with the
same dimensions and data type. The data type can be real (floating
point) or complex. The following example demonstrates several
features of <tt>sfmath</tt>.
<pre>
bash$ sfmath n1=629 d1=0.01 o1=0 n2=40 d2=1 o2=5 \
output="x2*(8+sin(6*x1+x2/10))" > rad.rsf
bash$ < rad.rsf sfrtoc | sfmath output="input*exp(I*x1)" > rose.rsf
bash$ < rose.rsf sfgraph title=Rose screenratio=1 wantaxis=n | sfpen
</pre>
The first line creates a 2-D dataset that consists of 40 traces 629
samples each. The values of the data are computed with the formula
<font color="#cd4b19">"x2*(8+sin(6*x1+x2/10))"</font>, where <tt>x1</tt> refers to the
coordinate on the first axis, and <tt>x2</tt> is the coordinate of the
second axis. In the second line, we convert the data from real to
complex using <tt>sfrtoc</tt> and produce a complex dataset using
formula <font color="#cd4b19">"input*exp(I*x1)"</font>, where <tt>input</tt> refers to the
input file. Finally, we plot the complex data as a collection of
parametric curves using <tt>sfgraph</tt> and display the result using
<tt>sfpen</tt>.  The plot appearing on your screen should look similar
to the figure.
[[Image:rose.png|frame|center|This figure was created with <tt>sfmath</tt>.]]
One possible alternative to the second line above is
<pre>
bash$ < rad.rsf sfmath output=x1 > ang.rsf
bash$ sfmath r=rad.rsf a=ang.rsf output="r*cos(a)" > cos.rsf
bash$ sfmath r=rad.rsf a=ang.rsf output="r*sin(a)" > sin.rsf
bash$ sfcmplx cos.rsf sin.rsf > rose.rsf
</pre>
Here we refer to input files by names (<tt>r</tt> and <tt>a</tt>) and combine the names in a formula.

Functions can be nested and combined, and variable names, as well as the <tt>input</tt> keyword, may be combined with the axes keywords <tt>x1</tt>, <tt>x2</tt>, etc. Here is an example that shifts a wavelet -0.4s to t=0, computes the frequency-domain complex square root (equivalent to the convolutional square root in time), and shifts the result back to t=0.4s:
<syntaxhighlight lang="bash">
pi=3.14159265
tshift=0.4

sfspike n1=256 k1=101 |\
sfbandpass flo=4 fhi=30 |\
sffft1 opt=n sym=y |\
sfmath output="sqrt(input*exp(2*$pi*$tshift*x1*I))*exp(-2*$pi*$tshift*x1*I)" |\
sffft1 opt=n sym=y inv=y |\
sfgraph |\
sfpen
</syntaxhighlight>

==sfpad==
{| class="wikitable" align="center" cellspacing="0" border="1"
! colspan="4" style="background:#ffdead;" | Pad a dataset with zeros.
|-
! colspan="4" | sfpad < in.rsf > out.rsf beg#=0 end#=0 n#=
|-
|  colspan="4" | <br>n#out is equivalent to n#, both of them overwrite end#.<br><br>Other parameters from the command line are passed to the output (similar to sfput).
|-
| ''int    '' || '''beg#=0''' ||   || 	the number of zeros to add before the beginning of #-th axis
|-
| ''int    '' || '''end#=0''' ||   || 	the number of zeros to add after the end of #-th axis
|-
| ''int    '' || '''n#=''' ||   || 	the output length of #-th axis - padding at the end
|}

<tt>sfpad</tt> increases the dimensions of the input dataset by padding
the data with zeroes. Here are some simple examples.
<pre>
bash$ sfspike n1=5 n2=3 > one.rsf
bash$ sfdisfil < one.rsf
   0:             1            1            1            1            1
   5:             1            1            1            1            1
  10:             1            1            1            1            1
bash$ < one.rsf sfpad n2=5 | sfdisfil
   0:             1            1            1            1            1
   5:             1            1            1            1            1
  10:             1            1            1            1            1
  15:             0            0            0            0            0
  20:             0            0            0            0            0
bash$ < one.rsf sfpad beg2=2 | sfdisfil
   0:             0            0            0            0            0
   5:             0            0            0            0            0
  10:             1            1            1            1            1
  15:             1            1            1            1            1
  20:             1            1            1            1            1
bash$ < one.rsf sfpad beg2=1 end2=1 | sfdisfil
   0:             0            0            0            0            0
   5:             1            1            1            1            1
  10:             1            1            1            1            1
  15:             1            1            1            1            1
  20:             0            0            0            0            0
bash$ < one.rsf sfwindow n1=3 | sfpad n1=5 n2=5 beg1=1 beg2=1 | sfdisfil
   0:             0            0            0            0            0
   5:             0            1            1            1            0
  10:             0            1            1            1            0
  15:             0            1            1            1            0
  20:             0            0            0            0            0
</pre>
You can use <tt>sfcat</tt> to pad data with values other than zeroes.

==sfput==
{| class="wikitable" align="center" cellspacing="0" border="1"
! colspan="4" style="background:#ffdead;" | Input parameters into a header. 
|-
! colspan="4" | sfput < in.rsf > out.rsf [parameter=value list]
|}


<tt>sfput</tt> is a very simple program. It simply appends parameters
from the command line to the output RSF file. One can achieve similar
results by editing the file by hand or with standard Unix utilities like
<tt>sed</tt> and <tt>echo</tt>. <tt>sfput</tt> is sometimes more
convenient because it handles input/output operations similarly to
other Madagascar programs.
<pre>
bash$ sfspike n1=10 > spike.rsf
bash$ sfin spike.rsf
spike.rsf:
    in="/tmp/spike.rsf@"
    esize=4 type=float form=native
    n1=10          d1=0.004       o1=0          label1="Time" unit1="s"
        10 elements 40 bytes
bash$ sfput < spike.rsf d1=25 label1=Depth unit1=m > spike2.rsf
bash$ sfin spike2.rsf
spike2.rsf:
    in="/tmp/spike2.rsf@"
    esize=4 type=float form=native
    n1=10          d1=25          o1=0          label1="Depth" unit1="m"
        10 elements 40 bytes
</pre>

==sfreal==
{| class="wikitable" align="center" cellspacing="0" border="1"
! colspan="4" style="background:#ffdead;" | Extract real (sfreal) or imaginary (sfimag) part of a complex dataset. 
|-
! colspan="4" | sfreal < cmplx.rsf > real.rsf
|}


<tt>sfreal</tt> extracts the real part of a complex type dataset. The
imaginary part can be extracted with <tt>sfimag</tt>, and the real
and imaginary part can be combined together with <tt>sfcmplx</tt>.
Here is a simple example. Let us first create a complex dataset 
with <tt>sfmath</tt>
<pre>
bash$ sfmath n1=10 type=complex output="(2+I)*x1" > cmplx.rsf
bash$ fdisfil < cmplx.rsf
   0:          0,         0i         2,         1i         4,         2i
   3:          6,         3i         8,         4i        10,         5i
   6:         12,         6i        14,         7i        16,         8i
   9:         18,         9i
</pre>
Extracting the real part with <tt>sfreal</tt>:
<pre>
bash$ sfreal < cmplx.rsf | sfdisfil
   0:             0            2            4            6            8
   5:            10           12           14           16           18
</pre>
Extracting the imaginary part with <tt>sfimag</tt>:
<pre>
bash$ sfimag < cmplx.rsf | sfdisfil
   0:             0            1            2            3            4
   5:             5            6            7            8            9
</pre>

==sfreverse==
{| class="wikitable" align="center" cellspacing="0" border="1"
! colspan="4" style="background:#ffdead;" | Reverse one or more axes in the data hypercube. 
|-
! colspan="4" | sfreverse < in.rsf > out.rsf which=-1 verb=n memsize=sf_memsize() opt=
|-
| ''int    '' || '''memsize=sf_memsize()''' ||   || 	Max amount of RAM (in Mb) to be used
|-
| ''string '' || '''opt=''' ||   || 	If y, change o and d parameters on the reversed axis;
:if i, don't change o and d
|-
| ''bool   '' || '''verb=n''' ||  [y/n] || 	Verbosity flag
|-
| ''int    '' || '''which=-1''' ||   || 	Which axis to reverse.
:To reverse a given axis, start with 0,
:add 1 to number to reverse n1 dimension,
:add 2 to number to reverse n2 dimension,
:add 4 to number to reverse n3 dimension, etc.
:Thus, which=7 would reverse the first three dimensions,
:which=5 just n1 and n3, etc.
:which=0 will just pass the input on through unchanged.
|}
Here is an example of using <tt>sfreverse</tt>. First, let us create a
2-D dataset.
<pre>
bash$ sfmath n1=5 d1=1 n2=3 d2=1 output=x1+x2 > test.rsf
bash$ < test.rsf sfdisfil
   0:             0            1            2            3            4
   5:             1            2            3            4            5
  10:             2            3            4            5            6
</pre>
Reversing the first axis:
<pre>
bash$ < test.rsf sfreverse which=1 | sfdisfil
   0:             4            3            2            1            0
   5:             5            4            3            2            1
  10:             6            5            4            3            2
</pre>
Reversing the second axis:
<pre>
bash$ < test.rsf sfreverse which=2 | sfdisfil
   0:             2            3            4            5            6
   5:             1            2            3            4            5
  10:             0            1            2            3            4
</pre>
Reversing both the first and the second axis:
<pre>
bash$ < test.rsf sfreverse which=3 | sfdisfil
   0:             2            3            4            5            6
   5:             1            2            3            4            5
  10:             0            1            2            3            4
</pre>
As you can see, the <tt>which=</tt> parameter controls the axes that are
being reversed by encoding them into one number.
When an axis is reversed, what happens with its axis origin and
sampling parameters? This behavior is controlled by <tt>opt=</tt>. In
our example,
<pre>
bash$ < test.rsf sfget n1 o1 d1
n1=5
o1=0
d1=1
bash$ < test.rsf sfreverse which=1 | sfget o1 d1
o1=4
d1=-1
</pre>
The default behavior (equivalent to <tt>opt=y</tt>) puts the origin
<tt>o1</tt> at the end of the axis and reverses the sampling parameter
<tt>d1</tt>.  Using <tt>opt=n</tt> preserves the sampling but reverses
the origin.
<pre>
bash$ < test.rsf sfreverse which=1 opt=n | sfget o1 d1
o1=-4
d1=1
</pre>
Using <tt>opt=i</tt> preserves both the sampling and the origin while
reversing the axis.
<pre>
bash$ < test.rsf sfreverse which=1 opt=i | sfget o1 d1
o1=0
d1=1
</pre>
One of the three possible behaviors may be desirable depending on the
application.

==sfrm==
{| class="wikitable" align="center" cellspacing="0" border="1"
! colspan="4" style="background:#ffdead;" | Remove RSF files together with their data.
|-
! colspan="4" | sfrm file1.rsf [file2.rsf ...] [-i] [-v] [-f] 
|-
|  colspan="4" | Mimics the standard Unix rm command.<br><br>See also: sfmv, sfcp.
|}

<tt>sfrm</tt> is a program for removing RSF files. Its arguments mimic
the arguments of the standard Unix <tt>rm</tt> utility: <tt>-v</tt>
for verbosity, <tt>-i</tt> for interactive inquiry, <tt>-f</tt> for
force removal of suspicious files. Unlike the Unix <tt>rm</tt>,
<tt>sfrm</tt> removes both the RSF header files and the binary files
that the headers point to. 
Example:
<pre>
bash&#36; sfspike n1=10 > spike.rsf datapath=./
bash&#36; sfget in < spike.rsf
in=./spike.rsf@
bash&#36; ls spike*
spike.rsf  spike.rsf@
bash&#36; sfrm -v spike.rsf
sfrm: sf_rm: Removing header spike.rsf
sfrm: sf_rm: Removing data ./spike.rsf@
bash&#36; ls spike*
ls: No match.
</pre>
==sfrotate==
{| class="wikitable" align="center" cellspacing="0" border="1"
! colspan="4" style="background:#ffdead;" | Rotate a portion of one or more axes in the data hypercube. 
|-
! colspan="4" | sfrotate < in.rsf > out.rsf verb=n memsize=sf_memsize() rot#=(0,0,...)
|-
| ''int    '' || '''memsize=sf_memsize()''' ||   || 	Max amount of RAM (in Mb) to be used
|-
| ''int    '' || '''rot#=(0,0,...)''' ||   || 	length of #-th axis that is moved to the end
|-
| ''bool   '' || '''verb=n''' ||  [y/n] || 	Verbosity flag
|}

<tt>sfrotate</tt> modifies the input dataset by splitting it into
parts and putting the parts back in a different order. Here is a quick example.
<pre>
bash&#36; sfmath n1=5 d1=1 n2=3 d2=1 output=x1+x2 > test.rsf
bash&#36; < test.rsf sfdisfil
   0:             0            1            2            3            4
   5:             1            2            3            4            5
  10:             2            3            4            5            6
</pre>
Rotating the first axis by putting the last two columns in front:
<pre>
bash&#36; < test.rsf sfrotate rot1=2 | sfdisfil
   0:             3            4            0            1            2
   5:             4            5            1            2            3
  10:             5            6            2            3            4
</pre>
Rotating the second axis by putting the last row in front:
<pre>
bash&#36; < test.rsf sfrotate rot2=1 | sfdisfil
   0:             2            3            4            5            6
   5:             0            1            2            3            4
  10:             1            2            3            4            5
</pre>
Rotating both the first and the second axis:
<pre>
bash&#36; < test.rsf sfrotate rot1=3 rot2=1 | sfdisfil
   0:             4            5            6            2            3
   5:             2            3            4            0            1
  10:             3            4            5            1            2
</pre>
The transformation is shown schematically in Figure~(fig:rotate).
[[Image:rotate.png|frame|center|Schematic transformation of data with <tt>sfrotate</tt>.]]

==sfrtoc==
{| class="wikitable" align="center" cellspacing="0" border="1"
! colspan="4" style="background:#ffdead;" | Convert real data to complex (by adding zero imaginary part).
|-
! colspan="4" | sfrtoc < real.rsf > cmplx.rsf pair=n
|-
|  colspan="4" | <br>See also: sfcmplx
|-
| ''bool   '' || '''pair=n''' ||  [y/n] || 	y - use odd elements for the real part and even ones for the imaginary part
|}

The input to <tt>sfrtoc</tt> can be any <tt>type=float</tt> dataset:
<pre>
bash$ sfspike n1=10 n2=20 n3=30 >real.rsf
bash$ sfin real.rsf
real.rsf:
    in="/var/tmp/real.rsf@"
    esize=4 type=float form=native 
    n1=10          d1=0.004       o1=0          label1="Time" unit1="s" 
    n2=20          d2=0.1         o2=0          label2="Distance" unit2="km" 
    n3=30          d3=0.1         o3=0          label3="Distance" unit3="km" 
        6000 elements 24000 bytes
</pre>
The output dataset will have <tt>type=complex</tt>, and its binary will be twice the size of the input:
<pre>
bash$ <real.rsf sfrtoc >complex.rsf
bash$ sfin complex.rsf 
complex.rsf:
    in="/var/tmp/complex.rsf@"
    esize=8 type=complex form=native 
    n1=10          d1=0.004       o1=0          label1="Time" unit1="s" 
    n2=20          d2=0.1         o2=0          label2="Distance" unit2="km" 
    n3=30          d3=0.1         o3=0          label3="Distance" unit3="km" 
        6000 elements 48000 bytes
</pre>

==sfscale==
{| class="wikitable" align="center" cellspacing="0" border="1"
! colspan="4" style="background:#ffdead;" | Scale data.
|-
! colspan="4" | sfscale < in.rsf > out.rsf axis=0 rscale=0. pclip=100. dscale=1.
|-
|  colspan="4" | <br>To scale by a constant factor, you can also use sfmath.
|-
| ''int    '' || '''axis=0''' ||   || 	Scale by maximum in the dimensions up to this axis.
|-
| ''float  '' || '''dscale=1.''' ||   || 	Scale by this factor (works if rscale=0)
|-
| ''float  '' || '''pclip=100.''' ||   || 	data clip percentile
|-
| ''float  '' || '''rscale=0.''' ||   || 	Scale by this factor.
|}

<tt>sfscale</tt> scales the input dataset by a factor. Here
are some simple examples. First, let us create a test dataset.
<pre>
bash&#36; sfmath n1=5 n2=3 o1=1 o2=1 output="x1*x2" > test.rsf
bash&#36; < test.rsf sfdisfil   
   0:             1            2            3            4            5
   5:             2            4            6            8           10
  10:             3            6            9           12           15
</pre>
Scale every data point by 2:
<pre>
bash&#36; < test.rsf sfscale dscale=2 | sfdisfil
   0:             2            4            6            8           10
   5:             4            8           12           16           20
  10:             6           12           18           24           30
</pre>
Divide every trace by its maximum value:
<pre>
bash&#36; < test.rsf sfscale axis=1 | sfdisfil
   0:           0.2          0.4          0.6          0.8            1
   5:           0.2          0.4          0.6          0.8            1
  10:           0.2          0.4          0.6          0.8            1
</pre>
Divide by the maximum value in the whole 2-D dataset:
<pre>
bash&#36; < test.rsf sfscale axis=2 | sfdisfil
   0:       0.06667       0.1333          0.2       0.2667       0.3333
   5:        0.1333       0.2667          0.4       0.5333       0.6667
  10:           0.2          0.4          0.6          0.8            1
</pre>
The <tt>rscale=</tt> parameter is synonymous with <tt>dscale=</tt> except
when it is equal to zero. With <tt>sfscale dscale=0</tt>, the dataset gets
multiplied by zero. If using <tt>rscale=0</tt>, the other parameters are
used to define scaling.  Thus, <tt>sfscale rscale=0 axis=1</tt> is
equivalent to <tt>sfscale axis=1</tt>, and <tt>sfscale rscale=0</tt>
is equivalent to <tt>sfscale dscale=1</tt>.

==sfspike==
{| class="wikitable" align="center" cellspacing="0" border="1"
! colspan="4" style="background:#ffdead;" | Generate simple data: spikes, boxes, planes, constants. 
|-
! colspan="4" | sfspike < in.rsf > spike.rsf mag= nsp=1 k#=[0,...] l#=[k1,k2,...] p#=[0,...] n#= o#=[0,0,...] d#=[0.004,0.1,0.1,...] label#=[Time,Distance,Distance,...] unit#=[s,km,km,...] title=
|-
|  colspan="4" | <br>Spike positioning is given in samples and starts with 1.
|-
| ''float  '' || '''d#=[0.004,0.1,0.1,...]''' ||   || 	sampling on #-th axis
|-
| ''ints   '' || '''k#=[0,...]''' ||   || 	spike starting position  [nsp]
|-
| ''ints   '' || '''l#=[k1,k2,...]''' ||   || 	spike ending position  [nsp]
|-
| ''string '' || '''label#=[Time,Distance,Distance,...]''' ||   || 	label on #-th axis
|-
| ''floats '' || '''mag=''' ||   || 	spike magnitudes  [nsp]
|-
| ''int    '' || '''n#=''' ||   || 	size of #-th axis
|-
| ''int    '' || '''nsp=1''' ||   || 	Number of spikes
|-
| ''float  '' || '''o#=[0,0,...]''' ||   || 	origin on #-th axis
|-
| ''floats '' || '''p#=[0,...]''' ||   || 	spike inclination (in samples)  [nsp]
|-
| ''string '' || '''title=''' ||   || 	title for plots
|-
| ''string '' || '''unit#=[s,km,km,...]''' ||   || 	unit on #-th axis
|}


<tt>sfspike</tt> takes no input and generates an output with
"spikes". It is an easy way to create data. Here is an example:
<pre>
bash&#36; sfspike n1=5 n2=3 k1=4 k2=1 | sfdisfil
   0:             0            0            0            1            0
   5:             0            0            0            0            0
  10:             0            0            0            0            0
</pre>
The spike location is specified by parameters <tt>k1=4</tt> and
<tt>k2=1</tt>. Note that the locations are numbered starting from 1. 
If one of the parameters is omitted or given the value of zero, the 
spike in the corresponding direction becomes a plane:
<pre>
bash&#36; sfspike n1=5 n2=3 k1=4 | sfdisfil   
   0:             0            0            0            1            0
   5:             0            0            0            1            0
  10:             0            0            0            1            0

</pre>
If no spike parameters are given, the whole dataset is filled with ones:
<pre>
bash&#36; sfspike n1=5 n2=3 | sfdisfil
   0:             1            1            1            1            1
   5:             1            1            1            1            1
  10:             1            1            1            1            1
</pre>
To create several spikes, use the <tt>nsp=</tt> parameter and give a comma-separated list of values to <tt>k#=</tt> arguments:
<pre>
bash&#36; sfspike n1=5 n2=3 nsp=3 k1=1,3,4 k2=1,2,3 | sfdisfil
   0:             1            0            0            0            0
   5:             0            0            1            0            0
  10:             0            0            0            1            0
</pre>
If the number of values in the list is smaller than <tt>nsp</tt>, the
last value gets repeated, and the spikes add on top of each other,
creating larger amplitudes:
<pre>
bash&#36; sfspike n1=5 n2=3 nsp=3 k1=1,3 k2=1,2 | sfdisfil
   0:             1            0            0            0            0
   5:             0            0            2            0            0
  10:             0            0            0            0            0
</pre>
The magnitude of the spikes can be controlled explicitly with the
<tt>mag=</tt> parameter:
<pre>
bash&#36; sfspike n1=5 n2=3 nsp=3 k1=1,3,4 k2=1,2,3 mag=1,4,2 | sfdisfil
   0:             1            0            0            0            0
   5:             0            0            4            0            0
  10:             0            0            0            2            0
</pre>
You can create boxes instead of spikes by using <tt>l#=</tt>
parameters:
<pre>
bash&#36; sfspike n1=5 n2=3 k1=2 l1=4 k2=2 mag=8 | sfdisfil
   0:             0            0            0            0            0
   5:             0            8            8            8            0
  10:             0            0            0            0            0
</pre>
In this case, <tt>k1=2</tt> specifies the box start, and <tt>l1=4</tt>
specifies the box end.
Finally, multi-dimensional planes can be given an inclination by using <tt>p#=</tt> parameters:
<pre>
bash&#36; sfspike n1=5 n2=3 k1=2 p2=1 | sfdisfil
   0:             0            1            0            0            0
   5:             0            0            1            0            0
  10:             0            0            0            1            0
</pre>
Note that <tt>sfspike</tt> interprets the <tt>p#=</tt> parameters in stepout from sample to sample, not in terms of axes units (i.e. s/m). 

More than one p parameter can be specified. In this case, a (hyper) plane will be created.

When the inclination value is not integer, simple linear interpolation is used:
<pre>
bash&#36; sfspike n1=5 n2=3 k1=2 p2=0.7 | sfdisfil
   0:             0            1            0            0            0
   5:             0          0.3          0.7            0            0
  10:             0            0          0.6          0.4            0
</pre>
<tt>sfspike</tt> supplies default dimensions and labels to all axes:
<pre>
bash&#36; sfspike n1=5 n2=3 n3=4 > spike.rsf
bash&#36; sfin spike.rsf
spike.rsf:
    in="/var/tmp/spike.rsf@"
    esize=4 type=float form=native 
    n1=5           d1=0.004       o1=0          label1="Time" unit1="s" 
    n2=3           d2=0.1         o2=0          label2="Distance" unit2="km" 
    n3=4           d3=0.1         o3=0          label3="Distance" unit3="km" 
	60 elements 240 bytes
</pre>
As you can see, the first axis is assumed to be time, with sampling of
<math>0.004</math> seconds. All other axes are assumed to be distance, with
sampling of <math>0.1</math> kilometers. All these parameters can be changed on
the command line.
<pre>
bash&#36; sfspike n1=5 n2=3 n3=4 label3=Offset unit3=ft d3=20 > spike.rsf
bash&#36; sfin spike.rsf
spike.rsf:
    in="/var/tmp/spike.rsf@"
    esize=4 type=float form=native 
    n1=5           d1=0.004       o1=0          label1="Time" unit1="s" 
    n2=3           d2=0.1         o2=0          label2="Distance" unit2="km" 
    n3=4           d3=20          o3=0          label3="Offset" unit3="ft" 
	60 elements 240 bytes
</pre>

==sfspray==
{| class="wikitable" align="center" cellspacing="0" border="1"
! colspan="4" style="background:#ffdead;" | Extend a dataset by duplicating in the specified axis dimension.
|-
! colspan="4" | sfspray < in.rsf > out.rsf axis=2 n= d= o= label= unit=
|-
|  colspan="4" |    This operation is adjoint to sfstack.
|-
| ''int    '' || '''axis=2''' ||   || 	which axis to spray
|-
| ''float  '' || '''d=''' ||   || 	Sampling of the newly created dimension
|-
| ''string '' || '''label=''' ||   || 	Label of the newly created dimension
|-
| ''int    '' || '''n=''' ||   || 	Size of the newly created dimension
|-
| ''float  '' || '''o=''' ||   || 	Origin of the newly created dimension
|-
| ''string '' || '''unit=''' ||   || 	Units of the newly created dimension
|}

<tt>sfspray</tt> extends the input hypercube by replicating the data
in one of the dimensions. The output dataset acquires one additional
dimension. Here is an example:
Start with a 2-D dataset
<pre>
bash&#36; sfmath n1=5 n2=2 output=x1+x2 > test.rsf
bash&#36; sfin test.rsf
test.rsf:
    in="/var/tmp/test.rsf@"
    esize=4 type=float form=native 
    n1=5           d1=1           o1=0          
    n2=2           d2=1           o2=0          
        10 elements 40 bytes
bash&#36; < test.rsf sfdisfil
   0:             0            1            2            3            4
   5:             1            2            3            4            5
</pre>
Extend the data in the second dimension
<pre>
bash&#36; < test.rsf sfspray axis=2 n=3 > test2.rsf
bash&#36; sfin test2.rsf
test2.rsf:
    in="/var/tmp/test2.rsf@"
    esize=4 type=float form=native 
    n1=5           d1=1           o1=0          
    n2=3           d2=1           o2=0          
    n3=2           d3=1           o3=0          
        30 elements 120 bytes
bash&#36; < test2.rsf sfdisfil
   0:             0            1            2            3            4
   5:             0            1            2            3            4
  10:             0            1            2            3            4
  15:             1            2            3            4            5
  20:             1            2            3            4            5
  25:             1            2            3            4            5
</pre>
The output is three-dimensional, with traces from the original
data duplicated along the second axis.
Extend the data in the third dimension
<pre>
bash&#36; < test.rsf sfspray axis=3 n=2 > test3.rsf
bash&#36; sfin test3.rsf
test3.rsf:
    in="/var/tmp/test3.rsf@"
    esize=4 type=float form=native 
    n1=5           d1=1           o1=0          
    n2=2           d2=1           o2=0          
    n3=2           d3=?           o3=?          
        20 elements 80 bytes
bash&#36; < test3.rsf sfdisfil
   0:             0            1            2            3            4
   5:             1            2            3            4            5
  10:             0            1            2            3            4
  15:             1            2            3            4            5
</pre>
The output is also three-dimensional, with the original data replicated
along the third axis.

==sfstack==
{| class="wikitable" align="center" cellspacing="0" border="1"
! colspan="4" style="background:#ffdead;" | Stack a dataset over one of the dimensions.
|-
! colspan="4" | sfstack < in.rsf > out.rsf scale= axis=2 rms=n norm=y min=n max=n prod=n
|-
|  colspan="4" | <br>This operation is adjoint to sfspray.
|-
| ''int    '' || '''axis=2''' ||   || 	which axis to stack. If axis=0, stack over all dimensions
|-
| ''bool   '' || '''max=n''' ||  [y/n] || 	If y, find maximum instead of stack.  Ignores rms and norm.
|-
| ''bool   '' || '''min=n''' ||  [y/n] || 	If y, find minimum instead of stack.  Ignores rms and norm.
|-
| ''bool   '' || '''norm=y''' ||  [y/n] || 	If y, normalize by fold.
|-
| ''bool   '' || '''prod=n''' ||  [y/n] || 	If y, find product instead of stack.  Ignores rms and norm.
|-
| ''bool   '' || '''rms=n''' ||  [y/n] || 	If y, compute the root-mean-square instead of stack.
|-
| ''floats '' || '''scale=''' ||   || 	optionally scale before stacking  [n2]
|}

While <tt>sfspray</tt> adds a dimension to a hypercube,
<tt>sfstack</tt> effectively removes one of the dimensions by stacking
over it. Here are some examples:
<pre>
bash&#36; sfmath n1=5 n2=3 output=x1+x2 > test.rsf
bash&#36; < test.rsf sfdisfil
   0:             0            1            2            3            4
   5:             1            2            3            4            5
  10:             2            3            4            5            6
bash&#36; < test.rsf sfstack axis=2 | sfdisfil
   0:           1.5            2            3            4            5
bash&#36; < test.rsf sfstack axis=1 | sfdisfil
   0:           2.5            3            4
</pre>
Why is the first value not 1 (in the first case) or 2 (in the second
case)? By default, <tt>sfstack</tt> normalizes the stack by the fold
(the number of non-zero entries). To avoid normalization, use
<tt>norm=n</tt>, as follows:
<pre>
bash&#36; < test.rsf sfstack norm=n | sfdisfil 
   0:             3            6            9           12           15
</pre>
<tt>sfstack</tt> can also compute root-mean-square values as 
well as minimum and maximum values.
<pre>
bash&#36; < test.rsf sfstack rms=y | sfdisfil
   0:         1.581         2.16        3.109        4.082        5.066
bash&#36; < test.rsf sfstack min=y | sfdisfil
   0:             0            1            2            3            4
bash&#36; < test.rsf sfstack axis=1 max=y | sfdisfil 
   0:             4            5            6
</pre>

==sftransp==
{| class="wikitable" align="center" cellspacing="0" border="1"
! colspan="4" style="background:#ffdead;" | Transpose two axes in a dataset. 
|-
! colspan="4" | sftransp < in.rsf > out.rsf memsize=sf_memsize() plane=
|-
|  colspan="4" | <br>If you get a "Cannot allocate memory" error, give the program a<br>memsize=1 command-line parameter to force out-of-core operation.
|-
| ''int    '' || '''memsize=sf_memsize()''' ||   || 	Max amount of RAM (in Mb) to be used
|-
| ''int    '' || '''plane=''' ||   || 	Two-digit number with axes to transpose. The default is 12
|}

The <tt>sftransp</tt> program transposes the input hypercube
exchanging the two axes specified by the <tt>plane=</tt> parameter. 
<pre>
bash&#36; sfspike n1=10 n2=20 n3=30 > orig123.rsf
bash&#36; sfin orig123.rsf 
orig123.rsf:
    in="/var/tmp/orig123.rsf@"
    esize=4 type=float form=native 
    n1=10          d1=0.004       o1=0          label1="Time" unit1="s" 
    n2=20          d2=0.1         o2=0          label2="Distance" unit2="km" 
    n3=30          d3=0.1         o3=0          label3="Distance" unit3="km" 
        6000 elements 24000 bytes
bash&#36; <orig123.rsf sftransp plane=23 >out132.rsf
bash&#36; sfin out132.rsf 
out132.rsf:
    in="/var/tmp/out132.rsf@"
    esize=4 type=float form=native 
    n1=10          d1=0.004       o1=0          label1="Time" unit1="s" 
    n2=30          d2=0.1         o2=0          label2="Distance" unit2="km" 
    n3=20          d3=0.1         o3=0          label3="Distance" unit3="km" 
        6000 elements 24000 bytes
bash&#36; <orig123.rsf sftransp plane=13 >out321.rsf
bash&#36; sfin out321.rsf 
out321.rsf:
    in="/var/tmp/out132.rsf@"
    esize=4 type=float form=native 
    n1=30          d1=0.1         o1=0          label1="Distance" unit1="km" 
    n2=20          d2=0.1         o2=0          label2="Distance" unit2="km" 
    n3=10          d3=0.004       o3=0          label3="Time" unit3="s" 
        6000 elements 24000 bytes
</pre>
<tt>sftransp</tt> tries to fit the dataset in memory to transpose it
there but, if not enough memory is available, it performs a slower
transpose out of core using disk operations. You can control the
amount of available memory using the <tt>memsize=</tt> parameter or
the <tt>RSFMEMSIZE</tt> environmental variable.

==sfwindow==
{| class="wikitable" align="center" cellspacing="0" border="1"
! colspan="4" style="background:#ffdead;" | Window a portion of a dataset. 
|-
! colspan="4" | sfwindow < in.rsf > out.rsf verb=n squeeze=y j#=(1,...) d#=(d1,d2,...) f#=(0,...) min#=(o1,o2,,...) n#=(0,...) max#=(o1+(n1-1)*d1,o2+(n1-1)*d2,,...)
|-
|  colspan="4" | <br>Other parameters from the command line are passed to the output (similar to sfput).
|-
| ''float  '' || '''d#=(d1,d2,...)''' ||   || 	sampling in #-th dimension
|-
| ''largeint'' || '''f#=(0,...)''' ||   || 	window start in #-th dimension
|-
| ''int    '' || '''j#=(1,...)''' ||   || 	jump in #-th dimension
|-
| ''float  '' || '''max#=(o1+(n1-1)*d1,o2+(n1-1)*d2,,...)''' ||   || 	maximum in #-th dimension
|-
| ''float  '' || '''min#=(o1,o2,,...)''' ||   || 	minimum in #-th dimension
|-
| ''largeint'' || '''n#=(0,...)''' ||   || 	window size in #-th dimension
|-
| ''bool   '' || '''squeeze=y''' ||  [y/n] || 	if y, squeeze dimensions equal to 1 to the end
|-
| ''bool   '' || '''verb=n''' ||  [y/n] || 	Verbosity flag
|}


<b>sfwindow</b> is used to window a portion of the dataset. Here is
a quick example: Start by creating some data.
<pre>
bash&#36; sfmath n1=5 n2=3 o1=1 o2=1 output="x1*x2" > test.rsf
bash&#36; < test.rsf sfdisfil
   0:             1            2            3            4            5
   5:             2            4            6            8           10
  10:             3            6            9           12           15
</pre>
Now window the first two rows:
<pre>
bash&#36; < test.rsf sfwindow n2=2 | sfdisfil
   0:             1            2            3            4            5
   5:             2            4            6            8           10
</pre>
Window the first three columns:
<pre>
bash&#36; < test.rsf sfwindow n1=3 | sfdisfil
   0:             1            2            3            2            4
   5:             6            3            6            9
</pre>
Window the middle row:
<pre>
bash&#36; < test.rsf sfwindow f2=1 n2=1 | sfdisfil
   0:             2            4            6            8           10
</pre>
You can interpret the <b>f#</b> and <b>n#</b> parameters as
meaning "skip that many rows/columns" and "select that many
rows/columns" correspondingly. Window the middle point in the dataset:
<pre>
bash&#36; < test.rsf sfwindow f1=2 n1=1 f2=1 n2=1 | sfdisfil
   0:             6
</pre>
Window every other column:
<pre>
bash&#36; < test.rsf sfwindow j1=2 | sfdisfil
   0:             1            3            5            2            6
   5:            10            3            9           15
</pre>
Window every third column:
<pre>
bash&#36; < test.rsf sfwindow j1=3 | sfdisfil
   0:             1            4            2            8            3
   5:            12
</pre>
Alternatively, <b>sfwindow</b> can select a window from the minimum and maximum
parameters. In the following example, we are
creating a dataset with <b>sfspike</b> and then windowing a portion of it
between 1 and 2 seconds in time and sampled at 8 miliseconds.
<pre>
bash&#36; sfspike n1=1000 n2=10 > spike.rsf         
bash&#36; sfin spike.rsf
spike.rsf:
    in="/var/tmp/spike.rsf@"
    esize=4 type=float form=native 
    n1=1000        d1=0.004       o1=0          label1="Time" unit1="s" 
    n2=10          d2=0.1         o2=0          label2="Distance" unit2="km" 
        10000 elements 40000 bytes
bash&#36; < spike.rsf sfwindow min1=1 max1=2 d1=0.008 > window.rsf
bash&#36; sfin window.rsf
window.rsf:
    in="/var/tmp/window.rsf@"
    esize=4 type=float form=native 
    n1=126         d1=0.008       o1=1          label1="Time" unit1="s" 
    n2=10          d2=0.1         o2=0          label2="Distance" unit2="km" 
        1260 elements 5040 bytes
</pre>
By default, <b>sfwindow</b> "squeezes" the hypercube dimensions
that are equal to one toward the end of the dataset. Here is an
example of taking a time slice:
<pre>
bash&#36; < spike.rsf sfwindow n1=1 min1=1 > slice.rsf
bash&#36; sfin slice.rsf 
slice.rsf:
    in="/var/tmp/slice.rsf@"
    esize=4 type=float form=native 
    n1=10          d1=0.1         o1=0          label1="Distance" unit1="km" 
    n2=1           d2=0.004       o2=1          label2="Time" unit2="s" 
        10 elements 40 bytes
</pre>
You can change this behavior by specifying <b>squeeze=n</b>.
<pre>
bash&#36; < spike.rsf sfwindow n1=1 min1=1 squeeze=n > slice.rsf
bash&#36; sfin slice.rsf 
slice.rsf:
    in="/var/tmp/slice.rsf@"
    esize=4 type=float form=native 
    n1=1           d1=0.004       o1=1          label1="Time" unit1="s" 
    n2=10          d2=0.1         o2=0          label2="Distance" unit2="km" 
        10 elements 40 bytes
</pre>