LEADER: fast estimates of asteroid shape elongation and spin latitude distributions from scarce photometry
Tampere University of Technology, Department of Mathematics, PO Box 553, 33101 Tampere, Finland
Received: 13 June 2017
Accepted: 25 September 2017
Context. Many asteroid databases with lightcurve brightness measurements (e.g. WISE, PanSTARRS1) contain enormous amounts of data for asteroid shape and spin modelling. While lightcurve inversion is not plausible for individual targets with scarce data, it is possible for large populations with thousands of asteroids, where the distributions of the shape and spin characteristics of the populations are obtainable.
Aims. We aim to introduce a software implementation of a method that computes the joint shape elongation p and spin latitude β distributions for a population, with the brightness observations given in an asteroid database. Other main goals are to include a method for performing validity checks of the algorithm, and a tool for a statistical comparison of populations.
Methods. The LEADER software package read the brightness measurement data for a userdefined subpopulation from a given database. The observations were used to compute estimates of the brightness variations of the population members. A cumulative distribution function (CDF) was constructed of these estimates. A superposition of known analytical basis functions yielded this CDF as a function of the (shape, spin) distribution. The joint distribution can be reconstructed by solving a linear constrained inverse problem. To test the validity of the method, the algorithm can be run with synthetic asteroid models, where the shape and spin characteristics are known, and by using the geometries taken from the examined database.
Results. LEADER is a fast and robust software package for solving shape and spin distributions for large populations. There are major differences in the quality and coverage of measurements depending on the database used, so synthetic simulations are always necessary before a database can be reliably used. We show examples of differences in the results when switching to another database.
Key words: methods: analytical / methods: numerical / methods: statistical / minor planets, asteroids: general / techniques: photometric
© ESO, 2017
1. Introduction
Photometric observations of total (diskintegrated) brightnesses are by far the most abundant source of information on asteroids in the sense of population coverage (Ďurech et al. 2015). There are currently some thousand shape and spin models from photometry for individual asteroids – for example, see the Database of Asteroid Models from Inversion Techniques (DAMIT^{1}) site – and tens of thousands more will be obtained from sparse photometry from various sky surveys such as Panoramic Survey Telescope and Rapid Response System (PanSTARRS), Gaia, and Large Synoptic Survey Telescope (LSST). Infrared photometry from, for example, the Widefield Infrared Survey Explorer (WISE) mission can also be added to visual data to obtain size and surface material parameters for thousands of targets (Ďurech et al. 2016).
Even larger populationlevel attributes can be obtained by using all available photometric measurements from the rich survey databases also when the data are not sufficient for individual models, as shown in Nortunen et al. (2017) and Cibulková et al. (2017). We refer to such data as scarce photometry to distinguish it from sparse photometry. The possibility of obtaining populationlevel attributes is based on the principle of using even a few data points per target to construct a set of brightness variation estimates for a given population (defined by the user). As shown in Nortunen et al. (2017), such estimates are sufficient for robustly reconstructing the shape elongation and (ecliptically symmetric) spin latitude distributions of the population when there are thousands of samples available. The determination of the spin distribution requires the population to be dominated by orbits concentrated near the ecliptic plane.
In this paper, we describe in detail the software package Latitudes and Elongations of Asteroid Distributions Estimated Rapidly (LEADER) for obtaining the shape and spin distributions. The mathematical principles are presented in Nortunen et al. (2017), but the computational aspects require further exposition given here. The main issues are the choice of the grid discretization level in solving the inverse problem, experimenting with various setup choices to establish the stability of the result, the “deconvolution” of the result distributions in the shapespin plot plane, and the mandatory use of simulations (from synthetic data for the same observing geometries as in the input database) to check the reliability of the inversion. We emphasize that, while the method itself is quite robust in the sense that the inaccuracy of the underlying ellipsoidal model is well tolerated since only the largescale elongation and spin estimates are obtained, the properties and coverage of the database entirely dictate the reliability of the results.
This paper is organized in the following manner. In Sect. 2 we present the main algorithm, some computational details, and the visualization by deconvolution, while in Sect. 3 we discuss the implementation of the necessary simulations from synthetic data and the comparison routines between populations. In Sect. 4 we illustrate the use of the LEADER package with examples. We sum up in Sect. 5, and in the appendix we describe some main components of the software.
2. Main algorithm for computing distributions
The main algorithm consisted of roughly three phases. First, we had a forward model for computing the brightness variation estimate η from observations. Then, we considered the inverse problem and determined the joint (p,β) distribution. Finally, we plotted the results and applied a deconvolution filter to smoothen the solution.
2.1. Forward model: brightness variations
We considered a population of N asteroids. For our analysis, we required a large population, with N ≳ 1000. In the forward model, our observable is the brightness deviation estimate η. Our objective is to obtain one or more ηs for each asteroid. With the observed brightness values L, we defined η as (1)where Δ(L^{2}) is a measure of variation for L^{2} as defined in Nortunen et al. (2017):
We used the squared brightness L^{2} for convenience; a more detailed explanation is given in Nortunen et al. (2017). We were then able to compute the amplitude A from each η, and from all the amplitude values of the population, we constructed the cumulative distribution function (CDF) C(A). To obtain the brightness deviation estimates, we analyzed the data file of asteroid i, where i = 1, ..., N, repeating the following steps for each file:

1.
From the data file, we read the Julian date, intensity, and thedirections of the Earth and Sun (these should be computed if notgiven) as seen from the asteroid translated into the origin. Wedenoted the direction vectors of the Earth and Sune_{earth} and e_{sun}, respectively. For each data point, we computed the phase angle α between the Sun and the Earth, and required that (2)In Nortunen et al. (2017), we used a limit of α_{tol} = 30° for the WISE database. For the PanSTARRS1 database, we used α_{tol} = 20°. We filtered out each brightness measurement where the phase angle exceeds this given tolerance. In addition, we required that a set of measurements has been done within a small enough change in geometry, and at least five brightness values are required for a valid η estimate. For PanSTARRS1, the phase angle α changes rapidly in time, so we could have used a condition that all measurements in a set are done within, for example, three days, to keep the change in the phase angle small. For the WISE database, more liberal rules can be used with the time span, as the phase angle changes at most 0.4° within a oneweeklong set of measurements. With this set of rules, we separated the measurements on a single data file into multiple sets. After that, we analyzed each set individually for computing η.

2.
Phase angle correction. Depending on the phase angle α, we should apply either an exponential or a linear correction to the brightness values L(α), as explained in Kaasalainen et al. (2001). Let us consider a single set of measurements. If min_{i}α_{i}< 8° in the set, we required that the phase angle changes no more than one degree, that is: (3)If the phase angle does change for more than one degree, we applied an exponential correction. In other words, we applied an exponential curve aexp(bα) into the (α,L) data, and used it to normalize the brightness data into the form (4)If min_{i}α_{i} ≥ 8° in the set, we required that the phase angle changes at most two degrees: (5)If the phase angle changes more, we applied a linear correction by using a least squares fit of the form (aα + b) into the (α,L) data, and use it to normalize the brightness data into (6)We note that with a small number of points, the correction may be unstable and can be omitted.

3.
For each measurement set, we computed the brightness variation η using Eq. (1). In Nortunen et al. (2017), we derived how η and the amplitude A are directly related: (7)We omitted any complexvalued or nonfinite amplitudes.
When all the brightness variations and amplitudes (η and A) have been computed, we sorted the amplitudes in an increasing order. Then the CDF of A is simply C(A_{i}) = i/k for i = 1, ..., k.
2.2. Inverse problem: obtaining the joint distribution for latitudes and elongations
In the inverse problem, our objective is to determine the distributions of two parameters, the shape elongation p ∈ [0,1] and the spin latitude β ∈ [0,π/ 2]. We modelled the asteroids with the shape of a triaxial ellipsoid, with semiaxes a, b and c, and a ≥ b = c = 1. With this model, the shape elongation is simply p = b/a. Here a small p value corresponds to an elongated body, while p = 1 describes a sphere. For the spin latitude, β = 0 represents a spin direction that is perpendicular to the ecliptic plane, while β = π/ 2 means the spin direction is in the ecliptic plane^{2}. Our model does not have any way of distinguishing whether the spin latitude is above or below the ecliptic plane. We assumed that the observations are concentrated near the ecliptic plane. This assumption does not usually hold entirely; the amount of variation in the ecliptic latitudes of the observations depends on the orbits of the population and the sampling epochs in database used. As a result, the computed β distribution tends to be less accurate than the p distribution. We tested the validity of the assumption, and for the databases we used in this paper, some 95% of the observations were concentrated within a ± 20° sector near the ecliptic plane.
To solve the inverse problem, we created a grid of N_{p} × N_{β} bins for our (p,β) values. For our algorithm, we typically used N_{p} = 20 and N_{β} = 29, so every bin is approximately 0.05 × 0.05 units in size, and equally spaced. We chose (p_{i},β_{j}) as a random point near the centre of each bin. Alternatively, as the values p< 0.4 are expected to have lower occupation numbers than the higher p values, we may lower the resolution for such values. Similarly, since the occupation numbers for β are expected to be somewhat proportional to sinβ (this means a uniform density on a sphere), low β values were expected to have lower occupation numbers, so we may lower the resolution of, for example, β values smaller than π/ 4.
In Nortunen et al. (2017), we derived how the CDF C(A) can be expressed in an analytical integral form when we have infinite observations available in every geometry. When we used a grid of points (p_{i},β_{j}), we were able to express C(A) as the following superposition: (8)Here w_{ij} are the occupation numbers (weights) of each bin (p_{i},β_{j}), and (9)are analytical basis functions, where Next, we constructed the data matrix M (with k rows and N_{p}·N_{β} columns) such that each column of M contains a basis function F_{ij}(A). When we write C(A) = :C ∈ R^{k}, the superposition of Eq. (8)can be written as a linear system, (10)where the unknown vector w ∈ R^{Np·Nβ} contains the occupation numbers w_{ij} of each bin. Before solving the system, the use of regularization is highly recommended, especially for the more unstable β. Let R_{p} be an ((N_{p}−1)·N_{β}) × (N_{p}·N_{β}) matrix that is meant to smooth the solution for p, and R_{β} be the respective (N_{p}·(N_{β} −1)) × (N_{p}·N_{β}) regularization matrix for β. For indices ij, we have: and similarly for R_{β}. The regularization matrices approximate the gradients at each w_{ij} in the p and βdirections. For the regularization parameters, we typically used values δ_{p} = 0.1 and δ_{β} = 1. Now, let us create an extended matrix and an extended vector : with our extended linear system being (11)To obtain the occupation numbers w_{ij}, we solved for w from Eq. (11)by using, for example, Matlab’s^{3} linear least squares method with a positivity constraint w_{ij} ≥ 0. The peak of the joint (p,β) distribution is simply the (p_{i},β_{j}) bin with the highest occupation number w_{ij}.
2.3. Visualization
To estimate the goodness of the fit C(A) = ∑ _{ij}w_{ij}F_{ij}(A) from Eq. (8), or the equivalent form C = Mw from Eq. (10), we may plot C and Mw in the same plot. The relative error ∥C−Mw∥/∥C∥ is usually less than 1% when the population contains at least 1000–2000 objects. For the actual joint distribution f(p,β), where and are the occupation numbers w_{ij} normalized such that ∑ _{i,j}f(p_{i},β_{j}) = 1, we may plot the solution (p,β,f(p,β)) either as a threedimensional surface plot, or alternatively as a contour plot. The marginal density functions (DFs) can also be computed for p and β: (12)As a postprocessing tool, we may apply deconvolution to correct for dispersion in the obtained solution. The deconvolution is used as a primarily visual tool, and it is applied only for the joint (p,β) distribution, not the marginal DFs. In order to know what kind of postsolution correction is needed, synthetic simulations (Sect. 3.1) should be performed on the database used. With the synthetic simulations, we gained understanding of the accuracy levels of the method, and were able to detect systematic errors associated with the database. As the solution tends to spread when moving away from the peak, it is a common procedure to introduce dampening to bins away from the peak. Let i^{∗} and j^{∗} be the indices for the statistical peak of the solution, that is, the bin with the highest occupation number. Then, the dampening we applied is (13)Usually it suffices to choose n = 1, but if heavier dampening is required, we may choose a larger n. As the solution of the shape elongation p is often shifted too much to the left (towards more elongated values), for example, by the amount , we may additionally shift the p values to the right (towards more spherical values): (14)For the WISE database, we chose in Nortunen et al. (2017). The more noisy the database is, the higher the required shift is. For β, the error behaviour is much harder to model and may lead to exaggerated correction, so typically we did not apply any correction in the β direction. However, we acknowledge that the solution tends to avoid extreme ends, so values near β = 0 (perpendicular to the ecliptic plane) and values near β = π/ 2 (in the ecliptic plane) have a tendency to shift away from the end points, moving towards the middle.
The computation times of the main algorithm are negligible. Depending on the size of the inspected population, reading the geometries from a database may take a few minutes, while the computation of the solution via the inverse problem is even faster. Therefore, it is easy to experiment with different grids for the solution of the inverse problem, or to test different populations. The latter means it is also fast to compare populations, which we will discuss in more detail in Sect. 3.2.
3. Main implementations of LEADER
In Sect. 3.1, we discuss an essential test when using the LEADER package: accuracy estimation by running simulations on synthetic data. Simulations are the only way to gain an understanding of how applicable the main algorithm is for a given database, and they should be performed on every database before those can be reliably used. In Sect. 3.2 we describe an extension of the main algorithm, an application for comparing shape elongation and spin latitude distributions of two populations.
3.1. Verifying the method using simulations based on synthetic data
The only way to test the correctness of the obtained solution is to run simulations based on synthetic data, where the (p,β) distribution of the artificial population is known, and see how accurately the solution is obtained. The test can additionally be used to detect systematic errors. The level of accuracy has a strong dependence on the database used (Nortunen et al. 2017), so whenever we start to use a new database, it is necessary to run simulations to see how well our method performs with the database.
The synthetic simulations begin by choosing a single peak for the (p,β) distribution. Let us denote this peak (p^{∗},β^{∗}). Now, let us assume we run the simulation for N asteroids in a population. Then, we repeat the following steps for each asteroid i, where i = 1, ..., N:

1.
We choose an asteroid model from DAMIT^{4}, with a shape elongation p = b/a (here a is the longest diameter in the equatorial xyplane, and b is the width in the corresponding orthogonal direction) that is close to the peak value p^{∗}. For example, we have set a criterion, (15)We may apply basic transformations, such as stretching, on DAMIT objects in order to get the intended shape elongation value for the asteroid. We compute the normal and area for each facet of the body.

2.
The next step is to construct the brightness data and our η estimate. We choose a β from Gaussian distribution, with β^{∗} as the mean value, and 0.05 as the standard deviation (with a restriction that β ∈ (0,π/ 2)). We also fix the longitude λ by choosing it from a uniform distribution, λ ∈ [0,2π].

3.
Next, we need geometries from the asteroid database we are studying. We read the data from a data file belonging to that database. We extract the direction vectors of the Sun and the Earth from the data file, and filter out the cases when the condition of Eq. (2)is violated. Next, we transform these vectors to the asteroidfixed frame using a coordinate transformation (Kaasalainen et al. 2001), with our fixed β and λ. We denote the direction vectors in the asteroid’s own frame e_{sun} and e_{earth}. Then we compute the total brightness L for the target. Finally, we add a minor Gaussian perturbation to L to simulate noise.

4.
From now on, we proceed as in Sect. 2.1. We require at least five L values for a valid η estimate, and if necessary, we only consider the measurements that have been done within a short time span to keep the changes in geometry small. If necessary, we apply phase correction on the sets, and finally, we compute η and A for each set as in Eqs. (1)and (7). We sort the amplitudes in an increasing order, and construct the CDF C(A) as in Sect. 2.1. As the actual p and β distributions are known with the synthetic data, it is often illustrative to plot the marginal DFs and the contour plot of the joint (p,β) distribution.
The next phase is to obtain the solution of the inverse problem. This is done in the same way as that described in Sect. 2.2. The graphical presentation of the results, as well as possibly applying a deconvolution filter, is done identically to Sect. 2.3.
3.2. A comparison of two populations
A comparison between userdetermined populations is typically desirable. Taking into account that the database used tends to cause biases in the “absolute” values of distributions, the relative differences between distributions can be expected to be more robust results from the database.
For populations S_{1} and S_{2}, we ran the algorithm described in Sect. 2. We collected the p, f(p), β and f(β) information from both algorithms, where p and β are the grid points, and f(p) and f(β) are the marginal density functions, computed from the normalized occupation numbers as in Eq. (12).
First, we computed the CDFs for the marginal DFs, denoted F_{p} and F_{β}: (16)With the CDFs F_{p}(S_{1}), F_{p}(S_{2}), F_{β}(S_{1}) and F_{β}(S_{2}) computed, we were able to compute the statistical differences as defined in Nortunen et al. (2017): (17)where we computed cases k = 1, k = 2 and k = ∞, with the scaling factors α_{1} = 1/4, α_{2} = 1 and α_{∞} = 2. As a general rule of thumb, the statistical difference between two populations can be considered significant if D ≳ 0.2. Naturally, one number does not tell everything about the quality of the statistical difference, which is why we used several different norms. Plotting the DFs (and CDFs) of both populations in the same figure is often more illustrative in terms of analyzing differences.
4. Tests with synthetic data and other examples
In this section, we demonstrate two applied examples of our method. First, we run a series of synthetic simulations to test the validity of the WISE database. In all of our examples, we used a customized version of the WISE database, which was compiled from the original data in Ďurech et al. (2016). The algorithm described in Sect. 3.1 was executed several times for different (p,β) peak values and population sizes. In the second example, we computed the (p,β) distributions for two populations using the algorithm from Sect. 2, and compared the distributions using the algorithm from Sect. 3.2. Both examples were executed using Matlab software^{3}.
We checked the accuracy of our method by observing how well we were able to reconstruct the peak of the joint (p,β) distribution. In our setup, we took a population of N asteroids from DAMIT, with the shape elongation and spin latitude and longitude known for each asteroid model. We used the geometries of the asteroid, the Sun, and the Earth that were computed for the WISE database in Ďurech et al. (2016). We observed how accurately the (p,β) peak was computed for the population, and repeated this 50 times, each time having a different, randomly generated (p,β) peak. This way, we gained a good understanding of how accurately the method computes the peaks for the given population size, no matter where the most occupied bin lies in the (p,β)plane. After this, we repeated the same simulation setup for another population size; we considered populations ranging from 100 to 5000 asteroids in order to see how the accuracy of the method improves with a growing number of asteroids in the population.
Plots of different population sizes are presented in Figs. 1 and 2. Each plot draws the actual p or β peak versus the computed p or β peak, respectively, also showing the ideal case when the solution is completely accurate. As we can see from the plots in Fig. 1, the variance in the p peak decreases noticeably when the population size increases. In addition, there is an obvious systematic error that the p peak has been shifted “down”, towards more elongated shapes. This shift is mainly due to the model (and data) noise. For example, a spheroidal shape with surface irregularities, estimated at p = 1, produces photometric variation interpreted as p< 1 by the smooth ellipsoidal model. For peaks with p ≳ 0.6, the shift is about 0.1 units, whereas for lower peaks, there is a bigger shift. In realistic populations, shapes with p< 0.5 are rare, so it is safe to assume that the value of the p peak is much higher than 0.5. Hence, it usually suffices to expect that the computed p peak is 0.1 units too low for WISE data.
In order to check whether the above result is a bias related to the asteroid database or our way of determining the shape elongation, we repeated the simulations with an alternative way to compute p. We considered the contours of the DAMITbased shapes and computed the ellipsoid (with semiaxes a and b, where a ≥ b) that best fits the contours, and then we computed the shape elongation: p = b/a. We observed that this produced no difference in the simulations, so we conclude that the bias is caused by the database; the bias caused by our definition of p is random rather than systematic. From now on, we will use our original definition of p (rather than finding the best elliptical fit for contours), as it is more simple and computationally faster. To correct the systematic error encountered with the WISE database, we shift the p values up by 0.1 units in the deconvolution phase.
Fig. 1 Synthetic simulations illustrating how the accuracy of the p solution increases when the population size increases from 100 to 5000 asteroids, with the geometries from WISE. The black, dashed “y = x” line presents the ideal case when the computed solution is completely accurate. 

Open with DEXTER 
Fig. 2 Synthetic simulations similar to the ones in Fig. 1, but for β. 

Open with DEXTER 
Fig. 3 Contour solution of the joint (p,β) distribution, computed from WISE (top) and combined WISE & PanSTARRS1 (bottom) databases. 

Open with DEXTER 
As we can see in Fig. 2, the location of the computed β peak is a coarse approximation of the actual position. The variance in the accuracy is large. While a bigger population does improve the accuracy of the peak, the variance always remains to some extent. Low β values perpendicular to the ecliptic plane have a lower variance, while the high β peaks in the ecliptic plane have a high variance, suggesting that the solution is moderately unstable if the actual β is high. Since the tail of the computed distribution tends to spread towards the ecliptic plane, as seen in Nortunen et al. (2017), it is not surprising that the computed peak tends to shift away from the ecliptic plane. Similarly for low β values, the computed peak shifts away from the low end of the β range. A systematic correction for β is complicated to implement, so we merely note that caution should be used with the obtained β solution, as our method yields coarse estimates of the spin latitude distribution.
The simulations presented in Figs. 1 and 2 were done with an equally spaced (p,β) grid. We tried simulations where we utilized a lower resolution for different p and β values, but we noticed that it reduced the accuracy of the method, especially for β, when using WISE database. The most accurate results were obtained by using an equally spaced grid. Alternatively, as it is extremely unlikely that p< 0.4, we could cut our grid by including only values of p ∈ [0.4,1]. However, we choose to include the whole interval [0,1] for completeness, since having high occupation numbers on low p values is usually a good indicator of some systematic error; that is, low p values in the grid are useful for error checking.
For future research, we are interested in testing our method with the PanSTARRS1 database (Cibulková et al. 2017). We tested the accuracy of the method with PanSTARRS1 using synthetic simulations, and found that the shape elongation p is highly accurate, whereas the solution of the spin latitude β is more unstable, as can be expected due to the scatter in the observation geometries. A more detailed analysis of using our method with the PanSTARRS1 database is presented in Cibulková et al. (2017). We are also interested in whether we can “combine” the databases by taking η estimates computed from measurements taken from both WISE and PanSTARRS1 databases. Therefore, we have performed a preliminary examination on whether our method works and is accurate with combined databases. When we are trying to determine the shape and spin distributions of a certain subpopulation, such as an asteroid family, we may not always have a sufficient number of targets in separate asteroid databases. In such cases, we may attempt to supplement the number of ηs by taking observations from multiple databases, provided our method remains accurate.
To estimate the bias caused by databases, we inverted a subpopulation of about 70 000 asteroids from the WISE database, and compared the results with those obtained from a joint WISE and PanSTARRS1 sample (containing the aforementioned WISE asteroids, and an additional subset of about 70 000 asteroids from the PanSTARRS1 database that is used in Cibulková et al. 2017). In the latter sample, we took the brightness variations η from each asteroid, and used the combined ηs to construct the CDF C(A). Here we only considered the asteroids that had brightness data available in both databases; this was to avoid selection effects, that is, to ensure that the differences between databases are not caused by observed targets being different. We plotted the computed (p,β) distributions in Fig. 3, and a comparison of the marginal DFs and their CDFs in Fig. 4 for p, and in Fig. 5 for β. No deconvolution has been used for the joint distributions of Fig. 3, as we wanted to preserve the information about multiple peaks as well as the spreading behaviour.
Fig. 4 Comparison of the marginal DFs (top) for WISE population and combined WISE & PanSTARRS1 databases, and of their marginal CDFs (bottom), for shape elongation p. 

Open with DEXTER 
Fig. 5 Comparison of the marginal DFs (top) for WISE population and combined WISE & PanSTARRS1 databases, and of their marginal CDFs (bottom), for spin latitude β. 

Open with DEXTER 
Even with half of the population taken from WISE, the combined WISE and PanSTARRS1 population provides noticeably different distributions compared to merely considering the WISE population. For the marginal distributions of the shape elongation, the differences are minor, despite the moderately high D(L^{i}) values in Fig. 4. The second double peak is so dominant in the combined database that it dampens the more elongated peak observed from the mere WISE population. In addition, the right tail of the distribution is somewhat heavy when using combined databases. Nevertheless, a visual inspection indicates that these differences in the p distributions are small, and should not be considered significant. Meanwhile, the differences are greater in the marginal distributions of the spin latitude in Fig. 5, as the β peak has been shifted towards the ecliptic plane in the combined population.
Fig. 6 Synthetic simulations illustrating the accuracy of the p (top) and β (bottom) solutions, using geometries from WISE and a combined WISE & PanSTARRS1 database. The black, dashed “y = x” lines present the ideal case. 

Open with DEXTER 
To determine whether the differences are caused by our method or the databases, we ran synthetic simulations, using geometries obtained from the abovementioned databases. Once again, we only took geometries from the targets that had been in observed in both databases. The results from the synthetic simulations have been plotted in Fig. 6, in a format similar to that in Figs. 1 and 2. The accuracy of the p solution is clearly unaffected by the combining of the databases, so the slight differences in the computed p distributions in Fig. 4 are likely to be caused by the differences in the databases. As β is more sensitive to the observing geometries, the accuracy of the β solution deteriorates faster than that of the p solution when we combine databases.
We conclude that combining databases may lead to instabilities, so combining databases may not be useful if they are constructed in a different manner. If a database contains enough targets in a subpopulation, combining it with other databases will not necessarily increase the accuracy of the method. Furthermore, this example shows how much bias the choice of database introduces, which further emphasizes that no database should be blindly used with our method. Instead, every time we are introduced to a new database, it is crucial to run synthetic simulations in order to test the validity and error levels of the method.
5. Discussion and conclusions
We have introduced the LEADER software package for the fast estimation of populationlevel distributions of shape elongations and spin latitudes of asteroids. The method itself is demonstrably robust and designed to extract all the available information from databases that contain at least a few photometric points for each target (and preferably thousands of targets), when the data are not sufficient for individual models of the targets. However, we cannot overemphasize the necessity of testing the usefulness of the database with simulations based on synthetic data at the same geometries as in the database. This is the only way to assess the reliability of the inversion result, which is essentially dictated by the database.
The synthetic simulations provide a good overview of the applicability of the database. If the geometries are undersampled or the number of reasonably good brightness variation estimates is low, the acquisition of useful populationlevel distribution estimates is not possible regardless of the method. We have also shown how different the results become when switching to another database. For example, when combining the η estimates from both WISE and PanSTARRS1, the bias introduced to the spin latitude is significant compared to using observations from merely the WISE database. As our examples show, some databases are more informative than others even if they all seem
to be extensive at a first glance. The simulations are necessary to determine this.
Acknowledgments
We would like to thank Matti Viikinkoski for his valuable comments and feedback with the software implementation. This research was supported by the Academy of Finland (Centre of Excellence in Inverse Problems). H.N. was supported by the grant of Jenny and Antti Wihuri Foundation. This publication also makes use of data products from NEOWISE, which is a project of the Jet Propulsion Laboratory/California Institute of Technology, funded by the Planetary Science Division of the National Aeronautics and Space Administration. In addition, this research made use of the NASA/IPAC Infrared Science Archive, which is operated by the Jet Propulsion Laboratory/California Institute of Technology, under contract with the National Aeronautics and Space Administration. We would also like to thank the reviewer for insightful comments that improved this paper. Last but not the least, we would like to thank Josef Ďurech and his research group for providing important feedback concerning the software package.
References
 Cibulková, H., Nortunen, H., Ďurech, J., & Kaasalainen, M. 2017, A&A, submitted [Google Scholar]
 Ďurech, J., Carry, B., Delbo, M., Kaasalainen M., & Viikinkoski, M. 2015, Asteroids IV, ed. P. Michel et al. (Tucson: U. Arizona), 183 [Google Scholar]
 Ďurech, J., Hanuš, J., AliLagoa, V., Delbo, M., & Oszkiewicz, D. 2016, Proceedings of the International Astronomical Union, S318, 170 [NASA ADS] [Google Scholar]
 Kaasalainen, M., Torppa, J., & Muinonen, K. 2001, Icarus, 153, 37 [NASA ADS] [CrossRef] [Google Scholar]
 Nortunen, H., Kaasalainen, M., Ďurech, J., et al. 2017, A&A, 601, A139 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Appendix A: Structure of the software package
In this appendix, we present the overall structure of LEADER as a roughscale pseudocode, as well as listing the key functions of the software package. First we discuss the main routine, then the variant used for synthetic simulations, and the application for comparing distributions. The software package presented here is designed for analyzing the WISE database.
The database is available for download in DAMIT, under the Software section^{5}.
Appendix A.1: Main routine
The principle of the main function leader_main_WISE has been presented on Table A.1. We assumed the population consists of N asteroids, and for each asteroid, we have a datafile available. For each datafile, we read the geometries, measurement times and brightness intensities using the function lcg_read_WISE. We split the observations into multiple sets according to the principles mentioned in Sect. 2.1: all measurements in the same set must be done within a small enough change in geometry, and at least five observations are required. Then for each set, we computed the brightness variation η and amplitude A using Eq. (7). After analyzing all datafiles, we had a list of amplitudes (the vector A_vec on Table A.1). We sorted the elements of the vector into an increasing order, and then we construct the CDF of A. In the inversion phase (function leader_invert), the matrix and the vector are generated as explained in Sect. 2.2. The Matlab function lsqnonneg computes the solution to the nonnegative leastsquares problem Finally, we used the function leader_plots to plot the solutions, and the function leader_postprocess_WISE to deconvolute the solution visually.
Pseudocode of the main function of LEADER.
Appendix A.2: Synthetic simulator
Pseudocode of the synthetic simulator implementation of LEADER.
We have presented the principle of the synthetic simulator leader_synth_main_WISE on Table A.2. We assumed a population of N asteroids. For each target, we fixed a desired (p_{wanted},β_{wanted}) value, and we chose the longitude λ_{wanted} from a random uniform distribution [0,2π]. Then, we kept selecting randomized and stretched asteroid models from the DAMIT database, until we had an asteroid with the desired p value (within a certain tolerance). Then we read measurement dates and geometries (direction of the Sun and Earth) from the database we were testing (for example, WISE). We fixed β = β_{wanted} and λ = λ_{wanted}, and used them to transform the direction vectors of the Sun and Earth into the asteroid’s own frame. Then, for each geometry, we used a scattering law to compute the brightness intensity L, adding a small Gaussian noise. The rest of the algorithm is identical to the main routine: we split the data into measurement sets, computed η and A for each of them, constructed the CDF C(A), used the subfunction leader_invert to compute the solution distribution, and finally visualized the solution with leader_plots and leader_postprocess_WISE subfunctions. It is recommended that some finetuning is done in the deconvolution function (leader_postprocess_WISE) to reshape the solution closer to the synthetic (p,β) distribution from the forward model.
Appendix A.3: Comparison of asteroid populations
Pseudocode of the statistical comparison of two asteroid populations.
The application for the statistical comparison of asteroid populations is called ast_comparison_WISE, and its principle has been presented on Table A.3. We ran the main routine leader_main_WISE first for population 1, saving the used p and β grids and their marginal DFs into variables p1, beta1, fp1 and fbeta1. We did the same for population 2, saving the grids and their marginal DFs respectively into variables p2, beta2, fp2 and fbeta2. Then we called the subfunction KS_comparison. The subfunction constructs the CDFs of each marginal DF, with Cp1 and Cp2 being the CDF of p for populations 1 and 2, and Cb1 and Cb2 being the CDF of β for populations 1 and 2, respectively. To compare the CDFs, we interpolated Cp2 and Cb2 at the grid points of population 1. Then we computed the statistical differences as defined by Eq. (17). Finally, we plotted the marginal DFs and their CDFs as in Figs. 4 and 5 to illustrate the differences of the distributions.
All Tables
All Figures
Fig. 1 Synthetic simulations illustrating how the accuracy of the p solution increases when the population size increases from 100 to 5000 asteroids, with the geometries from WISE. The black, dashed “y = x” line presents the ideal case when the computed solution is completely accurate. 

Open with DEXTER  
In the text 
Fig. 2 Synthetic simulations similar to the ones in Fig. 1, but for β. 

Open with DEXTER  
In the text 
Fig. 3 Contour solution of the joint (p,β) distribution, computed from WISE (top) and combined WISE & PanSTARRS1 (bottom) databases. 

Open with DEXTER  
In the text 
Fig. 4 Comparison of the marginal DFs (top) for WISE population and combined WISE & PanSTARRS1 databases, and of their marginal CDFs (bottom), for shape elongation p. 

Open with DEXTER  
In the text 
Fig. 5 Comparison of the marginal DFs (top) for WISE population and combined WISE & PanSTARRS1 databases, and of their marginal CDFs (bottom), for spin latitude β. 

Open with DEXTER  
In the text 
Fig. 6 Synthetic simulations illustrating the accuracy of the p (top) and β (bottom) solutions, using geometries from WISE and a combined WISE & PanSTARRS1 database. The black, dashed “y = x” lines present the ideal case. 

Open with DEXTER  
In the text 