Radiation pattern of a few moment tensor earthquake sources

Four type of earthquake sources:
1) explosion M=[1 0 0;0 1 0;0 0 1]
2) single force M=[1 0 0;0 0 0;0 0 0]
3) double couple M=[1 0 0;0 0 0;0 0 -1]
4) CLVD M=[2 0 0;0 -1 0;0 0 -1]
and their corresponding radiation patterns of both P- and S-wave in the far field

Explosion source and vector source

Explosion source and vector source

Double-couple and CLVD source

Double-couple and CLVD source

Other than that explosion generates no S-wave, all other sources generate both P- and S-wave energy in the far field including DC, CLVD.

Automatically generate N distinct colors – A python script

Inspired by discussions here, I coded up a python script to generate N (<=20) distinct colors in RGB.

"""
To generate N distinct colorcode, N < 20
First created by Junwei Huang, March 21 2013
"""
from numpy import *
import matplotlib.pyplot as plt

def hex_to_rgb(value):
	value = value.lstrip('#')
	lv = len(value)
	return tuple(int(value[i:i+lv/3], 16) for i in range(0, lv, lv/3))

N=60

ColourValues=["FF0000", "00FF00", "0000FF", "FFFF00", "FF00FF", "00FFFF", "000000",
        "800000", "008000", "000080", "808000", "800080", "008080", "808080",
        "C00000", "00C000", "0000C0", "C0C000", "C000C0", "00C0C0", "C0C0C0",
        "400000", "004000", "000040", "404000", "400040", "004040", "404040",
        "200000", "002000", "000020", "202000", "200020", "002020", "202020",
        "600000", "006000", "000060", "606000", "600060", "006060", "606060",
        "A00000", "00A000", "0000A0", "A0A000", "A000A0", "00A0A0", "A0A0A0",
        "E00000", "00E000", "0000E0", "E0E000", "E000E0", "00E0E0", "E0E0E0"]
M=len(ColourValues)

plt.figure()
plt.hold('on')
for i in range(N):
	plt.bar(i,1.0,color="#"+ColourValues[mod(i,M)])
	print hex_to_rgb(ColourValues[mod(i,M)])

plt.show()

As seen in the figure, it is distinct if N < 20.

RGB color code

RGB color code

3D visualization of DXF STL file using mayavi python script

extensive googling doesn’t find a good solution. Writing a python script from scratch actually takes much less time.

"""
To read a STL file and plot in mayavi
First created by Junwei Huang @ Toronto, Feb 26, 2013
"""

from numpy import *
from mayavi import mlab

STLfile="test.stl"
f=open(STLfile,'r')

x=[]
y=[]
z=[]

for line in f:
	strarray=line.split()
	if strarray[0]=='vertex':
		x=append(x,double(strarray[1]))
		y=append(y,double(strarray[2]))
		z=append(z,double(strarray[3]))

triangles=[(i, i+1, i+2) for i in range(0, len(x),3)]

mlab.triangular_mesh(x, y, z, triangles)
mlab.show()
tunnel geometry from STL file

tunnel geometry from STL file

Q: What about DXF file?
A: Convert DXF to STL file first before using this script:)

Why a c++ program based on openmpi 1.4.3 failed to execute the system() function

It wasted me a full day to find out why running a c++ program with a system() function call failed.

Background: to output moderate size of files by 30 processors and merge them to a single file.

MPI_Barrier( MPI_COMM_WORLD );
if(my_rank==0) {
   int sysflag;
   if (!system(NULL)) {printf("nError!...command processor is unavailablen");exit(1);}
   for (int i=0;i stackX.bin ",i);sysflag=system(cmd);//printf("system returned %dn", sysflag);
      } else {
         sprintf(cmd,"cat PE%d_loc_stackX.3d >> stackX.bin ",i);sysflag=system(cmd);//printf("system returned %dn", sysflag);
         sprintf(cmd,"rm PE%d_loc_stackX.3d",i);sysflag=system(cmd); //printf("system returned %dn", sysflag);
      }
   }
}
MPI_Barrier( MPI_COMM_WORLD );

That “Command processor is unavailable” is always the case, although this “save-individually-and-merge-later” works for another program. Extensive webpage search may indicate openmpi sometimes doesn’t support fork() which is used by system() call. However, this is not the reason in my case, since it works for other programs.

Options:
(a) all slave processors send data to the master which writes to a single file;
(b) all processors write their own file and merge to a single file afterwards;
(c) use MPI I/O.

My choice:
(b). (b) is faster than (a) and is easier to implement than (c).

My guess:
My understanding is the system() will spawn a sub-process costing the same amount of memory as other processors. Therefore in my case, one node running 5 processors already consumed 90% of memory and thus won’t have enough memory to run another sub-process. Consequently, function call of system() fails. In addition, system() does not throw exceptions upon failure. Thus the c++ program doesn’t execute the command passing to system() and appears nothing happened.

My solution:
reduce the number of node from 5 to 4 solves the problem.

Doxygen Intergration in Visual Studio 2010 Environment

When the source code package is huge including more than 2000 files and each may have > 10,000 lines of code, it is a good idea to generate some documents using doxygen. “Doxygen is a documentation system for C++, C, Java, Objective-C, Python, IDL (Corba and Microsoft flavors), Fortran, VHDL, PHP, C#, and to some extent D.” (webpage), it works out of the box on Linux/UNIX. For windows under Visual studio environment, some extra settings are needed. Here are the steps:

1. Download the Win32 binaries here and install them. For example, you installed doxygen in c:program filesdoxygen.
2. Create a doxygen configuration file and save it as default.doxygen under the same folder as the project. Notice that to generate graphical representation of class hierachy, you need to install graphivz (link). Here is the one I used.

# Doxyfile 1.8.2

#---------------------------------------------------------------------------
# Project related configuration options
#---------------------------------------------------------------------------
DOXYFILE_ENCODING      = UTF-8
PROJECT_NAME           = "InSite"
PROJECT_NUMBER         = 2.16
OUTPUT_DIRECTORY       =
CREATE_SUBDIRS         = NO
OUTPUT_LANGUAGE        = English
BRIEF_MEMBER_DESC      = YES
REPEAT_BRIEF           = YES
ABBREVIATE_BRIEF       = "The $name class" 
 "The $name widget" 
 "The $name file" 
 is 
 provides 
 specifies 
 contains 
 represents 
 a 
 an 
 the
ALWAYS_DETAILED_SEC    = NO
INLINE_INHERITED_MEMB  = NO
FULL_PATH_NAMES        = YES
STRIP_FROM_PATH        =
STRIP_FROM_INC_PATH    =
SHORT_NAMES            = NO
JAVADOC_AUTOBRIEF      = NO
QT_AUTOBRIEF           = NO
MULTILINE_CPP_IS_BRIEF = NO
INHERIT_DOCS           = YES
SEPARATE_MEMBER_PAGES  = NO
TAB_SIZE               = 8
ALIASES                =
OPTIMIZE_OUTPUT_FOR_C  = NO
OPTIMIZE_OUTPUT_JAVA   = NO
OPTIMIZE_FOR_FORTRAN   = NO
OPTIMIZE_OUTPUT_VHDL   = NO
BUILTIN_STL_SUPPORT    = NO
CPP_CLI_SUPPORT        = YES
SIP_SUPPORT            = NO
IDL_PROPERTY_SUPPORT   = YES
DISTRIBUTE_GROUP_DOC   = NO
SUBGROUPING            = YES
TYPEDEF_HIDES_STRUCT   = NO
SYMBOL_CACHE_SIZE      = 0
#---------------------------------------------------------------------------
# Build related configuration options
#---------------------------------------------------------------------------
EXTRACT_ALL            = YES
EXTRACT_PRIVATE        = YES
EXTRACT_STATIC         = YES
EXTRACT_LOCAL_CLASSES  = YES
EXTRACT_LOCAL_METHODS  = YES
EXTRACT_ANON_NSPACES   = YES
HIDE_UNDOC_MEMBERS     = NO
HIDE_UNDOC_CLASSES     = NO
HIDE_FRIEND_COMPOUNDS  = NO
HIDE_IN_BODY_DOCS      = NO
INTERNAL_DOCS          = NO
CASE_SENSE_NAMES       = YES
HIDE_SCOPE_NAMES       = NO
SHOW_INCLUDE_FILES     = YES
INLINE_INFO            = YES
SORT_MEMBER_DOCS       = YES
SORT_BRIEF_DOCS        = NO
SORT_GROUP_NAMES       = NO
SORT_BY_SCOPE_NAME     = NO
GENERATE_TODOLIST      = YES
GENERATE_TESTLIST      = YES
GENERATE_BUGLIST       = YES
GENERATE_DEPRECATEDLIST= YES
ENABLED_SECTIONS       =
MAX_INITIALIZER_LINES  = 30
SHOW_USED_FILES        = YES
SHOW_DIRECTORIES       = NO
SHOW_FILES             = YES
SHOW_NAMESPACES        = YES
FILE_VERSION_FILTER    =
LAYOUT_FILE            =
#---------------------------------------------------------------------------
# configuration options related to warning and progress messages
#---------------------------------------------------------------------------
QUIET                  = NO
WARNINGS               = YES
WARN_IF_UNDOCUMENTED   = YES
WARN_IF_DOC_ERROR      = YES
WARN_NO_PARAMDOC       = NO
WARN_FORMAT            = "$file:$line: $text"
WARN_LOGFILE           =
#---------------------------------------------------------------------------
# configuration options related to the input files
#---------------------------------------------------------------------------
INPUT                  = .
INPUT_ENCODING         = UTF-8
FILE_PATTERNS          = *.c 
 *.cc 
 *.cxx 
 *.cpp 
 *.c++ 
 *.d 
 *.java 
 *.ii 
 *.ixx 
 *.ipp 
 *.i++ 
 *.inl 
 *.h 
 *.hh 
 *.hxx 
 *.hpp 
 *.h++ 
 *.idl 
 *.odl 
 *.cs 
 *.php 
 *.php3 
 *.inc 
 *.m 
 *.mm 
 *.dox 
 *.py 
 *.f90 
 *.f 
 *.vhd 
 *.vhdl
RECURSIVE              = YES
EXCLUDE                =
EXCLUDE_SYMLINKS       = NO
EXCLUDE_PATTERNS       =
EXCLUDE_SYMBOLS        =
EXAMPLE_PATH           =
EXAMPLE_PATTERNS       = *
EXAMPLE_RECURSIVE      = NO
IMAGE_PATH             =
INPUT_FILTER           =
FILTER_PATTERNS        =
FILTER_SOURCE_FILES    = NO
#---------------------------------------------------------------------------
# configuration options related to source browsing
#---------------------------------------------------------------------------
SOURCE_BROWSER         = YES
INLINE_SOURCES         = YES
STRIP_CODE_COMMENTS    = YES
REFERENCED_BY_RELATION = YES
REFERENCES_RELATION    = YES
REFERENCES_LINK_SOURCE = YES
USE_HTAGS              = NO
VERBATIM_HEADERS       = YES
#---------------------------------------------------------------------------
# configuration options related to the alphabetical class index
#---------------------------------------------------------------------------
COLS_IN_ALPHA_INDEX    = 5
IGNORE_PREFIX          =
#---------------------------------------------------------------------------
# configuration options related to the HTML output
#---------------------------------------------------------------------------
GENERATE_HTML          = YES
HTML_OUTPUT            = html
HTML_FILE_EXTENSION    = .html
HTML_HEADER            =
HTML_FOOTER            =
HTML_STYLESHEET        =
HTML_ALIGN_MEMBERS     = YES
HTML_DYNAMIC_SECTIONS  = NO
GENERATE_DOCSET        = NO
DOCSET_FEEDNAME        = "Doxygen generated docs"
DOCSET_BUNDLE_ID       = org.doxygen.Project
GENERATE_HTMLHELP      = NO
CHM_FILE               =
HHC_LOCATION           =
GENERATE_CHI           = NO
CHM_INDEX_ENCODING     =
BINARY_TOC             = NO
TOC_EXPAND             = NO
GENERATE_QHP           = NO
QCH_FILE               =
QHP_NAMESPACE          = org.doxygen.Project
QHP_VIRTUAL_FOLDER     = doc
QHG_LOCATION           =
DISABLE_INDEX          = NO
ENUM_VALUES_PER_LINE   = 4
GENERATE_TREEVIEW      = ALL
TREEVIEW_WIDTH         = 250
FORMULA_FONTSIZE       = 10
#---------------------------------------------------------------------------
# configuration options related to the LaTeX output
#---------------------------------------------------------------------------
GENERATE_LATEX         = NO
LATEX_OUTPUT           = latex
LATEX_CMD_NAME         = latex
MAKEINDEX_CMD_NAME     = makeindex
COMPACT_LATEX          = NO
PAPER_TYPE             = a4wide
EXTRA_PACKAGES         =
LATEX_HEADER           =
PDF_HYPERLINKS         = YES
USE_PDFLATEX           = YES
LATEX_BATCHMODE        = NO
LATEX_HIDE_INDICES     = NO
#---------------------------------------------------------------------------
# configuration options related to the RTF output
#---------------------------------------------------------------------------
GENERATE_RTF           = NO
RTF_OUTPUT             = rtf
COMPACT_RTF            = NO
RTF_HYPERLINKS         = NO
RTF_STYLESHEET_FILE    =
RTF_EXTENSIONS_FILE    =
#---------------------------------------------------------------------------
# configuration options related to the man page output
#---------------------------------------------------------------------------
GENERATE_MAN           = NO
MAN_OUTPUT             = man
MAN_EXTENSION          = .3
MAN_LINKS              = NO
#---------------------------------------------------------------------------
# configuration options related to the XML output
#---------------------------------------------------------------------------
GENERATE_XML           = NO
XML_OUTPUT             = xml
XML_SCHEMA             =
XML_DTD                =
XML_PROGRAMLISTING     = YES
#---------------------------------------------------------------------------
# configuration options for the AutoGen Definitions output
#---------------------------------------------------------------------------
GENERATE_AUTOGEN_DEF   = NO
#---------------------------------------------------------------------------
# configuration options related to the Perl module output
#---------------------------------------------------------------------------
GENERATE_PERLMOD       = NO
PERLMOD_LATEX          = NO
PERLMOD_PRETTY         = YES
PERLMOD_MAKEVAR_PREFIX =
#---------------------------------------------------------------------------
# Configuration options related to the preprocessor
#---------------------------------------------------------------------------
ENABLE_PREPROCESSING   = YES
MACRO_EXPANSION        = NO
EXPAND_ONLY_PREDEF     = NO
SEARCH_INCLUDES        = YES
INCLUDE_PATH           =
INCLUDE_FILE_PATTERNS  =
PREDEFINED             =
EXPAND_AS_DEFINED      =
SKIP_FUNCTION_MACROS   = YES
#---------------------------------------------------------------------------
# Configuration::additions related to external references
#---------------------------------------------------------------------------
TAGFILES               =
GENERATE_TAGFILE       =
ALLEXTERNALS           = NO
EXTERNAL_GROUPS        = YES
PERL_PATH              = /usr/bin/perl
#---------------------------------------------------------------------------
# Configuration options related to the dot tool
#---------------------------------------------------------------------------
CLASS_DIAGRAMS         = YES
MSCGEN_PATH            =
HIDE_UNDOC_RELATIONS   = YES
HAVE_DOT               = YES
DOT_FONTNAME           = FreeSans
DOT_FONTPATH           =
CLASS_GRAPH            = YES
COLLABORATION_GRAPH    = YES
GROUP_GRAPHS           = YES
UML_LOOK               = YES
TEMPLATE_RELATIONS     = YES
INCLUDE_GRAPH          = YES
INCLUDED_BY_GRAPH      = YES
CALL_GRAPH             = YES
CALLER_GRAPH           = YES
GRAPHICAL_HIERARCHY    = YES
DIRECTORY_GRAPH        = YES
DOT_IMAGE_FORMAT       = png
DOT_PATH               = "C:Program Files (x86)Graphviz 2.28bindot.exe"
DOTFILE_DIRS           =
DOT_GRAPH_MAX_NODES    = 50
MAX_DOT_GRAPH_DEPTH    = 1000
DOT_TRANSPARENT        = YES
DOT_MULTI_TARGETS      = YES
GENERATE_LEGEND        = YES
DOT_CLEANUP            = YES
#---------------------------------------------------------------------------
# Configuration::additions related to the search engine
#---------------------------------------------------------------------------
SEARCHENGINE           = YES

3. Add a new custom tool, called “DoxyGen”, with the following parameters:
Command: c:program filesdoxygenbindoxygen.exe (or where you installed it)
Arguments: “$(ItemDir)default.doxygen” (the config file – including the quotes)
Initial Directory: $(ItemDir)
Check the “Use the output window” box
4. Add another tool, to view the results new “View DoxyDoc” tool, to view the results:
Command your favorite browser, e.g. C:Program Files (x86)Mozilla Firefoxfirefox.exe
Arguments: “$(ItemDir)htmlindex.html” (including the quotes)
Initial Directory: leave empty
5. Under the solution view, right click the project and select “Add”->”Existing Item…” and in the pop-up window select default.doxygen under your project folder. Open it for editing in VC++, and Change the PROJECT_NAME and the PROJECT_NUMBER to match your project.
6. Choose Tools/doxygen from the menu, and watch the output dump into the Output window of visual studio. Once finished, check for any error messages.
7. You are done! Use Tool->View DoxyDoc to open the index.html file.

Doxygen in Visual Studio 2010

Doxygen in Visual Studio 2010


This article was inspired by this blog (link).

Real time cloud based hydraulic fracture monitoring for real

Came across a ppt slide showing some company is doing real time hydraulic fracture monitoring using cloud service. I think it is a neat technique all service companies can have, not just Schlumberger or Halliburton.

Background information:
(1) hydraulic fracturing is a technique to pump high pressurized fluid into petroleum reservoir to enhance oil recovery and at the same time reduce cost
(2) cloud technique is the computing facility offering the horsepower for data crunching and delivering processed results to terminals such as a desktop in a remote office or a smartphone in a manager’s pocket.
(3) the combination of the two enables geoscientists and managers in all scale enterprise to make real time decisions to assist field operations. Time and money are significantly reduced. Thousand hours or days of CPU time is no longer the privilege of a few rich corporations.

Two slides are here:

cloud based realtime hydraulic fracture monitoring

cloud based realtime hydraulic fracture monitoring


The original ppt is here

Use seismic Interferometry, reverse time migration, reverse time source imaging, time reversal mirror, back-projection method, back-propagation method, diffraction stacking, and adjoint method, to locate micro-earthquake

Believe it or not, when you are locating microseismic events without picking the arrivals you may be using one of the above techniques. More importantly they are actually describing the same thing (more or less) in the application of locating micro-earthquakes. Here is why (non-math version):
(1) Seismic Interferometry: is to calculate a green’s function by correlating a complex conjugate green function with another green’s function. In time domain, the complex conjugate green’s function is to back-propagate the other green’s function into the medium. If the complex conjugate green’s function is a natural green’s function recorded from the real medium, you are doing seismic interferometry.
(2) Reverse time migration: if the complex conjugate green’s function is based on a velocity model and the other green’s function is related to the reflection response of the medium, you are doing reverse time migration started since McMechan 1983.
(3) Reverse time source imaging: if the complex conjugate green’s function is based on a velocity model and the other green’s function is related to the transmission response of the medium, you are doing reverse time source imaging. In other field, people call it time reversal method, time reversal acoustics, or call the device that can reverse time signal time reversal mirror. Someone even finds that adding a “sink” (the opposite of source) at the original source location, super-resolution can be achieved. In other words, in a dark night you see a car 1~2 km away approaches you with its light on, you only see one big light not two till it gets closer. Now if your eye and brain can do time reverse, you see two lights crystal clear even when the car is still 1~2km away from you. Isn’t that cool? But you need more than your eye and brain to achieve this super-resolution.
(4) Time reversal mirror: see (3)
(5) Back-projection method: you use ray-tracing method to calculate the complex conjugate green’s function. In other words, you are lightening up the medium, the wavelength is zero and likely you only use the first arrivals. Some specific ray-tracing methods may give you later arrivals. So you can select what late arrivals include. The more the better, and of course, the complex conjugate green’s function from the natural medium include all multiples.
(6) Back-propagation method: the function of the complex conjugate green’s function, see (1).
(7) Diffraction stacking: it works like this, (a) define a source location, (b) use ray-tracing to calculate its first arrivals for all sensors, (c) sum up the amplitude along the arrival time curve. Obviously you are using ray-tracing to do the back-projection, so see (5).
(8) Adjoint method: it uses a model based complex conjugate green’s function and back-propagate the difference between modeled data and observed data into the medium, till the difference is reduced to a certain level. At the beginning, you have no idea about the source, so you just back-propagate the observation. Obviously, its first iteration is (3). After the first iteration, you have some ideas about the source so you add a sink at the estimated location and repeat back-propagation, this is similar to (4) with a sink. If it converged, you just did reverse time source inversion. If you are happy with the first iteration, you just finished the reverse time source imaging. The rest is somewhere in between, and I guess we don’t have a name for each of them yet. But considering that we already coined 8 different key words for basically the same thing, we will coin more fancy words for them in no time. Just keep an eye on the literature.

In summary, you are using seismic Interferometry, reverse time migration, reverse time source imaging, time reversal mirror, back-projection method, back-propagation method, diffraction stacking, and adjoint method, to locate micro-earthquakes; you are just unaware of it.

Access virtualbox Linux guest from windows host and access windows host share folder from virtualbox Linux guest

I have to do a lot of googling every time installing virtualbox. So I’d better summarize the process here for future reference.
The situation:
Installed Linux guest in Virtualbox on Windows 7 host.
The goal:
Windows 7 (the host) can access the Linux (the guest) and the guest and access the host as well.
The way:
(1) Create a shared folder in Virtualbox and call it D_DRIVE
(2) add mount.vboxsf D_DRIVE /home/junwei/D_Drive to the first line of /etc/rc.local
(3) sudo reboot
Now the guest Linux can access files on the host Windows via D_Drive.
(4) install samba on the guest Linux
“Ubuntu Software Center”->search samba->install it
(5) In Samba Server Configuration, add share. For example, I share /home/junwei
(6) sudo poweroff
(7) In Virtualbox, use two network adapter, one is NAT to allow the guest Linux to access Internet, one is “host-only” for Samba.
(8) power up the guest Linux
(9) On Windows host network, look for the name of the Virtual Linux machine. Double click it and provide the name and the password. There you go!

The benefit of doing the above is
(a) calculation results on Linux can be directly accessed from Windows software without copying to the Linux mounted shared folder.
(b) calculation results on Windows can be directly accessed from Linux software without copying to the Windows network shared folder.

Linux_Windows


Enjoy Lindows/Winux programming!

Install Madagascar (svn 7866) under sage

As I mentioned in the last post, SageMath is a powerful platform for scientific computation integrating your favorite software packages. One of my favorite packages, for example, is Madagascar for seismic data processing and visualization. So I decided to integrate it with sage. I took the following steps:

0. Sage is installed at $SAGE_ROOT
1. unzip one version of Madagascar (in my case it is svn 7866 checked out in Nov 2011). Newer versions should work as well.
2. under the directory of madagascar source, type

sage -sh

to get into the sage shell, followed by
3.

./configure API=python

Sage will automatically configure it with Sage environmental variables and will install it to $SAGE_ROOT/local
4. exit sage shell, and type

sage --scons install

5. now you should find all Madagascar binary executables in $SAGE_ROOT/local/bin together with other binaries files:

x11pen             sfsort            sfplanemis3     sfitaupmo2         sfederiv2d          sfabalance                  rst2s5.py                   mcube
vpwhitepruf        sfsnrstack        sfplanemis2     sfitaupmo          sfdwt97             sfaastack                   rst2xml.py                  optimal
vpvrms             sfsnr             sfplane         sfisolr3           sfdwt               sfTestcdstep                rstpep2html.py              class.x
vpsg               sfsmstack         sfpick3         sfisolr2one        sfduwt              sfScanCoef                  ipbori                      cws.x
vpreflkine         sfsmoothreg2      sfpick          sfisolr2tau        sfdsr2              pspen                       PolyGUI                     mori.x
vpreflexpt         sfsmoothreg       sfphaserot      sfisolr2           sfdsr               ppmpen                      python2                     nef.x
vpreflector        sfsmoothderLS     sfpfactor2      sfisin2ang         sfdrays             pngpen                      python2-config              poly.x
vppen              sfsmoothder2      sfpetscawefd2d  sfiphase           sfdrayinte          pdfpen                      python-config               mori-11d.x
vpnmotraj          sfsmoothder1      sfpermlr2       sfinvrec1          sfdowmf             oglpen                      python                      nef-11d.x
vplot2png          sfsmoothder       sfpermlr3       sfinvbin1          sfdottest           latex2wiki                  python2.7                   cws-11d.x
vplot2gif          sfsmoothcur       sfpermlr1       sfinvbin           sfdots              jpegpen                     srptool                     class-11d.x
vplot2eps          sfsmooth          sfpen           sfinttest2         sfdonf              gdpen                       certtool                    poly-11d.x
vplot2avi          sfslice           xtpen           sfinttest1         sfdominantf         mayavi2                     gnutls-cli-debug            mori-6d.x
vphyplay           sfslant           sfpefdeburst    sfintshow          sfdomf              tvtk_doc                    psktool                     nef-6d.x
vpheadray          sfsizes           sfpef           sfintervalVTI      sfdoeps             editra                      gnutls-cli                  cws-6d.x
vpfrancis          sfsinc            sfpatch         sfinterpt          sfdmo               helpviewer                  gnutls-serv                 class-6d.x
vpellipse          sfsin             sfparcel        sfinterp2          sfdmeig             img2py                      libgnutls-extra-config      poly-6d.x
vpdcretard         sfsimilarity      sfpad           sfinterleave       sfdixshape          pyalamode                   libgnutls-config            mori-5d.x
vpdc               sfsimenv          sfovc           sfintegral1        sfdix               pycrust                     tconic                      nef-5d.x
vpcsp              sfsignoi          sfopsmigrk      sfintbin3          sfdivn              pywrap                      indep                       cws-5d.x
vpcroshyp          sfsigmoid         sfomp           sfintbin           sfdistmap           img2png                     allisog                     class-5d.x
vpconvert          sfsic3d           sfofsemb        sfint3d            sfdistance          img2xpm                     twist                       poly-5d.x
vpcmp              sfsic             sfofpwd2        sfinstattr         sfdisfil            pyalacarte                  torsion                     mori-4d.x
vp7ab              sfshotprop        sfofpwd         sfinmo3            sfdips              pyshell                     conductor                   nef-4d.x
tiffpen            sfshotholes       sfofilp         sfinmo             sfdiplet            pywxrc                      tate                        cws-4d.x
svgpen             sfshotconstkirch  sfoff2abs3      sfinitial          sfdipfilter         xrced                       findinf                     class-4d.x
sfztrace           sfshot2grid       sfoff2abs       sfinfill           sfdip2              wx-config                   ratpoint                    poly-4d.x
sfzomva            sfshot2cmp        sfofd2_7        sfin               sfdip               wxrc                        mwrank                      vsyasm
sfzomig3           sfshoot2          sfofd2_5        sfimray            sfdimag             wxrc-2.8                    opencdk-config              ytasm
sfzomig            sfshifts          sfofd2_25       sfimpl3            sfdijkstra          vtkpython                   fplll_micro                 yasm
sfzoeppritz        sfshift           sfofd2_13       sfimpl2            sfdiffraction       lproj                       llldiff                     freetype-config
sfzerocross        sfshearer         sfofd2_10       sfimpl1            sfdiffoc            vtkEncodeString             generate                    patch
sfzero             sfsharpsimi       sfofd1_5        sfimospray         sfdifference        vtkWrapPython               fplll_verbose               bzless
sfzcwt             sfsharpen         sfofd1          sfimage            sfdespike3          vtkWrapPythonInit           fplll                       bzcmp
sfxlagtoang2d      sfshapesigk       sfoctentwt      sfimag             sfdespike2          ctest                       sage_singular               bzegrep
sfxcor2d           sfshapebin1       sfocparcel      sfreal             sfdespike           cpack                       Singular                    bzfgrep
sfwmf              sfshapebin        sfoclet         sfigrad            sfderiv             cmake                       singular                    bzdiff
sfwkbjTI           sfshapeagc        sfnsmooth1      sfidempatch        sfdepth2time        ccmake                      LLL                         bzmore

and python files under

$SAGE_ROOT/local/lib/python2.7/site-packages/rsf$ ls
apibak.py    gui             sfbash.pyc      sfgee.py        sfjun.pyc       sfmain.py       sfplplot.pyc      sfseplib_compat.py   sftour.pyc       use.py         vplot2png.py
apibak.pyc   __init__.py     sfbrowaeys.py   sfgee.pyc       sfjyan.py       sfmain.pyc      sfpsava.py        sfseplib_compat.pyc  sftrip.py        use.pyc        vplot2png.pyc
api.py       __init__.pyc    sfbrowaeys.pyc  sfgeneric.py    sfjyan.pyc      sfmccowan.py    sfpsava.pyc       sfslim.py            sftrip.pyc       user           vpplot.py
book.py      latex2wiki.py   sfchen.py       sfgeneric.pyc   sfkourkina.py   sfmccowan.pyc   sfrickettj.py     sfslim.pyc           sftriptrace.py   version.py     vpplot.pyc
book.pyc     latex2wiki.pyc  sfchen.pyc      sfgodwinj.py    sfkourkina.pyc  sfparvaneh.py   sfrickettj.pyc    sfsongxl.py          sftriptrace.pyc  version.pyc
conf.py      path.py         sfcuda.py       sfgodwinj.pyc   sflcasasan.py   sfparvaneh.pyc  sfroman.py        sfsongxl.pyc         sftsai.py        vpconvert.py
conf.pyc     path.pyc        sfcuda.pyc      sfivlad.py      sflcasasan.pyc  sfpens.py       sfroman.pyc       sfsumain.py          sftsai.pyc       vpconvert.pyc
conjgrad.py  prog.py         sfeffsilva.py   sfivlad.pyc     sflexing.py     sfpens.pyc      sfsalah.py        sfsumain.pyc         sfxuxin.py       vplot2avi.py
doc.py       prog.pyc        sfeffsilva.pyc  sfjeff.py       sflexing.pyc    sfpetsc.py      sfsalah.pyc       sfsuplot.py          sfxuxin.pyc      vplot2avi.pyc
doc.pyc      proj.py         sffomels.py     sfjeff.pyc      sflibplot.py    sfpetsc.pyc     sfsaragiotis.py   sfsuplot.pyc         sfyliu.py        vplot2eps.py
dottest.py   proj.pyc        sffomels.pyc    sfjennings.py   sflibplot.pyc   sfplot.py       sfsaragiotis.pyc  sftariq.py           sfyliu.pyc       vplot2eps.pyc
flow.py      recipes         sfgchliu.py     sfjennings.pyc  sfllisiw.py     sfplot.pyc      sfseismic.py      sftariq.pyc          tex.py           vplot2gif.py
flow.pyc     sfbash.py       sfgchliu.pyc    sfjun.py        sfllisiw.pyc    sfplplot.py     sfseismic.pyc     sftour.py            tex.pyc          vplot2gif.pyc

Now you can start a sage notebook and text the following code:

import m8r as sf
import numpy, pylab
f = sf.spike(n1=1000,k1=300)[0]
# sf.spike is an operator
# f is an RSF file object
f.attr()
# Inspect the file with sfattr
b = sf.bandpass(fhi=2,phase=1)[f]
# Now f is filtered through sfbandpass
c = sf.spike(n1=1000,k1=300).bandpass(fhi=2,phase=1)[0]
# c is equivalent to b but created with a pipe
g = c.wiggle(clip=0.02,title='Welcome to Madagascar')
# g is a Vplot file object
g.image()
# Display it on the screen
d = b - c
# Elementary arithmetic operations on files are defined
g = g + d.wiggle(wanttitle=False)
# So are operations on plots
g.show()
# This shows a movie

The result looks like the following

Link to the Madagascar sage notebook page.

install MayaVi2 under sage (v5.0.1) for 2D/3D interactive visulization

According to Wikipedia, “MayaVi is a scientific data visualizer written in Python, which uses VTK and provides a GUI via Tkinter”, while “Sage (previously SAGE, System for Algebra and Geometry Experimentation[3]) is mathematical software with features covering many aspects of mathematics, including algebra, combinatorics, numerical mathematics, number theory, and calculus”. Under the framework of sage, multiple scientific libraries were combined in order to offer a “open source alternative to Magma, Maple, Mathematica, and MATLAB”. From my point of view, although matplotlib under sage can produce high quality figures, it’s 3D scientific visualization capability can be better improved if MayaVi is included in sage as a package. Unfortunately, till version 5.0.1, sage community does not offer an official spkg package to install MayaVi. Users should either install MayaVi outside of sage using system python or install MayaVi under sage using sage’s own python and other libraries. The former is easy and MayaVi community provides good guidance. The latter can integrates the power of both MayaVi and Sage, but is not straightforward. Information scatters over the world wide web and none provides a systematic treatment. In this blog, I summarize how I installed MayaVi2 (v3.4.0) under sage (v5.0.1). I am aware a newer version of MayaVi (v4.2.0) exists, however, I still don’t know why it doesn’t work under sage v5.0.1. Let me know if you manage to get it work. Now, let’s start to install MayaVi2 under Sage.

0. sage v5.0.1 should be already installed on your system. I use $SAGE_ROOT and $SAGE_LOCAL to denote the path of sage root directory and $SAGE_LOCAL=$SAGE_ROOT/local.
1. make sure cmake is installed in your system. I use cmake 2.8.3 under ubuntu 11.04. Cmake package for sage, such as the one on http://sage.math.washington.edu/home/jsp/SPKGS/MayaviOSMesa/ this page, is not useful in my case. Because cmake under sage can’t find paths for X11 related libraries which are located in both /usr/lib and /usr/local/lib. For this reason, I use cmake outside of sage.
2. download mesa-7.2.spkg, vtk-cvs-20090316osmesa.spkg, configobj-4.5.3.spkg, wxPython-2.8.11.0.spkg to $SAGE_ROOT.
3.

# install a local osmesa, this will work according to Ondrej and Prabhu:
./sage -i mesa-7.2.spkg

4.

# install the vtk-cvs-... (this is vtk-5.5.x) with modifications from Prabhu
# with osmesa enabled, see above.
# We use -f -m to keep the installation files in spkg/build,
# so we have VTKDATA and other stuff available
./sage -f -m vtk-cvs-20090316osmesa.spkg

5. The above installation will fail. No panic. This is the step to get out of sage. The error message is about python.h is not found. Go to directory: $SAGE_ROOTspkgbuildvtk-cvs-20090316osmesa, and edit file spkg-install. Basically just to change “2.5” to “2.7”. You may tell the age of this spkg. Under Linux shell, change directory to $SAGE_ROOTspkgbuildvtk-cvs-20090316osmesasrcVTKvtk-build and paste the following command:

export $SAGE_ROOT=/your/path/to/sage-5.0.1
export $SAGE_LOCAL=/your/path/to/sage-5.0.1/local
cmake -DBUILD_SHARED_LIBS:BOOL=ON
        -DVTK_WRAP_PYTHON:BOOL=ON
        -DPYTHON_EXECUTABLE:FILEPATH=$SAGE_LOCAL/bin/python
        -DPYTHON_INCLUDE_PATH:PATH=$SAGE_LOCAL/include/python2.7
        -DPYTHON_LIBRARY:FILEPATH=$SAGE_LOCAL/lib/python2.7/config/libpython2.7.a
        -DVTK_DATA_ROOT:FILEPATH=../VTKData
        -DVTK_USE_GUISUPPORT:BOOL=ON
        -DVTK_USE_DISPLAY:BOOL=ON
        -DVTK_USE_GL2PS:BOOL=ON
        -DVTK_USE_X:BOOL=ON
        -DBUILD_TESTING:BOOL=OFF
        -DCMAKE_INSTALL_PREFIX:PATH=$SAGE_LOCAL
        -DCMAKE_USE_RELATIVE_PATHS:BOOL=OFF
        -DVTK_OPENGL_HAS_OSMESA=ON 
        -DVTK_USE_OFFSCREEN=ON 
        -DOSMESA_INCLUDE_DIR=$SAGE_LOCAL 
        -DOSMESA_LIBRARY=$SAGE_LOCAL/lib/libOSMesa.so 
        -DOPENGL_INCLUDE_DIR=$SAGE_LOCAL/include 
        -DOPENGL_gl_LIBRARY=$SAGE_LOCAL/lib/libGL.so 
        -DOPENGL_glu_LIBRARY=$SAGE_LOCAL/lib/libGLU.so 
      ..

Afterwards, I like to use

ccmake ../

to check all parameters were set properly and all libraries are found. Then I type

make -j2

since I have two CPU cores in my Linux machine, and

make install
CWD=`pwd`
cd Wrapping/Python
sage --python setup.py install
cd $CWD
cp -a $SAGE_ROOT/local/lib/vtk-5.5/libvtk* $SAGE_ROOT/local/lib

You may tell that I was just following the steps in the spkg-install file, just not running it within sage or sage shell.

6.

# install configobj.spkg, needed by mayavi2 (and all other enthought packages)
./sage -i configobj-4.5.3.spkg

7.

# for now install wxPython
./sage -f -m wxPython.2.8.11.0.spkg

8.

# at last
./sage -f -m mayavi-3.1.1.rev23224.spk

The end of installation. Now go to sage notebook webpage, and evaluate the following command:

from enthought.mayavi import mlab as M
M.options.offscreen = True
M.test_mesh()
M.axes()
M.title("Mesh")
M.savefig("mesh.png")

If you should be able to see something like this:

mesh generated by Mayavi2 on sage notebook

Last but not least, you can now generate a x3d file on sage notebook and view it under any brower that has a viewer plugin. In my case I installed freeWRL. You can download the x3d file and view it on your own machine. Here is the source code:

#Test embedding x3d files
from enthought.mayavi import mlab
mlab.options.offscreen = True
mlab.clf()
import numpy
# A simple triangular mesh
mlab.test_triangular_mesh()
# Some nice vectors.
x, y, z = numpy.mgrid[-2:3,-2:3, -2:3]
r = numpy.sqrt(x**2 + y**2 + z**4)
u = y*numpy.sin(r)/(r+0.001)
v = -x*numpy.sin(r)/(r+0.001)
w = numpy.zeros_like(z)
mlab.quiver3d(x,y,z,u,v,w,scale_factor=1)
# Export a x3d file for online 3D viewing.
# FreeWRL (http://freewrl.sourceforge.net) works great!
mlab.savefig("a.x3d")

a.x3d viewed from freewrl standalone application


The worksheet of sage notebook is published at
https://junweihuang.info:8000/home/pub/1/