{'fraction_non_alphanumeric': 0.025395076260362144, 'fraction_numerical': 0.0083398031505923, 'mean_word_length': 4.739542906425183, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Handling Nested Parallelism and Extreme Load Imbalance in an Orbital Analysis Code

Nested parallelism exists in scientific codes that are searching multi-dimensional spaces. However, implementations of nested parallelism often have overhead and load balance issues. The Orbital Analysis code we present exhibits a sparse search space, significant load imbalances, and stopping when the first solution is reached. All these aspects of the algorithm exacerbate the problem of using nested parallelism effectively. In this paper, we present an inspector/executor strategy for chunking such computations into parallel wavefronts. The presented shared memory parallelization is no longer nested and exhibits significantly less load imbalance. We evaluate this approach on an Orbital analysis code, and we improve the execution time from the original implementation by an order of magnitude. As part of a Graduate Computer Science course in Parallel Programming models, we show how the approach can be implemented in parallel Perl, Python, Chapel, Pthreads, and OpenMP. Future work includes investigating how to automate and generalize the parallelization approach.


I. INTRODUCTION
Nested loop parallelism is natural to express in programming models such as OpenMP, but difficult to efficiently realize when sparse computation spaces with significant load imbalances and early termination criteria are involved. In this paper we present an approach to parallelizing such computations on shared memory machines.
The orbital analysis code we worked with in this case study consists of a six-deep nested loop structure including the outermost loop over particles (see Fig. 1). Each of the six nested loops can be executed in parallel, thus the computation experiences significant nested parallelism. One problem is that the iteration space is sparse: particles are checked for consistency and the (p,q) ratio is checked to avoid equivalent repeats. Specifically, there is a condition checked at line 10 before the computation for a particular (p,q) iteration executes. Most of the points in the parameter space fail this check. Another problem is early termination: this code will return upon finding parameters that satisfy the check in the innermost loop at line 16.
longer than the a short running particle. The execution time for a given particle cannot be predicted without running it fully. In the worse case, in which a single processor is given all of the long running particles, no gains are achieved from parallelization. A second alternative would be to use nested parallelism. However, the early termination check in the innermost loop makes nested parallelism impractical. If all size loops where specified as parallel, all iterations would execute even when early termination is possible, which is frequent. Additionally, there would need to be a reduction computation that determines the earliest values of (p,q,m,n,r) where the condition was satisfied, because that is the correct result for the program.
A third alternative is to use task-based parallelism. This would work by spawning off tasks for each call to check-Libration(). While implementing the task-based model two problems arise. First, the non-determinism introduced by the task parallelism loses the ordering information guaranteed by the serial code (i.e., early termination strikes again). This prevents us from knowing if the first value returned is the optimal solution. Second, there is too much task overhead. Each call to checkLibration is lightweight, but the amount of calls made is high. In an extreme case one particle spawned over 140,000 tasks.
To handle the load imbalance, sparse iteration space, and early termination issues, we developed a finer-grained parallelism internal to each particle. The approach consists of building parallel wavefronts of tuples in the search space at a particular nesting depth. Figure 2 shows the pseudocode for the algorithm implemented. In the new algorithm, for each (p,q) ratio that passes the check on Line 10, a subset of the search space is collected. The CheckLibration() function can then be called on that subset in parallel. The final loop at Line 23 will check if any of the tuples in the just executed wavefront passed the libration check and thus the computation should terminate early.
The original orbital analysis code was implemented in Perl and would have taken more than a month to analyze the monthly observations from a new telescope. Parallelizing the computation and porting the Perl analysis script to more efficient programming models makes the execution time practical (less than a week). As part of a Graduate Parallel Programming Models class, we evaluated the process of implementing the parallel analysis script in Perl, Python, Chapel, Pthreads, and OpenMP. We compare code snippets from the various programming models to exhibit different approaches for implementing the presented parallelization strategy Significant speed-up was achieved over the original Perl version. Using a 12-Core machine over 3x speed-up was achieved in every implementation versus their own serial version. The PThreads version was the overall fastest, bringing the worst case per particle down from 541 second to 5.2 seconds, a 4x speed-up versus its own serial version, and 103x speed-up versus the original Perl baseline. This improvement allows for the total analyses of a month's worth of data (approximately 40,000 particles) to be performed in several hours.
The Astronomy community has been developing more software in Python. The serial Python version performs comparably to the baseline Perl. With the described algorithm the total execution time for the most costly particle was still brought down to 78.8 seconds, or about a 6.5x speedup over the original Perl code. This brings the computation down into the realm of acceptability.
In this paper, we make the following contributions: • Provide a description of the application: details about the problems faced by the scientist that necessitate code performance improvement. • Parallelizing the code: Why simple solutions do not work, discussion of the inherent workload issues, and the final solution. • Description of how the parallelization approach can be mapped to the programming constructs in various pro-

II. ORBITAL ANALYSIS APPLICATION DETAILS
The Kuiper belt is a population of small bodies in the outer solar system. When new objects are discovered in the Kuiper belt, a common first step is to simulate the orbits of the observed objects forward in time and then analyze the results to distinguish between various types of orbital evolution. The problem is that aspects of this analysis experience significant load imbalance. There is additional interpretation overhead due to the code originally being written in Perl. Both issues can lead to analysis of a single particle taking in excess of 10 minutes. NASA's Large Synoptic Survey Telescope (LSST) will begin operations in the early to mid 2020s and is expected to discover and track about 40K Kuiper belt objects over a ten year survey (compared to the approximately 1000 currently tracked objects). Thus the need to do this analysis more efficiently is critical. In this section, we describe the analysis and its current performance bottleneck

A. Classifying Objects by Their Neptune Resonance
The objects in the Kuiper belt represent a record of the dynamical history of the solar system's giant planets. The distribution of Kuiper belt objects (KBOs) in orbital resonance with Neptune is of particular interest because it can serve as an observational test for theoretical models of the outer solar system's dynamical history [1]- [3]. An object is in orbital resonance with Neptune when there is an integer ratio of the number of times it orbits the sun and the the number of times Neptune orbits the sun. When KBOs are observed, their orbits must be analyzed to determine whether they are resonant. Such dynamical classification is important in prioritizing objects for continued scientific study. For example, some specific hypotheses for the history of the solar system can be tested by observationally determining the chemical compositions of resonant KBOs; such observations are costly because they can only be done on very large telescopes, and dynamical classification is necessary for efficiently planning them. The classification process, which is detailed in [4], entails the following basic steps: 1) The position of the object in the sky at a variety of epochs is determined from observations, 2) these positions are used to fit an orbit (a combination of position and velocity) for the observed object, 3) the orbit is numerically simulated forward in time under the gravitational influences of the solar system's planets, 4) and the simulated orbital history is analyzed to determine if the object is in resonance with Neptune.
If the object is in resonance with Neptune, then the analysis program determines the specific resonance the object is in and the amplitude of libration for the associated resonance angle. The resonance is labeled p:q according to the period ratio between the object and Neptune; a 3:2 (p = 3, q = 2) resonance is one where the Kuiper belt object's orbital period is 1.5 times Neptune's orbital period (see Figure 3). When an object's evolution is controlled by the resonance, an angle will librate around a central value, whereas objects not in resonance will have a resonance angle that freely circulates between 0 and 360 • as shown in the right side of Figure 3. For a given p and q, there are many possible angles determined by the integers m, n, r, and s that could librate; this is the origin of the nested loops shown in Figure 1.

B. Performance Bottleneck
The current performance bottleneck for this process is the last step where the simulated orbit is analyzed for libration of any relevant resonance angles. A large number of p : q resonances must be checked for each object: all pairings of the integers between 1 and pmax (on the order of 30 to 70), where pmax is a runtime parameter, which specifies the granularity of the angles checked. Objects that are near the correct period ratio for a large set of p, q values take a long time to analyze.
To our knowledge, there is no standard, open source code available to performs this analysis. Researchers generally report the results of such analyses in papers, but do not make their codes available or report specific details about the analysis methodology. We investigate the workload of a Perl script that planetary scientist Dr. Kat Volk wrote to perform the analysis.
A Serial Perl implementation of this algorithm takes approximately 10 minutes on a Xeon Westmere-EP Dual 8-core Processor node to fully analyze an object that is close to many potential resonance ratios (some of these are eventually categorized as non-resonant). It is not possible to determine up front whether a particle will pass the initial resonance checks in seconds or require most of the 10 minutes, thus causing one level of load imbalance in the workload. Typically, we expect ∼25% of all observed KBOs to be actually non-resonant and thus require the full analysis time [5].

C. Future Performance Demands
In the last 20 years, approximately 1000 objects have been observed and classified; because they were typically discovered in groups of 10 to 50 [6], an analysis time of up to 10 minutes per particle has not been an issue. However, when the Large Synoptic Survey Telescope (LSST) comes online in the mid to late 2020s, the expectation is that forty thousand KBOs will be discovered within the first few years of the survey, and they will be continuously tracked over the survey's ten year lifespan [7]. Each month the LSST will scan the full southern There is a large divide seen between particles that can be almost immediately rejected, and those that require a significant search time to confirm or reject sky multiple times and will provide new measurements of the position of these objects. Each set of new observations produces a more accurate orbit for the object, which requires re-analyzing the orbit for resonant behavior (until a sufficiently accurate orbit is determined). With 25% of the 40,000 objects requiring 10 minutes of wall-clock time each month, the analysis time for one month's worth of data is approximately 2 months. This represents a significant performance bottleneck for researchers. We also note that code that enables more efficient identification of resonant behavior would be useful for analyzing test particles in numerical simulations of the outer solar system; such simulations are used to produce the theoretical predictions the observations are meant to test [8].

III. PARALLELIZATION PROBLEMS DUE TO IRREGULAR WORKLOAD
The goal of the orbital analysis code is to classify particles as either resonant or not resonant with Neptune. Since there is no communication between the computations for individual particles, parallelization over the particles would be the simplest. Figure 4 shows the per particle contribution to execution time on our test set of particles. The unpredictable analysis times of the particles can easily lead to a small set of cores receiving all of the most computationally expensive particles, resulting in little to no improvement in the execution time. Therefore in the next section, we present a more effective strategy for parallelizing the computation. This section details our experimental methodology and the load imbalance.

A. Methodology
We ran experiments on the HPC system of the University of Arizona. The machine is an SGI Altix UV 1000 consisting of 170 Nodes, a given node consists of Xeon Westmere-EP Dual 8-core Processors running at 2.66 GHz. The machine was running Red Hat 6.0 Linux.
The original Perl code and the multithreaded version were run using Perl 5.10.1 and using Parallel::Loops from CPAN [9]. The OpenMP and Pthreads versions of the code were both compiled with g++ (GCC) 4.4.4 20100726 (Red Hat 4.4.4-13). Chapel was compiled using the Cray Chapel Compiler version 1.12.0. Python was run using Python 2.7.9, and the Multiprocessing module from the Python standard library.
Both of the datasets we used consist of sets of particle simulations. Each particle has its run discretized into a number of time steps, which records the relevant angular information at the given point in the simulation. The more time steps, the finer grained the analysis and the longer processing takes. The first, used for the bucketing stage, consists of 500 particles whose orbits were divided into 9629 time steps. This was used only for one step due to the original prohibitively long testing time of the 500 particles. The second, consists of a set of 82 particles, each of 9629 time steps, this was used for all cross-language execution time comparisons. The third consists of 100 particles consisting of 50,000 times steps, analysis of which is currently too time consuming to perform.

B. The Algorithm
The resonance check for a given particle consists of the following steps as seen in Figure 1: • If the isConsistent() call on Line 3 returns false the particle can be immediately rejected as non-librating. • Checks all unique ratios of p : q on Lines 8-10. • Calculates each possible angular variation of the particle based on the (p, q, m, n, r, s) values created in Lines 8-14 and calls checkLibration() to see if the particle librates with those angles. As mentioned, the analysis of each particle is fully independent from that of any other particle, and so the particle level was an obvious place for parallelization. In practice, this consists of parallelizing the for loop around the call to CheckRes() as seen in the pseudocode in Figure 1. Each particle is checked using the CheckRes() function, which does all the necessary work to determine if a particle is resonant or not.
Unfortunately, this did not improve performance much in the best case, and increased execution times in the worst case. Figure 5 shows the behavior of the this method in each implementation as the number of threads increases. The behavior is irregular, with several implementations showing fluctuations in execution times as the number of cores changed.

C. Workload Characterization
The particles can be placed into three categories: 1) the range of the particle's semi-major axis is too great and therefore can be rejected outright, 2) the particle is in resonance, or 3) the particle cannot be outright rejected, but is not in resonance. Each of these categories in turn tends to affect execution time in different ways. The first case is the fastest, contributing a negligible amount of time to overall execution time. The second case is highly variable, the search space is checked until a resonant angle is found. Though usually near the beginning of the search, it can theoretically require analyzing all possible ratios to find a single resonant angle. The third is the worst case, a particle that does not resonate but cannot be rejected outright. The particle at this point requires analyzing the entire search space before it can be rejected conclusively.
As the search space is quite large, these types tend to dominate run time. All three of these particle types are spread throughout the data and one cannot know which category it belongs to without actually running the code. The analysis times for individual particles were examined further. The original Perl code was used to classify a set of 500 particles, and analysis time for each particle was recorded. Figure 4 shows the execution time distribution for analyzing the resonance of these particles. The results are stark, showing a bimodal distribution of particles either finishing nearly instantly, or running for several minutes at a time. This high variability in the work required to classify a particle explains the poor performance of the naive parallelization. The distribution of the analysis times cannot be known ahead of time, and so any gain from this type of parallelization will happen randomly, depending on how the compiler chooses to distribute the particles.

IV. PARTICLE-INTERNAL PARALLELIZATION
Parallelizing the computation that occurs within each particle will help avoid the significant load imbalance issues that occur between particles. In this section, we present an approach to find wavefronts of parallel computation within the analysis for each particle.

A. Issues with Straight-Forward Parallelization
The most obvious means to parallelizing the analysis within each particle would be to parallelize work at the outermost p and q loops. This has two issues due to requirements of the algorithm.
One issue is early termination. In the serial code, as soon as a result is found the code returns. By doing this the code can often avoid much computation. As the code searches for a solution, it can often find a valid result on the first several checks. In the worst case, there can be over 270,000 such checks when pmax = 30. If we cannot return early, we are forcing all particles to be searched exhaustively. This increase in work leads to worse performance, even with parallelism added.
The other issue is one of ordering. As the algorithm progresses it checks possible values in a specific order. The earliest lexicographic iteration in this space is the least likely to be a false positive for resonance. This means that models that allow communication to end early, but cannot guarantee ordering of results are also invalid for our purposes.

B. Wavefront Parallelization
The load imbalance and deeply nested structure prevent consistent performance gains from being achieved across the different categories of particles. To overcome this, the inner most calls were factored out, and the nested structure was flattened into a single list that could be iterated over easily. The pseudocode for this is given in Figure 2. The internal parallelism algorithm can be broadly broken into the following steps: • An array of tuples is created, (p, q, m, n, r, s), instead of calling CheckLibration() directly, • CheckLibration() is mapped in parallel over the array of tuples, • a new array is returned of the results of each CheckLibration() call, and • these results are scanned serially to find the first occurrence of resonance, or lack thereof. • If a result is found, the tuple associated with that result is returned. • Otherwise, we continue searching until all possibilities are exhausted. This requires extra overhead in storing the tuples, but it fulfills all of the restrictions discussed. Returning the results as an array means that it does not matter how the individual processors run, the order of the results are the same as if it had run serially. We can then perform a serial scan of the results, finding the one that occurred first lexicographically. This allows us to reconstruct which result is the optimum in the case of false positives being returned.
As previously discussed, as soon as a result is found the code can stop. Often libration is found before searching the entire space. Building up the entire search space heavily penalizes those that could have terminated early.
To solve this a final addition to the algorithm was added. Instead of building up the whole search space at once, we build up a subset (or wavefront), where some common prefix of the (p, q, m, n, r, s) tuple is kept constant. We than check all tuples in this subset in parallel. If a result is found we can terminate early, if no result is found a new subset can be generated in check. This partial generation technique allows  for us to parallelize over the search space, but without forcing searching the entire space.

C. Placing the Parallelism
Determining the correct prefix to keep for each check is difficult to determine a priori. If you check on each iteration of p than each check is still checking a very large portion of the search space, and so the benefits of being able to return early are lessened. If you stop and check on each iteration of n loops you face the opposite problem, you end up making so many parallel calls the overhead costs hurts the results.
An empirical approach was used to best identify where to implement the parallel subset search. Figure 6 shows the execution times over the 82 particle set. This test was done in Python. On the ends we see the expected bimodal behavior, the shortest possible common prefix does poorly, as does the longest. As we slide towards the middle we see significant improvement in overall performance. The (p, q) prefix outperforms the others, and was chosen as the placement for parallelism for all further experimentation.

V. IMPLEMENTING THE PARALLELIZATION IN VARIOUS
PROGRAMMING MODELS The original code written by the planetary scientist was written in Perl. In the context of a graduate parallel programming models course held in the Spring of 2016, we implement the designed parallelization of the libration analysis code in the following parallel programming models: • Perl Thread Library, • C++ and OpenMP, • C++ and Pthreads, • Python Multiprocessing Library, and • Chapel. We experimented with a number of parallel programming models with the primary evaluation metric being the performance. However, the maintainability and probable evolution of 1 for(int p=0;p<pmax;p++) 2 for(int q=9;q-p<pmax && q>=0;q--) for(int n = p-q-m ; n >= 0 ; n--){ 7 for(int r = p-q-m-n ; r >= 0 ; r--){ 8 int s = p-q-m-n-r; 9 angles

A. Parallel Perl
Due to the algorithm being originally implemented in Perl work was done to parallelize in Perl. This was problematic though, and Perl was ill-suited to this task. The multithreading library of Perl is officially recommended against by the language designers. Each new thread spawns an entire new instance of the Perl interpreter, a heavyweight action. The results of this was that even with the full 12 cores of the system the execution time was over 10 times slower than serial implementation. The performance and source lines of code count were not recorded for Parallel Perl.

B. C++
C and C++ are popular languages for writing high performance code. For this comparison two of the most common parallel models were chosen: • OpenMP: A library consisting of a set of compiler directives that the user uses to indicate where parallelism will be inserted. • Pthreads: A library for manual thread managements, which enables the user to specify all details about how to divide up information and send it out to different threads for processing. First, the Perl code was translated to equivalent C++ code. From this new base code, the internal parallelization algorithm was implemented using each of the two models.

C. OpenMP
OpenMP is a parallelization library and compiler for C, C++, and Fortran [10]. It represents parallelism as a set of compiler directives known as pragma. The user must identify places where parallelism is to be inserted, and then inserts a pragma that the compiler uses to parallelize the program.
The compiler then does the necessary conversions to divide the iterations of the loop and assign it to individual threads, which will then be run in parallel. In many real world cases, more information is required for the compiler to perform its work correctly, and so there are a variety of possible arguments and different pragma that can be indicated by the user. In all cases though, OpenMP is dependent on the user identifying and specifying the type of parallelism to be used. Figure 7 shows the OpenMP implementation of our algorithm. For our code we first perform the build up of a list of tuples, as shown in Figure 2. Because data is being written to a shared memory variable, result we must be cautious of the shared variables since it may lead to correctness issues if not handled. This takes the form of specifying which variables are private to each thread, and thus can be copied, and which are shared, in which case the compiler must handle writing to a shared memory space. Finally, since OpenMP does not allow branching in or out of a structured block, we place the results in a global list of results. This list is scanned to see if resonance was found, and then to identify which values led to resonance.

D. Pthreads
Pthreads is a low level C/C++ library for creating shared memory multi-threaded programs. Programmers define and create the number of execution threads, partition the computation and the data, and explicitly define which thread is going to do what part of the computation. Figure 8 shows the internal parallelization using Pthreads. Just like other models, we create a subset of computation space that we want to check. Next we need to create a number of threads that are going to process part of this subset. In Pthreads model each thread's execution starts from a function that we pass in to pthread_create() constructor.

E. Python
Python has become one of the most popular languages for planetary scientists [11]. Due to the usage of a Global Interpreter Lock (GIL) in the reference interpreter for the language parallel options are limited. The GIL prevents multiple threads from executing Python bytecode at the same time. The Multiprocessing library in Python sidesteps the GIL by using subprocesses rather than threads. This is more heavyweight than threading in other languages, but is a necessity forced upon the users by the GIL. Figure 9 shows the snippet for the particle-internal parallelization in Python. In the Python Multiprocessing library, the standard form of data parallelism is represented as a map function called over a collection. In this case, the CheckLibration() function will be mapped over the subset of tuples the inspector has collected. The results are guaranteed to be ordered the same as the input tuples' orders, allowing us to guarantee our search for the optimum.
1 // ti is the array of thread ids 2 pthread_t * ti = new pthread_t[threads]; 3 4 for(p=1;p<pmax;p++) 5 for(int q=p;p-1<pmax && q>0;q--) 6 7 // inspecting for tuples of work 8 if(checkRatio(p,q,ratios)) 9 for(int m=p-q;m>=0;m--) 10 for(int n=p-q-m;n>= 0;n--) 11 for(int r=p-q-m-n;r>=0;r--) 12  There is a mismatch between the angles that change from each call, and the particle that remains fixed. This led to us having to create a partial function to combine the checkLibration() function and the particle data. This creates a new function which has the same execution as checkLibration, but with the value for particle held fixed across all calls.

F. Chapel
Chapel is the newest programming language of the three, a model designed from the ground-up to support parallelism [12]. Created by the super-computer manufacturer Cray, Chapel represents an attempt to integrate parallelism from the ground up into a language. Parallel features in this language are syntactic elements of the language, with keywords indicating when parallelism is to be used.
The language allows us to match closely to the pseudocode algorithm originally show. We build a set, we have a parfor over that set, and then we check the results for a positive result. The difference between the serial version, and the parallel version is as simple as replacing for with forall. In our case, because we have a collection of tuples, but a fixed particle for each set, we are forced to use a partial function fit into the map model.

G. Source Lines of Code
We also compare the total Source Lines of Code (SLOC) for the serial and parallel versions of each language.
Python requires the least amount of code, but required certain constructs, such as partial functions, which may be unfamiliar to users from a imperative or object oriented background. This is a case where the difference in SLOC can hide the actual complexity of the implementation from readers.
Though Chapel is more verbose than Python, it still has fewer lines of code than the serial Perl. The transition to the parallel algorithm was also simple, with the final parallelization change requiring only a keyword change from for to forall. Both the Python and Chapel code also had smaller differences between the serial and parallel versions than with the C++ versions.
For the OpenMP version the majority of the difference came from the refactoring necessary to fit the final parallel algorithm chosen for this task. The actual insertion the code related to OpenMP was much less: 3 lines overall including the pragma, the import, and a required initialization statement to control number of cores used. PThreads had the largest difference by a wide margin. Beyond the changes required for the parallel algorithm, significant changes were required to fit that within the framework of PThreads. Explicitly managing all of the threading logic leads to significant tangling with the underlying logic of the chosen algorithm.
All 4 versions required some amount of refactoring, both to implement the new algorithm and to add the parallel logic. The Chapel and OpenMP versions were the least invasive in the insertion of parallel logic. And only these two allow for switching between parallel and serial code without requiring any changes to the final code. In this regard, they have the much less tangling than the other two implementations.

VI. PERFORMANCE COMPARISON
We tested the implementations on two different datasets. The first was the same test set used for the Perl baseline data. This consists of 82 particles, represented as location and orientation at 9629 different time steps in the simulation. The long run-time at this granularity prevented analysis of the particles at a finer level of time. Each version showed at least 6x speed-up over the baseline code, with the OpenMP version achieving 103x speed-up, and a final runtime of 5.2 seconds versus the baseline run of 541 seconds.
Due to the success in improving performance on the dataset we did further testing on a second dataset of particles. These modeled each particle at 50,000 time steps and represented the granularity that the scientists would prefer to analyze. With the Execution times for parallelization of the internal particle check. In this model each particle is run one-by-one, but search for each particle is parallelized. The Pthreads and OpenMP versions overlap almost exactly.
baseline Perl version a single particle could take more than 45 minutes to analyze. Performing analysis at this level was even less feasible for the incoming amount of data. For these long running particles we were able to get significant improvement as well. In the worst case, the Python version had 8.4x speedup. In the best case, the OpenMP version once again performed best and gave us a 105x speed-up, bringing the worst particle down from 45 minutes to 25 seconds.
Finally, we give a Source Lines of Code (SLOC) comparison across language. This is done to give an approximation of the effort required to implement each version of the code. With this metric, both the Python and Chapel version outperformed the baseline Perl code.

A. Check-Internal Parallelization
The final version chosen for experimentation uses the Internal Parallelization model over the Naive model, because that version was able to handle the irregular workload. In this the checks over the search space for a given particle was parallelized, this has us handling each particle one by one, but handling each of them much quicker. These implementations were tested on the set of 82 short-length particles. Figure 12 shows the raw execution time changes for each implementation. Figure 13 shows the Speed-Up curves for each of the implementations. With this we are able to handle the unpredictable workload much better. The new implementation is entirely lacking in the sharp changes in performance found in the particle level parallelism. Each programming model shows a similar curved improvement as the number of cores increased. This shows that the algorithm chosen improves execution time of the longest running particles without too heavily penalizing those that can quickly return an  answer. Thus, the flattened tuple structure chosen overcomes the workload imbalance issues inherent to the dataset.
The Python version performs well for an interpreted language. As the number of cores increases we see that the execution time approaches the speed of the compiled Chapel language. This alone represents significant speed-up, enough to bring the total execution time into the realm of feasibility.
The Chapel version performs well with a low number of cores, but does not scale as well as other versions. Only barely outperforming the Python version, surprising considering that Chapel is a compiled language. Though it did not perform as well as one would hope, the language is still young and the team has specifically targeted performance as their major focus for improvement in the coming years.
Both of the C++ versions match each others performance. OpenMP and Pthreads each show the same curve, and in general seem to be functionally identical to one another. Overall these two versions perform the best, with the Pthreads version barely out-performing the OpenMP code. The final implementation brings then execution time down to 5.2 seconds from the 541 second baseline serial Perl execution, giving a final speed up of 103x.

B. Large Time Scale Particles
Testing was extended to the second set of longer time scale particles. These particles represent a level of detail that the researchers would prefer for their analysis, but were unable because of the prohibitive execution times. With the original Perl implementation each particle was taking more than 10 minutes to analyze on the low end. On the high end a single particle could took up to 45 minutes. This quickly led to an unworkable situation. Figure 14 shows a comparison of the execution times of the 5 worst particles from the long-time scale data set. The three versions shown perform well, bringing the execution times down much closer to the realm of feasibility. The OpenMP version is the fastest. but even in the worst case there is 8.2x speed-up over the 45 minute long per particle time for the baseline serial Perl implementation. Even with the slowest version we can show that the run-time can be brought down into a more feasible realm. This problem of handling a larger dataset is interesting in that it is a common problem seen in the Scientific Computing community. So often are they limited by the speed at which they can do computation that they often work with simplified workloads to bring execution time down into the realm of feasibility. And so speed-up for them often represents a means to continue on to more complex computation, or finer grained analyses of their data. In this case, execution time was able to be improved so significantly that it brings nor only the current time scale level into a more feasible time frame, but to allow for better analysis to be done in the future.

VII. RELATED WORK
There has been much work on improving the performance of nested parallelism. Blikberg and Sørevik [13] argued for flattening nested parallelism into a single level of parallelism. They then have an approach for load balancing when there is some model of how much work each task is doing. In the orbital analysis application the amount of work that will be needed per particle cannot be modeled ahead of time. Thibault et al. [14] are dealing with the problem of expressing data affinity to the underlying runtime system. The orbital analysis program has severe load imbalance issues, not data locality issues. Dimakopoulos et al. [15] created a nested parallelism benchmark and then showed that many nested parallelism implementations of OpenMP have a lot of overhead. A significant amount of research investigates the advantages and disadvantages of various programming languages in the context of scientific computing. Caie et al. [16] describe various libraries and capabilities in Python such as NumPY and the ease of calling Fortran and C and how those impact the performance of stencil computations that occur when solving partial differential equation solvers. Many others have compared various parallel programming languages in terms of their performance and programmability with various benchmarks and applications [17]- [22]. This study focuses on characterizing the workload and performance alternatives for implementing a specific analysis needed for the Large Synoptic Survey Telescope (LSST) project.
Shen et al. [23] describe a parallel algorithm implemented in MPI for the analysis of objects that are close to the earth. The algorithms in question are different than those we study in this paper.

VIII. CONCLUSION
In this paper, we discuss problems that arise for the analysis of orbital particle data. We provide analysis of the deeply nested structure, and the sparseness of the search space, which lead to significant load imbalance in naive attempts at parallelization. We describe a solution to these issues, presenting an algorithm for flattening the nested structure to allow better parallelization.
We then discuss the issues that arise during parallelization in Perl, Python, OpenMP, PThreads, and Chapel. We include code snippets which show the implementation details for our algorithm, and the language specific differences that arise because of this.
In the end, we show significant speed-up in all languages, excepting Perl. Special attention was given to improvements in Python, a language which has seen widespread adoption in scientific computing. We achieved 9.1x speed-up in Python, a language whose serial version performed comparably to the baseline. In the best case we show the Pthreads version achieved a 103x speed-up over the original performance. These improvements allow for analysis of finer-grained data that had previously been considered infeasible due to the execution time bottleneck.