This page describes the package NA-sampler. It consists of a set of user callable subroutines which implement the NA multi-dimensional parameter space search algorithm. The package contains source code for the NA routines and an example driver program. The latter applies the NA algorithm to the estimation of crustal seismic velocity profiles from inversion of seismic receiver functions (RFI). Several independent `utility' programs are also provided which manipulate the output files of the main program and plot some results. Some of these programs are specific to the RFI problem, and others may be more generally useful.These pages are meant only as a guide. For more information look at the various README files in the src sub-directories, the comments in the source code for the RFI example, and the input files for the utility plot programs. Finally if you are that adventurous you can look at the code itself.
This package performs a derivative free parameter space search using the Neighourhood Algorithm. See references below for details. The latest details of NA-sampler and the related ensemble inference algorithm NA-Bayes can be found at the NA homepage (http://rses.anu.edu.au/~malcolm/na/na.html)
InstallationWhen the tar file is unpacked it produces four directories:
data - contains data input files and example output files created by RFI example problem. Also contains data files for running splot graphics program, and RFI plot programs. (See utility programs for more details.)
bin - contains executable for RFI example program and all utility programs.
docs - contains these HTML pages.
> cp macros/makena.macros_sun makena.macros
To compile the complete set of programs from the src directory
If compilation is successful, all executables are placed in the bin directory and all examples demonstrating the utility programs can be run from the data directory (See below). The reverse of make all is,
which will remove all .o files and all executables.
Note that the Demo directory, the Makefile and the file rfi_na.f can serve as a template for new user applications. If you are not interested in running, or looking at output from the seismic example program then you can remove all RFI specific source code and programs by,
> \rm -fr rfi_subs
> cd ../data
> \rm -fr demo_results rfi_files
Parallelization with MPI
In May 2002 message passing interface (MPI) calls were added to the NA source code. This allows the possibility of running the search algorithm across a cluster of machines simultaneously (e.g. a linux Beowulf cluster). The advantage of this is that the forward modelling for each model in an ensemble can be performed on separate machines, giving much faster overall run times. Note that this does not mean that the user's own forward modelling code needs to be parallelized with MPI. All MPI calls appear in the NA source code and are used to farm the forward modelling out to different processors.
The MPI code is activated by a switch in the installation file `makena.macros'. If you have say a linux cluster running MPI then you need to set the TYPE=mpi and the MF77 variable to the location of the MPI fortran 77 compiler on your system. See the example makena.macros files for examples. (If you want the usual single processor version then set TYPE=serial.) If installation occurs with these variables then you must execute the program using the mpirun command provided by your MPI installation in the usual way. For example one might run the seismic example executable with a command like,
which will execute the program rfi_na with the forward modelling distributed across the first 20 machines that your MPI is set up to run on.
Note that when the number of processors (e.g. 20 in this example) divides equally into the number of models generated at each iteration (see parameter Ns2 below) then one gets perfect parallelization in the forward modelling part of the code, and this usually means a significant overall speed up, e.g. a factor of Ns1/20 in the forward modelling part of the code. Since for many applications most time will be spent on forward modelling then this can represent a significant speed up of the overall code. When the number of processes is not a factor of Ns2 then one gets a smaller speed up factor. Note that there is no point in asking for more processors than Ns2 because there will be nothing for the extra processors to do.
Compiling problems and platform dependency
So far I have succesfully compiled and run all codes in this package on a SUN workstation (running Solaris 2.6 or later) using both the native compiler and the GNU compilers g77 and gcc; on a Compaq Alpha (running OSF1 V4.0) using the Digital Fortran and C compilers; on a linux box with g77, or intel fortran/C compilers ifc and icc. In all cases exactly the same results were obtained. Using the intel fortran and C compilers on Linux I obtained different answers. Possibly due to rounding error changing the Monte Carlo sequence. The default makena.macros file is for a linux box running the intel compiler. When installing on other platforms SUN you will need to edit the file src/makena.macros appropriately.
If you are using a Dec/Compaq alpha/compaq platform then you will need to use the compiling option -assume byterecl which makes the Fortran open statement for direct access files the same as on a sun, i.e. with record length specified in bytes rather than words. This can be done by editing the file src/makena.macros, which is used to define compilers and options during the make. An example file makena.macros_alpha is provided as an example.
Note that for the example problem direct access output (nad) files are provided with the package (in directory data/demo_results). To allow platform independence both little endian and big endian binary versions are included. The make file uses your choice of platform (specified in file makena.macros) to decide which type of binary file to link to the examples. This is a change over earlier distributions and means that you will now be able to run all of the utility programs, e.g. splot (showing the models generated by the NA ) on the example output files provided. The downside of this is that the distribution tar file is larger.
Some programs (from T. Shibutani, RCEP, DPRI, Kyoto Univ., Japan) are provided in the directory src/utl/byteconv for converting SAC files between a sun and a Compaq Alpha machine. This may be useful for seismologists with many sac files. Note that SAC files are the format used to store the input and output seismic receiver functions used in the RFI example. They are specific to this example and are read in by the user routines. They are not a generic NA input file.
In January 2004 the source codes of the NA_sampler package were changed to incorporate dynamic memory allocation using Fortran 90. This means that array bounds in file na_param.inc only need to be set if using Fortran 77 compatibility mode, which is itself set by -DNA_F77_COMPATIBILITY=1 in the file makena.macros. Try the following
> diff makena.macros_alpha_f90 makena.macros_alpha_f77
Code structure and interfacing
The figure above shows the structure of the NA code and how it is interfaced with the user code. The green routines must be provided by the user. The main routine passes control to the NA code by calling the subroutine `na' with no arguments. This then calls three user written subroutines.
forward - calculates the misfit/objective function for a given input model.
writemodels - allows the user to write out, in their own format, all of the models, misfit values etc, from the ensemble produced by the NA.
The file rfi_na.f contains the main driver program and the three subroutines which implement the RFI example problem. In this case the driver routine merely calls na and exits. All data I/O and initialization is performed inside the routine user_init. The three RFI interface routines make use of many other problem specific subroutines and libraries. Source files for all of these are in the directory src/rfi_files. They are compiled and linked into the main program by the Makefile.
To use the NA subroutines you need to do the following,
Note that the only essential information that the NA routines need from the user routines is the definition of the parameter space (from user_init) and the misfit for a given model (from forward). The routine writemodels is called once after the NA has finished and can be null.
Output filesThe two generic output files are na.sum and na.nad The first is an ascii file with summary information on the performance, the second is a direct access file binary file containing all of the models generated by the NA during the search. The format of this file is described in Appendix C. The NAD file is used as a compact storage format which allows users to add their own information in a header. NAD files are read in by certain utility programs like the plot program splot (See appendix A). NAD files are also read in by the appraisal algorithm NA-Bayes, which is a separate package that can be downloaded. See the NA homepage for details.
Code developmentThe development of this code is an ongoing research effort and I expect that the NA subroutines will be updated from time to time, however, it is intended that the interface routines will remain unchanged. (The input file na.in will probably change) In this way the incorporation of updated NA routines into an existing application should be quite straightforward, i.e. compiling the new na.f to obtain na.o and linking with the existing driver code. Note that the NA source code in sub-directory NA_src should not normally need editing ! (Edit at your own risk !)
The code was developed on a sun network under Solaris 2.6. Note that there are some issues to consider when transporting the code across to new platforms and compilers systems. The ones that I am aware of are described in the section on Compiling problems.
Controlling the Neighbourhood algorithmThe NA control parameters and other options are read in from the file na.in.
Each line is described below.
# Neighbourhood Algorithm input options file
ialg : Algorithm type (1=Uniform search, 2=Neighbourhood algorithm)
maxit : Maximum number of iterations
Ns1 : Number of models for initial sample
Ns2 : Number of models generated at each iteration
Nr : Number of Neighbourhoods re-sampled at each iteration
YorN, iseed : Use quasi random number generator ? (y/n); random seed
istype : Type of initial sample (0=random;+-1=read in a NAD file)
olevel : Output info level (0=silent;1=summary info;2=1+write models)
YorN : Turn timing mode on ? (y/n)
YorN : Turn debug mode on ? (y/n)
The parameter istype controls the type of initial ensemble of models used. Several options are available.
If a quasi random number generator is not chosen then a pseudo-random number generator is used. (See the papers cited below or Press et al 1992 Numerical Recipes, C.U.P., for a discussion of the difference.)
Tuning the control parameters
The two control parameters usually need to be tuned for the problem in hand and the style of sampling that you want. There are no well defined rules but experience has proven to be a useful guide. Most of my understanding is based on empirical evidence from applications at RSES. The usual objective with a direct sampling algorithm is to try and achieve a balance between exploration of the parameter space and exploitation of the previous samples. The more explorative an algorithm is, the more likely it will be robust against local minima, and the starting random seeds, the more exploitative it is the faster it will reduce the objective function. The amount of exploration and exploitation performed by the NA is largely controlled by the two control parameters ns (ns2) and nr.
nr = Number of Voronoi cells re-sampled at each iteration.
See the papers below for more details. As ns and nr are increased it will be more explorative and as they are decreased it will be more exploitative.
Our experience has shown that reasonable results can be obtained after enough iterations when ns is of the order of twice the number of dimensions, and nr is in the range 2 - ns/2. The number of iterations is more difficult to judge. I have used used 10 - 100 times ns. It is clear that the larger values of ns and nr (i.e. more explorative) require more iterations to concentrate sampling (especially in high dimensions). In our experiments we have found in all cases that the sampling does eventually concentrate quite well, even for larger nr values, although you might not want to wait that long !
My strategy is often that it is better to get good sampling slowly than poor sampling quickly. However, small values of nr do seem to concentrate sampling quickly, especially in small dimensional (e.g. d ~ 4) problems. These issues are also discussed in Sambridge (1999a). I would be interested to hear of other peoples experiences (see contacting the author ).
Utility programsWhen the package is compiled the following general utility programs are put in the bin directory (source in src/utl).
nad2asc - reads in a nad file and writes out all models to an ascii file
asc2nad - goes the other way
splot - reads in a nad file and plots a scatter diagram
of all models (PS or X)
plot_model - a density plot of all velocity profiles produced by the NA search.
synrf - calculates a synthetic receiver function from an input velocity profile.
synrf_noise - Same as synrf but adds noise to the result.
nad2rfi - converts direct access NAD file to an ascii file containing
rfi2nad - performs opposite conversion to nad2rfi.
Further details can be found at the RFI utility programs page.
Running the receiver function inversion example
The example driver program provided with the NA-sampler package is for the inversion of receiver functions for crustal seismic structure. This is described further in the two NA papers cited below. Note that the driver program is not intended to be a general purpose receiver function code, but merely an example of using the NA routines on a realistic seismic problem. Changes to the code may be required if it is applied to other data sets (For example, the maximum number of layers, and waveforms may need adjusting etc.)
To run the RFI examples after compilation of the package (see above),
> cp na.in_examples/na.in_9_2 na.in
This will run a short (3 iteration) test run with Ns=9 and Nr=2. The output should be similar to that in demo_results/na.sum_3_9_2. A similar set of steps is required to run one of the longer example inputs in directory na.in_examples. You could try editing na.in to adjust input parameters and repeating the test run.
The RFI subroutines use several input and output files. All of which are in sub-directories of data/rfi_files. The main one is rfi.in, which contains parameters associated with the input data and the names of the other I/O files.
Two other input files are used. In the default setup these are
rec_func - in sub-directory ORF, which contains the one or more input (stacked) receiver functions.
rec_models - in sub-directory NA_MDL, which contains all models generated at each iteration of the NA in a convenient format.
rec_func - in sub-directory NA_SRF, which contains the predicted receiver function from the best fit model.
Conditions of use
Interested parties should contact the author directly.
2. Due acknowledgment must be made of the use of the NA routines in research reports or publications. Whenever such reports are released for public access, a copy should be forwarded to the author.
3. The NA routines may only be used for not-for-profit research and development, unless it has been agreed in advance with the author in writing.
Contacting the author
The purpose of this documentation and the associated README files, is to reduce the likelihood of running into problems. I would be grateful if you do not contact me with queries on the installation, or running, of the code unless you have thoroughly looked for the answer yourself. (Note that either these HTML pages, the README files, example files and the source code itself may provide an answer). Note that for detailed queries on the RFI example it may be best to contact Dr. T. Shibutani, (email@example.com) who wrote the original versions of the RFI code and plot programs. His contribution is much appreciated.
I am interested in obtaining feedback on people's experiences with this code. Especially to hear what type of problems it has been applied to. I expect that the code will be updated from time to time. Information on new developments will be posted on the NA homepage. If you develop your own web pages showing work where NA has been used, or publish papers then I would be happy to add a reference on the NA homepage.
Geophysical Inversion with a Neighbourhood Algorithm I: searching
a parameter space,
A related paper describes the Neighbourhood re-sampling algorithm (implemented in program NA-bayes) for calculating Bayesian integrals from a finite set of samples produced by the NA algorithm or any other search method (e.g. GA or SA etc.).
Geophysical Inversion with a Neighbourhood Algorithm II: appraising
Postscript and PDF versions of these papers can be found on the NA homepage.
Enquires to Malcolm Sambridge: firstname.lastname@example.org
Last modified July 2002.