R Reuillon
email: romain.reuillon@isima.fr

D R C Hill

C Gouinaud

Z El Bitar

V Breton
email: breton@clermont.in2p3.fr

Monte Carlo Simulation With The GATE Software Using Grid Computing

Keywords: G.3. [Probability and Statistics]: Statistical Computing. J.3 [Life and Medical Sciences] Health. I.6.3 [Simulation and Modeling] Application Algorithms, Performance, Design, Reliability, Experimentation Monte Carlo, Grid computing, GATE simulation

Monte Carlo simulations needing many replicates to obtain good statistical results can be easily executed in parallel using the "Multiple Replications In Parallel" approach. However, several precautions have to be taken in the generation of the parallel streams of pseudo-random numbers. In this paper, we present the distribution of Monte Carlo simulations performed with the GATE software using local clusters and grid computing. We obtained very convincing results with this large medical application, thanks to the EGEE Grid (Enabling Grid for E-sciencE), achieving in one week computations that could have taken more than 3 years of processing on a single computer. This work has been achieved thanks to a generic object-oriented toolbox called DistMe which we designed to automate this kind of parallelization for Monte Carlo simulations. This toolbox, written in Java is freely available on SourceForge and helped to ensure a rigorous distribution of pseudo-random number streams. It is based on the use of a documented XML format for random numbers generators statuses.

INTRODUCTION

Monte Carlo simulations (MCS) are widely used in emission tomography; for protocol optimization, design of processing or data analysis methods, tomographic reconstruction, or tomograph design optimization. GATE [START_REF] Jan | GATE: a simulation toolkit for PET and SPECT[END_REF] is a Monte Carlo simulation tool based on the Geant4 package and dedicated to Single Photon Emission Computed Tomography and Positron Emission Tomography simulations. It was designed to be flexible and precise, thus GATE simulations are computer intensive and cannot be used in a clinical context. This work presents a distributing method and a tool for the parallelization of MCS. This method is then applied to a practical application in image reconstruction using GATE and execution times are given for clusters and the EGEE European grid environment.

MATERIAL AND METHODS

MCS are commonly considered to be naturally parallel [START_REF] Mascagni | Parameterizing parallel multiplicative lagged-Fibonacci generators[END_REF]. It is widely assumed that with N processors executing N replicates of a Monte Carlo calculation, the pooled result will achieve a variance N times smaller than a single instance of calculation in the same time [START_REF] Mascagni | SPRNG: a scalable library for pseudorandom number generation[END_REF]. In the next sections we discuss why we changed the default Pseudo-Random Number Generator (PRNG) of the GATE software and will also present how we separate experiments to avoid correlations that could slow down the convergence.

A. Using GATE with a better Pseudo-Random Number Generator (PRNG)

GATE simulations were initially based on the "James Random" algorithm [START_REF] James | A review of pseudorandom number generators[END_REF], [START_REF] Marsaglia | Toward a Universal Random Number Generator[END_REF] as implemented in the "Class Library for High Energy Physics" (CLHEP) [START_REF] Lönnblad | CLHEPa project for designing a C++ class library for high energy physics[END_REF]. This generator is 21 years old and has been shown to have poor statistical properties. We checked that it succeeded in only 36 tests out of 122 using the recent and already well-known statistical test battery "TestU01" of L'Ecuyer [START_REF] Ecuyer | TESTU01: a software library in ANSI C for empirical testing of random number generators[END_REF]. We therefore modified GATE to use the Mersenne Twister 19937 [START_REF] Matsumoto | Mersenne Twister: A 623-dimensionally equidistributed uniform pseudorandom number generator[END_REF] as implemented in CLHEP. This generator is recent, has a huge period of 2 19937 and is equidistributed in 623 dimensions. It passes almost all the tests of the test battery TestU01 and it is fast.

Parallelization of PRNG

For quantitative Monte Carlo simulations the MRIP or "Multiple Replication In Parallel" parallelization approach ( [START_REF] Hill | Object-oriented pattern for distributed simulation of large scale ecosystems[END_REF], [START_REF] Pawlikowski | Towards credible and fast quantitative stochastic simulation[END_REF]), allows a maximum speed up if many replications of the same experiment have to be made in order to obtain a good approximation of the result. However, when parallelizing the underling pseudo-random number generator (RNG), correlations within and between the random numbers streams generated in each processor have to be avoided [START_REF] Mascagni | Parameterizing parallel multiplicative lagged-Fibonacci generators[END_REF]. Different parallel generation techniques of pseudo random numbers can be found in the following documentation: [START_REF] Traore | The use of random number generation for stochastic distributed simulation: application to ecological modeling[END_REF]. In the "central server approach", a central RNG generator provides numbers for all simulation jobs. This approach is the natural one but doesn't fulfill the requirements for a good parallel RNG [START_REF] Coddington | Random number generator for parallel computers[END_REF] and creates a bottleneck that slows down the distributed simulation. The "sequence splitting" or "blocking" consists in splitting the RNG cycle into non-overlapping contiguous sections [START_REF] Coddington | JAPARA -A java parallel random number library for high-performance computing[END_REF]. This technique must be used with caution because long range correlation in the parallelized generator might become short range inter-sequences correlations. Instead of unrolling the generator, one might consider to randomly generate states of the pseudorandom number generator (a status is archiving a precise state). The average minimal distance between n statuses should be in this case 1/n 2 times the distributed generator period [START_REF] Wu | Parallel use of multiplicative congruential random number generators[END_REF]. This is possible using a cryptographic generator or a hash function to generate the statuses. The distribution of the Mersenne Twister 19937 algorithm is achieved that way in the library SPRNG [START_REF] Mascagni | SPRNG: a scalable library for pseudorandom number generation[END_REF]. The "leap frog" technique distributes the sequences to the processor like a deck of cards to card players. Each process of the distributed simulation uses 1 number out of n in the original sequence. This last technique requires a generator that allows cycle division [START_REF] Jan | GATE: a simulation toolkit for PET and SPECT[END_REF]. With this technique long range correlations in the original sequence might also become short range intersequences correlations if the interval of sampling in the original sequence is not chosen carefully [START_REF] Wu | Parallel use of multiplicative congruential random number generators[END_REF]. The "independent sequences" technique produces different cycles of numbers depending on the initial seed. This technique is available for a few generators like some lagged Fibonacci pseudo-random number generators. [START_REF] Coddington | JAPARA -A java parallel random number library for high-performance computing[END_REF]. This last technique is close to the "parameterization" technique that might be used with some RNG like the Mersenne Twister [START_REF] Matsumoto | Dynamic creation of pseudorandom number generators[END_REF] or linear congruential generator with Mersenne or Sophie-Germain prime moduli [START_REF] Mascagni | Parallel linear congruential generators with Sophie-Germain moduli[END_REF]. It generates algorithm parameters leading to the generation of highly independent random number streams. Within the current state of the art, we are not able to provide a theoretical proof of independence between pseudo-random number streams. However, various approaches can be tested empirically, implying heavy computation that can be achieved once for many applications under a precise experimental framework. The Mersenne Twister 19937 has a very long period of 2 19937 drawings. It is an already parameterized version of the generic Mersenne Twister algorithm and it has no efficient cycle division technique available. The Mersenne Twister 19937 generator is well suited to a parallelization using the "sequence splitting" technique. To achieve the non-overlapping condition of the "sequence splitting" parallelization technique, we first estimated the number of random numbers drawn using a simulation job designed to have an average execution time of 12 hours (on an average working node of the European grid). This estimation led to 12 billion drawings per job. Then we generated over 6000 statuses for the Mersenne Twister spaced by 15 billion numbers each. We used a similar approach in [START_REF] Maigne | Parallelization of Monte Carlo simulations and submission to a grid environment[END_REF]. The generated status were archived and converted into a documented XML format, in order to be reused with different implementations of the Mersenne Twister 19937 algorithm.

Creation of a generic parallelizing tool for MCS

We designed and implemented an open source software tool in Java called "DistMe", which is dedicated to parallelize stochastic simulations. This tool contains a status database and is able to create jobs for various distributed environments independently from the random number generation library. It is based on the intensive use of a documented XML generic format for the pseudo-random number generator statuses [START_REF] Reuillon | Rigorous Distribution of Stochastic Simulations Using the DistMe Toolkit[END_REF]. DistMe is fully usable and its sources can be found on the Internet (http://sourceforge.net/projects/distme). With this tool, we could generate GATE jobs for any distributed execution system: basic scripts (using ssh for instance), bags of work for the European Grid using JDL descriptors (the European grid Job Description Language) and more specially "OpenPBS" (using Portable Batch System scripts). A tutorial is available (www.isima.fr/~reuillon) as well as 9000 statuses for the Mersenne Twister 19937 algorithm, spaced of 15 billion drawings each.

Hardware

We could access 650 worker nodes of the EGEE European computing grid (Enabling Grid for E-sciencE, www.eu-egee.org) mainly in France, United Kingdom, Netherlands and Poland. We also had at our disposal two clusters hosted by local research laboratories (the LIMOS/ISIMA cluster composed of 14 biprocessors and the LAMI/IFMA cluster composed of 28 biprocessors managed by an "OpenPBS" system. Each processor of the cluster is an Intel Xeon 3 GHz with hyper-threading.

Merging of simulation results

Each simulation generates two binary output files requiring about 10 megabytes of storage space. Simulation output files produced on the grid were automatically registered and copied on a Storage Element (SE). When all the simulations were completed, a script using grid commands retrieved these files from the SE into a local machine. When the computing was performed on a local cluster, the retrieving of the simulation output files was achieved using a regular and local FTP commands (File Transfer Protocol). The merging of all files was performed using a simple C code and required less than 5 minutes for less than 30 Gigabytes (on a local desktop computer -Xeon 3 GHz with simple SATA disk).

RESULTS

PRNG Parallelization

The computation of the 6000 PRNG status is not possible in parallel and took around 80 days on a single node of the ISIMA cluster running at 2.4 GHz. Once the list of status is generated, it is important to test the resulting random number series. Indeed, a good parallel PRNG must behave like several good sequential PRNG. Each sequence was then tested using the statistical tests battery for sequential PRNG TestU01. As shown in Table I, only 2% of sequences failed in more than 5 tests and no sequence failed in more than 10 tests out of the 122 tests of the battery. This calculation has been made on the ISIMA cluster and took 35 days at full cluster load. If the work has been made on a single machine it would have taken around 3 years.

Simulation execution

The two local clusters achieved 600 jobs, which all succeeded and 2300 jobs were executed on the EGEE European grid, (with 1811 usable results). Fig. 1 shows how the jobs were executed on the different calculation units. The IFMA cluster hosted 400 jobs, 200 jobs ran on the ISIMA cluster, 200 on the Polish worker nodes, 499 on the Dutch ones, 922 in England and 190 in France. A variance study on the final results showed that after the execution of 2000 jobs a convergence was reached. The curve on Fig. 2 shows an asymptotic behavior around 2000 jobs. Executing more jobs was then un-necessary.

Distribution of the computing time

The time required for the total execution of the simulation on a single sequential computation unit (Intel Xeon 3 GHz) is 906 days/CPU (Central Processing Unit).

The execution time of the distributed simulation on clusters is inversely proportional to the number of processors from which they are composed since the migration time of the jobs is neglected compared to the total execution time. Hence, the gain factor was 84 since the number of local bi-processors is 42 (14 biprocessors on ISIMA cluster; 28 bi-processors on IFMA cluster) resulting in 84 execution units running in parallel. For the grid the problem is a bit trickier. The execution power of the grid is virtually unlimited, supposing that the number of processors available on the grid is always greater than the number of submitted jobs. Furthermore, we may consider that there is no latency time in the grid architecture: no job migration time, no data migration time, and that execution units are homogenous and fast. Hence, the ideal proportional gain factor in time execution is proportional to the number of jobs in which the sequential simulation is distributed. In our case the number of jobs has been arbitrarily set to 1813 jobs. Each job should run during 12 hours on a fast local computer for completion. This mean the total execution time of the simulation is 22356 hours / CPU (Central Processing Unit). To compute the gains in the next paragraph we compute virtual execution times of jobs under certain conditions:

taking into account the average execution times of all jobs or only the longest execution time of a job taking into consideration that we have access to a limited part of the grid or consider the grid as able to execute all our jobs in parallel

In each case, we compute virtual execution times for our simulation jobs. After that we consider that the virtual jobs are executed perfectly in parallel. By consequence the resulting gain is the total execution time of the simulation divided by the duration of one virtual job.

On a real grid, the gain factor in term of execution time is penalized by the latency time affected by the performances of the targeted worker nodes and the migration time in which a job passes through the following states: submitted, waiting, ready and scheduled. Taking into account the latency and the average power of the execution units, we have computed from the execution log files the average execution time among all our simulation jobs on the grid. The average execution time for our simulation is 24.675 hours. This lead to what we have called the average theoretical gain. In our case we obtain a gain of 906. Unfortunately, the number of worker nodes available in our real execution environment was 650. It is greatly inferior to the number of jobs we had to execute. This means, that we could only execute 650 jobs concurrently. This impact negatively on the virtual execution time of our global simulation jobs, increasing it to 70,722 hours by job in average, thus the average practical gain factor is 316.

The end-user might be interested in the global execution time of the simulation corresponding to the time between the submission of the jobs and the return of the results from the last job. Supposing a concurrent submission the simulation ends when the last part of the results is returned. From the log files, the longest job with the longest execution time over all jobs is 36 hours. This leads to a minimal theoretical gain of 621.

CONCLUSION AND DISCUSSION

By distributing the calculation on many execution units our nuclear medicine simulation was achieved in a few days. It would have taken more than three years on a single powerful computer without distributing the simulation using the MRIP approach. We have not repeated this simulation to study the grid and cluster overhead, since we may obtain different execution times with different grid/cluster loads. The simulation results were directly used by scientists working in nuclear medicine [START_REF] Lazaro | Feasibility and value of fully 3D Monte-Carlo reconstruction in single photon emission computed tomography[END_REF], [START_REF] Lazaro | Fully 3D Monte Carlo reconstruction in SPECT: a feasibility study[END_REF].

Improvements can be made in the following directions: the different gain factors might be improved using more worker nodes of the grid and optimization techniques for stochastic simulations distribution like the "N out of M" strategy presented in [START_REF] Li | Improving Performance via Computational Replication on a Large-Scale Computational Grid[END_REF]. Furthermore, each random number sequence has been tested individually with the best test battery presently available, but tests have to be done to check that the correlation between the sequences is acceptable using the parallel PRNG tests described in [START_REF] Jan | GATE: a simulation toolkit for PET and SPECT[END_REF] and implemented in SPRNG. It represents a huge amount of calculation and it will be achieved using the internet computing platform BOINC [START_REF] Anderson | BOINC: A System for Public-Resource Computing and Storage[END_REF]. The use of the DistMe toolbox requires the downloading of the statuses from the internet and a manual operation to insert them in a database on the local computer running DistMe. To simplify this task for the end user, the statuses and the tests results will be published via a central web service and DistMe will gain a transparent access to this web service. Last but not least, it might be interesting to optimize the status generation phase, by combining the sequence splitting technique with highly independent random numbers sequences obtained using a "parameterization" technique and then generating the statuses for each sequence in parallel. A pseudorandom number generation library, "DistRNG", is being implemented and already allow the use of cutting edge parallelization techniques.

Table 1 .Figure 1 .Figure 2 .

 112 Figure 1. Repartition of the jobs

ACKNOWLEDGMENTS

The authors would like to thank the Auvergne Regional council (France) for its research sponsoring, the EGEE consortium and also ISIMA/LIMOS & IFMA/LAMI laboratories which provided an easy access to their local clusters.