dataesr/scientific-paragraphs-categorization huggingface dataset:
Overview of the data pipeline HAP-Fr is implemented in the historical analysis chain of H. E. S. S. , the H. E. S. S. Analysis Package (HAP). Its architecture relies on a custom framework developed for both the on-line and off-line analysis, SASH (Storage & analysis software at H. E. S. S. ), allowing a seamless integration of the data analysis code for the on-line data quick view during the data acquisition  [2]. Its implementation is made in C++, using the ROOT data analysis framework on which the raw H. E. S. S. data format is based. ROOT offers a serialisation mechanism for the transport and storage of objects and classes. MySQL was chosen to design databases of monitoring, calibration and astrophysical data. The SASH framework and data analysis classes are handled by the C-shell and python scripting languages, allowing to implement the analysis steps for massive datasets.
To derive high level products, custom tools have been developped. The package START  [12]  allows to reconstruct spectral shape of γ-ray sources. It relies on libraries provided by the ROOT framework, MINUIT2, and the algorithms for numerical integration and interpolation from the GNU Scientific Library (GSL). A 'forward-folding' method is used to adjust an assumed spectral shape on data, taking into account the instrument functions  [13]. In order to test the consistency of the fitted spectral parameters, a MC tool was developed within the software, providing a simulated number of 'measured' ON and OFF events in each energy bin for a given observation and a source spectral shape. This allows to test the robustness of the spectral START convergence  [14]. The START package can also convert H.
E. S. S. spectra and responses in OGIP-compliant format (pha, arf and rmf files). We have developed a python package for the Sherpa environment, 'HSPEC', that allows to perform full spectral analyses for large number of files with a profile likelihood technique  [13]. HSPEC is part of GammaPy. A 2D-morphology Maximum-Likelihood fitting package was also developed and has already been used for a number of scientific results (e. g. HESS J1018  [15] ).
Reconstruction and γ/background discrimination This section describes the algorithms used to produce the DL3 data. Events are recorded either in the stereoscopic mode, i. e. when any pair of telescopes among 5 triggers the system, or in a monoscopic mode when only the large dish CT5 triggers. The mono mode is designed to lower the energy threshold of the H. E. S. S. II phase down to 20 GeV for specific scientific topics such as γ-ray pulsars  [5] , γ-ray bursts or high redshift Active Galactic Nuclei (AGN). The overall scheme for event reconstruction and discrimination against background is similar for the mono and stereo modes, albeit some differences as outlined below. The shower images are extracted from the calibrated raw data (DL1) after application of a tail-cuts algorithm with two thresholds ((5,7) photo-electrons or p.
9% and 4. 4% respectively), showing that the cuts efficiencies are under control. For PG 1553+113, the analysis gives about 5. 5σ / √ h, confirming the expected performance derived by simulations. While satisfactory, the results on PKS 2155-304 could be further improved by using an analysis configuration dedicated to brighter sources. For these datasets, the spectra measured with START give very satisfactory results and will be presented in future communications.
Conclusions HAP-Fr is a software pipeline allowing to analyse the H. E. S. S. data taken with the full array composed of 5 telescopes. It processes the raw data by using any H. E. S. S. calibration or simulated products to derive high level products, such as spectra or phasograms. This pipeline has been adapted to account for the running of our array of different telescope types and new algorithms have been set up to analysis data in mono mode (with CT5 only) and in stereo (CT1-5). A special care has been carried on the γ/background discrimination. A multi-variate analysis using boosted decision trees is used together with robust discriminant parameters built with the usual Hillas analysis and Model3D fit of the shower. One forest of decision trees is trained for each interval of the parameter space of discriminant variables and observation parameters.
Their operating points are derived by optimising the signal significance for a fixed γ-ray flux and spectral shape. The performances of the stereo analysis for a Faint source configuration reach the expected ones by comparing them with the preliminary results on AGN datasets. The mono analysis is also quite promising despite modest performances of the energy estimator and the difficulty of the γ/background separation. It allows to treat robustly data on strong and very soft sources, such as the Vela pulsar.
This work describes the challenges presented by porting parts of the gysela code to the Intel Xeon Phi coprocessor, as well as techniques used for optimization, vectorization and tuning that can be applied to other applications. We evaluate the performance of some generic micro-benchmark on Phi versus Intel Sandy Bridge. Several interpolation kernels useful for the gysela application are analyzed and the performance are shown. Some memory-bound and compute-bound kernels are accelerated by a factor 2 on the Phi device compared to Sandy architecture. Nevertheless, it is hard, if not impossible, to reach a large fraction of the peak performance on the Phi device, especially for real-life applications as gysela. A collateral benefit of this optimization and tuning work is that the execution time of Gysela (using 4D advections) has decreased on a standard architecture such as Intel Sandy Bridge.
This work describes the challenges presented by porting parts of the gysela code to the Intel Xeon Phi coprocessor, as well as opportunities for optimization, vectorization and tuning that can be applied to other applications. Section 2 is devoted to the evaluation of Xeon Phi performance using generic micro-benchmarks. This overview of the observed performance and the encountered difficulties on Xeon Phi will help us to understand the bottlenecks and benefits of this architecture. Using these information, we will able in following sections to port the gysela application on this device. To simplify our study, we have extracted a mini-application from the gysela code with a reduced code line count, and with less physics submodels. The mini-application is much easier to modify in order to optimize and to tune the code.
Xeon Phi architecture The benchmarks presented in this document have been mostly obtained on the Helios machine on IFERC computing facility (Rokkasho, Japan). The hardware we used consists of Xeon Phi cards with 60 cores (1. 052 GHz), where each core is capable of executing four concurrent threads. For comparison, the host processor (Sandy Bridge, Intel Xeon E5-2600) consists of two 2. 7 GHz processors with eight cores each and with two admissible threads per core. The types of problems the Xeon Phi is well suited for are intensive numerical computations. Additionally, for better performance the computations can use one of the highly optimized vector math libraries that were implemented using assembly language constructs tuned specifically for the Xeon Phi architecture, or very well vectorized code.
Programming models The developer has currently two different programming models to implement an application on such architecture. In the following, these two programming models will be compared regarding three criteria: impact on the code of an existing MPI application, concurrent usage of host and coprocessor resources and availability in production. The first strategy is the so-called "offload mode", similar to GPUs programming models. It consists in introducing pragma statements in the code around the computing kernels that should be executed on the coprocessor. At execution time, when a process running on the CPUs reaches this portion of the code, the execution on the host stops, data are transfered on the coprocessor and the execution resumes on the coprocessor. Once the execution reaches the end of the pragma section, the execution stops on the coprocessor, data are transfered back on the CPU and the execution resumes on the CPU.
Asynchronous pragmas are available and enable concurrent execution on both the CPU and the coprocessor. MPI is then used to exploit the resources of several nodes. The MPI processes running on the CPUs (typically one per socket) implement these pragmas and offload data and computations on the coprocessors. On the one hand, the advantage of this approach is that it is ready for production today on all the Xeon Phi clusters. On the other hand, to implement such offload scheme efficiently, the developer has to minimize the data transfer between the host and the coprocessor as it can very easily become a strong bottleneck for the performance. Additionally, to distribute computations on both the host and the coprocessor at the same time, the usage of asynchronous offload pragmas is required. This makes the data synchronization between both computation units particularly complex.
At the same time a single Sandy Bridge core is much more powerful than a single Xeon Phi core, so serial parts of the code will execute faster on a Sandy Bridge. So depending on the computation and for a given load, the execution time will likely be different on the host and on the coprocessor. This translates differently in the two programming models. For the native mode, as the large majority of MPI applications distributes equally the same computations with the same load on the different nodes, some tasks will finish earlier than others then leaving their computing unit idle. A strategy to alleviate this would consist in modifying the MPI application such that different groups of MPI tasks can have different loads. But this kind of modifications can be very difficult to implement and is highly application dependent.
But having the hardware does not seem to be enough. Indeed, Supermic is equipped with it but still demonstrates similar performance as Robin which has only one IB port connected to CPU1. In short, the network performance is very inhomogeneous between the different computing elements of such MIC cluster. If one considers the native mode use case where a standard MPI application is deployed, the communication performance will likely be limited by the slowest path. So one should expect latencies five times larger and bandwidth 20 times lower than on the network of a "standard cluster" connecting only CPUs. The reasons for this network performance inhomogeneity are not yet completely understood. Setting the appropriate software stack and operating system configurations for that type of hardware seem to be complex and still an issue.
gysela application The gysela code is a non-linear 5D global gyrokinetic full-f code which performs flux-driven simulations of ion temperature gradient driven turbulence (ITG) in the electrostatic limit with adiabatic electrons. It solves the standard gyrokinetic equation for the full-f distribution function, i. e. no assumption on scale separation between equilibrium and perturbations is done. This 5D equation is self-consistently coupled to a 3D quasineutrality equation. The code also includes other features not described here (ion-ion collisions, several kind of heat sources). The code has the originality to be based on a semi-Lagrangian scheme  [SRBG99]  and it is parallelized using an hybrid OpenMP/MPI paradigm  [LCGS07, CLS09]. The coordinate systems we consider is as follows. The spatial coordinates consist in (r, θ) the polar coordinates in the poloidal plane (the origin for the radius r is the magnetic axis and θ is the angle), ϕ the angle in the toroidal direction.
Splitting scheme versus 4D advection, benefits expected The usual way to perform a single advection (transport mechanism caused by the flow) in the gysela code [GBB + 06, LCGS07] consists of a series of directional advections: ( v /2, φ/2, rθ, φ/2, v /2). Each directional advection is performed with the semi-Lagrangian scheme. This Strang-splitting converges in O(∆t 2 ). It decomposes one time step into four 1D advections and one central 2D advection (in the poloidal plane (r, θ)). Algorithm 1 shows a single 1D advection in the ϕ direction. Let us now consider an alternative method that avoids the Strang splitting. Algorithm 2 sketches the corre- ). The difference with the previous scheme is twofold. The feet of the characteristics has to be computed in four dimensions and the spline coefficients come from a tensor product of splines in four dimensions.
gysela mini-app In practice, we have used the OpenMP paradigm on a shared memory node throughout all the codes presented hereafter. In this study, the focus was made on the optimization/vectorization issue inside one single node (2 sockets, Sandy Bridge, 16 cores) and one Xeon Phi (60 cores). This enables to look at optimization and vectorization issues, and compare raw performances. The gysela mini-app is a 4D drift-kinetic code (µ = 0) in which we have removed a lot of diagnostics and physics submodels from the original gysela code. Furthermore, the mini-app contains the 4D advection capability which is not yet used in the production code. A small run with this mini-app on one single Sandy Bridge node takes 10 hours for the following domain size: N r = 128, N θ = 256, N ϕ = 32, N v = 64 with 4000 time steps.
This run fully used 16 cores thanks to OpenMP parallelization. The code was compiled with -xAVX compiler flag in order to have some automatic vectorization whenever it is possible. A profiling of the application using a single process with 16 threads has been performed on the Sandy Bridge partition with the Scalasca toolkit. As a result 98% of the run time is concentrated in one single routine and can be split into: • Spline construction 10% • Vlasov Solver (4D advection) 89% split into: -Computation of the feet of the characteristics 41% -Interpolations with 4D stencil using spline coefficients and feet 59% The spline construction represents 10% of the execution time which is mostly spent in memory copies. In principle, it shows some space for improvements as it has not been completely parallelized.
However, as it represents only 10%, we will focus on the parallel region which consumes 89% of the dominant subroutine and mainly on the 4D interpolation routine (51% of the total run time). We first tried a top-down approach to improve the performance of the 4D interpolation kernel. We extracted the most computation intensive part of the kernel in a single subroutine of less than 200 lines of code. Then, a lot of optimization techniques were investigated on this subroutine to get a vectorization of parts of the kernel and cache-friendly behaviour. We have speed up the initial version by a factor 2 on medium test cases on the Xeon Phi (domain size 128 × 128 × 64 × 32, 45s for one advection step with the initial version and 25s after the improvements have been done). Nevertheless, on a Sandy node with 16 cores and a similar test case, the code is eight times faster (3.
Choosing interpolation kernels The choice of interpolation method in a semi-Lagrangian scheme is crucial. It determines the numerical quality of the scheme and its computational cost. The tensor product is employed in this work to achieve multi-dimensional interpolations. Therefore, we need to fix first the method for the 1D interpolation. We will investigate the cubic spline scheme, which is known to be a well-balanced compromise between cost and quality for plasma simulations. This is the scheme used in the gysela code in production. Nevertheless, cubic splines lead to rather complex computational kernel. Thus, we have chosen to first look at Lagrange polynomials 3. The Lagrange interpolator is simpler and thus gives us the opportunity to easily investigate several optimization alternatives.
This simplifies the code and allows to isolate the code that performs the interpolation. This unrealistic assumption will be relaxed in the last section of this paper. The kernel in Figure  6  performs the same calculation as Figure  5  but using dedicated intrinsics  [INT] , that enables to use Xeon Phi SIMD instruction set. Figure  7  shows the interpolation with the 2D tensorial product of Lagrange polynomial (order 4). The number of FLOP grows from 7 in the inner loop of Figure  5  to 35 in Figure  7.
Memory-bound kernels For Lagrange 1D and 2D, the outer loops were parallelized through OpenMP directives and there is no dependency between loop iterations. The most inner loop is vectorized using SIMD features. Several versions of the code have been written using Intel intrinsics for Xeon as well as standard C or Fortran code plus Intel SIMD directives. The intrinsic versions require to write explicit low-level algorithms. This ends up with a non-portable code (specific to Phi architecture), and also some troubles to maintain the code. In the case of kernels instrumented with directives, it is quite difficult to find the correct location for directives, and also which directive to use. Indeed, in its report, the compiler can say that a loop or a block is vectorized, even if the quality of the vectorization will not give large benefits for example by masking all elements of the vector but one.
A look at the generated assembly code 6  helps a lot to control if the compiler is doing a good choice of SIMD instructions for vectorization. At the end of the day, both versions (directives and intrinsics) give approximately the same performance and we have checked that the most inner loop contains nearly the same floating-point SIMD instructions. For these memory-bound kernels, the prefetch mechanism is important. In practice, it means loading data in advance thanks to ad-hoc directives or intrinsics. We have observed a 20% reduction of execution time by good choice of prefetch parameters. To determine the best choice for the prefetch parameters, a script that scans an extensive set of possible combinations has been made. It is also interesting to tune the kernels in order to work on aligned data.
Compute-bound kernels Cache blocking (loop tiling) is crucial to save computation time. As caches are small on Phi, designing algorithms that are cache friendly is quite a hard task. Nevertheless, we manage to save execution time with proper loop tiling on the 3D kernels (50% gain). We have compared C and Fortran kernels that are similar line-by-line (with SIMD directives for the vectorization). The language has some influence on the performance (at the percent level), but on some kernels C is faster, on some others the Fortran version is faster. So no clear trend could be identified here. If the language does not seem to be an issue, the way algorithms are written triggers different kind of optimizations performed by compiler. These have a large impact on execution time. For example, splitting the body of a 10-40 lines loop into multiple loops of less than 10 lines of code each (if it is semantically correct) can lead to effective speedups.
One has to exhibit much more parallelism to get the performance out of the Phi device. It also means that for some domain sizes or parameter sets, it is difficult to get enough threads to get the best performance on the Phi. These rather simple compute-bound kernels are vectorized with difficulty by the Intel compiler despite the usage of the available SIMD directives. A small modification in the code can lead to bad vectorization at some location. When this kind of event happens, the code can slow down by a factor of 4 suddenly. It is extremely difficult to identify this kind of vectorization problem (compiler does not report it) and to set up a work around. As a practical method, we detect such weak link by commenting a group of code lines and looking at the impact on performance. Of course, the numerical results are wrong, but iterating this process on other groups of lines enables to localize the weak link and sometimes to find a cure.
Getting efficient kernels A partial conclusion of this study is that achieving performance even on a simple kernel is a challenge on Phi. Compared to the relative easiness to get a reasonably efficient code on Sandy, the programmer has a lot of constraints to fulfill on Phi. The developer has to interact finely with the compiler and profiling tools to improve the performance to an acceptable level. If it is possible to take care of some of the important points on a simplified kernel code, it can be very difficult to satisfy some constraints throughout a long program. The non exhaustive list of theses points could be: good vectorization, fine grain parallelism (threads), data alignment, cache blocking (adapted to small caches), prefetching, loop splitting. Applying these hints in a long code, in which the computation time is distributed over several routines, is not simple.
Porting mini-app on Xeon Phi In the previous section, we have studied computation kernels closer and closer to the ones used in gysela so as to identify guidelines and good practices to reach acceptable performance on the Xeon Phi. This section builds up on this experience to target a more realistic kernel that uses a 4D tensor product of splines instead of Lagrange polynomials and a variable displacement field for the feet of characteristics. In a first phase, we keep focusing on the interpolation part similarly to previous sections while in a second phase, we analyze the complete advection code. To do so, we first consider the complete advection code but we comment out the parts we want to ignore for the first phase. This approach enables us to easily extend the previous interpolation kernels and to integrate it without bothering with other issues coming from other parts of the application.
Advection kernel design The kernel we design in this section is not a synthetic kernel anymore but is instead a complete implementation of Algorithm 2 intended to be included in the gysela code. As we have observed in previous paragraphs, there is only a small impact on performance if one choose SIMD directives instead of intrinsics for vectorization. For maintenance issue and because gysela is a Fortran code, we made the choice to write the kernels in Fortran. Designing a complete kernel increases the amount of code that has to be optimized compared to the small (typically less than 150 lines of code) kernels of the previous sections. This is the first source of increased complexity. Another complexity comes from the use of a variable displacement field for the feet of the characteristics.
Also, this computation of the foot mixes floating-point (for the values) and integer (for the location) calculations that depend on each other. Furthermore, we did not found a way to vectorize both floating-point and integer operations efficiently. Then we focus on the vectorization of costly floating-point computations. This leeds to an innermost loop of the algorithm that is not completely vectorized with SIMD instructions since some integer operations remain. Taking all these aspects into account, we have applied strategies to reach good performance on the Xeon Phi. This results in the code presented in Figure  8  for the outer loops and Figure  9  for the kernel called from these. In order to come up with all these optimizations, we had to test several ways of writing the algorithms and closely look at the generated assembly code as it is always difficult to predict the impact of a change on the generated code.