extfieldsofstudy:
Physics
Computer Science

text:
Parallel Programming with Matrix Distributed Processing

Matrix Distributed Processing (MDP) is a C++ library for fast development of efficient parallel algorithms. It constitues the core of FermiQCD. MDP enables programmers to focus on algorithms, while parallelization is dealt with automatically and transparently. Here we present a brief overview of MDP and examples of applications in Computer Science (Cellular Automata), Engineering (PDE Solver) and Physics (Ising Model).

Introduction
Matrix Distributed Processing (MDP) [1][2] is a collection of classes and functions written in C++ for fast development of parallel algorithms such as solvers for partial differential equations, mesh-like algorithms, and various types of graph-based problems. These algorithms find frequent application in many sectors of physics, engineering, electronics and computational finance.
MDP includes: • a natural syntax for the algorithms that is transparent to the underlying parallelization; • parallelization based optimization algorithms implemented in MPI; • functions for linear algebra computations with support of a Maple-like syntax; • statistical functions; • fitting functions; • a Parallel SIMulator (PSIM) which allows MDP algorithms to run on single processor machines without MPI (uses fork for creating the processes and socket pairs for communication).
MDP was originally developed for and still constitutes the core of FermiQCD [3][4], a library for Lattice QCD computations. Lattice QCD is a Monte Carlo based numerical approach to the physics of composite particles made of quarks (e.g., protons and neutrons) and it is considered one of the most computationally intensive projects of modern physics.
FermiQCD, was developed by the author at the University of Southampton (UK) and Fermilab (Department of Energy) and it is currently used by other physics departments around the world [5][6] [7][8] [9].
While Lattice QCD is currently one of the main applications of MDP, its range of applicability is not limited to it as we show in the following examples.

Example: Parallel Game of Life
As a first example we show here how to write, in a few lines, a parallel program to find stable configurations for the Game Of Life. The game works as follows: • It is defined on a board (N × N cells) with periodic boundary conditions.
• The game consists of a series of iterations and, at each iteration, if a cell is dead and the number of beighbor cells who are alive is greater or equal 3, the cell dies; if a cell is alive and the number of beighbor cells who are alive is 2 or 3, the cell remains alive; otherwise the cell will be dead.
• Given a random starting configuration for the board, we want to determine its evolution and whether it reaches a stable configuration.
• Line 14 opens the parallel communication channels, line 37 closes them.
• Line 15 declares a 2D array containing the board size, L.
• Line 16 declares a lattice object (board) that represents the set of board sites and its topology. By default a mesh topology with toroid boundary conditions. Two arguments are passed to the constructor: the dimension, and size of the board. For each board site, it is possible to specify on which parallel process it is going to be allocated as well as alternative topologies.
• Line 17 declares a field of cells (C) on the board (the cells that live on it).
• Line 18 declares an auxiliary field that will be necessary for the computation.
• Line 19 declares a x a variable that represents a generic site of the board and will be used for looping over the sites and therefore the cells that live on them.
• Line 21 sets the initial configuration of the board by looping over all sites x and setting the corresponding cell C(x) to a random alive or dead. Each process loops only over the sites that are allocated locally by the process. Notice how MDP provides a local random number generator for each site of the board, board.random(x), which is vital in order to be able to reproduce results of stochastic algorithms.
• Lines 23-35 loop over the iterations of the game.
• Line 24 is performs a critical operation; it informs MDP that the value of the field C has changed and parallel communication is needed to synchronize those buffers that contain copies of non-local sites (synchronization).
• Lines 25-28 apply the rules of the game at each site. The new values of the states for the cells are stored in the auxiliary field newC. Each process loops over the local sites only. For each local site (x), x+0 represents the neighbor site when coordinate 0 is incremented by one, x-0 represents the neighbor site when coordinate 0 is decremented by one, x+1 represents the neighbor site when coordinate 1 is incremented by one, etc. 1 .
• Lines 29-33 perform two operations in a single parallel loop: the new states for the cells (newC) are copied back into C; the number of cells that have changed state are counted and the number is stored in diff. 1 This notation may appear bizzarre but it is nothing more than a sum ( x is a vector that represents a site on the board and the integer − → i represents a unit verson in direction i.
• Line 34 terminates the iteration if no cell has changed its state.
• Line 36 saves the board configuration in a file cells.dat.
Most of the parallel work is done by the constructor (mdp lattice) of the board which, from the board topology, determines how to partition it and determines the optimal communication patterns; by the constructor (mdp field) of the field of cells (C) which allocates the local cells and the buffers to store copies of non-local cells; and by the method update which performs communication to copy remote cell values into the local buffers.
We'll show in the next few examples that is equally easy to implement fields of any class of objects. In the next example we'll consider a field of matrices.

Example: Parallel PDE Solver
Consider here, as a different example that presents many similarities with the one above, the following Laplace equation: The initial conditions are ϕ initial (x) = 0. We will also assume that x i + L i = x i (torus topology). Solution: In order to solve eq. (1) we first discretize the Laplacian (∇ 2 = ∂ 2 0 + ∂ 2 1 + ∂ 2 2 ) and rewrite it as where µ is a unit vector in the discretized space in direction µ. Hence we solve for ϕ(x) and obtain the following recurrence relation The following is a typical MDP program that solves eq. (1) by recursively iterating eq. (4). Notice how the program is parallel but there are no explicit calls to communication functions: phi.save("field_phi.mdp"); // save field 24 mdp.close_wormholes(); // close communications 25 }

Notes:
• Line 04 declares the size of the box used to approximate the space L = {L 0 , L 1 , L 2 }.
• Line 05 declares a 3-dimensional lattice, called space, on the box L. MDP supports up to 10dimensional lattices. By default, a lattice object is a mesh with torus topology.
• Line 06 declares a site variable site x that will be used to loop over the lattice.
• Line 07 declares a field of 2 × 2 matrices, called phi, over the lattice space.
• Lines 11-12 initialize the field phi. Notice that phi is distributed over the parallel processes and forallsites is a parallel loop.
• Line 16 loops over all sites in parallel.
• Line 23 saves the field. Notice that all fields, including user-defined ones, inherit save and load methods from the basic mdp field class.

Example: Parallel Ising Model
As one more example of usage of MDP we report here a simple program for the Ising model. • Lines 3-4 declare a 1D lattice of 100 points (line).
• Line 5 declares a field of integers (spin) on the lattice.
• Line 7 sets the total magnetization M for this initial spin configuration.
• Line 9 sets the initial configuration: all field variables equal to 1.
• Line 14 computes the energy variation (dE) of each Markov Chain Monte Carlo (MCMC) step.
• Lines 15-16 perform the Monte Carlo accept-reject. If a change is accepted the spin at site x is flipped and the total magnetization M changes (line 16).
Note how at each MCMC step, first the code tries to flip the spins at even locations then, after it updates the lattice sites, it tries to flip the spins at odd locations (line 13). This guarantees computation results are independent on parallelization of the lattice line.
Since this even-odd distinction is common in many lattice algorithms, MDP stores all even lattice sites and all odd lattice sites close together. This speeds up loops over one of the two subsets and also speeds up communication. In fact, in this example, we are able to limit the synchronization (update) to the site of a given parity (line 18).

Linear Algebra and Other Tools
MDP includes a Linear Algebra package and other tools. Some of the most important classes are: • class mdp real, that should be used in place of float or double; • class mdp complex, for complex numbers; • class mdp array, for vectors and/or multidimensional tensors; • class mdp matrix, for any kind of complex rectangular matrix; • class mdp measure, for error propagation; • class mdp jackboot, a container for jackknife and bootstrap algorithms.
The most notable difference between our linear algebra package and other existing packages is its natural syntax.
For example: A=Random.SU (7); B=exp(A)+inv(A)*hermitian(A)+5; reads like A and B are matrices A is a random SU (7) matrix B = e A + A −1 A H + 5 · 1 Note that each matrix can be resized at will and is resized automatically when a value is assigned. MDP includes functions for fitting such as the Levenberger-Marquardt algorithm.

Lattice, Site, and Field
An mdp lattice is the class that describes the space on which fields are defined; it stores the topology of the space (by default that of a torus in d dimensions) and information about partioning of the space over the parallel processes. In MDP, a lattice is a graph, defined as a collection of points (lattice sites) connected by links (they specify the topology). Each site is uniquely mapped to one of the parallel processes.
The only restriction is that the graph must have a degree less than 20. From now on we will assume the default topology of a torus; therefore the lattice should be thought of as a mesh in d ≤ 10 dimensions.
The constructor class mdp lattice determines on which process to store each site, determines the neighbors of each site, and determines the sizes of the buffers where each process keeps copies of those sites that are non-local but are neighbors of the local sites.
The constructor also allocates a parallel random number generator so that each site of the lattice has its own independent random number generator. This is important for parallel Monte Carlo applications of MDP and ensures reproducibility of computations on different architectures. On each lattice it is possible to allocate fields. Some fields are built-in, for example This code declares a 10 × 10 × 10 lattice (cube) and a field (psi), that lives on the cube. The site variables of psi, psi(x) belong to class S (assuming x is an mdp site on the cube).
User-defined fields can be saved, loaded, and synchronized psi.save("filename"); psi.load("filename"); psi.update(); Synchronization means that all processes will perform MPI communications to make sure all buffers that contain copies of non-local site variables are updated with the proper values. The method update should be called immediately after the local site variables of a field have been changed.
Once a field object is declared in the field constructor, each process dynamically allocates memory for the buffers that store the copies of those sites that are non-local, but are neighbors of the local sites. These buffers are created in such as way to ensure optimal communication patterns.
Every time a field changes, for example in a parallel loop such as the program notifies the field that its values have been changed by calling phi.update(); The method update performs all required communication to copy site variables that need to be synchronized between each couple of overlapping processes.
Lattice sites are represented by objects of class mdp site. Site objects can be looped over in parallel loops (such as forallsites) but it also possible to explicitly address a specific site by specifying the site coordinates. Obviously, only the process that stores a site locally should address a specific site. Class mdp site has methods to check if a site is local, if it is non-local but a local copy is present, and which process stores the site locally.

Optimal Communication Patterns
In MDP, the lattice objects, according to the lattice topology and the parallel partitioning, determine the optimal way to store site variables in memory and performing parallel communication. This information is then used by the field method update that performs the synchronization of the field variables.
Note that MDP does not attempt to overlap computation and communication. By "optimal communication" pattern we mean that, under the assumptions below, the method update minimizes network traffic and data copies.
Current optimizations are based on the assumption that each processing node has one and only one network card and that the network is isotropic (latency and bandwidth for each couple of nodes is the same). This assumption is generally true for Ethernet and Myrinet clusters.
Communications are optimal in the sense that: • Each process retrieves all non-local site variables in a single send/recv for each process that contains sites which are neighbors of local sites.
• Two processes that do not store neighbor sites do not communicate with each other.
• No process is involved in more than a single send and a single receive at one time.
• Each process stores close in memory those copies of non-local site variables which are local to the same process. In this way synchronization does not require the use of additional buffers receiving buffers.
This technque is particularly efficient for algorithms that only require next-neighbor communication and run on a all-to-all network topology such as Ethernet or Myrinet.
It is possible, in principle, to change the above communication patterns to optimize communication for other network topologies.
Although communications are currently based on MPI, they do not make use of communication tags. It is therefore possible, in principle, to speed-up communication by using a faster tagless and bufferless protocol such as Myricom GM.
Our communication patterns have the effect of making communication almost insensitive to network latency, and communication speed is dominated by network bandwidth. Benchmarks are very much application dependent since parallel efficiency is greatly affected by the lattice size, by the amount of computation performed per site, processor speed and type of interconnection. In many typical applications, like the one described in the preceding example, the drop in efficiency is less than 10% up to 8 nodes (processes) and less than 20% up to 32 (our tests are usually performed on a cluster of Pentium 4 PCs (2.2GHz) running Linux and connected by Myrinet).

MDP and PSIM
For portability reasons, MDP is based on MPI. Nevertheless it is desirable to be able to run, test and debug MDP programs on a single node (with single or multi processor architecture) without having to install MPI. The latest version of MDP includes a Parallel SIMulator (PSIM). Despite the name, this is not quite a simulator but an emulator, i.e., a message passing library that uses local (unix/posix) socket pairs. PSIM is Objected Oriented and is not based on MPI.
When compiling with PSIM, the parallel processes are created at start-up by forking. The number of parallel processes is specified at runtime by passing the following command line argument to any MDP executable program, -PSIM_NPROCS=4 (this makes 4 parallel processes). PSIM also creates a communication log that can be used for debugging. MDP with PSIM has been tested on Linux, Mac and Windows (with cygwin).
For single processor node, using PSIM does not introduce any speed-up, but, for a small number of processes (2-16) it does not slow down the code either. For multi-process shared memory architecture, parallelization of PSIM should produce a speed-up comparable with MPI. We have not been yet performed such tests.
Moreover, PSIM should perform well on openMosix clusters as soon as openMosix starts supporting migratable sockets since all communication between the parallel processes will be done by the operating system. Unfortunately openMosix does not support migratable sockets yet.

Conclusions
MDP is an easy, powerful, and reliable tool for developing efficient parallel numerical applications. We have shown here examples in Computer Science (Cellular Automata), Mathematics (PDE solver) and Physics (Ising model).
MDP enables the programmer to focus on algorithm design while parallelization issues are dealt with automatically in a way transparent to the programmer.
The underlying communication functions are written in MPI but it is possible also to compile it without MPI using the built-in PSIM emulator which enables one to run parallel processes on a single node and/or single processor architectures, such as a PC. This is useful for testing and debugging purposes.
All of the features here described are fully functional and have been tested in real-life applications such as FermiQCD, developed by the University of Southampton (UK) and Fermilab (Department of Energy).


Design and Implementation of Components in the Earth System Modeling Framework

The Earth System Modeling Framework is a component-based architecture for developing and assembling climate and related models. A virtual machine underlies the component-level constructs in ESMF, providing both a foundation for performance portability and mechanisms for resource allocation and component sequencing.

Introduction
Component-based design is a natural fit for climate modeling. At its simplest, a software component is a code that has a standard calling interface and behavior and a coherent function (see, for example, http://www.corba.org and http://www.cca-forum.org). Components are ideally suited for the representation of a system comprised of a set of substantial, distinct and interacting domains, such as atmosphere, land, sea ice, and ocean. Further, since Earth system domains are often studied and modeled as collections of processes (radiation and chemistry in an atmosphere, for example), it is convenient to model climate applications as a hierarchy of nested components.
Component-based software is also well suited for the manner in which climate models are developed and used. The multiple domains and processes in a model are usually developed as separate codes by specialists. The creation of a viable climate application requires the integration, testing, and tuning of the pieces, a scientifically and technically formidable task. When each piece is represented as a component with a standard interface and behavior, that integration, at least at the technical level, is more straightforward. Similarly, standard interfaces help to foster interoperability of components, and the use of components in different contexts. This is a primary concern for modelers, since they are motivated to explore and maintain alternative versions of algorithms (such as dynamical cores), whole physical domains (such as oceans), parametrizations (such as convection schemes), and configurations (such as "standalone" versions of physical domains).
The Earth System Modeling Framework (ESMF; Hill et al., 2004; see http://www.esmf.ucar.edu) project is a large, multi-agency collaboration whose goal is to develop common modeling infrastructure and deploy it in climate, weather, and data assimilation applications. It is currently in its third year of development. The climate applications that are evaluating and adopting ESMF include the Community Climate System Model (CCSM; Boville and Gent, 1998; see http://www.ccsm.ucar.edu), the National Oceanic and Atmospheric Administration (NOAA) Geophysical Fluid Dynamics Laboratory (GFDL) models (see http://www.gfdl.noaa.gov/~fms), the new NASA Goddard Earth Observing System (GEOS-5) model, and the MITgcm (see http://www.mitgcm.org). The Strategic Plan for the U.S. Climate Change Science Program (CCSP), a report by the CCSP and the Subcommittee on Global Change Research (July 2003), recognizes ESMF as the realization of common modeling infrastructure for the climate community.
The ESMF software consists of an infrastructure of utilities and data structures for building model components, and a superstructure for combining them into applications. In this paper we concern ourselves with the superstructure, and how ESMF components and related constructs are designed and implemented.

The ESMF Superstructure
The ESMF superstructure defines an architecture for assembling Earth system applications from modeling components. A component may be defined in terms of the physical domain or process that it represents, such as an atmosphere or sea ice model. It may also be defined in terms of a computational or scientific function, such as an I/O or diagnostic system. Climate models often require that such components be coupled together to create an application. By coupling we mean the data transformations and, on parallel computing systems, data transfers, that are necessary to allow data from one component to be utilized by another. ESMF offers a suite of regridding methods and other tools to simplify the organization and execution of intercomponent data exchanges.
All ESMF methods can be called from Fortran 90, and some, including those at the component level, can be called from C++. The framework is implemented in a combination of these languages. Since most climate modelers work in Fortran, in this paper we focus on Fortran bindings. ESMF currently runs only in single program multiple datastream (SPMD) mode, although execution of components can be either concurrent or sequential. It is supported on IBM, SGI, Compaq, Mac OS X, Cray X1, and many Linux platforms.
The ESMF is implemented using an object-oriented design strategy. Superstructure and infrastructure are organized as sets of classes, which are data structures with associated methods. The main classes in the superstructure are component, state, and application driver.

COMPONENT
An ESMF component has two parts: one that is supplied by the ESMF and one that is supplied by the user.
The part that is supplied by the framework is in the form of a predefined Fortran derived type. ESMF provides two kinds of components: a gridded component (Grid-Comp) and a coupler component (CplComp). A gridded component represents a physical domain in which data are associated with one or more grids; for example, a dynamical core. A coupler component arranges and executes data transformations and transfers between one or more gridded components. All gridded components and coupler components possess initialize, run, and finalize methods with standard interfaces. These methods can be multiphase. Components also have create and destroy methods.
The second part of an ESMF component is the user code that will perform the computational work. The user code must be split up so that it too contains clear initialize, run, and finalize subroutines that match standard interfaces. The arguments must be ESMF data structures, although these data structures need not be used in the body of the code. Users set entry points within their code so that their initialize, run, and finalize subroutines are callable by the framework. In practice, setting entry points means that a special method called set services 1 is written by the user for every component. Within set services there are calls to ESMF set entry point methods that associate the name of a user's Fortran subroutine with the corresponding standard component operation. For example, a user might create an ESMF gridded component representing an ocean model known as the Parallel Ocean Program (POP) by making the following calls from a driver or parent component: If the Fortran subroutine names of the user's initialize, run, and finalize methods were popInit, popRun, and popFinal, respectively, the set services method would contain the following fragment:

STATE
ESMF components exchange information with other components only through states. A State is a Fortran derived type that can contain other ESMF types representing fields, bundles of fields on the same grid, and arrays. It can also contain other states. A gridded component is associated with two states: an import state and an export state. Its import state holds the data that it receives from other gridded components. Its export state contains data that it can make available to other gridded components.

APPLICATION DRIVER
The application driver (AppDriver) is a small, generic driver program that contains the "main" routine for an ESMF application. It can be thought of as the "cap" for a hierarchical application. Infrastructure classes closely associated with the ESMF superstructure include the Clock class and the virtual machine (VM) class. An ESMF clock contains information about the start time, stop time, and time-step of a component, based on a variety of supported calendars. The time-step may be changed during component execution. The clock may be associated with multiple alarms that "ring" for periodic or unique events. The VM controls computational resource allocation and is described in detail in Section The ESMF Virtual Machine and Component Parallelism.
A very simple ESMF coupled application might involve an AppDriver, a parent gridded component, two or more child gridded components that require an intercomponent data exchange, and a coupler component. Calls cascade so that when, for example, the initialize routine of a parent component is called, it in turn calls the initialize routines of all its children. Figure 1 illustrates the structure of such an application.

Application Example
The newly developed NASA GEOS-5 atmospheric general circulation model (AGCM) provides a more compre-hensive example of how ESMF is being used to structure climate and related models. GEOS-5 is the atmospheric component of a variety of applications at the NASA Global Modeling and Assimilation Office (GMAO). It will be used for research in satellite data utilization; as part of the GMAO atmospheric data assimilation system; for weather, subseasonal, and seasonal to interannual forecasting; for atmospheric chemistry studies; carbon cycle research; and research on ocean-atmosphere and atmosphere-land interactions. In 2006, GEOS-5 is planned for use in a satellite-era reanalysis in support of the CCSP. Figure 2 shows GEOS-5 science components connected for a standalone configuration of the coupled atmosphere-land system. Each box in the GEOS-5 diagram is an ESMF gridded component. GEOS-5 adopts the hierarchical topology that is natural under ESMF. The main computational modules are located in the leaves of the hierarchy (finite volume dynamical core, catchment, infrared, etc.). These are connected through composite components (radiation, dynamics, physics, etc.), which implement higher levels of integration. A single composite component (AGCM) integrates the entire standalone GCM. At the same level is an I/O component, which handles the diagnostic history interface. The whole system is capped by a GEOS-5 gridded component and an application driver.
One of the advantages of the hierarchical structure is that branches of the tree can be easily pruned and capped to form more limited standalone applications. For example, the GEOS-5 single column model is constructed by simply connecting physics to a new application "cap" and discarding the dynamics branch, but retaining history diagnostics. Another advantage of the architecture is that standard interfaces at each level make it technically straightforward to swap in new components; for example, a new radiation module.

Component Implementation
The connection of specific subroutines to the initialize, run, and finalize entry points is implemented using a virtual function table. The function table code is implemented in C++ but is callable from standard Fortran 90. For each entry this table holds the name of the entry, an enum that defines the types of arguments to be passed to the subroutine, an array of argument pointers declared as type (void *), and a pointer to the subroutine to be called.
We can store a pointer to the subroutine to be called because Fortran 90 allows subroutines or function names to be passed as arguments. While the Fortran 90 language does not mandate how subroutines are passed when calling into C or C++ code, all compilers we have run on pass the subroutine as a pointer to the start of the execut-able code. This is stored in the virtual function table as void *func, and dispatched as (*func)(arg1, arg2, ...). The framework uses this feature to populate the virtual function table and to dispatch subroutines from the table. Connections are defined at run-time and not compile time. In the future this will allow components to be dynamically replaced based on the evolution of a simulation, for example.
When a component method such as ESMF_Grid-CompRun() is called, control is transferred to the framework before user code is executed. The framework determines which execution threads should be active in the component, may optionally validate arguments, may optionally supply default values for any missing arguments, and then calls through the virtual function table to execute the code for the requested method. This isolates the components from any language dependence; components written in Fortran 90, C, or C++ can call subcomponents which are implemented in Fortran 90, C, or C++ without change.
A pointer to the virtual function table is stored as an entry in standard ESMF component derived types as an opaque 4-or 8-byte integer. All subroutines involving the function table simply require that the Fortran 90 interface has the component type as one of the arguments. The opaque entry is then passed to a C++ subroutine which uses it as the this pointer in the virtual function table class method.
We find that ESMF components work well with the module concept in Fortran 90, especially for controlling the scope of symbols. Module data and methods can be private to the module. The parent component includes the child code with the Fortran use statement. Only the set services subroutine name must be a public symbol, so it can be referenced by the parent component. All other module methods, even those which implement initialize, run, and finalize, can be private. Applications are linked into a single executable, which simplifies the interaction with the various batch systems common on high-end supercomputers, and requires no run-time dynamic loading of user code.

The ESMF Virtual Machine and Component Parallelism
Components are the building blocks of any ESMF application. We extended the component concept to include parallelism by making ESMF components the very units of parallel execution. In ESMF, components offer both data and task parallelism, and provide powerful concepts that aid the computational scientist to write highly efficient and scalable code.
The management of the available computational resources maps quite naturally onto the component hierarchy. When a component is created it obtains a set of resources from its parent component. In turn, when the component creates children of its own, it divides up its resources and provides them to its children. Critical in the design of the parallel aspects of components was to guarantee that components remain self-contained and could be written without detailed knowledge of the par-ent component code. This has been achieved by running each component in a separate parallel execution context. 2 In order to provide the desired transparent interface to the parallel execution environment, ESMF components utilize the ESMF VM class. VMs are not new; the idea of a generic machine representation dates back to the 1960s and has appeared in a wide diversity of software packages and languages, including early IBM systems (Adair et al., 1966), Java (see http://java.sun.com), and, in the high performance computing realm, Parallel Virtual Machine (PVM; Geist et al., 1994).
The ESMF VM's interaction with the rest of the framework can be divided into two parts, as follows.
1. The execution engine provides separate parallel execution contexts or VMs for every component.  The VM abstracts away many details of the underlying execution environment. Each VM can be viewed as a generic representation of hardware and system software resources available to each individual ESMF component. There is exactly one VM object associated with each component in an ESMF application. The VM handles resource management tasks and provides a topological description of the underlying configuration of compute resources accessible to the component. The communications interface provided by the VM utilizes this information to offer the best possible performance.

USING THE VIRTUAL MACHINE
At the beginning of every ESMF application, before making any framework calls, the user must make a call to an ESMF_Initialize method; likewise, at the end, the user must call ESMF_Finalize. These calls set up and shut down the whole framework, including the VMs. By default, ESMF is based on MPI-1 (see http://www.mpiforum.org) and, as such, the main program is launched like any other SPMD message passing interface (MPI) application. ESMF does not directly interact with batch or queue management systems, and obtains its computational resources through the system-dependent MPI startup facility, such as mpirun, prun, poe, etc. During the initialization phase, a default global VM is created. The user code can gain access to this VM by either providing an optional VM argument to ESMF_ Initialize() or by calling ESMF_VMGetGlobal() anywhere in the ESMF application. The default global VM is equivalent to MPI_COMM_WORLD and provides the parallel execution context for the main program. Although there is no component associated with the main program, the default global VM can be seen as the parent of all VMs of the ESMF application.
The basic elements of a VM are persistent execution threads (PETs). PETs are equivalent to POSIX threads and have a lifetime of at least that of the associated component. By their nature, PETs provide an elementary OS unit of execution, associated with the virtual address space of the POSIX process in which the PET is running. Furthermore, each VM contains a mapping of its PETs onto the set of unique hardware processing elements (PEs) that were discovered during the framework initialization phase.
The PEs in a VM correspond to the smallest processing units that are recognized by the operating system. The VM layer keeps internal information about the topology of the PEs by noting whether they are part of a single CPU, lie within a single system image, or are connected via an interconnection fabric. All VMs in an ESMF application share this hardware map.
ESMF-level threading and user-level threading are two special features of the ESMF's VM implementation.
ESMF-level threading provides coarse-grained parallelism by allowing multiple PETs to run within the same POSIX process. A component must be completely threadsafe to utilize this type of threading but can expect major reduction in its communication cost between threaded PETs. The VM communication API is completely transparent with respect to ESMF-level threading. User-level threading, on the other hand, does not require the entire component to be thread-safe. A typical user-level threading approach would employ OpenMP for fine-grained loop-parallelization. The VM assists user-level threading by controlling resource management, and allowing single PETs to be associated with multiple PEs. When during the execution of a component a region is encountered that profits from user-level threading, the user code can inquire how many PEs are associated with a particular PET and spawn this number of temporary threads without risking oversubscription of the available PEs. The decision to use ESMF-level and/or user-level threading is made on a per component basis, allowing a component writer to choose the best approach for the particular problem.

COMPONENT-VIRTUAL MACHINE INTERACTIONS
There are four distinct phases in the life cycle of an ESMF component during which it interacts with the execution engine of the VM.
1. Component Creation. When a parent component or driver creates a child component, by either calling ESMF_GridCompCreate() or ESMF_Cpl-CompCreate(), it needs to provide a reference to its own VM. It is in the component creation call that a parent indicates which of its compute resources it wants to give to its child. This is done by adding the optional parameter petList to the arguments of the creation call. A petList is a list of parent PETs that are to be given to the child. The default, when no petList argument is provided, is to give all parent PETs to the child. 2. Component SetServices. As described in Section 2, The ESMF Superstructure, for each component the user writes a set services routine containing calls to set entry point methods. These link the initialize, run, and finalize subroutines of the user's code with the associated ESMF component derived type. In addition to setting entry points, during the child's set services method, the child component can also set certain properties of its VM. For example, when the call ESMF_GridComp-SetVMMaxPEs(pop, rc=rc)is executed within set services it will associate as many PEs as possible for each child component PET. It is important to note that the set entry point code is executed within the context of the parent VM, not the child VM (which does not yet exist). This means that users should not add additional operations such as data allocations to the SetServices method because they will not be associated with the correct VM. On return from the set services call the framework starts up the child VM and puts it on hold. 3. Call of standard component method. The child VM becomes active and executes the initialize, run, or finalize component method. When the end of the routine is reached the VM is placed back on hold, waiting for the next invocation, and control is handed back to the parent VM. 4. Component Destroy. The child VM is shut down and all resources are released.

CONCURRENT AND SEQUENTIAL EXECUTION OF COMPONENTS
By default, ESMF components run sequentially in the order in which they are called by their parent component or driver. ESMF also provides a very easy to use concurrent component model. The approach is an extension of the non-blocking concept well known from MPI's communication interface. In this model the run sequence (item 3 in the previous list) is split into two separate phases: (a) Non-blocking call of registered component method. The child VM becomes active and executes the registered component routine within the child VM. When the end of the registered routine is reached the VM is placed back on hold, waiting for the next invocation. The parent VM does not wait until the child method has reached the end but continues execution in its own context immediately. (b) Parent component calls wait. In order to synchronize its children, the parent component uses ESMF_GridCompWait() or ESMF_CplComp-Wait() to wait for the completion of a previously invoked child component method. The wait call will block all parent PETs until the child component's method has reached its end and the child VM has been placed on hold.
Within ESMF's concurrent component model, it is the parent's prerogative to decide if and which child components it will run concurrently. The child component code is unaffected by this. As with many other aspects, ESMF does not ensure that it is semantically correct to run certain components concurrently. It is up to the application writer to decide what makes sense and what does not.

Performance Overhead
At its relatively early stage of development, ESMF has just begun to perform exhaustive performance tests and optimizations for low-level communications, regridding methods, halo operations, and other parts of the framework. However, it has been a high priority from the start of the ESMF project to ensure that the basic component architecture of the framework is not associated with substantial performance overhead. In this section we describe results of initial performance evaluations. Table 1 shows the performance overhead associated with calling component methods through a virtual function table. Tests were run on the Compaq SC45 at NASA Goddard Space Flight Center. The "w/o ESMF threads" columns show the overhead for the case where ESMF threading is not enabled (although the user may choose to implement threading using OpenMP on their own). Timings are shown with and without a barrier synchronization upon completion of the component call. 3 We have also examined the impact of ESMF on the performance of a real code, albeit one that, unlike climate models, emphasizes efficiency of initialization as well as efficiency of run-time. The National Centers for Environmental Prediction (NCEP) Spectral Statistical Interpolation (SSI) package (Parrish and Derber, 1992) is a three-dimensional variational analysis of observations used by the National Weather Service in their data assimilation system to initialize their global atmospheric model. The physical domain of SSI is the global atmosphere from the surface to the stratopause. Instead of physics parametrizations, the SSI is comprised of forward model elements and their adjoints. Some major forward model elements are radiative transfer algorithms for each satellite instrument, convective and large-scale precipitation, spectral transform, grid interpolation, the balance equation, and the divergence tendency equation. The analysis minimizes a combination of fits to observations, fits to model background, and a set of dynamical constraints. The minimization is performed in spectral space. Computation of forward models and their adjoints are required at every iteration, and both initialization and run routines are executed at every cycle. One cycle of initialize and run was performed for this experiment.
The version of the SSI code tested with ESMF used only the superstructure classes. It was not coupled with any other ESMF components nor did it use any part of the ESMF infrastructure layer.
In addition to the general component overhead, the SSI code with ESMF also had algorithmic overhead, due to the lack of ESMF support for spectral fields at the time of the experiment. The result was that there were three extra data conversions in the ESMF code. The performance analysis was run on an IBM Power 3 cluster at the National Center for Atmospheric Research (NCAR).  Figure 3 shows the total ESMF overhead including the extra data conversion from the spectral fields into Gaussian grids and vice versa. The overhead was about 4 seconds (slightly increasing with increasing number of processors), or 0.26% for eight processors to 1.82% for 128 processors. For the SSI application, this timing difference is insignificant.

Conclusions
The ESMF hierarchical, component-based architecture is natural for the construction of climate and related applications. ESMF components are linked to a powerful VM construct which offers integrated parallelization and a generic representation of the high performance computing hardware and software environment. Initial performance tests indicate that the architecture will not impose a prohibitive performance overhead on applications. The ESMF project will continue to evolve and optimize its software in close collaboration with the Earth system modeling community.

ACKNOWLEDGMENTS
The ESMF Joint Specification Team defined the framework architecture. Thanks to the many members of this multi-agency team for their thoughtful advice, patience, and persistence. The NASA Earth Science Technology Office (under Cooperative Agreement CAN-00OES-01), the National Science Foundation, and the Department of Defense High Performance Computing Modernization Program support the research and development of ESMF.

AUTHOR BIOGRAPHIES
Nancy Collins is one of the lead architects of the Earth System Modeling Framework. She designed and implemented the component layer in ESMF, along with many other aspects of the framework. Previously she worked at IBM Research in Yorktown Heights, NY, in the Computer Science Department, where she specialized in the development of software for scientific visualization on parallel systems. She holds a Bachelor of Science in physics and chemistry from Whitworth College and a Masters of Science in computer science from the University of Colorado, Boulder.
Gerhard Theurich is an employee of Silicon Graphics Incorporated (Professional Services) and is currently part of the ESMF development team. His contributions to date focus on the low-level communications and execution layer of the ESMF. Before joining the ESMF development team, he was part of the applications team at the NASA Center for Computational Sciences (NCCS) at the Goddard Space Flight Center. His current interests include performance portability of parallel applications and multiparadigm parallel computing. He received a degree in Diplom Physik from the University of Erlangen, Germany, a Ph.D. in computational material science from the University of California, Santa Barbara, and worked as a postdoctoral associate at the University of Pennsylvania. During this time he developed and contributed to several electronic structure and quantum chemistry codes.
Cecelia DeLuca is currently the manager of the Earth System Modeling Framework implementation team at the NCAR. In this position, she is responsible for translating the needs of a broad scientific community and innovations in computer science into a robust infrastructure for climate and weather simulation. Previously she worked as a lead developer of a framework for real time signal processing at Massachusetts Institute of Technology Lincoln Laboratory. She is experienced with architecting large geophysical codes, assessing and selecting computational advances for use in production models, and using a systematic software process, including requirements collection. Ms. DeLuca holds a Masters of Science in engineering from Boston University and a Masters of Science in meteorology from Massachusetts Institute of Technology.
Max Suarez is a Meteorologist at Goddard Space Flight Center. He received a B.S. and M.E. in engineering science from the University of Florida, and an M.A. and Ph.D. in geophysical fluid dynamics from Princeton University. His interests include large-scale atmosphere/ ocean interactions, climate modeling, numerical methods, parametrization of subgrid-scale processes in atmospheric models, maintenance of the atmospheric general circulation, and climate sensitivity.
Atanas Trayanov is senior computational analyst at the Global Modeling and Assimilation Office, NASA Goddard Space Flight Center. He received an B.S and M.S. in physics in 1979 from Sofia University, Sofia, Bulgaria, and a Ph.D. in physics in 1988 from the Bulgarian Academy of Sciences, Sofia, Bulgaria. His interests include climate modeling, high performance computing, and parallel algorithms.
V. Balaji is the Head of the Modeling Systems Group at the NOAA Geophysical Fluid Dynamics Laboratory. He received a M.Sc. in physics from the Indian Institute of Technology and a Ph.D. in physics from Ohio State University. His interests include climate simulation, high performance computational systems and algorithms, and collaborative, standardized approaches to modeling.
Peggy Li received a B.S. and an M.S. in electrical engineering from the National Taiwan University, Taiwan, and a Ph.D. in computer science from the California Institute of Technology in 1986. She was the technical group supervisor for the Parallel and Distributed Application Technology Group (formally Automation and Scheduling Technology) in the Applications Development Section from 1996 to 1999. She has led several R&D tasks in the parallel and distributed computing area since she joined the Jet Propulsion Laboratory in 1989, including the development of a distributed discrete event simulation framework, a remote interactive visualization and analysis (RIVA) system for planetary data visualization. She is currently leading a parallel volume rendering project (ParVox) and a Fault Detection, Isolation, and Recovery (FDIR) Server task for TMOD Technology. Her research interests include scientific visualization, parallel rendering, and distributed computing.
Weiyu Yang is a research scientist in the NOAA NCEP. He received a B.S. and M.S. in atmospheric science from Peking University, China, and a Ph.D. in atmospheric science from the Institute of Atmospheric Physics, Chinese Academy of Sciences in 1989. He is the lead developer of the ESMF project at NCEP, developing the ESMF versions of the NCEP operational data assimilation systems, the global forecast model, and the ensemble forecast system. He designed the special parallel structure of the current NCEP operational data analysis system to obtain the best performance to fit the restrictive operational time requirement. Previously he worked at the Supercomputer Computational Research Institute, Florida State University, where he specialized in the development of the adjoint models of the NASA GEOS1 model and the NASA SISL model and researched on related data assimilation issues.
Chris Hill is a researcher in the Department of Earth, Atmospheric and Planetary Sciences at the Massachusetts Institute of Technology. His research is directed at applying advanced computing technology to the study of planetary scale circulation. He is particularly interested in how cutting edge technologies become effective research tools and how algorithms, hardware and software combine to make this possible and profitable. He received a B.S. in physics from Imperial College, London.
Arlindo da Silva is a Meteorologist in the Global Modeling and Assimilation Office at NASA Goddard Space Flight Center. He received a B.S. and M.S. in physics from Catholic University of Rio de Janeiro, Brazil and a Ph.D. in meteorology from the Massachusetts Institute of Technology. His current research interests include techniques for global atmospheric data assimilation, physicalspace analysis systems, error covariance modeling, bias estimation and correction, quality control, land-surface, precipitation and aerosol data assimilation, and efficient methods for assimilation of remotely sensed data. Other research interests not in the area of data assimilation include aerosol forcing of climate, hydrological cycle of the subtropics, estimation of fluxes of heat, momentum and fresh water over the global oceans for observational studies and forcing ocean models.


Upgrading Pulse Detection with Time Shift Properties Using Wavelets and Support Vector Machines

Current approaches in pulse detection use domain transformations so as to concentrate frequency related information that can be distinguishable from noise. In real cases we do not know when the pulse will begin, so we need a time search process in which time windows are scheduled and analysed. Each window can contain the pulsed signal (either complete or incomplete) and / or noise. In this paper a simple search process will be introduced, allowing the algorithm to process more information, upgrading the capabilities in terms of probability of detection (Pd) and probability of false alarm (Pfa).

INTRODUCTION
In previous work ( [5] and [3] ) we presented a new algorithm for pulse detection using a window where either a complete pulse or noise alone is found. After the input data is sampled, its discrete wavelet transform is computed to obtain wavelet detail coefficients for different scales (see [2], [4] and [6] ). This representation is then introduced into an optimum linear function (a Support Vector Machine) using all available information to make the decision about the presence or absence of such a pulse (see figure 1). This algorithm had a process gain 15 dB better than previous work using wavelet filters alone, i.e., without an optimum information processing decision function, described in [7] and [8]. Note that the output of the linear SVM ( [1] and [9] ) can be either "hard" (if it has two possible values only, 1 and -1) or smooth (when it is a floating point unconstrained number).

Figure 1. Pulse detection scheme
Usually an SVM gives a hard output, but here we have separated the last step: if the smooth output is greater than some threshold the data is positive-class, otherwise it is negative-class. This step loses information, which can be used for further processing, as, for example, class probability. Only in case the decision function is final should the threshold be applied. But in order to make the algorithm useful in the real world, this approach requires a mechanism to inspect the time domain such that whenever a pulse is emitted you can guarantee that at least one window will be filled with the complete pulse alone.

TIME-SEARCH PROCEDURE
Our time-search procedure is defined as follows. Given some time interval [t 0 , t f ], where amplitude versus time function is defined. Let T s be the time elapsed between two samples (usually defined by the discretization hardware resources). Let S be the number of samples needed by one complete known pulse. Let S t be the number of samples (of size T s ) between t 0 and t f . See figure 2.

Figure 2. Data windows in time-search procedure
Because of the definition of our approach described in [5], windows W must be of the same size as one complete pulse, that is, each window contains S consecutive samples. On interval [t 0 , t f ] we will obtain S t -S + 1 windows, such that each window will be defined as the interval [S k ,S k+S ]. Therefore, two consecutive windows will have the same data except for one sample, and all common samples in the second window will be shifted one position to the left with respect to the first window.

OUTPUT ANALYSIS
One of the main observations we did was how the window had to fit precisely on the pulse so as to let the SVM perform a correct positive classification. If we execute the complete-pulse optimum SVM over all windows we can see that only a very small subset (less than 7 elements) around the complete-pulse window will show some positive classification degree. All other time-windows behave much as a noise-alone window would about the positive/negative classification rate, regardless of the pulse energy. Therefore, given one single window, we define positive class when at least 99% of the pulse is inside the window, and negative class otherwise. On the other hand, we can train more SVMs to detect different incomplete pulses, that is, windows were only a portion of the pulse is present plus some noise (to fill the size of the window). Obviously, each one of these SVMs alone will produce poorer results than the complete-pulse SVM, as the detection procedure is applied on less information. Nevertheless, each one processes a different projection (through wavelet transform) of the pulsed signal, and so, they bear different information, as the experiments show. In figure 3 a three-fold multiple window processing is shown. For each window three linear detectors are executed to detect different portions of the pulse. For one pulse, three detectors in three different known windows should detect this event.

EXPERIMENTS
Our experiments had the following setup: chirp pulse (see [3]), 1024 samples size; mother wavelet Daubechies 5, using d 4 wavelet coefficients; white Gaussian noise with zero mean and deviation equals one. Three linear detectors (SVMs) were trained against noise alone, such as to detect the complete pulse (named 0-shift), 11 noise samples plus 1013 pulse samples (named 11-shift), and a similar approach for 23 noise samples (named 23-shift) and having P fa = 10 -3 , as established in [3]. Obviously, incomplete-pulse classifiers have less information and will have worse classification rates than its complete-pulse counterpart, but it must still be considered as a source of discriminatory information when used on its own. For each generated pulse observation, we computed 3 similar Wavelet + Decisionfunction schemes applied to the complete pulse, incomplete 11-samples-shift pulse, and incomplete 23-samples-shift pulse respectively. For each one, the corresponding above-defined SVM was executed. In figure 2 we described how one processing window contains a subset of sampled input data. In that example, subset w k+2 would be defined as the complete-pulse window, and subsets w k+1 and w k would be defined as n-samples-shift pulse, having the first one lower n value than the second one. Using the three variables as previously defined, the covariance matrix is calculated. Table 1 shows the covariance matrix using 5000 observations on positive data, while table 2 shows similar features for noise-alone observations. As wavelets are linear filters, these matrices shall have equal coefficients using an infinite number of observations, regardless of the signal-to-noise ratio on the pulse or even the presence of signal.  Table 2. Covariance matrix of 5000 noise observations and three variables (0-shift, 11-shift and 23-shift).
Note that although the three variables are clearly not independent, they are not completely correlated, that is, all three variables hold some small portion of individual, unshared information. Therefore, as we are able to process this kind of data together efficiently using machine learning algorithms, we can be sure to obtain better results than any of the three sources alone. Moreover, the more incomplete pulse information sources we use, the better results we will obtain. As more variables are introduced in the model, new uncorrelated information will be harder to find, but nevertheless there will always be some degree of upgrade.

CONCLUSIONS
In this paper we have introduced a new time-search algorithm to allow the linear detector (described in previous papers) to work under real world conditions. This algorithm solves the detector constraints, but it also upgrades the whole system results. In previous work we saw how the application of machine learning models was limitless on the known signal pulse detection problem using wavelets. In this paper we have seen a new unlimited line of information processing: the use of additional Wavelet + Decision-function scheme applied to determined incomplete pulsed signals. This new set of features provides the final model with uncorrelated information, upgrading the classification rates. Our next step will be to generate a multiple source model using machine learning algorithms (most probably SVMs), to confirm the hints provided by the experiments shown in this paper. Further analysis is needed to determine how an increased number of processing units in one multiple-source decision function will upgrade probability of detection and / or probability of false alarm.


extfieldsofstudy:
Mathematics
Computer Science
Physics

text:
An ant-based algorithm for annular sorting

In this paper we describe a minimal model for annular sorting by Leptothorax ants. Simulation results are consistent with the structures observed in actual ant colonies.

Introduction
The ability of social insects to collectively solve problems has been well-studied and documented [1]. The behaviour of foraging ants, for example, has been abstracted to provide algorithmic solutions that are robust, distributed and flexible [3,4]. The particular behaviour that we will focus on is the clustering or sorting of ant corpses or larvae [2]. Abstract models of these behaviours have been successfully applied to, amongst other problems, numerical data analysis, data mining and graph partitioning [6]. In this paper, we focus on the task of brood sorting. This behaviour, when observed in Leptothorax unifasciatus [5], leads to the formation of a single cluster of offspring made up of concentric rings of brood items, with the youngest items (eggs and micro-larvae) being tightly packed at the centre, and successively larger larvae arranged in increasingly wider-spaced bands moving out from the centre of the cluster.
In this paper we investigate an intriguing hypothesis by Franks et al. concerning the biological mechanisms underlying annular sorting. In their article, the authors state that "The mechanism that the ants use to re-create these brood patterns when they move to a new nest is not fully known. Part of the mechanism may involve conditional probabilities of picking up and putting down each item which depend on each item's neighbours . . . The mechanisms that set the distance to an item's neighbour are unknown. They may be pheromones that the brood produce and which tend to diffuse over rather predicable distances . . . " [5] We have constructed a corresponding minimal model, using only stochastic ant behaviour and pheremone concentrations, that appears to account for the Figure 1: Annular sorting in a real Leptothorax colony (taken from [5], with permission). emergence of annular clusters of objects in simulated Leptothorax colonies. In Section 2 we present the background to the problem, before describing our model in Section 3. We conclude with a summary of our results and discussion of their implications.

Annular sorting
Franks et al. [5] carried out an observational biological study of the brood sorting behaviour of Leptothorax ants. The nesting behaviour of this species made them particularly suitable for study as they nest in single clusters in flat rock crevices, a situation that is easy to replicate and monitor in a laboratory. Photographs were taken of the ants' brood cluster and individual items in the cluster were classified, before a tessellation was applied to determine density and distance from the cluster centroid. This study showed that Leptothorax ants placed smaller brood items at the centre of the cluster with a greater density, forming an annular cluster. This arrangement reasserted itself when the ants were forced to migrate to a new nesting site, and proved to be an ongoing process. This structure is illustrated in Figure 1, where three different types of brood item are arranged in roughly-sorted concentric rings.
Wilson et al. have recently proposed the first model of "ant-like annular sorting" to simulate the behaviour of Leptothorax ants using minimalist robot and computer simulations [9]. Three models for annular sorting were presented: "Object clustering using objects of different size", "Extended differential pull-back" and " Combined leaky integrator. The first was run exclusively as a computer simulation, since modifying robots to allow them to move objects of different sizes proved to be too complex. Despite this, the computer simulation modelled physical robot behaviour faithfully, preserving the limitations of movement inherent in simple robots, and even going so far as to build in a 1% sensor error that matched the rate seen in the machines.
The first model explored the theory that annular sorting occurs solely due to the different sizes of the objects involved, Wilson et al. compared this to the manner in which muesli settles in transit, with smaller clusters falling to the bottom, leaving the larger ones on top. The simulation modelled agents who picked up the first object they encountered while unladen and deposited it at the next object they encountered. The results displayed a slight increase in the quality of clusters when there was a greater number of different object sizes; however, clusters tended to form at the edge of the area and were often "inside out", i.e. with larger objects at the centre surrounded by smaller objects. This lead Wilson et al. to conclude that merely using objects of different sizes did "not create a sufficient muesli effect".
The second and third models attempted to recreate ant brood clustering behaviour by assigning more complex behaviour to the ants. Wilson et al. hypothesised that ants were able to recognise inherent differences in larval growth, and created annular clusters by depositing different objects at different differences away from each other depending on size. Due to the limitations of the robots, this was implemented by having the agents reverse a distance depending on the kind of object carried before depositing it. Initial results with this approach were not good. As a result Wilson et al. proposed a third model, calling it the "Leaky integrator". This allowed the agents to have an adaptive amount of "pullback" that varied according to how many objects of the same kind had been encountered in the last n seconds. Initial tests with this system produced poor results but results improved when a genetic algorithm was used to select parameter values.

Our proposed model
We now propose an alternative model to explain the phenomenon of annular clustering. Ants are represented as "agents", and brood items of different sizes are represented by "objects". Agents and objects are spatially-distributed randomly on a two dimensional "board" of fixed dimensions.
Each object has a placement score; agents move randomly across the board, and when they collide with an object they calculate its placement score. This score is then used to probabilistically determine whether the agent should pick up the object and become laden. Laden agents carry objects around the board, and at every timestep they evaluate what placement score the carried object would have if it were to be deposited at the current point. This score is then used to probabilistically determine whether the object should be deposited.
Each category of object to be clustered posses three attributes: size, min- imum perimeter and maximum perimeter. Size is important, as objects may not overlap, the other two attributes are both functions of the object's size.
Minimum and maximum perimeter are used in the calculation of the object's placement score: When evaluating the placement score of an object, the agent counts how many nearby objects fall within the minimum perimeter (these count towards a weighted penalty), and how many fall within the maximum perimeter (these count towards a weighted bonus).
It is important to note that the maximum and minimum perimeters are only calculated for the object whose placement score is currently being evaluated. Agents take no account of whether a placement results in a good or bad score for neighbouring objects, thus small objects with smaller minimum perimeters will frequently be placed quite close to larger objects, resulting in a penalty score for those larger objects. As a result groups of smaller objects will tend to "force out" larger objects. It could be argued that this process implements a stronger version of the size-based muesli analogy proposed by Wilson et al. [9] As stated earlier, the minimum and maximum perimeter values are a function of an object's size (we describe the exact function shortly). The minimum perimeter of an object acts as a "repulsive" force, whereas the maximum perimeter acts as an "attractive" force. These coupled forces ensure that objects maintain a proportional "exclusion zone", around themselves, whilst ensuring the coherence of a single cluster around some centroid ("centre of gravity"). We now illustrate this with an example, depicted in Figure 2.
This example shows how larger objects are "forced" out of clusters by their proximity to other objects. The large object in Figure 2 has its minimum perimeter depicted as a dotted line, and its maximum perimeter as a solid line. In Figure 2(a), the object has a high-scoring placement score, as it has no objects within its minimum perimeter (which would attract a penalty), and several objects in its maximum perimeter (attracting a bonus score). However, if several objects are later deposited in the central object's immediate vicinity ( Figure  2(b), its placement score changes for the worse, as these objects contribute a Figure 3: Poor clustering caused by zero delay significant penalty. As a result, the next time this object is encountered by an agent it will be carried until it once again has a beneficial placement score (Figure 2(c), where it will be deposited. As we can see, this calculation of placement scores has the effect of moving larger objects away from clusters of smaller objects, whilst maintaining relative proximity.
Before embarking on a full implementation of our algorithm, we carried out test runs using a prototype system. This highlighted two, previously unforseen, effects that we deal with in the full algorithm. Early investigations showed that ants had a tendency to construct good annular clusters, but would then deconstruct them at the outer edge, as there was no termination criterion defined for the algorithm. We solved this problem by introducing the notion of "energy"; each ant starts with a fixed amount of energy, represented as an integer value, which is decremented every time the ant picks up an object. Once the agent's energy level reaches zero is is removed from the system, and the simulation terminates when there are no more agents left. This prevents the simulation from ending when objects are still being carried, but also has the beneficial effect of smoothly and gradually reducing the overall activity within the system as the simulation approaches termination. The biological validity of this approach is unclear, but it appears to make the clustering of outer objects much more realistic.
Another problem that was highlighted by the prototype algorithm was that of "chaos"; ants would frequently deposit an object and then immediately pick it up again, as they were still in contact with it. Conversely, ants would also frequently deposit objects after taking a single "step", particularly if the placement score was moderately good. This behaviour creates clusters in which overall placement was quite poor, as large objects were rarely removed from the centre of the cluster (Figure 3). In order to deal with this problem, a "cooling down" delay period was introduced; agents must wait a set number of steps after picking up or depositing an object before they carrying out any further placement calculations. This appeared to solve the problem; different (non-zero) values for the delay variable impacted only on the run-time of the simulation, and did not affect the quality of clusters generated. A delay value of 4 was chosen for what follows.

The algorithm
We now describe in detail our algorithm for annular sorting. Ants are modelled by agents, each of which has the following attributes: • Location (application representation) • Laden (true/false)

• Object (application representation of object carried if Laden == true)
• Energy (integer) • Delay (integer) Objects are modelled as spheres, and have a single attribute, from which their minimum and maximum perimeters are calculated: • Size (real from the set [0. 5, 1.5, 3], corresponding to small, medium or large) The following variables are defined to deal with objects (all distances are measured from the edge of objects, rather than from their centre): • BASESCORE (real) • BONUS (real) • PENALTY (integer) • BONUSMULTIPLIER (real) (to calculate maximum perimeter) • PENALTYMULTIPLIER (real) (to calculate minimum perimeter) There are n objects and m ants initially distributed at random on a twodimensional "board" of fixed size. The number of ants and numbers of objects of each size may be specified in advance. Ants may move over other ants or over objects; there are no spatial restrictions other than those imposed by collisions between unladen ants and objects.
The pseudo-code expression of the algorithm is as follows:

Implementation
The pseudo-code above was implemented using the Breve 1 multi-agent system environment [7]. This package facilitates the simulation of decentralised systems, in a similar fashion to the well-known Swarm 2 system. The benefits of using such a framework are derived from the fact that it automatically handles issues such as collision detection, construction of object neighbourhoods, and the user interface. This allows the programmer to concentrate effort on the important aspects of the model's implementation, rather than on "house keeping" issues. All simulations were run on an Apple Macintosh under MacOS, although the code runs equally well on both Linux and Windows machines. An additional Java application was written to deal with post-processing of object coordinates. All code is available, on request, from the corresponding author (MA).

Results
We found that the simulation consistently generated final cluster patterns that are strongly reminiscent of those observed in real colonies by Franks et al. [5]. A typical final pattern is depicted in Figure 4; note the densely-packed core of smaller items at the centre of the cluster, followed by successive "rings" of increasingly larger and more widely-spaced items, as observed in Figure 1. In Experiment 2, we carried out an additional 50 runs, replacing the following parameter values to reflect the object distribution of Colony A in [5]: • Number of small brood items = 85 For each experiment, we calculated the mean distance from centroid for each object type. Direct comparisons with the biological study were not possible, due to the different spatial units used; however, dividing the results for distance from centroid by the value of the smallest object gives the ratios described in Table  1. Values for Franks 1 refer to Colony A before migration, and Franks 2 refers to Colony A after migration.

Discussion
Comparing the values, we see that the model achieves values for distance from centroid that are very close to those observed in the biological study. The difference between the simulation and reality is about the same as the difference due to migration in the two biological studies.
We also carried out tile size analysis of our results, in order to assess the distribution of objects of different types. However, the only Voronoi tesselation software available to us (VpPlants 3 ) did not allow the imposition of a bounding box, so we were unable to obtain meaningful comparisons with the tiling results described by Franks et al. [5].
This study has demonstrated a minimalistic model that can consistently sort objects into annular clusters. There are real world examples where this might be useful, particularly when one has objects that may change over time and require different amounts of separation. Some of the aspects that set this study apart from others such as Wilson et al. [9] also make it harder to create a physical implementation. In particular the model relies on the fact that agents are able to pull poorly placed objects out of the centre of a cluster without disturbing other objects or risking collision, something that is hard to do with simple robotics. These difficulties will become less significant as the state of the art advances.
The model proposed in this study works well as an explanation of the sorting behaviour observed in Leptothorax ants. Unfortunately it is not possible to compare the results with those given by Wilson et al. [9] directly as a different performance metric was used. However the model described in this paper does not make use of a complex "leaky integrator" system, or have the unrealistic constraints on movement that were introduced by the robotic implementation used by Wilson et al. [9].
Towards the end of this study an additional paper by Sendova-Franks et al. [8] became available. It presented a further biological study of brood sorting by Leptothorax Ants. As with the previous study [5], ants were induced into nest migration and observed. This time the focus was on the behaviour of the ants and the authors were able to observe two distinct phases. During the "clustering" phase, ants carried brood from the old hive to the new one. During this phase carrying distance was determined by brood size. Ants formed a small cluster inside the new hive with small objects furthest from the entrance and large objects closest. In the second phase the ants spread out the pre-clustered brood items to achieve an annular cluster. This work naturally has some impact on the model proposed in this study but it does not necessarily invalidate it. The behaviour described in this paper could well explain what occurs during the clustering phase seen in [8], the only difference being that the ants start the second phase with the brood already partially clustered. At present it is not clear whether the first phase is actually essential to the process.

Conclusions
In this paper we have described a minimal model of annular brood sorting in Leptothorax ants. Simulation results are consistent with the structures observed in actual ant colonies. Our model uses only stochastic ant behaviour and a brood item "repulsion" and "attraction" mechanism, which supports the previous argument that this emergent behaviour may well be linked to pheromone deposition. Further work in this area will include a more detailed analysis of the results in terms of object distribution, as well as continued experimental investigations into the underlying biological mechanisms.