Debugging data transfers in CMS

The CMS experiment at CERN is preparing for LHC data taking in several computing preparation activities. In early 2007 a traffic load generator infrastructure for distributed data transfer tests was designed and deployed to equip the WLCG tiers which support the CMS virtual organization with a means for debugging, load-testing and commissioning data transfer routes among CMS computing centres. The LoadTest is based upon PhEDEx as a reliable, scalable data set replication system. The Debugging Data Transfers (DDT) task force was created to coordinate the debugging of the data transfer links. The task force aimed to commission most crucial transfer routes among CMS tiers by designing and enforcing a clear procedure to debug problematic links. Such procedure aimed to move a link from a debugging phase in a separate and independent environment to a production environment when a set of agreed conditions are achieved for that link. The goal was to deliver one by one working transfer routes to the CMS data operations team. The preparation, activities and experience of the DDT task force within the CMS experiment are discussed. Common technical problems and challenges encountered during the lifetime of the taskforce in debugging data transfer links in CMS are explained and summarized.

Introduction
The CMS experiment is a large particle physics experiment located at CERN, Geneva, Switzerland that is presently being commissioned for beginning of data taking. CMS depends on a world wide distributed data grid of about 50 computing and storage clusters to archive and analyze its data. Individual clusters vary both in size (10TB to a few PB) as well as expertise of their operations teams. The several hundred end-to-end data transfer links between these sites needed to be commissioned. In July 2007 CMS created a "Debugging Data Transfers" (DDT) task force to coordinate the debugging of data transfer links in the preparation period and during the CSA07 data transfer test [1]. The CSA07 service challenge was a data challenge in 2007 designed to test the transfer system at 50% of the design goal for 2008. The goal of the DDT task force was to deliver fully debugged and operational end-to-end links to the CMS data operations team. We aimed to commission the most crucial transfer routes among CMS tiers by designing and enforcing a clear procedure to debug problematic links. The procedure aimed to move a link from a debugging phase in a separate and independent environment to a production environment when a set of agreed conditions were achieved for that link.
This note details the activity of this task force before, during and since CSA07. Section 2 describes the task force charge and scope. Section 3 describes some of the details of the CMS Computing Model and the CSA07 data challenge relevant to this task force. Section 4 describes briefly the system components used to transfer data across the wide area network. Section 5 details the metric used to commission links and progress in link commissioning. Section 6 concludes with a description of the current status of link commissioning in CMS.

Task Force Charge within CMS
The DDT task force was focused on the status of data transfer links, defined as unidirectional end-toend data transfer between site A and site B. The responsibilities of the task force were set out to be: • To define details on how the metrics are measured to put links in/out of production status.
• To define a procedure, including a set of steps or stages to pass that gets a link from a decommissioned state to production. • Definition of the procedure to commission a link, including documentation of the kinds of tests, and tools to use. This includes helping sites to resolve their problems by pointing them to storage element (SE) support channels for the SE they have chosen to deploy, for example. The task force is the first point of contact for the site administrators. The task force thus facilitates information exchange. • Documentation and creation of a list of known problems encountered, and instructions for solving them. • Creating a table that keeps track of the matrix of status of all links.
• Reporting weekly on the status of this matrix.

The CMS Computing Model
The CMS computing model [2] has three tiers of computing facilities. These sites are interconnected by high-speed networks of 1-10Gbps. Data flows between and within each of these tiers. These include the Tier 0 at CERN (T0), used for data export from CMS and archival to tape, and 8 Tier 1 (T1) centres, including one at CERN, used for the tape backup and large-scale reprocessing of CMS data, and distribution of data products to the Tier 2 centres. The T1 centres are typically at national laboratories with large computing facilities and archival storage systems. There are ~50 Tier 2 (T2) facilities, where data analysis and Monte Carlo production are primarily carried out. These centres are typically at universities and do not have tape backup systems, only disk storage. The CMS computing model envisions commissioning all links between the CERN to T1 sites, and T1 sites to CERN (14 links), all other T1-T1 cross-links (42 links), all T1 to T2 downlinks (~400 links). and all T2 to "regional" T1 uplinks (~50 links). Therefore, the total number of links to be commissioned in the computing model is ~500. This number will increase with the addition of new T2 sites. At the beginning of the task force only about 30 transfer links were sufficiently stable to be considered commissioned. T2 to non-regional T1 uplinks were not a priority in 2007 but were commissioned if the sites wished. Each T2 is associated to a T1 (called the "regional" or "associated" T1), although in some cases this T1 is not geographically near the T2. The T2-to-T2 cross-links are not part of the computing model, but in fact are used especially within the same country as in the United States, Germany and Belgium. These ~2500 links are not included in the scope of the computing model but were also considered by the DDT task force if the sites wanted to commission them. Likewise, links that begin or end at a Tier 3 (T3) site are not in the computing model. A T3 site is typically a small or mediumsized computing facility associated with a T1 or a T2, usually at a university or research institute. These links were commissioned on request of the T3 site. There were only six T3 sites active within CMS that attempted data transfers during the CSA07 period.

CSA07 Workflows
The CSA07 data transfer challenge preparation began with the production of Monte Carlo datasets, primarily at T2 sites. These datasets were transferred to CERN via the regional T1 sites at which point the CSA07 workflows would begin at CERN. For this activity to take place, commissioned data transfer uplinks from T2 sites to at least one T1, and from those T1 sites to T1_CERN needed to exist.
The subsequent CSA07 workflows included the transfer of data from the T0_CERN to the 7 T1 sites, data reprocessing at the T1's, and data transfer of the resulting skimmed datasets to the T2 sites, either directly or via a T1-T1 cross-link. The CSA07 project is described in [1], including the various metrics measured to determine success.
After several iterations of discussion it was decided that the priority of debugging effort was to be first the CERN-T1 and T1-CERN links, then all other T1-T1 links, next the T1 to and from associated T2 links, establishing at least one link per T2 in each direction so that the T2 site is useful to data operations for data analysis and Monte Carlo production activities, and lastly the links for a T1 to and from non-regional T2 sites.

PhEDEx, the Data Transfer Middleware
PhEDEx [3,4] is the data transfer middleware of the CMS experiment. Within PhEDEx there are several "instances", which generally means separate databases, accounting, etc. The "Production" instance is for commissioned links only, and carries out the CSA07 workflows and Monte Carlo production transfers. A "Debug" instance was used from early August 2007 to handle test transfers and in September became used exclusively for the test transfers.
The PhEDEx LoadTest [5,6] is the main way that data transfer links are tested within CMS. All DDT traffic is in the context of the PhEDEx LoadTest. The procedure is to inject files at a certain rate into the database and queue them for transfer over the various links. Injection rates are now tuneable on the Web. The injected logical files map to 256 real files at each site, so that files are transferred multiple times without the need to constantly create new files. As seen in Figure 1, LoadTest transfers began in the Debug instance of PhEDEx in early August and ramped up before the start of CSA07 to about 1GB/s. This rate was maintained more or less constantly over the course of CSA07 despite the addition of new links throughout September and October. This indicates that there are real limitations in the capacities of sites to transfer data in and out of their storage elements. These limitations are well below the physical network capacities of the links between sites.

Storage Elements
The CMS T1 and T2 sites use various storage elements (SE) systems, which are briefly described.
The dCache storage system [7] was developed at DESY and FNAL. The system gets support from both US OSG storage support group as well as dCache User Forum in Europe. It is used by FNAL, PIC, FZK and IN2P3 as the T1 storage backend and by a majority of the active T2 centres that have deployed distributed storage solutions. As the system is developed in Java, it can be deployed on different platforms and different operating systems and a number of these have been used also in the CMS associated sites.
The CASTOR [8] storage system was developed at CERN and is used by the T1 centres at CERN, CNAF, RAL and ASGC. There are no T2 sites that use this storage system. Developed as part of the gLite middleware, DPM [9] is mostly used by smaller T2 centres. About 10 T2 sites use DPM as their storage solution.
The T1 sites also have tape back ends to their storage systems, which were not exercised in this project.

Transfer Protocols
Transfers between some sites are made with 3-rd party SRM transfers [9] using GridFTP. However, where load is an issue, transfers are scheduled by FTS [10]. FTS is part of the gLite middleware and therefore used mostly by the EGEE sites. Its main features include submission of data transfer jobs, which are scheduled by an FTS server based on the settings of the channel utilized for that specific source/destination combination. It allows sites to set limitations on the number of files in transfer, number of streams used etc. Transfers to or from a T1 site use the FTS server at the T1 site. Regional T2 sites have dedicated FTS channels at the servers at their associated T1, while "non-regional" transfers use T1-STAR channels.

Networking
No major issues with networking handicapped the debugging of transfers. However, it is apparent that majority of links are not performing anywhere close to the speeds per stream that they should be able to achieve, and no coordinated effort has been done to identify fully the reasons. Only a handful of sites have performed testing to understand the network path between them and to try to tune their storage accordingly. The majority of storage nodes are running un-tuned default kernel configurations that do not favour high-speed long distance transfers, but are sometimes tuned for the requirements of the storage elements at the sites.
Most T1 sites are interconnected through 10Gbps networks, although they serve T2 sites that are often connected through national network infrastructures with a more limited capacity (1-2 Gbps). This imposes limitations on some regions. Some hardware issues that caused network outages in Brazil occurred during CSA07. Network bandwidth was rarely a limiting problem otherwise.

Link Commissioning Metric
The first activity of the DDT task force was to define and implement a metric by which links can become commissioned and subsequently handed over to data operations. There are several stages through which a link passes from "NOT-TESTED" to "COMMISSIONED": • NOT-TESTED: links never actually tested, i.e. links showing no successful transfer attempts within PhEDEx. • PENDING-COMMISSIONING: links that have transferred successfully at least one file in PhEDEx, but have not yet passed the requirements below for link commissioning. • COMMISSIONED: links that are demonstrated to work, and can be delivered to data operations. Note that this commissioning does not imply that the link or the site has met the requirements of the computing model or the current service challenge, but simply that the link has passed some minimum requirements to be considered usable for data operations. To be COMMISSIONED during 2007, a link was required to: • Transfer 300GB/day for 6 out of 7 consecutive days, and transfer a total of 2.3TB during that same 7-day period. • For links involving an endpoint at a T2, this requirement was relaxed to 4 out of 5 days, and a total transfer volume of 1.7TB. This is to match the service requirement of business hours only support committed to by the T2 sites. • PROBLEM-RATE, for links that were working but whose rate has dropped off. To remain COMMISSIONED, a link had to transfer at least 300GB/day for a single day at least once every 7 days. Otherwise, the link had to be re-commissioned by following the procedure above.
These requirements were developed with the idea of having a higher threshold to commission than to decommission the link. These thresholds can be increased in time as networks and sites develop, since the rates implied by 300GB/day are of the order of 3-4MB/s per link, far below the commitments envisioned in the computing model which envision T2 sites being able to download a total of up to 5TB/day from T1 sites, or over 60MB/s sustained downloads. However, the computing model also envisions that transfers will occur in bursts, not as a continuous rate over several days. The metric used during CSA07 deviated from this model to prove the stability of data transfer links.
A member of the DDT team presented a list of changes to the COMMISSIONED or PROBLEM-RATE links to the data operations team at daily and weekly meetings. Any exceptions to the commissioning or decommissioning metric were decided in this meeting and were only rarely approved. Full enforcement of these procedures began in early September 2007 in a phased manner.
In 2008-2009 these commissioning metrics were revised to more closely match the Computing Model requirements for higher rates and transfers in bursts. Links were required to transfer at least 1.65 TB (>20MB/s) in a 24 hour period, or 422 GB (>5MB/s) for T2-T1 uplinks only.
To remain commissioned, all links were periodically exercised in a random order and were required to meet half of the above metric goals within 3 calendar days. During link exercising about 1% of transfer links showed problems that resulted in decommissioning. Almost all such problem links were re-commissioned in a short period of time.

Monitoring of Link Status
A tool was developed by B. Bockelman and Sander Sõnajalg, a CERN summer student, to extract transfer volume data from PhEDEx and apply the DDT commissioning criteria. This tool takes data transfers from both the Production and Debug instances of PhEDEx into consideration. Figure 2 shows an example of this DDT Matrix. Green links are those that are COMMISSIONED, red are PROBLEM-RATE, light blue are PENDING and white links are NOT-TESTED.

Link Commissioning Progress
During the CSA07 T1-T1 link debugging was identified as an urgent priority, as only a few T2 sites were able to commission links from multiple T1s and keep them commissioned over a long time period. To allow T2 sites to participate in the physics transfers phase of CSA07, the data would have to be copied from one T1 to another T1 before the T2 could download the data from the associated T1. As can be seen on Figure 3, around the beginning of October a major increase in commissioned links shows the time when DDT effort was prioritized for T1-T1 matrix commissioning. In a short time period over half of the links were commissioned. No major technical obstacles were overcome during this period. This was in part legacy from previous efforts to bring up the transfer infrastructure and in part a matter of focus and attention of the T1 administrators to make this activity a priority, mainly in setting up LoadTest samples, creating and approving subscriptions, etc. Although there were ongoing technical challenges before and after the links became commissioned, there was no major technical obstacle to commissioning these links in October.
The largest overall increase in commissioned links came in the period leading to the start of CSA07. At this time, we began to enforce the policy to use only commissioned links for production data transfers, which clearly encouraged sites to put in the effort to commission their links and participate in more than just the test exercises. Most of the progress in commissioning the T2 data transfer links involved configuration of FTS channels or configuration of the PhEDEx transfer software at the T2.

Current Status and Conclusions
During the task force, an improvement in the number and quality of data transfer links was achieved, through the hard efforts of site administrators, PhEDEx developers, data operations and networking experts, etc. The initial mandate of the DDT task force concluded with the end of the CSA07 service challenge in November 2007. However, the effort was considered useful and continued in 2008 in a modified form. The DDT task force continues to aid sites in their data transfer link commissioning efforts. The task force undertook dedicated campaigns in 2008 that included helping sites complete the commissioning of all of their downlinks from the T1 sites, for example, which are documented in the CMS twiki [11]. Most sites were also encouraged to commission a second T2-to-T1 uplink, so that Monte Carlo production data could be uploaded to several T1 sites.
Firstly, the metric was modified to more closely match the CMS computing model requirements. To commission links from 2008 onwards, a data link must transfer at a rate of at least 20MB/s over 24 hours. Recognizing that uplinks from T2 to T1 sites have a lower requirement in the computing model, they are only required to transfer 5MB/s. Secondly, it was recognized by the CMS computing management that continual exercising of transfer links placed an unnecessarily large burden on site administrators and storage systems. Within the computing model data links are foreseen to be transferring data in bursts with periods of inactivity. To more closely match the model, the requirements for links to stay commissioned was changed. A testing program was organized in 2008 in which the DDT task force attempted to exercise each link in rotation for 12 hours, trying to meet the metric goal of 20MB/s or 5MB/s for a T2-T1 uplink. Links were only decommissioned if they failed this transfer exercise for three days, or developed an obvious problem and the data operations team requested the disabling of the transfer link. The status of commissioned links as of March 2009 is as follows. Thirty-five T2 sites have all of their downlinks from T1 sites COMMISSIONED, and 2 more have seven out of eight links COMMISSIONED. Forty-three T2 sites have at least two COMMISSIONED uplinks to T1 sites. All T1-T1 links are COMMISSIONED. Transfer quality continues to improve, as seen in Fig. 5.
The experience of the DDT task force showed that dedicated computing campaigns result in rapid progress in commissioning data transfer links and improving permanently the quality of data transfers, for example. The procedures and metrics developed by the DDT task force are now part of the commissioned links overview by the CMS Site Commissioning project [12,13].
In conclusion, a focused effort to debug transfer links within CMS proved to be useful in helping to maintain a working system for data transfers, documenting common problems and solutions, and alerting site administrators to problems. This effort is continuing with requirements and testing exercises that more closely match the CMS computing model and expected data transfer patterns at the start of data taking later this year.


Computational advances in gravitational microlensing: a comparison of CPU, GPU, and parallel, large data codes

To assess how future progress in gravitational microlensing computation at high optical depth will rely on both hardware and software solutions, we compare a direct inverse ray-shooting code implemented on a graphics processing unit (GPU) with both a widely-used hierarchical tree code on a single-core CPU, and a recent implementation of a parallel tree code suitable for a CPU-based cluster supercomputer. We examine the accuracy of the tree codes through comparison with a direct code over a much wider range of parameter space than has been feasible before. We demonstrate that all three codes present comparable accuracy, and choice of approach depends on considerations relating to the scale and nature of the microlensing problem under investigation. On current hardware, there is little difference in the processing speed of the single-core CPU tree code and the GPU direct code, however the recent plateau in single-core CPU speeds means the existing tree code is no longer able to take advantage of Moore's law-like increases in processing speed. Instead, we anticipate a rapid increase in GPU capabilities in the next few years, which is advantageous to the direct code. We suggest that progress in other areas of astrophysical computation may benefit from a transition to GPUs through the use of"brute force"algorithms, rather than attempting to port the current best solution directly to a GPU language -- for certain classes of problems, the simple implementation on GPUs may already be no worse than an optimised single-core CPU version.

Introduction
Gravitational microlensing is the deflection of light by matter in the regime where multiple-imaging and high magnification events occur, but the angular separation of images is too small to resolve. The observable effect is a change in brightness of a source over a period of time, with the details of the light curve depending on the complexity of the lens distribution and the surface brightness profile of the source. Gravitational microlensing has been used to study high magnification events on local scales, viz. the searches for dark matter compact objects in the Galactic bulge and halo (see for example references in and citations to Alcock et al. 1993;Aubourg et al. 1993;Udalski et al. 1993;Calchi Novati et al. 2002;Wyrzykowski et al. 2009) and extrasolar planets (e.g. Gould & Loeb 1992;Bond et al. 2004;Sumi et al. 2010), and quasar microlensing on cosmological scales (e.g. Vanderriest et al. 1989;Irwin et al. 1989;Kochanek 2004;Bate et al. 2008;Dai et al. 2010). For recent reviews of microlensing on all scales, see Wambsganss (2006), Kochanek et al. (2007), Gould (2009), Mao (2008) and Gaudi (2010). 6 Research undertaken as part of the Commonwealth Cosmology Initiative (CCI: www.thecci.org), an international collaboration supported by the Australian Research Council Email addresses: nbate@astro.swin.edu.au (N.F. Bate), cfluke@astro.swin.edu.au (C.J. Fluke) Microlensing on local scales occurs in the low optical depth regime, where the source is affected by at most a few microlenses along the line of sight. Large numbers of stars must therefore be monitored in order to observe a single microlensing event. Conversely, the microlensing optical depth at the location of lensed quasar images is likely to be ∼ 1. This microlensing signal is therefore more complex, as it is the result of a large number of microlenses acting on the source essentially all the time (Wambsganss, 2006). In the remainder of this work, we focus on high optical depth microlensing.
A key element in many microlensing calculations is the determination of the magnification map; a typical example is shown in Figure 1. This map provides an estimate of the pointsource brightening over a finite region of a two-dimensional source plane. The magnification of an extended source is found by convolving a source profile (intensity and geometry) with the magnification map. This approach has been used to probe the structure of quasar emission regions, including accretion discs (e.g. Kochanek 2004;Eigenbrod et al. 2008;Floyd et al. 2009) and broad emission line regions (e.g. Lewis & Ibata 2004;Abajas et al. 2007).
The limitation in many microlensing calculations is the computational overhead in producing magnification maps with sufficient resolution and high (statistical) accuracy. Fitting observed light-curves obtained over long periods of time requires large maps with high pixel resolution to resolve small sources (such as quasar X-ray emission regions). Similarly, fitting Figure 1: A sample microlensing magnification map for a microlensing model with the following parameters: total convergence σ = 0.4, shear γ = 0.0, N pix = 1024 2 pixels in the source plane and N rays = 2000 light rays per source pixel (on average). These parameters correspond to N * = 558 lenses. The source plane has a side length of 20 Einstein radii. The greyscale covers the range from low magnification (black) to high magnification (white). multi-wavelength observations requires both large maps and high pixel resolution in order to simultaneously model emission from small optical emission regions and larger infrared emission regions. The situation is even worse when attempts are made to simultaneously model microlensing of accretion discs and the much larger broad emission line regions, a difference in scale of 2-3 orders of magnitude. Finally, many maps must be generated in order to ensure statistical independence of simulated measurements. It is only reasonable to sample any one magnification map as many times as you can place a source on it without overlapping.
In this paper we compare three techniques for computing magnification maps: a hierarchical tree code for single-core CPUs; a parallel, large data tree code for use on a CPU-based cluster supercomputer; and a direct, inverse ray-shooting code for use on graphics processing units (GPUs). We consider the accuracy of the codes, to determine whether the approximations inherent in the tree code are valid over a wider range of parameter space than has been feasible to test before, and compare run-times between the single-core tree code and the GPU code over a range of astrophysically-motivated parameters. While the "brute force" approach to ray shooting has not been considered a feasible option on CPU architectures, we show that the GPU implementation is highly competitive, and is indeed faster than the tree code in many circumstances. The expected increase in GPU speeds in the years ahead bodes well for advances in gravitational microlensing simulations. Our work hints at the potential benefit to other fields of astrophysical com-putation from the implementation of simple, "brute force" algorithms as a means of accelerating the adoption of GPUs in astronomy for time-consuming computations.
The remainder of this paper is set out as follows. In Section 2, we describe the process of generating magnification maps, and summarise key features of three gravitational microlensing codes which can be used for this task. In Section 3, we present a series of code comparisons (accuracy and speed). In Section 4, we discuss the relative merits of hierarchical and direct rayshooting approaches, and suggest the classes of problems to which each is best-suited. We present our conclusions in Section 5.

Generating magnification maps
In general, magnification maps are calculated using the (twodimensional) gravitational lens equation which relates the locations of the source, y, and an image, x, via the deflection angle, α(x). The magnification is: where A = ∂y/∂x is the Jacobian matrix of equation (1). Equation (1) gives a one-to-one relationship between an image location and the corresponding source location (x → y), but a one-to-many relationship for the alternative case of mapping a source location to its multiple images (y → x). While analytic solutions exist for a small number of simple lens models (e.g. Schneider et al. 1992), in general a numerical approach is used to calculate µ over a finite region of the source plane. This can be achieved most easily using inverse ray-shooting (Kayser et al. 1986;Schneider & Weiss 1986;Schneider & Weiss 1987). Here, light rays are propagated backwards from the observer, through the lens plane, where they are deflected using equation (1), and then mapped to a grid of pixels in the source plane. This approach avoids the need to determine the location of every image. By counting the number of light rays reaching each source plane pixel, N i j , compared to the average per pixel if there was no lensing, N rays , an estimate of the per-pixel magnification, µ i j , is obtained: For cosmological microlensing problems, the appropriate model for the deflection angle describes N * compact lens masses in a smooth background mass distribution, acted on by an external (macro-model) shear γ. This model, an extension of the one proposed by Chang & Refsdal (1979), is suitable for considering the impact of lensing due to individual stars where the external shear is due to the large-scale galactic mass distribution. The lens equation in this case is: where the convergence σ = σ c + σ * combines contributions from smooth matter σ c , and compact objects σ * . This is the model for the deflection of light by a screen of mass that we will use throughout this paper. The relationship between σ * and the number of lenses N * is: If the deflecting mass distribution were smoothed out, a circle in the image plane would map onto an ellipse with major-to-minor axis ratio |1 − σ + γ|/|1 − σ − γ| in the lens plane. The area of this 'minimum' ellipse is set by the size of the image plane. In practice, the graininess of the matter distribution can cause rays to be scattered outside the receiving area. Similarly, rays from a long way outside the shooting region can be scattered into the receiving square in the source plane. For that reason, the lenses in the codes considered here are distributed in a circle of area A with a diameter larger than the long axis of the minimum ellipse.
The number of calculations, Φ, required to obtain a mangification map with N pix pixels scales as: where N FLOP is the number of floating point operations for each deflection calculation, and N rays is chosen to achieve a desired level of accuracy. It was quickly realised that the direct rayshooting approach was not practical on single CPUs: a scenario, Φ fiducial ∼ 3×10 16 , with N pix = 2500 2 pixels, N rays = 500, N * = 10 6 lenses, 1 and N FLOP = 10 would take months to years to calculate (Wambsganss, 1999) on hardware that was available at the time. Wambsganss (1990;1999) introduced a method based on a hierarchical tree code (Barnes & Hut, 1986). Here, the contribution of each lens depends on its distance from the light ray -lenses that are further away are grouped together into pseudo-lenses with a higher mass. This approach reduces N * , thus increasing the speed of the overall calculation, but at the cost of a slight error in the magnification map due to the approximation. This approach has been used to constrain both the size and temperature profile of quasar accretion discs (e.g. Wambsganss, Paczynski, & Schneider 1990;Rauch & Blandford 1991;Morgan et al. 2006;Pooley et al. 2007;Keeton et al. 2006;Anguita et al. 2008;Eigenbrod et al. 2008;Bate et al. 2008;Floyd et al. 2009), dark matter fractions in the lensing galaxies (e.g. Pooley et al. 2009;Dai et al. 2010), and possibly even the orientation of the quasar with respect to the line of sight (Poindexter & Kochanek, 2010).
Despite the much-improved speed compared to direct rayshooting, the hierarchical method is complex to implement, and requires a certain amount of expert knowledge. These facts make slight modifications or additions to the code quite challenging. Parameters such as the required accuracy (which determines whether to use a particular lens cell, or its four subcells, based on the distance between the ray and the cell's centre of mass) must be determined experimentally by repeatedly 1 Relevant for the case σ ∼ 1 running the code until acceptable results are obtained. Additionally, the generation of the tree data structure imposes a memory overhead which exceeds typical single CPU memory sizes for large N * . A parallel hierarchical tree code, suitable for running billion-lens calculations, has now been developed by Garsden & Lewis (2010). This approach overcomes the single-CPU memory limitation by distributing the computation over a number of nodes in a CPU-based cluster supercomputer.
Alternatives to the hierarchical method include the use of polygonal cells (Mediavilla et al., 2006), optimised to map areas of the image plane to the source plane thus reducing N rays required for a given accuracy, and the contouring method for point sources (Lewis et al. 1993;Witt 1993), which maps an infinite line plus N * closed loops from the image to source plane. The polygonal cell approach results in its own computational overhead in determining the correct polygonal lattice, and suffers from the need for a substantial number of polygons in the vicinity of caustics, which results in lengthy computation times. The contouring method has limited applicability to extended sources. Wyithe & Webster (1999) presented an extension of the contouring method to cover extended sources, however the additional computations required slowed the code to the point where it could no longer compete with the simplicity of inverse ray-shooting.
A recent alternative is the appearance of the GPU as a mathematical co-processor.
Originally developed to enhance the rate of generating graphics for the computer game industry, GPUs are revolutionising scientific computing (Fournier & Fussell 1988;Tomov et al. 2003;Venkatasubramanian 2003;Owens et al. 2005). Their highly parallel processing architecture, accessed through flexibile programming libraries such as NVIDIA's CUDA 2 or the crossplatform OpenCL 3 standard created by the Khronos Group, is enabling speed-ups of O(10) to O(100) times over a broad range of algorithms. Indeed, the notion of general purpose computing on graphics processing units (GPGPU) has motivated a reinvestigation of time-consuming processing tasks.
From a conceptual and implementational point of view, direct ray-shooting is a very simple algorithm with a high degree of parallelism: the deflection of each light ray is independent of every other light ray, and the deflection of a light ray due to a single lens is independent of all other lenses. The existence of such fine-grained parallelism means that ray-shooting can be considered an "embarassingly parallel" algorithm -ideal for GPGPU. Thompson et al. (2010) reported on a CUDA implementation of direct ray-shooting, and demonstrated that their code could calculate a magnification map for the Φ fiducial scenario in ∼ 7 hours using a four-GPU NVIDIA S1070 Tesla unit. Moreover, they were able to calculate a magnification map with one billion lenses (N pix = 512 2 and N rays = 42) in 24 hours, with a peak sustained processing rate of 1.28 Tflop/sec. While the CUDA implementation of direct ray-shooting gave an O(100) processing speed-up compared to using the same method on a single CPU, a critical comparison omitted from Thompson et al. (2010) was between the tree code and the GPU code.
In this paper, we compare: • The "industry standard" single-core tree code by Wambsganss (1990;1999), which we refer to as CPU-T; • The parallel, large-data tree code by Garsden & Lewis (2010), which we refer to as CPU-P; and • The CUDA direct ray-shooting code by Thompson et al. (2010), which we refer to as GPU-D.
We now describe salient points of these three approaches in more detail, including implementation issues that are relevant for the code comparison in Section 3.

CPU-T
The most computationally intensive part of the ray-shooting algorithm is the calculation of the deflection of individual rays by each lens. The hierarchical tree code attempts to minimise the number of calculations necessary by grouping together lenses that are sufficiently distant from a given light ray. A group of lenses is treated as a single pseudo-lens with a total mass equal to the sum of the individual lens masses, located at their centre of mass.
The first step in the calculation is to divide up the lenses into a cell structure, starting with the root cell, which is the size of the lens plane. Each cell that contains more than one lens is subdivided into four cells with half the side length of the parent cell. This process is continued iteratively until only cells containing zero or one lens remain. This calculation need only be carried out once per magnification map, and typically takes only a small fraction of the code's total run-time. A determination must be made for each light ray regarding which lenses need to be included individually and which can be grouped together into pseudo-lenses. The lens equation then becomes: where N L is the number of lenses to be included directly, and N PL is the number of pseudo-lenses with total mass m PL j and centre of mass position x PL j . We note that the actual implementation of the tree code uses higher order multipoles of the mass distribution in order to increase accuracy. This calculation scales on average as O(logN * ).
Additional techniques are used in the standard hierarchical tree code in order to improve the accuracy and speed of the calculation. For example, "test rays" are used to determine the cell structure, and then the same cells are used for a number of actual rays surrounding the test ray. Also, the deflection due to pseudo-lenses is only calculated directly for the test rays; for the actual rays, pseudo-lens deflections are interpolated between those calculated for the test rays. For full details, see Wambsganss (1999).

CPU-P
A tree code requires the following data structures to be stored in memory: 2-D magnification map array; lens array (mass and 2-D location for each lens); and the cell tree (stored as an array of indices), which provides a significant memory overhead as N * increases. In order to increase the scale of microlensing simulations beyond those that were achievable with CPU-T, where typical single-core memory permits a maximum of 20 million lenses, Garsden & Lewis (2010) implemented a distributed, parallel version of the tree code (CPU-P). They take advantage of the inherent parallelism of the ray deflections by distributing rays between multiprocessors of a CPU-based cluster supercomputer, using a simple but effective approach to static load balancing between computing nodes. To handle the large memory requirements for billion-lens simulations (∼ 200 GB), they use a combination of local memory, distribution of data between nodes (thus sharing the memory requirements), and disk-based storage (although this limits processing speeds somewhat by the requirement of costly file accesss).
The CPU-P code implements the same cell-tree and microlensing algorithms as CPU-T, but surrounded by a quite different environment within a running simulation. The generation of lens masses at the beginning of the run, and the ray-shooting itself, use as much parallelism as is available. The generation of the cell tree has not been parallelised, but is executed concurrently with another phase in the simulation. All these operations are logically the same as in the CPU-T version but execute without knowledge that they are running concurrently, or that their data may be contained within huge disk files. Disk files naturally degrade performance and parallel speed-up, but there are various caching schemes which can make up for that. CPU-P has been built to have the same accuracy as CPU-T, because the same numerical algorithms are used, however due to the reimplementation in a different programming language (FORTRAN for CPU-T, C++ for CPU-P), with different random number generators, and a different precision of integers (arbitrary precision has been implemented), there will naturally be minor differences.

GPU-D
The GPU parallel implementation of direct ray-shooting is achieved by splitting the total number of light rays into N B batches, which are deflected in parallel. For convenience, the deflection calculation of equation (4) is recast as: where T t=1 N t = N * for N T parallel threads, and each processing step considers only a single light ray. The contributions from γ and σ c are determined at the end of the deflection calculation.
The GPU-D code, described in Thompson et al. (2010), was implemented in CUDA according to the following algorithm: 1. N * lens coordinates are obtained (see below) and stored in CPU memory. Coordinates are loaded into GPU device memory as required; 2. N B = 2 17 light ray coordinates are pseudo-randomly generated on the CPU and loaded into GPU device memory; 3. GPU computation is initialised; computation is split into groups of N T = 128 threads; 4. Each thread group loads N T lenses and calculates deflection on N T rays; this is repeated until all lenses and rays are exhausted; 5. Once computation on the GPU is complete, results are copied back from the GPU device memory to system memory; 6. CPU maps the ray locations onto the source pixel grid in order to obtain the magnification map; 7. Steps 2-6 repeated until an average of N rays rays per source pixel is reached.
The only change in approach from that of Thompson et al. (2010) is the way the N * lenses coordinates are generated: for this work, we first run a tree code and output lens positions to a data file. These coordinates are then read from the data file by the GPU-D code, and subsequent processing is carried out as described above. Without this modification, we would not be able to make pixel-by-pixel comparisons of the magnification maps, and would be left with just statistical comparisons based on, for example, magnification probability distributions. For full details on choice of optimisation of GPU parameters, especially N B and N T , see Thompson et al. (2010), as these are somewhat hardware dependent.

Comparing the methods
In this section, we consider the relative accuracies of CPU-T, CPU-P and GPU-D over a much wider range of parameter space than has been considered previously, and we compare processing times between CPU-T and GPU-D, noting regions of parameter space where each code might currently be preferred. The former aspect is important in assessing the validity of all microlensing results that have relied on the tree code to date, the latter in motivating a choice of code for future applications. The choice of model parameters for comparison is somewhat arbitrary. We attempt to cover the majority of σ and γ parameter space typically probed by realistic lens models, as well as typical map sizes and pixel scales used in those analyses (see for example: Witt et al. 1995;Kochanek 2004;Morgan et al. 2006;Bate et al. 2008;Mediavilla et al. 2009).

Accuracy
Since the tree code uses an approximation to the lens locations, in some sense the direct approach produces a magnification map that is more accurate for a given configuration of lenses. An unanswered question is whether the hierarchical approach maintains sufficient accuracy as the number of calculations, Φ, is increased beyond approximately 2 × 10 14 . This restriction is largely historical, due to the impracticality of using the direct ray-shooting approach on a single-core CPU beyond this limit.
We conducted simulations for a range of parameters, summarised below: • N pix = 1024 2 , 2048 2 and 4096 2 ; • N rays = 100, 200, 500, 1000 and 2000; and These parameter choices allow us to probe the range: where we assume that N FLOP = 10. We choose all lenses to have the same lens mass as this does not influence the run-time.
The CPU-T/CPU-P accuracy parameter was set to 0.60 for all of our simulations, as this was found to reliably produce acceptable magnification maps. CPU-P maps were generated using a minimum of 2 compute nodes, to ensure that the parallelisation was not producing inaccurate results. There is a subtle difference between the way N rays is defined for CPU-T/CPU-P and GPU-D. In both codes, the same number of rays are shot across the entire shooting grid. The CPU-T/CPU-P rays are shot on a regular grid, thus in the absence of lensing exactly N rays are received at each CPU-T/CPU-P pixel. The GPU-D rays are shot from random positions in the shooting grid, and so on average N rays are received at each GPU-D pixel. This will introduce Poisson noise into the GPU-D magnification maps, however the total number of rays shot is so large (10 8 to 4 × 10 10 in our simulations) that any error introduced by shooting rays randomly will be insignificant.
We use two methods for assessing accuracy: the root mean square difference between corresponding magnification maps, and a Kolmogorov-Smirnov (K-S) test on probability distributions for change in magnitude ∆m. The K-S test results in a measurement of the significance with which we can reject the null hypothesis, that the two samples were drawn from the same distribution. Change in magnitude is calculated from the (per pixel) map magnifications using the following relationship: Throughout, we define "excellent agreement" as a root mean square difference between maps of less than 0.2 magnitudes (often much less than 0.2 magnitudes), and a K-S test probability that the null hypothesis is not rejected of greater than 99.7%.

No shear
To begin with, we set both the external shear, γ, and the smooth matter fraction, s = σ c /σ, to zero. This allows us to consider just the effect of compact objects, which are directly related to the number of lenses via equation (5). Later, we allow both γ and s to vary.
As indicative results, we present the calculated RMS values and K-S probabilities, P(KS ), for the N pix = 4096 2 case in Table 1. For σ ≥ 0.3, we find excellent agreement between the magnification maps from CPU-T and GPU-D when N rays ≥ 500. The σ = 0.1 and 0.2 cases differ slightly. Here, a low number of lenses (N * = 65 and 163 respectively) mean that there is very little variation in magnification across the maps. The K-S test probabilities are therefore strongly affected by small differences in magnification probability distributions. For N pix = 1024 2 and 2048 2 , N rays = 1000 are sufficient to obtain agreement between CPU-T and GPU-D. For N pix = 4096 2 even N rays = 2000 are not enough. In other words, a larger number of rays per pixel are required for agreement between CPU-T and GPU-D when σ is low (≤ 0.2).
We note also that for a given σ and N rays , the RMS error tends to increase with increasing N pix . For example, when σ = 0.4 and N rays = 500, we obtain the following RMS differences: 0.039 mag (N pix = 1024 2 ), 0.081 mag (N pix = 2048 2 ) and 0.168 mag (N pix = 4096 2 ). In other words, for a fixed accuracy we must increase N rays by a factor of ∼ 2 if we increase N pix by a factor of 4.
For the comparisons between CPU-P and GPU-D, we find the same results. In fact, the CPU-T and CPU-P maps are virtually indistinguishable down to the accuracy of our comparison. They therefore behave identically when compared with GPU-D.

Shear and convergence
We then fixed the convergence at σ = 0.45, and examined the impact of varying shear over the range 0.0 ≤ γ ≤ 1.0 in increments of ∆γ = 0.1. The number of pixels was fixed at N pix = 2048 2 . Once again, the smooth matter fraction s was set to zero. We find very similar results to the no-shear case discussed above: for N rays ≥ 500, the results of CPU-T and GPU-D are indistinguishable (RMS difference < 0.2 magnitudes and K-S probability > 99.7%), and similarly for comparisons between CPU-P and GPU-D.

Smooth matter fraction
Finally, we test the accuracy of CPU-T and CPU-P when the convergence is split into a continuously distributed component σ c and a compact stellar component σ * . We set γ = 0.5 and σ = σ c + σ * = 0.45, and vary s in the range 0.0 ≤ s ≤ 0.9 in increments of ∆s = 0.1. Since there is no microlensing if s = 1, we also simulated the s = 0.99 case. For convenience, we choose N pix = 2048 2 . Table 1: Comparison between the CPU-T, CPU-P and GPU-D codes for magnification maps with N pix = 4096 2 , γ = 0.0 and s = 0.0 (corresponding to σ c = 0.0). RMS is the root mean square difference between magnification maps (compared pixel by pixel), and P(KS ) is the K-S test probability. The CPU-T and CPU-P comparisons are virtually indistinguishable.
CPU-T and GPU-D CPU-P and GPU-D N rays σ RMS P For s ≤ 0.8, we find the same results as the two cases discussed above: for N rays ≥ 500 the results of CPU-T and GPU-D are in excellent agreement. The s = 0.9 case has a slightly higher RMS difference at N rays = 500 of 0.4, although the K-S probability is still > 99.7%. For the s = 0.99 case, however, no agreement is reached between CPU-T and GPU-D. Visual inspection of the maps shows that too few rays have been shot in the CPU-T case. This suggests that more rays should be shot in CPU-T for higher smooth matter fractions; N rays = 1000 should be sufficient for s = 0.90, however more than N rays > 2000 is required in the s = 0.99 case. Once again, CPU-P performs identically to CPU-T.
By comparing the accuracy of the two tree codes, which rely on a series of approximations to the lens distribution in order to speed up calculation time, with GPU-D, which conducts the entire brute-force calculation, we have verified the accuracy of the standard tool for high optical depth microlensing simulations. Indeed, we have considered a wider range of parameters, as represented through the quantity Φ, than previous work (Φ > 2 × 10 14 ). While not unexpected, our results confirm that the tree code can indeed be used with confidence, provided a sufficiently large number of rays per pixel are shot (i.e. typically N rays ≥ 1000 for the parameter space explored here).

Timing Tests
Having established accuracy of the three codes under consideration, we now turn to timing tests. Specifically, we determine the time taken to generate magnification maps with a given set of model parameters, allowing us to determine the regions of (γ, σ, N pix ) parameter space where one of the implementations may be favoured.
We do not consider timing comparisons between CPU-T and CPU-P, as these are reported in Garsden & Lewis (2010). Since CPU-P run-times scale with the number of nodes available, it is difficult to establish which configuration provides the most appropriate comparison between CPU-P and GPU-D. Using the hardware cost as a factor in choosing between software alternatives, we suggest that the cost of purchasing a suitable high-end graphics card plus computer is much lower than the set-up costs for a modest CPU-based computing cluster, once networking overheads and cluster management are considered. We reflect further on the implications this has for the future of microlensing computations in Section 4.
Care must be taken when quoting absolute speeds for different hardware -we use our results as a guide to the relative processing speeds of CPU-T and GPU-D, but note that alternative hardware can change timing by factors of a few. Our benchmarks were conducted on the following equipment: • CPU-T Intel Xeon Dual Core 3060 CPU (2.4 GHz) with 1 GB RAM.
• GPU-D NVIDIA S1070 Tesla unit connected via dual PCIe x8, or higher, buses to an Intel Q6600 Quad Core CPU (2.4 GHz) with GB RAM. The Tesla S1070 comprises four GPUs, each with 240 stream processors, running at 1.296 GHz. Each GPU is able to access 4 GB of RAM. The quoted peak processing rate is 2.488 TFLOP/s, however, for GPU-D, a peak of 1.28 TFLOP/s was obtained using all four GPUs (Thompson et al. 2010). We note that the majority of our timing benchmarks in the present work were conducted using only two of the four GPUs available on the Tesla S1070 4 .
Since we have determined that the magnification maps produced by CPU-T are accurate for N rays ≥ 1000 across the majority of parameter space, we conduct our timing tests by comparing the total time taken by CPU-T and GPU-D to generate identical maps with N rays = 1000, for each combination of parameters. The results of these tests are discussed below.

No shear
Timing results for N pix = 1024 2 , 2048 2 and 4096 2 , with γ = 0 and s = 0, are plotted in Figure 2. Each panel shows both σ (along the top axis), and the logarithm of N * (along the bottom axis). Note that as s = 0 for these simulations, σ = σ * .
The performance of CPU-T is essentially constant as σ increases. The same is not true for GPU-D; it performs better for lower σ (lower N * ) values. This result is as expected; as the number of lenses is increased, the direct code is required to perform correspondingly more calculations, whereas the approximations in the tree code (specifically the depth to which the cell tree is accessed) mean it performs roughly the same number of calculations at each step.
The shape of the GPU-D timing curves is the same as those presented in Figure 2 of Thompson et al. (2010). For low numbers of lenses (N * 10 4 ), overheads induced by the need to transfer data between GPU and CPU affect the performance of GPU-D. For N * 10 4 , near the limit of our comparison, the performance of the code scales linearly with N * . While this linear growth suggests that there is a cross-over in run-time between CPU-T and GPU-D at around N * = 10 5 , it must be recalled that this is approaching the maximum lens limit of CPU-T, due to the memory requirements for storing the cell tree.
We note that although the performance of GPU-D is decreasing as σ (or equivalently, N * ) is increased, it is still faster than CPU-T for all of our γ = 0.0 simulations. Using all four GPUs on the Tesla S1070 unit, instead of the two we use for this detailed timing comparison, results in factor of two speed-up.

Shear and convergence
Real lensing galaxies are unlikely to be perfectly spherically symmetric, and so realistic lens models do not contain zero external shear. We therefore also conducted timing tests for varying γ. Equation 4 shows that the actual calculation of the effect of external shear should not significantly add to computation time -it simply adds one extra operation to the calculation of the deflection of each ray. However, the side lengths L ± of the rectangular shooting grid that is projected on to a square magnification map with side length L map depends on the shear in the following way: When 1 − σ − γ ≈ 0 or 1 − σ + γ ≈ 0, one side of the shooting grid can be very large. This leads to a large area for the lens plane, and hence a correspondingly large number of lenses (see Equation 5 and the discussion following it). We conducted timing tests for the case where N rays = 1000, N pix = 2048 2 , s = 0 and σ = σ * = 0.45, with 0 ≤ γ ≤ 1.0 and ∆γ = 0.1. The results of these tests are presented in Figure  3. The figure shows γ versus the logarithm of the time taken to generate the map. We also plot time against the theoretically predicted average magnification µ av in that figure, where µ av is given by: A negative magnification is interpreted as a parity flip in the image. For low magnifications (when 1 − σ − γ and 1 − σ + γ are very different from zero), GPU-D is faster than CPU-T. There is a crossover point at a magnification of ∼ 5 where CPU-T begins to out-perform GPU-D. For high magnifications, CPU-T is significantly faster than GPU-D. We emphasise again that these timing tests were conducted using only two of the available four GPUs. Running on all four GPUs should improve the performance of GPU-D by approximately a factor of two.
The timing behaviour of GPU-D with the addition of a nonzero shear is exactly the same as the zero shear case. Initially, overheads dominate the code's performance, but as the number of lenses is increased to ∼ 10 4 , the code begins to scale linearly with N * .
The same is not true for CPU-T -the timing behaviour in the zero shear and varying shear cases seems to be quite different. This is due to the way in which the shooting grid is stretched by the addition of a non-zero shear. Equation 11 shows that for certain combinations of σ and γ the shooting grid will be significantly stretched in one direction. A ray shot from one corner of this stretched grid is close to fewer lenses than a ray shot from a corner of a square with the same area and thus a greater number of lenses can be grouped together into pseudolenses. Although N * goes up for these cases, the number of lenses that need to be included directly in the calculation goes down for the tree code.
We expect this behaviour to be generic as σ is varied -CPU-T will be faster than GPU-D in the high magnification regime on current-generation hardware, whereas GPU-D will out-perform CPU-T in the low magnification regime. However, as we discuss in Section 4, the former situation may only be a limit in the short-term.

Smooth matter fraction
Finally, we present the results of the timing tests for varying smooth matter fraction in Figure 4. Tests were conducted for N pix = 2048 2 , σ = 0.45, γ = 0.5, and N rays = 1000. This combination of σ and γ corresponds to the "worst case" magnification for GPU-D timing, and the "best case" magnification for CPU-T timing (see the previous section).
Increasing the smooth matter fraction for fixed σ corresponds to decreasing N * , as mass is transferred from the compact stellar component to the continuously distributed component. We therefore see exactly the timing behaviour we expect: CPU-T takes roughly the same amount of time for each s, whereas GPU-D is faster for fewer lenses.
At the distances from the centre of lensing galaxies where lensed images are typically observed (∼ 2 − 10 kpc), we expect the smooth matter fraction to be quite high. Indeed, a few rough microlensing measurements of smooth matter percentage at lensed image positions exist: 80% to 85% at the location of the A and D images in SDSS J0924+0219 (Keeton et al., 2006); ∼ 80% in HE 1104+080 (Chartas et al., 2009); ∼ 90% in PG 1115+080 (Pooley et al., 2009); and ∼ 70% in RXJ 1131-1231 (Dai et al., 2010). For smooth matter fractions in this range, we find that GPU-D already performs comparably to CPU-T in high magnification systems, and out-performs it in low magnification systems.
We note that our accuracy tests in Section 3.1.3 showed that even N rays = 2000 were not sufficient for CPU-T to obtain results indistinguishable from GPU-D at s = 0.99, at least with the CPU-T accuracy parameter used throughout this paper. The s = 0.99 point in the Figure 4 CPU-T results is therefore likely to be lower than the actual run-time required to produce a useable map.

Discussion
Our intention in this work is not to imply that the tree code is excessively inaccurate. However, its reliance on an approximation to the lens distribution to reduce the number of calculations per magnification map warrants further investigation that has not been feasible until now. The standard tree code (Wambsganss 1990;1999) has been in widespread use throughout the lensing community since the early 1990s. Though it is unsurprising, it is nevertheless reassuring to know that CPU-T remains accurate over an extended range of Φ-space. Our comparison here demonstrates that researchers can, and should, continue to use the tree code for microlensing simulations with confidence. Additionally, we have demonstrated that this accuracy is also maintained by the CPU-P code, which enables exploration of problems orders of magnitude beyond those that can be achieved on a single-core CPU.
However, we stress that the direct inverse ray-shooting code is more accurate than its tree code counterparts by simple virtue that it includes contributions from all lenses without approximations -the entire calculation is done by brute force, taking advantage of the massively parallel GPU architecture to complete the computation in a reasonable timeframe. While CPU-T can indeed be used for a wide range of astrophysically relevant scenarios, GPU-D is now a highly competitive alternative.
Given that we have verified the accuracy of the tree codes, though, it is reasonable to ask in which areas of parameter space each code should be used. The larger the number of lenses, the better the tree codes perform relative to GPU-D. For the values of N pix and N rays we examined in our benchmarking, the crossover point where CPU-T begins to out-perform its GPUbased counterpart is N * ∼ 50, 000 lenses. This can occur either when (σ − γ) or (σ + γ) are close to 1. For scenarios with N pix > 4096 2 and N * > 10 5 , there are advantages in moving to a distributed CPU-based computing cluster, although this can be a financially costly approach. One of the perceived benefits of GPU computing is the heavily reduced dollar per Tflop/s they offer, which makes them an attractive alternative as "desktop" supercomputers (see Thompson et al. 2010). The downside of GPUs, at this point in time, is the paucity of astronomy codes that have been adapted to run on them, whereas CPUbased computing clusters can support a much wider range of existing codes.
Since about 2005, single-core CPU processing speeds have plateaued.
Additional computing performance has been achieved, and indeed Moore's Law growth in processing power has been maintained, by moving to multi-core CPUs. A growing range of applications across a variety of scientific disciplines are now taking advantage of the "many-core" processing architecture of GPUs, which are providing a preview of the future of CPU architectures. Yet if GPUs do indeed point towards a future where many-core processing is mainstream, the need to convert existing trusted, tested, accepted codes will become more pressing. A question that many researchers will have to answer is which code should be adapted?
We suggest that decisions to move other standard astronomy computations to GPU may benefit from a similar approach to Thompson et al. (2010). Rather than attempting to directly port an existing advanced code (i.e. CPU-T) to GPU, they considered a brute force approach that was not feasible for singlecore CPU, but which exhibited the type of fine-grained, massive parallelism that is ideal for GPGPU. This resulted in much faster code development time, despite the learning curve associated with adopting CUDA and the specifics of GPU parallel programming. Moreover, we have shown that in most astrophysically-motivated cases, we are already no worse off using GPU-D compared with CPU-T, with the added benefit of greater intrinsic accuracy.
The N * -crossover in timing tests that we have identified does depend on the hardwared used. For example, our benchmarking was performed using two of the four GPUs available on the Tesla S1070, while the CPUs used for CPU-T were a few years old. The key point, however, is that the GPU-D is already parallel -GPU performance will increase, and more GPU cores will be placed in individual machines, yet no further code development is required. The single-core tree code, however, is more heavily restricted -CPU-T can no longer make significant processing gains by simply relying on Moore's Law to make single-core CPU hardware faster.

Conclusion
In this work, we have examined the accuracy of three highly competive alternatives for generating microlensing magnification maps: a widely-used hierachical tree code for single-core CPUs; a parallel, large data implementation for a distributed CPU-based computing cluster; and a brute force solution that utilises the massive parallelism of modern graphics processing units. We have demonstrated that these three codes do indeed present comparable accuracy, and choice of approach will depend on considerations relating to the scale and nature of the problem to be investigated.
While considered impractical as a calculating approach even a few years ago, our timing tests show that direct, inverse rayshooting is now a viable, and indeed more accurate, option. Moving a simple algorithm to GPU, rather than attempting to directly implement a more advanced solution -originally developed to overcome a hardware restriction of single-core CPUsto achieve run-times that are already similar, suggests an approach that other early adopters of GPUs for astronomy computation may wish to consider.


extfieldsofstudy:
Computer Science
Medicine

text:
Insights from unifying modern approximations to infections on networks

Networks are increasingly central to modern science owing to their ability to conceptualize multiple interacting components of a complex system. As a specific example of this, understanding the implications of contact network structure for the transmission of infectious diseases remains a key issue in epidemiology. Three broad approaches to this problem exist: explicit simulation; derivation of exact results for special networks; and dynamical approximations. This paper focuses on the last of these approaches, and makes two main contributions. Firstly, formal mathematical links are demonstrated between several prima facie unrelated dynamical approximations. And secondly, these links are used to derive two novel dynamical models for network epidemiology, which are compared against explicit stochastic simulation. The success of these new models provides improved understanding about the interaction of network structure and transmission dynamics.

INTRODUCTION
The vast majority of infectious diseases can be considered as spread through a network of contacts between individuals or groups. These network concepts have been used, to great effect, in scenarios where the network of relevant contacts can be readily ascertained: for sexually transmitted infections, the network of human sexual contacts [1 -3]; and for bovine diseases, the network of animal movements as captured by the Cattle Tracing System [4 -6]. More recently, attention has focused on the implications of network structure for human infections transmitted through close contact (such as influenza, SARS and smallpox) using diarybased social encounter information [7], contact tracing [8] or predicted patterns of movements [9] to infer the appropriate network structure. In all cases, it is important to acknowledge the effects of network structure, both because of its impact on the uncontrolled epidemic [10,11] but also because of its utility in targeting controls and tracing the spread of infection [12]. It is, therefore, important that we develop methods of modelling infection dynamics on complex structured networks. Stochastic simulations obviously provide the most accurate and versatile models, but at the expense of tractability. There has, therefore, been considerable focus on developing less computationally intensive approaches that can help us interpret the difference between network-based predictions and those of mean-field (random mixing) models where local structure is ignored.
A recent review paper [13] counterposed three dynamical approaches to the study of epidemic dynamics on networks. These are: pairwise approaches [10], dynamic probability generating function (PGF) formalism [11,14] and heterogeneous mixing [15,16]. Following this approach, we demonstrate here that each of these three approaches is an approximation to the more general pairwise model of Eames & Keeling [17]. We then use the relationships between approaches to derive two new moment-closure-based models for network epidemics: a 'clustered PGF' model, which is capable of capturing epidemic dynamics on clustered networks of heterogeneous link distribution using a relatively small number of ordinary differential equations (ODEs), and a 'heterogeneous susceptible-infectious -susceptible (SIS)' model, which makes a significant but less dramatic reduction of dimensionality where acquired immunity is not long-lasting. Finally, we compare these two models to simulation on exemplar networks similar to those considered in more applied contexts.

DERIVING OTHER APPROACHES FROM PAIRWISE MODELS
We start by considering a general heterogeneous contact network with N nodes. Using the notation developed in Keeling [10] and Eames & Keeling [17], we use square brackets [] to represent the expected numbers of nodes, pairs or triples of any particular type; all notation used is summarized in table 1. The local structure of this network can be defined in terms of three main measures: the degree of distribution d k (where d k is the proportion of nodes having k network contacts), the clustering of contacts f that measures the ratio of triangles (groups of three nodes all connected to each other) in the network to triples of all types (lines of three nodes with or without a transitive link) and the assortativity within the network as captured by the matrix where n is the mean node degree, which compares the true number of pairs of a given degree with the expected number if half-links connected at random. The matrix C represents the extent to which the number of [kl] pairs is over-or under-represented with respect to a random process of pairing half-links [18]. We wish to explore how these three forms of local structure influence the types of model formulation that can be successfully applied.
The full pairwise equations of Eames & Keeling [17], from which we begin our analysis, are given for SIR-type infections by and ð2:1Þ where [A k ] refers to the number of nodes of type A with degree k; a disease state without a subscript implicitly contains the sum over all possible degrees (e.g. ; t is the transmission rate; and g is the recovery rate. (The equations for infections that obey SIS-type dynamics are derived by modifying how recovery acts.) As with many moment-based methods, these equations are exact but unclosed. If we continue to write higher-order equations for triples, quads, etc., then these equations will close only once we reach the network size N, destroying the main motivation for the use of ODEs to describe epidemics on networks. We, therefore, must seek approximations that will allow us to close these equations at a lower dimension (approximating the number of triples in terms of pairs and singles), which will aid computation and analytic understanding.

Dynamical and network assumptions
We present now a series of different assumptions about network structure and epidemiological dynamics that allow other approaches to be derived from this general pairwise model, discussing the reasoning behind each assumption and its range of validity. We pay particular attention to the number of equations required to simulate SIR dynamics in terms of the maximum node degree M.
2.1.1. Triple closure. Closure schemes for pairwise models, which have become essentially standard, have been presented for networks of significant clustering but with homogeneous degree [10] and also for unclustered networks with heterogeneous degree [17]. The number of pairs with one member having degree k, and with the other having degree l (equal to correlation matrix between nodes of degree k and degree l f clustering coefficient of the network (equal to the number of triangles divided by the number of triples) number of pairs with one member in state A and with degree k, and with the other member in state B and with degree l number of pairs with one member in state A and with degree k, and with the other number of pairs with one member in state A, and with the other member in state B (equal to

number of triples with one edge member in state
A and with degree k, with the middle member in state B and with degree l, and with the other edge member in state C and with degree m u(t) the fraction of degree one nodes that remain susceptible at time t Y(t) auxiliary variable used in clustered PGF model (equal to PGF for the network degree distribution (equal to P d k x k ) t rate of transmission of infection across a network link g rate of recovery from infection natural combination of these two schemes is where A, B and C stand for arbitrary disease states. While this closure is the natural extension of existing approximations, it does pose a question of interpretation compared to existing closures, since the term proportional to f cannot be rigorously interpreted as a prevalence of triangles owing to the lack of symmetry between the k, l and m nodes. An extended discussion of how to interpret asymmetric clustered closure appears in House & Keeling [19]. It is also worth noting that the unclustered pairwise equations can be reinterpreted in terms of neighbourhoods [20], although this approach has yet to be extended to the clustered case. The closure scheme (2.2) acting on the full pairwise equations (2.1) produces a closed ODE system of (5M þ 2)M independent equations for SIR dynamics, which can quickly outstrip even modern computational resources for graphs with 'fat tailed' degree distributions such as scale-free networks, where M is often very large. This motivates the investigation of further approximations that can allow us to reduce the system size.
One additional, commonly made assumption that will be used in deriving other approaches is the absence of triangles in the network, that is, The assumption introduced by Eames & Keeling [17] was to express the joint probability of a fully described pair as a product involving pair types and network structure: Under this assumption, pairs with a single indexing degree are the main variables, which significantly reduces the number of independent ODEs needed in the pairwise equations to 10M for SIR dynamics. A priori the accuracy of this assumption is not assured, since it depends on the precise dynamical process taking place on the network; however, for all epidemiological scenarios so far considered in the literature, the results using this approximation are in extremely close quantitative agreement with the full pairwise model. Most importantly, this assumption allows the consideration of assortative mixing and heterogeneous degree, in a system where a number of differential equations is linear rather than quadratic in the maximum node degree.
2.1.3. Detailed balance. The most general prevalence of pairs of connected nodes with degrees k and l is given by that is, if C k, l can take any value for combinations of k and l, then all possible distributions for [kl] can be enumerated. However, an assumption is often made that the correlation matrix C obeys C k;l ¼ 1 8k; l: ð2:5Þ In mean-field (non-network) models, such an assumption has been used, to great effect, to study the spread of sexually transmitted infections in risk-structured populations defined by individual-level data on sexual contacts [21]; however, it is recognized that, in general, most populations are assortative such that C k, l . 1 for similar k, l and C k,l , 1 for dissimilar k, l. In network models, the assumption that C k, l ¼ 1 means that link 'stubs' from each node are connected randomly, which is typically called detailed balance with respect to swapping randomly picked edges in physical science [18], although terminology and definitions concerning assortativity can be different in other subject areas, for example, the sexually transmitted infection literature (e.g. [22] which allows us to keep assortativity but loses the effects of network structure-the correlation between the states of connected nodes is lost-although some of these can be maintained through extra factors, as in Kiss et al. [23]. Application of this closure reduces the number of equations needed to 2 M for SIR dynamics.

Deconvolution of individuals.
The final approximation we consider has not been previously explicitly stated in this form, but can be used to derive the PGF approach from a general pairwise model. It can be viewed as the main assumption that allows the PGF formalism to approximate disease dynamics on heterogeneous networks within a low dimensional framework. This involves writing the joint probability [A k B] (defined above) as a product, and is analogous to the concepts used to generate equation (2.4), in that its validity depends on the independence of dynamics and network structure.
The underlying approximation is that given a node (of type A k ), the type of connected node (taken here to be B) is independent of k and therefore independent of the local network structure. For SIR dynamics, this assumption creates a set of equations whose dimension does not depend on the maximum node degree M, although the exact number of equations required depends on whether clustering is present, and which quantities one wishes to calculate over the course of an epidemic.

Relating ordinary differential equation-based approaches
Starting with the general pairwise model (2.1), the assumptions above can be used to derive other Approximations to infections on networks T. House and M. J. Keeling 69 approaches. Taking the assumptions (2.2) -(2.5) and (2.7) allows us to derive the PGF equations originally formulated from first principles in [14] (full equations given in electronic supplementary material). This is an interesting and unexpected result since the underlying arguments used originally to derive the pairwise and PGF models are quite different. The standard heterogeneous-mixing models associated with riskstructured populations (given in electronic supplementary material) are derived by putting the assumption (2.6) into the pairwise equations and ignoring the pair-level variables. We also note that for networks with no clustering and degenerate degree distribution (often called either regular graphs or homogeneous random networks), the PGF and pairwise approaches are formally identical, since assumption (2.7) is trivially satisfied.

CONSTRUCTION OF NOVEL MODELS
While the formal links between prima facie distinct epidemic models are intrinsically interesting, the main motivation for our work is to derive novel, parsimonious models for epidemics on a range of complex networks without having to argue the formulation from first principles. The fact that full pairwise models are simply written down, together with the new explicit form of the assumption underlying PGF models (2.7), means that it is possible to produce systematically much simpler ODE-based models for network-based epidemics. We now consider two such models, together with the potential limits on this methodology.

Incorporation of clustering
The relationships between different ODE approaches to network epidemics presented above open up the intriguing possibility of extending the PGF approach to include clustering so that both heterogeneity in link distribution and clustering can be analysed using a low-dimensional model with a small number of dynamical variables. This is indeed possible: using (2.4) and (2.7) together with the standard pairwise closure (2.2), the two triples that appear in the unclosed pairwise SIR equations can be closed through  determines the epidemic behaviour: and _ ½II ¼ 2tð½ISI þ ½SI Þ À 2g½II : ð3:3Þ The number of susceptibles can be tracked nondynamically through the relationship S ½ ¼ NgðuÞ: ð3:4Þ The primary significance of this result is that, regardless of the maximum node degree on a network, a system of six ODEs together with the non-dynamical relationship between [S] and u can be used to calculate an expected prevalence curve for heterogeneous, clustered networks.

Incorporation of other disease natural histories
In addition to the clustered PGF model derived above, which assumes SIR (susceptible -infectious -recovered) dynamics, other extensions of the PGF approach will clearly be possible, for example, the inclusion of a latent class or multiple infectious classes as in Kamp [24]. While this requires more variables, the system will still be of far lower dimension than the equivalent full pairwise model. However, other commonly studied disease behaviours, such as the SIS paradigm, which is ideal for sexually transmitted infections and has been extensively studied using pairwise approaches, are not straightforwardly reduced to the PGF formulation. This is because, under SIS dynamics, the set of ODEs governing the evolution of quantities such as P k k a [A k ] is not closed. This does not, however, preclude the use of approximation (2.7). Inserting this assumption into the standard pairwise SIS model gives the following new model: ð3:5Þ Unlike the clustered PGF model above, the number of these equations that need to be manipulated numerically is M þ 1, that is, linear in the maximum node degree; however, their dimensionality is still significantly lower than standard pairwise models and, perhaps more importantly, there are only two equations that need to be manipulated in analytic work on this system.

Assortativity
Deviation from null assortativity is a complexity that simply cannot be naturally incorporated in the PGF Approximations to infections on networks T. House and M. J. Keeling 71 equations without vastly increasing the dimensionality of the system. We provide a numerical demonstration of the importance of this observation in the electronic supplementary material. Nevertheless, as shown in Eames & Keeling [17], by using networks (2.2) and (2.4), an assortative model can be created that is linear in maximum node degree.

Comparison with simulation
We now test each of the two new models above (clustered PGF and heterogeneous SIS) against simulation and compare to other relevant ODE approaches.
For the clustered PGF, we start with a random graph of approximately 10 4 nodes and with Poisson parameter l ¼ 6. We then introduce a clustering coefficient of f ¼ 0.2 using the 'big V' rewiring [19,25], which does not change the degree distribution. Following [26], we bias node selection by k(k 2 1), so individual-level clustering is constant; while this will make C kl slightly deviate from unity, at the level of clustering we consider this should have negligible dynamical impact. We choose these network parameters because they separate the ODE approaches without being so large that concerns about global network properties like giant component size are posed. While it is probable that realistic networks for respiratory infection have more variable degree distributions, larger numbers of mean nodes and higher clustering coefficients, our approach is to start with a system with few underlying parameters and a transparent method for the introduction of clustering. Epidemic rates were t ¼ 0.8, g ¼ 1. Figure 1 shows the comparison of the clustered PGF model, which is in good agreement with simulation, and two other ODE approaches. The homogeneous pairwise model underestimates early growth in cases and, conversely, the unclustered PGF model overestimates early growth. This is as would be expected from general arguments about the effects of clustering [10] and population heterogeneity [15]. The electronic supplementary material shows how these other, existing, ODE approaches remain in good agreement with simulations on appropriate networks and also considers the impact of assortativity.
For the heterogeneous SIS model, we generate a scale-free network of approximately 10 4 nodes by using the standard [27] method with parameters m 0 ¼ 20, m ¼ 2, and removing nodes of degree 0. Since the scale-free property is consistent with observed sexual contact networks [28], and for sexually transmitted infections recovered individuals often fail to acquire long-lasting immunity, pairing scale-free networks with SIS dynamics is natural. Epidemic rates were t ¼ 1, g ¼ 1. Figure 2 shows that models incorporating heterogeneity both fit to simulation much better than the homogeneous pairwise model and that the addition of just one extra equation for the heterogeneous pairwise when compared with heterogeneous mixing significantly improves the fit to early growth behaviour.

DISCUSSION
We have analysed the conditions under which PGF and heterogeneous-mixing models can be derived from a general pairwise approach. These split into assumptions about the network itself, such as zero clustering (2.3) and non-assortativity (2.5), and assumptions about the interaction of dynamics with network structures such as (2.2), (2.4), (2.6) and (2.7). We have used these conditions to derive a clustered PGF model and a low-dimensional heterogeneous SIS model, which are likely to be of significant utility in the study of epidemics on networks owing to their relatively low dimensionality.
In general, the starting point for analysis of disease transmission on a network has to be the available data. Where the full network is known, or can be imputed with confidence, then explicit stochastic simulation of the epidemic process is likely to be the best approach. When only statistical properties of the network, such as degree distribution d k , correlation between degrees C kl or clustering coefficient f, are known, then there are two complementary approaches: either generate exemplar networks for simulation with the appropriate statistics or make use of ODE-based models of the kind considered here, which can be parametrized directly from the network statistics, are numerically tractable and mathematically transparent.
Within this second approach, each ODE-based model used to study network epidemics will have its own domain of validity. While the most general pairwise model can be applied to all compartmental paradigms, with clustering, degree heterogeneity and assortativity included, this comes with increased and potentially unnecessary computational overhead. The more tractable models, on the other hand, may suffer from important inaccuracies when the implicit assumptions they make are not justified.
A question is also posed, however, about the ultimate justification of even the general pairwise model. At present, it is widely believed that in an appropriate regime, this model will be 'exact' in the same way that the standard SIR equations for an epidemic tend to the mean behaviour of an exact stochastic epidemic model in the appropriately constructed large-population limit. Whether such a regime can be defined remains an important open problem, and one that will hopefully be resolved in the near future.
In conclusion, we hope that this study has clarified the relationship between diverse ODE-based models and extended the repertoire of models available for use.