Editing Guide to madagascar API (section)

==Python interface==
Examples that use the python interface are in the directory $RSFSRC/api/python/test.

The Python script clip.py is:.

<syntaxhighlight lang="python">
#!/usr/bin/env python

import numpy
import m8r

par = m8r.Par()
inp  = m8r.Input()
output = m8r.Output()
assert 'float' == inp.type

n1 = inp.int("n1")
n2 = inp.size(1)
assert n1

clip = par.float("clip")
assert clip

trace = numpy.zeros(n1,'f')

for i2 in xrange(n2): # loop over traces
    inp.read(trace)
    trace = numpy.clip(trace,-clip,clip)
    output.write(trace)

</syntaxhighlight>

Let us examine it in detail. 
<syntaxhighlight lang="python">
import numpy
import m8r
</syntaxhighlight>

The script starts with importing the <tt>numpy</tt> module and the
<tt>m8r</tt> API.
<syntaxhighlight lang="python">
par = m8r.Par()
inp  = m8r.Input()
output = m8r.Output()
assert 'float' == inp.type
</syntaxhighlight>

Next, we initialize the command line interface and the standard input and
output files. We also make sure that the input file type is floating point.
<syntaxhighlight lang="python">
n1 = input.int("n1")
n2 = input.size(1)
assert n1
</syntaxhighlight>

We extract the "<tt>n1</tt>" parameter from the input file.
Conceptually, the RSF data model is a multidimensional hypercube.  The
<tt>n1</tt> parameter refers to the fastest axis. If the input dataset
is a collection of traces, <tt>n1</tt> corresponds to the trace
length. We could proceed in a similar fashion, extracting <tt>n2</tt>,
<tt>n3</tt>, etc. If we are interested in the total number of traces,
like in the clip example, a shortcut is to use the <tt>size</tt>
method of the <tt>Input</tt> class.  Calling <tt>size(0)</tt> returns
the total number of elements in the hypercube (the product of
<tt>n1</tt>, <tt>n2</tt>, etc.), calling <tt>size(1)</tt> returns the
number of traces (the product of <tt>n2</tt>, <tt>n3</tt>, etc.),
calling <tt>size(2)</tt> returns the product of <tt>n3</tt>,
<tt>n4</tt>, etc.
<syntaxhighlight lang="python">
clip = par.float("clip")
assert clip
</syntaxhighlight>

The clip parameter is read from the command line, where it can be specified,
for example, as <tt>clip=10</tt>.
<syntaxhighlight lang="python">
trace = numpy.zeros(n1,'f')
for i2 in xrange(n2): # loop over traces
    inp.read(trace)
    trace = numpy.clip(trace,-clip,clip)
    output.write(trace)
</syntaxhighlight>

Finally, we do the actual work: allocate an array to hold a trace and  loop 
over input traces, reading, clipping, and writing out each trace.

Alternative code use inp.trace to allocate an array and read 
the whole file into memory. The loop is no longer needed:
<syntaxhighlight lang="python">
alltraces=inp.read()
alltraces = numpy.clip(alltraces,-clip,clip)
output.write(alltraces)
</syntaxhighlight>


===Compiling===
The python script does not require compilation. Simply make sure that
<tt>&#36;RSFROOT/lib</tt> is in <tt>PYTHONPATH</tt> and <tt>LD_LIBRARY_PATH</tt>.

===Using the Python API for interactive development in Jupyter notebooks===
Jupyter notebooks are a good way to prototype python code, explore data, and to 
integrate rich text to describe your code.  There are four example notebooks in 
$RSFSRC/api/python/test that demonstrate how to use Jupyter notebooks.  These are:
{|class="wikitable"
 |-
 | clip.ipynb || clip.py converted to a Jupyter notebook
 |-
 |simple_m8r_create_write_filter_read.ipynb || numpy, write rsf, Madagascar commands, read rsf, plot
 |-
 | tle_edge_preserve.ipynb || reproduces a paper from TLE with numpy and Matplotlib
 |-
 | file_filter_scons_example.ipynb || advanced m8r options;  File, filter, and scons from python  
|} 
They can be started form the command line with the commands:
<pre>
cd $RSFSRC/api/python/test
jupyter notebook clip.ipynb
</pre>


===Interactive mode usage without graphics===
Madagascar's Python API can be used interactively too. Create an input dataset with

<pre>
sfmath  n1=10 n2=9 output=x1+x2  > test.rsf
</pre>

Then, start the python interpreter and paste the following to its command line:

<syntaxhighlight lang="python">
import numpy, m8r

inp = m8r.Input('test.rsf')
n1 = inp.int("n1")
n2 = inp.int("n2")

data =inp.read(shape=(n2,n1))
data = data.transpose() # Example of numpy in action

print data
</syntaxhighlight>

You will get

<pre>
[[ 0.  1.  2.  3.  4.  5.  6.  7.  8.]
 [ 1.  2.  3.  4.  5.  6.  7.  8.  9.]
 [ 2.  3.  4.  5.  6.  7.  8.  9. 10.]
 [ 3.  4.  5.  6.  7.  8.  9. 10. 11.]
 [ 4.  5.  6.  7.  8.  9. 10. 11. 12.]
 [ 5.  6.  7.  8.  9. 10. 11. 12. 13.]
 [ 6.  7.  8.  9. 10. 11. 12. 13. 14.]
 [ 7.  8.  9. 10. 11. 12. 13. 14. 15.]
 [ 8.  9. 10. 11. 12. 13. 14. 15. 16.]
 [ 9. 10. 11. 12. 13. 14. 15. 16. 17.]]
</pre>

This code will also work in batch mode in a Python script, not only pasted to the interpreter's command line. 

===Graphics with Matplotlib===
Python can plot arrays directly from memory, without having to write a file to disk first. [http://matplotlib.sourceforge.net/ Matplotlib] is one of the [http://en.wikipedia.org/wiki/Category:Free_plotting_software several] packages that accomplish this. To create a figure, execute the code in the previous section, followed by:

<syntaxhighlight lang="python">
from pylab import *
imshow(data)
xlabel('X (m)')
ylabel('Y (m)')
title('Matplotlib example')
</syntaxhighlight>

If you want to pop up a figure in an interactive session, after pasting to a Python command line the code shown before, also paste:

<syntaxhighlight lang="python">
show()
</syntaxhighlight>

You will get Figure 1. The figure will pop up if you run the code in a script too, and the script will stop until the figure is manually closed. You must press the floppy disk button in order to save it. To have the image written to disk automatically, instead of <tt>show()</tt> use:

<syntaxhighlight lang="python">
savefig('myfile.png')
</syntaxhighlight>

[[Image:matplotlib_imshow.png]]

Putting it all together, here is a sample script reading a RSF file from stdin and printing out a figure:

<syntaxhighlight lang="python">
#!/usr/bin/env python
import m8r, numpy, sys, pylab

inp = m8r.Input('test.rsf')
n1 = inp.int("n1")
n2 = inp.int("n2")

data=input.read()

pylab.imshow(data)
pylab.savefig('out.png')
</syntaxhighlight>
===Python interfaces to the standalone programs===
The <tt>m8r</tt> module offers a way to call Madagascar standalone programs (including graphics) elegantly from inside a Python program, interactively or in batch mode. The blog contains examples of running the <tt>m8r</tt> module [http://www.ahay.org/rsflog/index.php?/archives/173-Extending-Python-interface.html from inside a SAGE notebook] or [http://www.ahay.org/rsflog/index.php?/archives/264-Running-Madagascar-in-an-interactive-console.html from inside an iPython shell].