A Deep Learning Neural Network Algorithm for Classification of Eclipsing
  Binary Light Curves
We present an image classification algorithm using deep learning
convolutional neural network architecture, which classifies the morphologies of
eclipsing binary systems based on their light curves. The algorithm trains the
machine with light curve images generated from the observational data of
eclipsing binary stars in contact, detached and semi-detached morphologies,
whose light curves are provided by Kepler, ASAS and CALEB catalogues. The
structure of the architecture is explained, the parameters of the network
layers and the resulting metrics are discussed. Our results show that the
algorithm, which is selected among 132 neural network architectures, estimates
the morphological classes of an independent validation dataset, 705 real data,
with an accuracy of 92%.
Burak Ulas
# A Deep Learning Neural Network Algorithm for Classification of Eclipsing Binary Light Curves

###### Abstract

We present an image classification algorithm using deep learning convolutional neural network architecture, which classifies the morphologies of eclipsing binary systems based on their light curves. The algorithm trains the machine with light curve images generated from the observational data of eclipsing binary stars in contact, detached and semi-detached morphologies, whose light curves are provided by Kepler, ASAS and CALEB catalogues. The structure of the architecture is explained, the parameters of the network layers and the resulting metrics are discussed. Our results show that the algorithm, which is selected among 132 neural network architectures, estimates the morphological classes of an independent validation dataset, 705 real data, with an accuracy of 92%.

## 1 Introduction

Deep learning techniques strengthen their solid ground in various areas from art to science every passing day. In the present, countless machine and deep learning applications let researchers achieve faster and more precise results in their studies, as well as changing the daily life. Convolutional neural networks, a specialized architecture in deep learning algorithms using neural networks, give an opportunity to use powerful methods in some processes such as image recognition and classification. The prototype of these networks, recognitron, was proposed by Fukushima (1980). Lecun et al. (1998) introduced the convolutional networks by remarking their performance on variability in 2D shapes. They also noted the advantages of fast learning in their handwriting experiment. The improvement in both hardware and software technology allows taking giant leaps in the usage and development of convolutional neural networks. For instance, a famous architecture AlexNet (Krizhevsky et al., 2012) classified 1.2 million images with an accuracy value of about 85% in ImageNet computer vision challenge. CoAtNet (Dai et al., 2021) also reached 91% accuracy by improving the model capacity and introducing the hybrid models.

Eclipsing binary stars are basically stellar systems showing light variations in their light curves due to occultations of the companions' light. Their importance arises from being tools for deriving the crucial stellar parameters precisely, and therefore, allowing the researchers to determine the structure of the stars in realistic estimations.
Guinan (1993) remarked that the analyses of their light curve enable us to estimate important parameters like mass, radius, luminosity, effective temperature as well as atmospheric properties. The systems appeared in several morphological types, mainly contact, detached and semi-detached, (Kopal, 1955; Bradstreet, 2005) that can be corresponded to various phases of their evolution. Thus, determining the morphological classes of these stars relieve information about related light parameters as well as the stellar evolution in different circumstances, which opens the doors to understanding the universe better as we know the binary and multiple systems are very common in our galaxy.

Researchers made efforts in detecting, fitting and classifying the light curves of binary systems using machine and deep learning algorithms. Wyrzykowski et al. (2003) proposed an algorithm using artificial neural networks and detected 2580 binary systems in Large Magellanic Cloud based on the OGLE data (Udalski et al., 1998). Prsa et al. (2008) presented an artificial neural network trained with data points of 33235 light curve samples of detached eclipsing binaries for deriving some physical parameters of eclipsing binary stars selected from several databases. The authors remarked that the success rate of the algorithm is more than 90% for OGLE and CALEB data sets, respectively. Kochoska et al. (2020) evaluate different fitting methods and concluded that machine learning techniques are useful tools for estimating the initial parameters of the binaries. The preliminary results of a systematic classification for the light curve morphologies of eclipsing binaries from TESS (Ricker et al., 2015) performing machine learning technique were published by Birky et al. (2020). Ulas (2020) suggested a deep learning image classification algorithm for the classification of light curve morphologies of ASAS-SN eclipsing binaries with an accuracy value of 92%. Lately, Cokina et al. (2021b) introduced a two-class (detached and over-contact) classification based on the 491425 synthetic light curve data generated by ELISa software (Cokina et al., 2021a). The authors accomplished 98% accuracy with their combined deep learning architecture.

Applications of machine learning techniques on astrophysical data are developing and promise novel results in the area. The potential of the subject motivated us to apply the method to the eclipsing binary light curves. In the following section, we introduce the details and structure of the light curve data used in the study. Sec. 3 deals with the details of our code and the architecture of the convolutional neural network. The results are discussed and the concluding remarks are given in the last section.

## 2 Light Curve Data

The algorithm needs light curve images corresponding to certain morphological classes of binary stars to train the machine and perform the classification. Therefore, we collected real light curve data to construct light curve images of eclipsing binary stars. The data for eclipsing binary stars with known morphological types are provided from three main data sources in this study; Kepler Eclipsing Binary Catalog (Kirk et al., 2016), All Sky Automated Survey (ASAS, Pojmanski, 1997) and Catalog and Atlas of Eclipsing Binaries (CALEB, former EBOLA, Bradstreet et al., 2004).

Kepler light curves were accessed through the Kepler Eclipsing Binary Catalog
(Kirk et al., 2016). The authors catalogued some basic properties of 2920 binary systems and indicated a parameter for their morphological classes. The parameter (\(c\)), introduced by Matijevic et al. (2012) using locally linear embedding method, is a classification criterion for contact (\(0.7<c<0.8\)), detached (\(c<0.5\)) and semi-detached (\(0.5<c<0.7\)) binary systems. We collect 1913 binary systems (239 contact, 1253 detached and 421 semi-detached) with corresponding \(c\) parameters and used their phase and detrended flux values from the catalogue to construct the light curve images.

ASAS (Pojmanski, 1997) variable star database was also used to gather light curve data and morphological classes of the eclipsing binary stars. The variability class of the targets in the catalogue was determined by Pojmanski (2002) using an approach based on multidimensional parametric space as well as an extended method using certain Fourier coefficients. We were able to collect data and morphological classes of 5907 binary systems through database query service1. The phases for the light curves were calculated by adopting the times of minimum and orbital period values from the ACVS (ASAS Catalog of Variable Stars) list given by the author. The magnitudes were also converted to normalized fluxes by deriving the maximum magnitudes for corresponding light curves for the systems.

Footnote 1: [http://www.astrouw.edu.pl/asas/](http://www.astrouw.edu.pl/asas/)

The CALEB data were achieved via the catalogue's web page2. The author catalogued light curves and observational properties of 305 individual stars with their morphological classes. Since the catalogue contains light curves in several filters for many stars, the actual number of the data exceeds the above-mentioned value. 1632 light curves from the database were included in our study. The light curve images were constructed by using the phase and flux values given by the catalogue.

Footnote 2: [http://caleb.eastern.edu](http://caleb.eastern.edu)

The \(256\times 256\) pixels light curve images were generated by plotting the data in \(0.25-1.25\) phase interval from the mentioned databases. The total number of data was decreased after eliminating the light curves which (\(i\)) show very large scattering, (\(ii\)) have very few data points and (\(iii\)) do not resemble the light curve of an eclipsing binary system. The incorrect orbital period values, especially in ASAS Catalog, were also responsible for the decrement in the number of light curves. Additionally, the entire dataset from three databases was checked by eye to prevent misclassification. The data in each class were also balanced. Namely, we limit the number of light curve images in each morphological class is to be equal, thus, it is one-third of the total number of data in a given database (e.g. 2286 light curves from ASAS contains 762 images from each individual class; contact, detached and semi-detached). We randomly chose 657, 2286 and 585 images from the final datasets of Kepler, ASAS and CALEB, respectively. Therefore, a total of 3528 light curves from all databases were selected to use for image classification. The number of the light curves in the training set is 2823 while the validation set covers 705 data, about 20% of the training set, following the Pareto principle (Moore, 1897; Juran and Godfrey, 1999). Fig. 1 is a Sankey diagram showing the relation among the morphologies, databases and datasets based on the number of light curves. Nine samples of data with different morphologies in the training set from three databases are illustrated in Fig 2.
Figure 1: Sankey diagram showing the distribution of 3528 light curve image data in three nodes; morphology, database and type of the dataset. C, D and SD refer to contact, detached and semi-detached morphologies, respectively. **tr** indicates training set, while **val** remarks the validation data. See text for details. Diagram created using SankeyMATIC\({}^{3}\).
Figure 2: Selected images from the training set generated by using the light curve data of KIC 12458133 (a), ASAS 065227-5524.6 (b), V572 Cen (c), KIC 06545018 (d), ASAS 075602-4454.8 (e), QX Car (f), KIC 03954798 (g), ASAS 101553-6012.9 (h) and TZ Lyr from three different databases. C, D and SD refer to contact, detached and semi-detached morphologies, respectively. Figure created using gnuplot4.
## 3 Architecture of the Neural Network

A Python (Van Rossum and Drake, 2009) code6 was written to set a deep learning neural network algorithm and thus train the machine to classify the light curve images generated from the light curve data. The code contains import procedure of the main and axillary sources in the preamble, that are _NumPy_ package (Harris et al., 2020), _os_ module (Van Rossum, 2020), _pandas_ library (Wes McKinney, 2010; The Pandas Development Team, 2020), _Random_ module (Van Rossum, 2020), _TensorFlow_ platform (Abadi et al., 2015) and _Keras_ API (Chollet et al., 2015) which are needed the processes in the algorithm work. It proceeds with seed fixing for NumPy, Python and TensorFlow to avoid randomness and make the results reproducible, and yet randomness that may arise from the calculations on the Graphics Processing Unit (GPU) still remains. It must be noted that when running the code on a GPU, randomness may alter the results slightly from one run to another due to the parallel operations, as remarked by the Keras team. The problem can be solved by conducting the calculations on a Central Processing Unit (CPU), however, neural networks are computationally expensive, and it takes an extremely long time to achieve the results on a CPU.

Footnote 6: [https://github.com/tensorflow/tensorflow/tensorflow](https://github.com/tensorflow/tensorflow/tensorflow)

The image data in three folders (C: contact, D: detached and SD: semi-detached) were indicated in the code and the total number of the data files inside those directories was commanded to display on the output. The sizes of input images were also defined. We applied data augmentation by adding random Gaussian blur to images in the training dataset. Augmentation enriches information about data by applying certain operations to a given dataset, and it helps prevent overfitting (Shorten and Khoshgoftaar, 2019). As some augmentation methods (e.g. flip and rotation) may cause changes in the shape of the light curve, we avoided employing further augmentation to our data. Training and validation directories were defined to direct the algorithm to the targeted data which was intended to be dealt with. The data generation process resizes the images to \(128\times 128\) pixels to save computing time and converts to greyscale since colour is not a distinctive feature of our data. In data generation, class_mode argument was categorical which defines 2D _one-hot_ encoded labels that contain one _hot_ (1) among all other _cold_ (\(\mathfrak{\theta}\)) values (Harris and Harris, 2007).

The backbone of the algorithm is the sequential model, a stack of layers, which includes several hyperparameters forming the convolutional neural network architecture which consists of convolutional, pooling, fully connected and output layers. Convolutional layers use kernels with the size of \((3,3)\), referring to the size of the convolutional window (Chollet et al., 2015). Rectified Linear Unit (ReLU) function, \(f(x)=max(0,x)\), was selected as the activation function for the layers, which assigns \(0\) for the values smaller than or equal to zero (Goodfellow et al., 2016). The training is done by using the stochastic gradient descent algorithm (Ruder, 2016) to achieve the converged result and then ReLU provides relatively more effortless optimization and calculation since it represents the nonlinearities with two linear functions. The convolution operation was done by applying a \(L2\) regularization penalty (Cortes et al., 2009). The regularization term is:
\[\lambda\sum_{i=1}^{N}\omega_{i}^{2} \tag{1}\]

where \(\lambda\) (=0.001 in our case), \(\omega\) and \(N\) are the regularization parameter, weight and the number of features, respectively. The term adds the squared weights to the loss function and controls the weights to be relatively small values, thus, preventing the model from overfitting and structural complexity. The padding hyperparameter was adjusted to same, which guarantees that the feature map is the same size as the input (Chollet et al., 2015). The stride value was left default, \((1,1)\) which corresponds that the filter moves one pixel at a time. We also applied max pooling operation (Christlein et al., 2019) with a pool size of \((2,2)\) between convolutional layers. It basically downsamples the input data by taking the maximum values within the pool size. The pooling also helps avoid overfitting and lowers the computation time. The above processes lead to the feature extraction and the next stage, the flattening operation, is necessary since the final convolutional layer does not cover the entire dimension of the input image (Basha et al., 2020). Flattening converts the data into a 1-dimensional array, the shape which is mandatory to make the algorithm be able to perform the classification. A Dropout layer with a rate of 0.5 follows the flattening, which avoids the model from overfitting, as mentioned by Srivastava et al. (2014). The last steps of the convolutional neural network include fully connected layers where the classification takes place. All the input neurons are connected to the neurons in the present layer at this stage (Geron, 2017), therefore, the dimensionality of the upper Dense layer is equal to the filter number of the last convolutional layer. The dimension was set to 3 in the output layer of the network since we have three classes (contact, detached and semi-detached). Probabilistic distribution was determined using softmax activation function as it is appropriate for multiclass classifications using categorical cross-entropy loss function, which is:

\[L=-\sum_{i=1}^{n}p(x_{i})\log_{e}(q(x_{i})) \tag{2}\]

given by Zhou et al. (2021), where \(p(x_{i})\) and \(q(x_{i})\) denote real and predicted distributions, and \(n\) is the number of classes. Layers from the first convolutional to the second last Dense layer, generally called hidden layers, consist of 5 trainable and 6 nontrainable layers.

The Adam algorithm (Kingma and Ba, 2015) was chosen as the optimizer, which uses the stochastic gradient descent method. Adam is appropriate for multiclass problems and can be adjusted with the learning rate hyperparameter. The learning rate is basically referring to the step size (Murphy, 2012) in the convergence of the learning process. Tuning this parameter plays an important role in obtaining reliable results during calculations. Large values can result in straying from convergence, while small values may cause taking long times of training. We control the learning process by monitoring the validation loss through EarlyStopping callback, which is known to boost the performance of algorithms (Yao et al., 2007). The arguments of early stopping were arranged to stop training when no decrement in validation loss is observed in 20 consecutive epochs, and therefore, training was prevented to be overfit. Another callback, ModelCheckpoint, was also included in the code to save the best model having
the maximum validation accuracy in a model file. We compiled our model using cross-entropy loss, as mentioned before, based on accuracy evaluation. The final operation, fitting the model, was done by specifying the number of training samples per iteration (batch_size=32), generators, and the callbacks remarked above. Additionally, in our code, we stored the number of filters in convolutional layers and learning rate values in variables (11, 12, 13, 14 and \(\mathtt{lrate}\)) to be able to test various architectures quicker, only by changing the set of variables.

The aim of the neural network is to minimize the cross-entropy type loss function and reach the maximum accuracy value. Accuracy is a measure of how model predictions are close to the ones from the real model in all classes, while loss, a cost function, measures whether the predictions give the correct values (Sammut and Webb, 2017). Specifically, in zero-one loss, 0 and 1 refer to correct and incorrect classifications, respectively. To achieve the best result based on the corresponding values, we employ a total of 132 different convolutional neural network architectures (Fig. 3) with three different learning rate values (\(10^{-3}\), \(10^{-4}\), \(10^{-5}\)) which were run using NVidia T4 GPU accelerator provided by Kaggle5 platform. Although the higher accuracy values were reached in some other architectures, the optimum result was achieved in the network of 5 convolutional layers having 32, 32, 64, 128, 256 filters, respectively (see Sec. 4). The final accuracy and loss values for the models with validation accuracies larger than 0.9 are represented in Fig. 4. Models with the learning rate value of \(10^{-3}\) are not included in the figure, since their validation accuracy never exceeded 0.9.

Footnote 5: [https://www.kaggle.com](https://www.kaggle.com)

## 4 Results and Conclusion

As mentioned in the previous section, an accuracy of 92% was achieved in architecture with 5 convolutional layers having 32, 32, 64,128 and 256 filters. The learning rate was \(10^{-4}\), and it takes 142 epochs for the architecture to achieve the result, whose training loss (0.233) is slightly lower than validation loss (0.257), and the training and validation accuracies (0.937 and 0.936) are close. Thus, the model can be considered prevented from overfitting and underfitting compared to other models shown in Fig. 4. A visualization of the final architecture is shown in Fig. 5. The learning curves, accuracy and loss values for both training and validation versus epoch, are plotted in Fig. 6. The trends of the curves imply a typical good fit. The KERAS model file containing the final architecture is provided through GitHub repository6.

Footnote 6: [https://github.com/burakulas/ebclass](https://github.com/burakulas/ebclass)

A plot of the confusion matrix (Fig. 7) for the validation dataset reveals the details of the classification result. 217 of contact, 228 of detached and 201 of semi-detached systems out of total 705 were correctly classified. The maximum misclassification was seen in semi-detached binaries; 22 of them are classified as detached systems. Selected light curves among the best and the worst classified data for each of the morphological classes are given in Fig. 8. True positive (TP), true negative (TN), false positive (FP) and false negative (FN) values calculated based on the confusion matrix are listed in Table 1.
\begin{table}
\begin{tabular}{c c c c c} \hline \hline  & _TP_ & _TN_ & _FP_ & _FN_ \\ \hline C & 217 & 458 & 12 & 18 \\ D & 228 & 448 & 22 & 7 \\ SD & 201 & 445 & 25 & 34 \\ \hline \hline \end{tabular}
\end{table}
Table 1: True positive (_TP_), true negative (_TN_), false positive (_FP_) and false negative (_FN_) values of classification for validation dataset covering 705 image data.

Figure 3: Schematic representation of 44 different architectures with their filter numbers which are employed with three specific learning rate values (\(10^{-3}\), \(10^{-4}\), \(10^{-5}\)) separately, and correspond to 132 different networks.
Figure 4: Final Accuracy (a and b) and Loss (c and d) values for different architectures whose validation accuracies are larger than 0.9. Filter orders refer to the filter numbers of convolutional layers. Blue and light blue refer to the training and validation datasets with the learning rate value of \(10^{-4}\). Red and orange bars represent the training and validation results when the learning rate was set to \(10^{-5}\).

Figure 5: Visualization of the final neural network architecture. Orange, red, green, blue and black colors refer to Convolutional, Max pooling, Flatten, Dropout and Dense layers. Figure created using visualkeras for Keras/TensorFlow (Gavrikov, 2020).
The classification report, indicating metrics for the classification of the validation dataset, is shown in Table 2. Precision is the ratio of true positives to the total number of true and false positives (\(TP/(TP+FP)\)), a measure of how trustable is the model in predicting positive samples (Ting, 2010). Recall is defined as the ratio of the number of correctly classified positives to the total number of positives, \(TP/(TP+FN)\), and it focuses on positive samples. F\({}_{1}\) score is the harmonic mean of precision and recall. In addition to these metrics, the subset accuracy of the classification, calculated using accuracy_score function of scikitlearn library (Pedregosa et al., 2011), is 92%. This is simply the percentage of correctly classified samples (Tsoumakas and Vlahavas, 2007):

\[\frac{1}{|D|}\sum_{i=1}^{|D|}I(Z_{i}=Y_{i}) \tag{3}\]

where \(Y_{i}\) and \(Z_{i}\) are actual and predicted labels, while \(|D|\) is the number of multilabel examples and \(I\) takes the value of 0 or 1 for false or true statements, respectively.

Furthermore, it is worth looking up the filters and the output of convolutional layers (feature maps) of the final architecture in the way of how the machine sees and processes the light curves through the network. As an example, we demonstrate the feature maps for the light curve of KIC 03954798 as proceeded along with the filters of the first and the fourth convolutional layers in Fig. 9 and Fig. 10. For the human eye, the deeper the

\begin{table}
\begin{tabular}{c c c c c} \hline \hline  & Precision & Recall & F\({}_{1}\) score & number of data \\ \hline C & 0.948 & 0.923 & 0.935 & 235 \\ D & 0.912 & 0.970 & 0.940 & 235 \\ SD & 0.889 & 0.855 & 0.872 & 235 \\ \hline Average & 0.916 & 0.916 & 0.916 & 705 \\ \hline \hline \end{tabular}
\end{table}
Table 2: Classification report. C, D and SD refer to contact, detached and semi-detached systems, respectively.

Figure 6: Learning curves, variation of training (black) and validation (red) accuracy and loss with epoch.
Figure 7: Confusion matrix for validation dataset (705 images) as obtained by using _metrics_ module of Scikit-learn (Pedregosa et al., 2011) based on Keras model file containing the final network architecture. C, D and SD refer to contact, detached and semi-detached morphologies, respectively.
Figure 8: Selection from the best and the worst classified light curves with their estimated morphological classes. The algorithm correctly classified KIC 10723143 (a), BW Aqr (b) and KIC 06852488 (c) with the probability of 99.9% as contact, detached and semi-detached, respectively. The semidetached binaries ASAS 125523-7322.2 (d) and KIC 10191056 (e) were estimated as contact and detached systems with 99.9% probability. KIC 03530668 (e) was also misclassified as 93.1% semidetached, while its actual class is detached. C, D and SD refer to contact, detached and semi-detached morphologies, respectively.
layer is, the harder the light curve perception is.

Beside fixing the seeds in our code for reproducibility, following the 2.0 version of the reproducibility checklist for machine learning given by Pineau et al. (2020) we addressed the details of our model in the previous section. The algorithm was explained in detail with necessary mathematical descriptions. The application platform and infrastructure used were also denoted. The sample size of the data was given, the number of examples in the train and validation set was remarked, and the data preparation process was denoted (Sec. 2). The dataset and the code executing the classification are downloadable6. We defined and give the metrics of the classification and indicated the classification report, which refers to the quality of the classification.

Footnote 6: [https://github.com/faceface/face/face/face/face](https://github.com/faceface/face/face/face/face).

The scientific importance of our neural network algorithm arises from its providing a way to distinguish the morphological types of eclipsing binary systems with high accuracy only using their light curve images. It is also a lot faster than other conventional methods conducting the same process, such as a workflow covering the trials of at least two morphologies in a widely known light curve analysis software and comparing the results to choose the best one. The determination of the morphological class is vital in the analysis of an eclipsing binary light curve in order to yield physically meaningful results, therefore, our algorithm can be applied to a light curve image before its analysis to establish a rapid and reliable morphological assumption for the light curve solution.

When it comes to comparing our results to other studies using machine learning algorithms related to the morphologies, our accuracy is found to be close to that of investigations in the literature. Although they were not to deal with the morphology alone, the accuracy in the three-layer artificial neural network by Prsa et al. (2008), which focused on detached morphological classification, was higher than 90%. An image classification algorithm proposed by Ulas (2020) was also reached an accuracy value of 91%. Cokina et al. (2021b) achieved 98% accuracy through their combined classifier, which was trained with synthetic light curve data constructed using ELISa software for detached and overcontact morphologies. We did not run our code using the images generated from the light curve of the above-mentioned studies, since a classification owes its resulting accuracy to properly collected training data as well as the architecture. A complete change in the training set most probably requires modification in the network architecture and hyperparameters to achieve the same accuracy, over 90%.

The accurate information on the classes of training samples plays a vital role in the quality of the results. Therefore, as a future of the study, we plan to improve our algorithm by collecting light curve images with more accurate information on their types, namely the light curves of the systems having the morphological classes determined by analyses through human-controlled software, since hands-on modelling is the finest approach as Kochoska et al. (2020) concluded. This is projected to be done by a detailed survey of the literature for individual analyses of eclipsing binary light curves. Additionally, we publish a data collection platform6 to where the researchers from the community can upload morphological information and light curve data of human-confirmed eclipsing binary stars. Thuswise, the number of training and validation samples, another crucial parameter, is also aimed to be increased. The increasing number of space telescope data of binary stars will boost the number of samples without a doubt, as long as morphological classes are precisely determined.
Figure 9: 32 (\(3x3\)) filters of the first convolutional layer (upper panel) and 32 feature maps of the light curve image of KIC 03954798 as output from the first convolutional layer with corresponding filters (lower panel). Figure created using Matplotlib(Hunter, 2007).
Figure 10: Same as the lower panel of Fig. 9, but for the fourth convolutional layer with 128 filters. Note that human perception for the light curve is almost lost.
Finally, our code and collected data are public, therefore, it is open to be improved by tuning the hyperparameters or altering the architecture by the researchers from the area.

## Acknowledgements

This paper includes data collected by the Kepler mission and obtained from the MAST data archive at the Space Telescope Science Institute (STScI). Funding for the Kepler mission is provided by the NASA Science Mission Directorate. STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555.

## Data Availability

The datasets generated and analyzed during the current study are available in the GitHub repository, [https://github.com/burakulas/ebclass](https://github.com/burakulas/ebclass).
