Analysis of close encounters with Ganymede and Callisto using a genetic n-body algorithm

In this work we describe a genetic algorithm which is used in order to study orbits of minor bodies in the frames of close encounters. We find that the algorithm in combination with standard orbital numerical integrators can be used as a good proxy for finding typical orbits of minor bodies in close encounters with planets and even their moons, saving a lot of computational time compared to long-term orbital numerical integrations. Here, we study close encounters of Centaurs with Callisto and Ganymede in particular. We also perform n-body numerical simulations for comparison. We find typical impact velocities to be between $v_{rel}=20[v_{esc}]$ and $v_{rel}=30[v_{esc}]$ for Ganymede and between $v_{rel}=25[v_{esc}]$ and $v_{rel}=35[v_{esc}]$ for Callisto.


Introduction
Jupiter's large icy moons such as Ganymede and Callisto show countless impact craters across their surface. Studying these craters gives deep insights into the impactors as well as the moons themselves. This is the first approach in the frame of future works of studying collisions with the outermost moon Callisto. We are especially interested in the Valhalla crater system, Callisto's biggest crater. This impact structure measures several hundred of kilometers in diameter and shows some extraordinary features such as an extensive ring system in the outskirts of the crater (Greeley et al., 2000;Stewart and Allen, 2002). It has been shown that studying the formation of the Valhalla crater reveals new insights regarding Callisto's subsurface composition (Winter et al., 2017).
In this first work we focus on the use of a genetic algorithm 1 (which is described in Section 2 and in the Appendix for further details) as a proxy tool to find preliminary orbits of possible close encounters for bodies in the Solar System. In particular, we selected the Centaurs' population. Currently, 423 Centaurs are known 2 . It is estimated 1 herein 'GA' 2 Data taken from JPL Small-Body Database Search Engine https://ssd.jpl.nasa.gov/sbdb_query.cgi 1 that the number of Centaurs with a diameter larger than 1 km lies between n ∼ 10 7 (Volk and Malhotra, 2008) and about n ∼ 8 · 10 9 (Di Sisto and Brunini, 2007; Fernández and Sosa, 2015; Napier et al., 2015). These bodies mainly origin from the Trans-Neptunian Objects and have very chaotic orbits. With semi-major axes between a = 5.5AU and a = 30.1AU, they lie between the giant plantes, from which they are frequently ejected out of the Solar System via close encounters or even impact on one of the planets or their moons. Their lifetime is of the order of 10 Myr (see e.g., Horner et al., 2004a,b;Bailey and Malhotra, 2009;. A GA in combination with standard orbital numerical integrators can be used as a good proxy for typical orbits of minor bodies in close encounters or impacts with planets and even their moons, saving a lot of computational time compared to full orbital numerical integrations (Section 3.1). We also perform n-body numerical simulations (Section 3.2). We extrapolated the typical orbits in close encounters with their osculating elements and velocities (Section 3.1 and 3.2). This kind of orbital analysis can also be very useful for studying the origin of impactors in the last million years, i.e. the Bosumtwi impactor (Galiazzo et al., 2013). We particularly design the genetic algorithm to investigate close encounters with Callisto -and for comparison -Ganymede. Close encounters are the first natural approach to actual impact scenarios which we will study in following works.

Methods
Measuring close encounters or even collisions between minor and major bodies in the context of n-body simulations is computationally demanding. Typically one has to constrain the parameter space of the minor bodies to selected regions within the Solar System (e.g. Kuiper belt objects or specific families of objects). We use a genetic algorithm to find asteroids of the Centaur type family which are likely to have a close encounter with the Jovian moons Ganymede and Callisto. Centaurs mainly origin from TNOs and in particular from the Kuiper Belt . The method is used to encounter the problem of a large parameter space of initial orbital elements (see Table 1). The genetic algorithm boosts the performance compared to classical searching grids by some orders of magnitudes within the given computation time. This allows for measuring a reasonable amount of datapoints to do a statistical analysis of typical close encounter families which marks the first step towards studying actual impact scenarios.

Genetic Algorithms
A genetic algorithm (Turing, 1950) is an iterative searching algorithm to find solutions for highly complex problems which can have large parameter spaces. GAs origin from the field of biology, therefore we may also use the corresponding terms. Further information is given in the Appendix. GAs are efficient optimization methods for highly complicated functions. The population can overcome local minima quite easily either by the means of crossover or mutation. Moreover, GAs tend to find all possible families of solutions, even if they are unrealistic or uncommon. However, GAs also have nega-tive aspects. The implementation can be quite tricky and there is no general rule how to implement them efficiently because the functions (for fitness, crossover and mutation) and parameters (number of generations, number of individuals) have to be adapted to the problem as well as to the hardware. If the minimum fitness limit which is needed to solve the problem is not met (due to poor convergence or no learning process), the GA will not find any solutions.

Genetic N-Body Algorithm
genbody is a GA which is being developed to find close approach-and collisional orbits of particles in the context of n-body simulations. Each generation of the GA corresponds to a distinct n-body simulation. During the simulation, the fitness of the population is measured. Afterwards, a new population is created by the means of crossover and mutation. A pseudo code of the genetic n-body algorithm is given in Algorithm 1.
Algorithm 1: genetic n-body algorithm initialize population for g in generations do while t < t end do take n-body step (t → t + h) measure fitness crossover mutation We use the code to find objects which are likely to collide with either Ganymede or Callisto within a certain time interval. The following list gives an overview how the GA and the n-body code are associated with each other. In analogy to section 2.1, we link the terms between the genetic algorithm and the actual problem: • population: Centaur type asteroids with random initial conditions • DNA: initial Keplerian orbital elements (a, e, i, ω, Ω, M ) • generation: numerical simulation of the system via an n-body method • evolution: iterative process of consecutive simulations • fitness: score of each test particle at each simulation given by the fitness function • crossover: combination of initial orbital elements of parent particles given by the crossover function • child: a test particle with a new set of initial orbital elements • mutation: small, random variation of initial orbital elements The so-called fitness function is the function to be optimized. We use a fitness function of f = 1/d rel 2 with d rel being the minimum relative distance between a particle and the corresponding moon during each generation. The squared distance was found to be useful if the fitness is used as a probability distribution to sample quite fit parents for crossover, ensuring a high fraction of fit parents and thus increasing the overall performance of the GA. Note that one can use other fitness functions to study completely different problems. One only has to find a quantity to score objects which show a specific behaviour to obtain a population which develops that behaviour. The fitness of each body is measured throughout the simulation. If a close approach with the moon occurs, a measurement is taken and the algorithm starts a completely new evolution process. If no close approach happens during the simulation, unfit objects are replaced by new children which are created by objects with a high fitness score. The initial orbital elements y ∈ {a, e, i, ω, Ω, M } of the new children are exact 1:1 copies of their corresponding parent. Note that one could use other crossover functions which include two or more parents to create a new child. Next, a small mutation is applied to the whole population in order to avoid to get stuck in local maxima in the learning curve (generation vs. mean fitness). Since our n-body system turns out to be highly sensitive to the mutation rate χ, it is crucial to find suitable values for each generation. A too small χ leads to a solution which gets stuck in a local maximum while a too large χ on the other hand leads to destruction of genetic information. This corresponds to random guessing without any learning process of subsequent generations of populations. Since the trajectories of small bodies in the Solar System are chaotic, the mutation has to be low enough to allow for similar trajectories of subsequent populations (except for the newborn child) within the simulation time. We therefore use an adaptive mutation rate χ which depends on the overall fitness evolution. A scaling factor ζ is introduced to control the learning curve by setting the amplitude of the mutation rate χ: If the mean fitness between two successive generations varies too much, ζ is decreased to prevent a loss of DNA information. If the mean fitness between two successive generations varies too little (by less than 1% in this work) or if the inidividuals of the population are getting too similar to each other, ζ is increased to ensure a healthy population. The mutation rate itself is given by the standard deviation σ of each initial orbital element throughout the population (e.g. χ e = ζ · σ e ). The initial mutation scaling is set to ζ = 0.1 for all evolutions. At this point it should be noted that the functions for fitness, crossover and mutation are empirical functions found to yield good results (high performance due to learning curves with steep slopes) for this specific problem. Mathematical formulations of the functions we use for fitness, crossover and mutation are given in the Appendix.
The size of the population is n pop = 30 and the total number of bodies in each simulation is n tot = 38, including the Sun, Venus, Earth, Mars, Jupiter, Saturn, Uranus and Neptune. The mass of Mercury was added to the Sun. The orbital elements for the massive bodies are obtained from the JPL HORIZONS system. We perform the simulations in the Jovian-centric system, as we expect only negligible changes of the results because the systematic error produced by the GA is significantly larger than the error made by not using a barycentric frame of reference. We use the Lie-Series n-body integrator as described in Hanslmeier and Dvorak (1984) with a numerical accuracy of ε = 10 −11 . Note that we do not explicitly include Ganymede and Callisto in the simulations due to a performance gain. In contrast to symplectic n-body integrators, the stepsize of the Lie-integrator is limited by the minimum relative distance between any of the objects, where h is the stepsize in days, k G is the Gaussian gravitational constant and ε is the numerical accuracy. The quantities ν and D are obtained from the integrator itself during each step, where ν is associated with the number of used Lie-terms and D depends on the masses and relative distances between the objects (Eggl and Dvorak, 2010). If we include the moons in the calculations, the stepsize would be limited by the distances within the Jupiter system and much higher computational capacities would be needed in order to get the same amount of results. Therefore, we use pseudo-moons which are characterized only by the distance to Jupiter without gravitation in order to be able to calculate distances. A measurement is taken if a particle has its closest approach within the Hill sphere of a moon. The stepsize of the integrator during a Hill sphere crossing has to be low enough to ensure that the maximum spatial distance between two consecutive frames does not exceed the size of the Hill sphere. Otherwise, the particle would overleap the Hill sphere and no measurement would be detected. Therefore, an additional restriction for the upper bound of the stepsize is used. Moreover, if the moons are included in the simulations, their very short orbital periods would lead to significant stability issues with their orbits even at very low stepsizes. Note that since no actual positions of the moons are given during the simulations, the Hill spheres are represented by torus-shaped regions around Jupiter with their first radius being the relative distance to Jupiter and the second radius being the Hill radius of the respective moon. This approach can be used because of the statistical nature of the study, expecting a large number of measurements. Since we are doing statistical studies of typical intersection velocities we found this simple approach to be highly effective. Due to the nature of the GA, we are able to use a large parameter space for the random initial conditions of the test particles (see Table 1). We set the maximum simulation time t max = 165 yr in order to allow at least a full orbital period for each individual. The maximum number of generations g is set to g max = 250. A new iteration is started either if there is no measurement (i.e., close encounter with either Ganymede or Callisto) within these g max iterations or as soon as such a close encounter is recorded. Finally, for a quick check of the results with the GA method, we also check its output with the ones from full integrations of Centaurs (through all their lifetime) taking data kindly provided by , see section 3.2 for more descriptions. We compare these data with the orbits found for Ganymede, as far as close encounters with Jupiter were taken in account up to d = 0.01 AU.

Analysis of GA measurements
After about 3300 CPU hours, total number of 531 and 625 measurements was obtained for Ganymede and Callisto, respectively. Each measurement corresponds to an individual evolution process of the GA and represents an individual close encounter scenario which contains all relevant information about positions, velocities, orbital elements, etc. about the bodies at both the time of closest approach and at the beginning of the respective simulation.  The terms pro-pro, pro-ret, ret-pro and ret-ret refer to the corresponding geometries, where the first part describes the prograde or retrograde movement with respect to the Jovian System (with other words: i < 90 • is pro and i > 90 • is ret) and the second part describes the direction of flight with respect to the corresponding moon within the Jovian System, moving parallel (pro) or antiparallel (ret) to the respective moon. More detailed statistics of the four classes are given in Table 2  Further observations can be obtained from the results: • The overall form of the relative velocity histograms can be reproduced by overlapping Gaussian distributions which are represented by the four classes.
• The classes are overlapping stronger for Ganymede than for Callisto, even swapping places (in velocity) when comparing class 2 and class 3 for Ganymede.
• There is a clear trend favouring retrograde encounters (for both Jupiter and the respective moon), with most close encounters being ret-ret (36.3% for Ganymede, 33.8% for Callisto).
The first row in Figure 2 shows the final semi-major axes at the time of the closest approach. While Ganymede shows two distinct peaks at a ≃ 5 AU and a ≃ 20 AU, Callisto only has one single, large peak at a ≃ 9 AU. The other panels shows the respective quantity at starting time of the last generation versus the same quantity at the time of the respective measurement. A 1:1 relation corresponds to no change of the orbit during the simulation, while a large scattering means significant changes.
As we can see for the eccentricities, there is a significant scattering towards larger values. The plots imply that the scattering for Ganymede close encounters is higher than for Callisto. More than half of the particles even undergo the transition from elliptic (e < 1) to hyperbolic (e > 1) trajectories. The strong scattering is also apparent in the semi-major axes. The inclination data shows only minor scattering between i = 50 • and i = 150 • , which implies quite stable configurations within the simulation time. However, at low inclinations particles are scattered over the full parameter space and prograde ones tend to get scattered stronger than their retrograde counterparts.
Several other obervations can be obtained from the datasets: 1. The needed number of generations for taking the measurement grows steadily until about g = 100 before falling again rapidly.
2. The most dramatic changes of initial orbits happen during the first few tens of generations.
3. It is easier for the GA to find intersection orbits after a short simulation time. Therefore, the time of measurement peaks towards low values with a mean simulation time of t mean = 45.5 yr.

The intersection probability is approximately twice as high for Callisto because
her Hill radius is larger than Ganymede's. For statistical reasons, we therefore assigned more computation time to Ganymede.
5. Two actual collisions are measured for Ganymede (with impact velocities of v rel = 34.6 km/s and v rel = 42.2 km/s, respectively), none for Callisto.

Preliminary Parent Bodies' Orbits
We take the orbital evolution of Centaurs to do a comparison between the predicted Centaurs' orbits at close encounters and the Jovian moons. We integrate forward a sub-sample of the Centaurs for 30 Myrs and check all the close encounters with Jupiter, using the Lie-integrator. This study considers only close encounters up to a distance of d = 0.01 AU from Jupiter similarly to , but with an encounter radius 4 times smaller at the cost of computational time. Thus, the comparison is limited to Ganymede only (since he has a distance of d ≃ 0.0072 AU to Jupiter. Callisto has a distance of about d ≃ 0.0126 AU to Jupiter.) 3 . We take 319 Centaurs with 15 clones in each interval ranging over 5 AU in semimajor axis (for a total of 5104 bodies) 4 in order to quickly get a statistical sample which covers the entire Centaur region from 5 AU to 30 AU in semi-major axis. This approach is sufficient to give an idea of these kinds of bodies approaching Jupiter and its moons and to have a compareable sample to the orbits produced by the GA.
From the evolution of 5104 objects, a total number of 292 measurements was obtained for Ganymede close encounters. From our sample of Centaurs we find that ∼ 22.6% can have close encounters with Jupiter. As the percentage of Centaurs which can cross Ganymede orbits is about 20.1% (8.7% with e < 1), almost all the Centaurs' close encounters with Jupiter (∼ 89%) can reach Ganymede's orbit in the range of its Hill sphere. Figure 3 shows the comparison between GA and Lie integrations for the respective orbital sample output in the Hill sphere of Ganymede. Figure 3 suggests that the datasets are in quite good agreement, with the GA and the Lie datapoints covering a similar area in the a − e phase space. However, a significant sample at high inclinations is apparent in the GA case, which makes the GA measurements not fully consistent with the measurements obtained by the Lie-integrations. An explanation for this effect is given in section 4. It turns out that these high-inclination datapoints are also responsible for the dense cloud in the a − e phase space between a = 15 − 30 AU and e = 0.7 − 0.8. Otherwise, the distributions in the a − i space are reasonably in agreement, especially at semi-major axes below a = 15 AU. Even the major cloud between a = 2.5 − 15 AU and i = 50 − 140 • is reflected well. In order to quantify the resemblence of the two distributions obtained via the GA and Lie-integrations, we deploy a two-dimensional Kolmogorov-Smirnov test (Press and Teukolsky, 1988). Due to the reasons mentioned in section 4 we exclude the nonphysical high-inclination datapoints with i > 160 • . For the a − e distributions we get a p-value of p = 0.27, for the more diverse a − i distributions p = 0.058. While the a − e distributions do not differ significantly, the a − i distributions show significantly less correlation, which indicates the need for a refined GA method as laid out in section 5.

Conclusions
The velocities of the four classes as described in section 3.1 can be explained by geometrical considerations: Class 1 encounters happen if the particles experience deceleration from the Jupiter flyby and additionally 'lose' relative velocity due to the parallel direction of flight compared to the respective moon. The intermediate classes 2 and 3 either experience deceleration or acceleration from the Jupiter flyby and additionally 'gain' or 'lose' relative velocity due to the antiparallel or parallel direction of flight compared to the respective moon. Both acceleration by Jupiter and a 'gain' of relative velocity is true for class 4.
Comparing the numbercounts of the classes, a correlation exists with higher fractions of both types of retrograde encounters being more probable. The higher fractions of retrograde encounters can be described by a simple geometrical effect: While a particle moves within the torus-shaped Hill region (the complete area which is accessed by the moon's Hill sphere during a full orbital period), the probability for a retrograde measurement is much higher than a prograde one since the typical relative velocity between the particle and the moon is high. This finding supports the observation of the heavily cratered front-side of the rotationally bound moons. For example, many large crater systems on Callisto (namely Valhalla, Asgard, Adlinda, Utgard, etc.) are located at the front-side.
Comparing class 2 and class 3, it seems that the influence caused by the flyby at Jupiter and the parallel or antiparallel direction of flight are similar in strength, produc-ing similar numbercounts for these classes. More discreet or weaker effects like the selection of the parameter space and the behaviour of the GA may influence the shape of the curves as well as the numbercounts for the classes. However, the overall shape of the individual components can be interpreted as Gaussian distributions, covering a large interval for possible close encounter velocities.
Recalling the distributions of semi-major axes in Figure 2, the second peak for Ganymede is caused by the stronger scattering of the classes ret-pro and ret-ret. This means that Jupiter scatters retrograde orbits stronger than prograde orbits at smaller distances to Jupiter. Many particles end up with high eccentricities, regardless of their initial values. Note that the simulation time of t max = 165 yr indicates that these drastic changes are caused by a single or only a handful of close encounters with Jupiter. The stronger scattering towards higher values of the semi-major axes and eccentricities for Ganymede is intuitive, because particles experience a stronger acceleration by Jupiter during the flyby. The data implies that more than half of all particles even undergo the transition from eccentric to hyperbolic trajectories during the simulations (taking into account we only considered elliptic orbits initially). The comparison between the final semi-major axes indicates that the large peak (at Callisto) vanishes at closer distances to Jupiter (towards Ganymede) and gets scattered over a wider range at typically higher values.
For the inclination, the GA has a selection effect for the sector a > 15 AU and i > 160 • , although the initial setup for the first generation resembles a realistic inclination distribution. This is clearly a selection bias of the GA, as it forces the population to settle either in very low or very high inclination orbits. This settling is due to the fact that statistically the relative distances between the test particles and the Jupiter system are smaller if their orbital planes are aligned. This leads to higher fitness values at very low/high inclinations. However, as a first approximation the GA works well as a predictor for close encounters for prograde orbits. Therefore, we recommend using the GA only for prograde orbits at current state.
The comparison with the numerical integations (see Figure 3) reveals that the GA represents the overall close encounter situation quite well in the a − e space and pretty well in the a − i space when excluding retrograde orbits with very high inclinations. Both density distributions cover similar areas in the parameter spaces with almost all datapoints from the Lie measurements having their counterparts in the GA measurements. The statistical analysis via a Kolmogorov-Smirnov test even suggests a correlation between the density distributions for the a − e phase space. Note that the simulations in the GA have a very short integration time of t max = 165 yr, which do not allow for modelling the long-term behaviour of chaotic orbits. However, the GA still yields a good coverage of parameter space, even with this short simulation time. Interestingly, the GA also finds orbits which lie far beyond the initial parameter space of semi-major axes given in Table 1. This indicates, that the implementation of the mutation function of the GA is flexible enough to explore the parameter space beyond its limits. The a − i parameter space is reproduced for inclinations up to i = 160 • , although the statistical analysis suggests no significant correlation in density. As already stated, the high number of retrograde measurements is a selection bias of the GA.
For future work with the GA we may include the simulation time and further properties of the orbits in the fitness function in order to avoid too short simulations and make high inclinations less likely, for example. We may also include the TNO region as an important extension to the parameter space.
The GA can easily be adapted in order to efficiently measure actual collisions with any given object. In this work we measure only two collisions with Ganymede because as soon as the first particle has its closest approach within the Hill sphere, a completely new evolution process is started.
In summary, the results from the GA are not fully consistent compared to the classical approach because the underlying principles are different. GAs in general tend to find all possible solutions to a given problem rather than only the realistic, physical ones. However, in our implementation this effect can be overcome by optimizing the GA either to avoid clearly non-physical solutions or to enhance realistic solutions by increasing the simulation time, enlarging the parameter space and refining the functions for fitness, crossover and mutation.

Discussion
The comparison with the Lie integrations reveal that the genetic n-body algorithm yield both a high number of physical as well as non-physical results. For example, an unrealistic high number of retrograde orbits are found while using a realistic probability distribution for the random initial inclinations. This is clearly a selection effect of the GA, as it finds that the probability for measuring retrograde encounters is significantly higher than for prograde ones. Several reasons such as a too short simulation time, too powerful crossover-and mutation functions or the choice of hyperparameters (such as the population size, number of generations, etc.) can be responsible for causing this high fraction non-physical results. However, it is expected that this high fraction can be significantly reduced by applying one or more of the following improvements: • A higher simulation time for each generation enables more dynamical effects in general. Therefore, the GA will also tend to produce more physically motivated results.
• The fitness function can be refined to avoid solutions which are clearly nonphysical. For example, one may introduce additional terms which depend on the inclination, the time of measurement, etc.
• Since crossover tends to find orbits within the initial random parameter space, this parameter space should be large enough, e.g. 50% larger compared to the parameter space of interest. Regions outside the inital random parameter space are only accessible via mutation.
• We found the mutation function to have a significant effect on the behaviour of the GA's learning curve (generation vs. mean fitness). An optimized mutation function can boost the overall performance of the GA drastically, enabling the use of a higher simulation time, a larger population, obtaining more measurements, etc. with the same computation resources.
• A larger population is able to cover the parameter space more homogeneously and may reveal further close encounter families with low probabilities.
• Like in this work, the non-physical results can be efficiently filtered by comparing the results with classical approaches.
However, the GA also yields useable results, especially on the small-scale close encounter dynamics. We find it to be an efficient tool to get a rough idea of the underlying dynamics of the problem and the expected families of solutions before investigating into more detailed analysis with classical approaches. The GA supports the use of existing approaches rather than replacing them. In this work, the GA efficiently finds all possible close encounter geometries even beyond the initial parameter space with low computational effort. The measurements cover all major areas of the parameter spaces in semi-major axis, eccentricity and inclination. Even a weak correlation in the density distributions for the a − e space is apparent when comparing the results of the GA with the long-term Lie integrations. This is quite interesting, as the used methods represent completely different approaches to the underlying problem of close encounters. This encourages more detailed studies with optimized algorithms. In theory, the GA can be used to study a variety of different problems in celestial mechanics, given the appropriate fitness function and adapted functions for crossover and mutation. Apart from the behaviour of the GA itself, relevant information is obtained from the measurements: • The four classes, which are motivated by geometrical considerations, can be distinguished well in the datasets. The classes allow for a more detailed and structured way for analysing close encounter events in the Jovian System.
• One may distinguish between different impact scenarios depending on the impact velocity. For example, there is no need for analysing retrograde collisions if the particle is classified as class 1 and vice-versa for prograde collisions and class 4.
• There are significantly more retrograde than prograde encounters for both moons. This fact is supported by the heavily cratered front-side of the rotationally bound moons.
• As shown in Figure 2, most of the particles get scattered quite drastically in eccentricity during the close encounter. A high fraction may even end up in hyperbolic trajectories.
• Moreover, the distributions of semi-major axes reveal a double-peak structure for Ganymede in contrast to a single-peak structure for Callisto. This can be explained by a stronger scattering of the classes ret-pro and ret-ret at smaller distances to Jupiter.

Author Contributions
PW wrote the code for the genetic algorithm and the transformations for the osculating elements. He wrote about 70% of the paper, especially the part inherent to the genetic algorithm, the core of this paper. MG has contributed performing the full orbital numerical part section, selecting the initial populations for this study, he suggested how to consider the close encounters with the Jupiter moons and he edited several parts of the paper for quality improvements in any section of the paper, apart the description of the genetic algorithm. He wrote about 15% of the paper. TM helped in editing the quality of the paper and wrote about 15% of the paper. He also helped in the description of the genetic algorithm and performing the Kolmogorov-Smirnov test.