{'fraction_non_alphanumeric': 0.0329697429574125, 'fraction_numerical': 0.02257421986151949, 'mean_word_length': 4.449245400041348, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
MAGIC: Microlensing Analysis Guided by Intelligent Computation

The modeling of binary microlensing light curves via the standard sampling-based method can be challenging, because of the time-consuming light-curve computation and the pathological likelihood landscape in the high-dimensional parameter space. In this work, we present MAGIC, which is a machine-learning framework to efficiently and accurately infer the microlensing parameters of binary events with realistic data quality. In MAGIC, binary microlensing parameters are divided into two groups and inferred separately with different neural networks. The key feature of MAGIC is the introduction of a neural controlled differential equation, which provides the capability to handle light curves with irregular sampling and large data gaps. Based on simulated light curves, we show that MAGIC can achieve fractional uncertainties of a few percent on the binary mass ratio and separation. We also test MAGIC on a real microlensing event. MAGIC is able to locate degenerate solutions even when large data gaps are introduced. As irregular samplings are common in astronomical surveys, our method also has implications for other studies that involve time series.


INTRODUCTION
When a distant star (called the source) gets sufficiently aligned with a massive foreground object (called the lens), the gravitational field of the lens focuses the light out of the distant star, thus making the distant star appear brighter (Einstein 1936;Paczynski 1986). For a typical source star inside the Milky Way, one can observe the time evolution of their brightness (i.e., light curves) and infer the existence and properties of companion objects to the lens by monitoring the deviations in the light curve from the single lens scenario (e.g., Mao & Paczynski 1991;Gould & Loeb 1992). This so-called gravitational microlensing technique has been frequently used to detect exoplanets and stellar binaries and are complementary to other techniques (see reviews by Gaudi 2012 andDong 2021).
Corresponding author: Wei Zhu (祝伟) weizhu@mail.tsinghua.edu.cn For a microlensing event with multiple lenses, the interpretation of the light curve can be challenging. First of all, the computation of the multiple-lens microlensing light curve can be time-consuming due to the finitesource effect (e.g., Dong et al. 2006;Bozza 2010). This is especially true when the microlens system consists of three or more objects (e.g., Gaudi et al. 2008;Kuang et al. 2021). Additionally, the likelihood landscape of the high-dimensional parameter space can be so pathological that traditional sampling-based methods may have a hard time searching for the correct solution (or solutions). This remains to be true even when the brute force search on a fine grid that is defined by a subset of model parameters is conducted. As a result, the current analysis of multiple-lens microlensing events is still case-by-case, with each event requiring hundreds of (or more) CPU hours as well as (usually) the supervision of domain experts (e.g., Dominik 1999;Gould & Han 2000;An 2005;Song et al. 2014;Zang et al. 2022;Zhang et al. 2022).
Several methods have been proposed to automate the parameter estimation problem in microlensing events with two lenses (Vermaak 2003;Khakpash et al. 2019;Zhang et al. 2021;Kennedy et al. 2021). These binarylens events are the overwhelming majority of all anomalous events in current microlensing surveys, as events with more lenses are expected to be rare (Zhu et al. 2014). Of the proposed methods, machine learning is a promising approach that may efficiently and accurately infer parameters out of binary microlensing light curves (Vermaak 2003;Zhang et al. 2021), in addition to its application to the identification (or classification) of microlensing events (Wyrzykowski et al. 2015;Godines et al. 2019;Mróz 2020). However, the existing machine learning methods for parameter inferences are not readily applicable to microlensing data from the ongoing surveys. In particular, Zhang et al. (2021) trained and evaluated their deep learning network on simulated light curves from the Roman microlensing survey ). While their model can successfully predict the posterior distributions of microlensing model parameters and identify degenerate solutions (Zhang et al. 2022), their simulated light curves contain ∼ 10 4 time stamps with equal steps and high photometric precision. Ongoing microlensing surveys from the ground, such as the Korea Microlensing Telescope Network (KMTNet, Kim et al. 2016), see light curves with lower signal-tonoise (S/N) ratios and irregular samplings (including data gaps due to seasonal, diurnal, and instrumental issues). These realistic features cannot be handled in the framework of Zhang et al. (2021).
Irregular time series are common in astronomical observations from the ground (and sometimes from space as well). These irregular time steps pose challenges to the family of recurrent neural networks (RNNs, Goodfellow et al. 2016), which is widely considered the standard machine learning method to deal with time series data. Variants of RNNs have been invented to handle irregular time series (Che et al. 2018). For example, Charnock & Moss (2017) imputed/interpolated on the originally unevenly sampled light curves and obtained evenly sampled ones, whereas Naul et al. (2018) added the time steps between consecutive observations as additional information into the standard RNN. Neither of these approaches would work well in the present situation. As the source enters or exits the caustic curve of the binary lens, microlensing light curves may change behaviours within brief time intervals. If these intervals fall into the data gap, it will be difficult for the standard interpolation methods to make up the missing features without introducing false signals (Rubanova et al. 2019).
In this paper, we present Microlensing Analysis Guided by Intelligent Computation, or MAGIC for short. MAGIC is a machine learning framework that can efficiently and accurately perform parameter inference in binary microlensing events with realistic sampling conditions. A schematic overview of the proposed framework is shown in Figure 1. The binary microlensing parameters are explained in Section 2. We divide these parameters into two groups and develop different machine learning methods in Sections 3 and 4. Section 5 describes the joint pipeline and its application to a real microlensing event. Finally, we discuss our results in Section 6.

BINARY MICROLENSING PARAMETERS
We describe the parameters that are required to model the binary microlensing events. These parameters are grouped and inferred separately.
The standard binary microlensing light curve can be defined by the following parameters (see Gaudi 2012 and references therein) (1) Here t 0 is the time of the closest approach between the source and lens under the rectilinear trajectory, u 0 is the separation between the source and the lens at t 0 , ρ is the size of the source star, and t E is the event timescale where µ rel is the relative proper motion between the lens and the source. The Einstein ring radius θ E is given by where κ ≡ 4G/c 2 au ≈ 8.14 mas/M , π rel ≡ au(D −1 L − D −1 S ), with M L and D L the mass of and distance to the lens object/system, respectively, and D S the distance to the source. Both u 0 and ρ are normalized to the Einstein radius θ E .
Three of the remaining parameters, (q, s, α), characterize the binary component. Here q is the mass ratio between the secondary and the primary components, s is the projected separation between the two binary components normalized to θ E , and α measures the angle from the binary axis to the direction of the lens-source relative motion. We follow the convention of Skowron et al. (2011) in the definition of the microlensing parameters except for α, for which our definition is offset by 180 • .

1
The baseline magnitude of the microlensing light curve is given by m base , and the fraction of the total baseline Neural CDE

Training
Marginal Density The schematic view of MAGIC. The binary microlensing parameters are divided into two groups. We use locator, which is a U-Net architecture, to infer the time-dependent parameters, namely the peak time t0 and the event timescale tE. The output is used to transform the light curve into a function of τ ≡ (t − t0)/tE, and estimator, which combines neural controlled differential equation (neural CDE) and mixture density network (MDN), is used to infer the other microlensing parameters. We refer to Sections 3, 4, and 5 for the detailed description of our method.
flux coming from the microlensing source is given by f S . That is, where A(u x , u y ) is the light curve in microlensing magnifications. The quantity, 1−f S , represents the blending fraction, namely the fraction of total flux coming from stars/objects other than the source. The parameters u x and u y measure the displacement from the center of mass of the binary (i.e., the chosen origin) to the source at time t. In the chosen reference frame, (5) In principle, since microlensing observations are usually taken by telescopes at different sites or in different filters, each event may have multiple sets of m base and f S parameters. In practice, as long as the different datasets cover a broad range of microlensing magnifications, they can all be aligned to a single dataset to relatively high precision without reference to specific microlensing models (e.g., Yoo et al. 2004;Gould et al. 2010a). We therefore only include one set of m base and f S in the present work.
The binary microlensing parameters, given by Equation (1), are treated/inferred differently in the present work. First of all, it is reasonable to assume that the baseline magnitude m base can be immediately read out of the ground-based light curve, which usually monitors (or can monitor) the source star for a time much longer than the event timescale. We therefore do not include m base in our method. The scaled source size, ρ, is constrained only when the source is very close to or crosses the caustic curve of the binary lens. For the chosen (effective) observational cadence of this work, the finite-source effect is not expected to be widely detected. We therefore fix ρ = 10 −3 and do not infer it. As will be shown later, our inference of the other parameters is not affected by whether or not the finite-source effect is present.
Of the remaining parameters, t 0 and t E define the position and extent of the light curve in time, whereas the other parameters, (u 0 , q, s, α, f S ), concern primarily the shape of the binary microlensing light curve. We show in Sections 3 and 4 that these two groups of parameters can be inferred separately. It is also worth noting that, once the other parameters are chosen, the source flux parameter, f S , can be derived analytically via the maximum likelihood method (Gould 1995). In Section 4, we make use of this property in the further optimization of the microlensing parameters.
Throughout this work, the VBBinaryLensing algorithm (VBBL, Bozza 2010; Bozza et al. 2018) inside the MulensModel (Poleski & Yee 2019) package is used to generate binary microlensing light curves. We describe the generation of the training data set in the two different networks separately.

Method
Together with u 0 , parameters t 0 and t E are the socalled Paczynski parameters (Paczynski 1986). Unlike the parameters that define the binary properties, these Paczynski parameters are constrained by the global properties of the light curve. In many low-q events, especially planetary events, one can accurately constrain the Paczynski parameters by masking out the relatively brief anomalous features (e.g., Gould & Loeb 1992). Of all three Paczynski parameters, the time-dependent ones, namely t 0 and t E , enter the microlensing phenomenon via the parameter τ (see Equation 5). For these reasons, we develop a stand-alone machine learning method to infer the time-dependent parameters, which is dubbed locator.
The basis of locator is the U-Net architecture. U-Net is a fully convolutional neural network that is originally designed to perform semantic segmentation of images (Ronneberger et al. 2015). U-Net takes the original image as the input, and outputs a mask of the same shape. The mask can be interpreted as a pixel-wise classification result, representing, for example, the probability of each pixel belonging to some specific segment. U-Net has been ubiquitously used in biomedical image processing (Siddique et al. 2021) since it was proposed in 2015. It has also seen a growing number of applications in astronomy (e.g., Van Oort et al. 2019;Makinen et al. 2021).
The architecture of our locator closely follows that of (Ronneberger et al. 2015). See their Figure 1 for an illustration. In details, locator differs in three aspects. First, all convolutional layers are one dimensional, because locator deals with time series data rather than image data. Second, the pure convolutional blocks are replaced by residual blocks (He et al. 2016) with Parametric Rectified Linear Unit (PReLU, He et al. 2015) activation functions to ease gradient back-propagation. Third, we set the channel and kernel sizes to be 128 and 7, respectively, to better match the shape of the input data. In practice, the performance of locator is not sensitive to the values of these hyper-parameters.
We design locator to identify the light curve segment that falls within a predefined time window ∆T ≡ [t 0 − kt E , t 0 + kt E ]. Here k is a parameter that can be tuned to optimize the output performance, and its choice will be discussed in Section 3.3. For each light curve, which is irregularly sampled (see Section 3.2), we first interpolate it with natural cubic splines and extract the values at n t = 4000 2 equally stepped time stamps {t j } nt j=1 . 3 These time stamps are assigned labels depending on whether or not they are inside ∆T , At run-time, locator takes the time series and attempts to produce a mask {p j } nt j=1 that best approximates the ground truth labels.
We train locator on simulated light curves that match the key properties of real microlensing light curves (see Section 3.2). The loss function that locator minimizes is given by Here L BCE is the binary cross entropy (BCE) loss function which is the standard loss function for binary classification (Goodfellow et al. 2016). The expectation runs over all time stamps of all light curves in the training set. The second term in Equation (7) is the Dice loss function (Sudre et al. 2017), which is introduced to assuage the class imbalance problem in events with small t E : The third term in Equation (7) is the regularization term to prevent the mask from wiggling up and down for more than twice The expectations in the last two terms are performed on all light curves in the training set. The output mask from locator contains information about the event peak time, t 0 , and the timescale, t E . Specifically, the centroid of the mask corresponds to t 0 and the total enclosed area is directly related to t E . In the continuum form, they are given by These integrals are evaluated via the trapezoidal rule in the discrete time series.

Data generation
To simulate light curves for the training of locator, we randomly sample 2000 time stamps out of a total duration of 270 days. This corresponds to an averaged cadence of ∼ 3 hr, with the lower and upper 10th percentiles corresponding to 20 min and 7.2 hr. For each light curve, we randomly draw t E from a truncated log-normal distribution with mean of 10 1.15 days, standard deviation of 10 0.45 days, and boundaries at 5 and 100 days. This log-normal distribution matches approximately the observed distribution of t E (e.g., Mróz et al. 2017), and the truncations are introduced to eliminate relatively short-timescale and long-timescale events, both of which may be subject to additional issues (e.g., parallax effect for long events) and are not the focus of the present work. As for t 0 , even though it can be any value during (and possibly outside) the sampled time interval, it is reasonable to believe that a very rough by-eye estimate can limit t 0 to within ∼ ±t E of the truth. We therefore randomly draw t 0 out of a uniform distribution between −t E and t E . To summarize, t E /day ∼ TruncLogNorm(5, 100, µ = 10 1.15 , σ = 10 0.45 ), Although the goal of locator is to infer t 0 and t E and not the other parameters, the training data set should be representative of the target-the binary light curves. We therefore sample the remaining parameters in the following way and ρ is fixed at 10 −3 , given that the finite-source effect is not expected to be detected for the adopted (effective) cadence. Note that the mass ratio q is limited to > 10 −3 simply because we want to focus on relatively massive binaries. With increased sampling cadence locator can in principle be extended down to lower q values. We defer a detailed exploration to future works.
With the sampled binary parameters and time stamps, we calculate the microlensing light curve using the VBBL algorithm (Bozza 2010;Bozza et al. 2018). To introduce noises to the simulated light curves, we adopt a simplified noise model. Gaussian noises with signal-tonoise ratios (S/N) of 33 in flux are assumed for all data points, and the flux values are converted to magnitudes for a baseline object with m base = 18.
Simulated light curves without significant binary features, defined as χ 2 pspl /d.o.f. < 1.25, 4 are rejected. In the end, two data sets each containing 10 5 light curves are generated, one for training and the other for validation and test.

Performance of locator
For our default run, we set k = 1/3 so that ∆T with the longest timescale (100 days) can still be fully covered in the simulation interval [−135, 135]. We train locator with the ADAM optimizer (Kingma & Ba 2014). The learning rate is set to 4 × 10 −3 initially and drops by ten percent after each epoch. Training samples are divided into mini-batches of size 128. After 13 epochs (∼2 hours) of training on one NVIDIA Tesla V100 GPU, the average loss drops to 0.16 on the validation set of size 1024. The optimized model is then used to predict the mask and estimate t 0 and t E according to Equation (11).
We evaluate locator on 16,384 light curves from the test set. The output masks are illustrated in Figure 2 and Figure 3, and the comparisons between the predicted and ground truth values of lg t E and t 0 /t E are shown in Figure 4. The accuracy of the prediction is measured by the root-mean-squared error (RMSE) between the ground truth and the prediction. They are indicated in the upper left corner respectively. We note that locator has very different performances on close (s < 1) binaries and wide (s > 1) binaries. This is especially true for the t 0 /t E parameter. We take the results of close binaries as a more reliable evaluation of the overall performance and defer to the end of this section a more detailed explanation of this phenomenon.
In the default run, locator achieves an accuracy in lg t E with RMSE= 0.049, which translates to a fractional uncertainty of 11.3% on t E . The accuracy on the ratio t 0 /t E is comparable at RMSE(s < 1)=0.062, which corresponds to an uncertainty of 1.9 days on t 0 for a typical event with t E = 30 days. An example event is shown in Figure 5, in which the light curve substantially deviates from the Paczynski (1986) model due to the caustic crossing feature and yet locator successfully recovers t 0 and t E at high accuracy.
As shown in the left panel of Figure 4, the scattering in lg t E increases for shorter events. This is probably because, as timescale decreases, fewer data points are available to locator in the parameter inference. For example, with k = 1/3 there are only ∼ 25 time stamps fall within ∆T for the shortest (t E = 5 days) event. To improve the performance of locator on short events, we select light curves with t E < 40 days and train another model with k = 2. After 8 epochs of training with the same learning rate schedule, the loss drops to 0.083 on the validation set. We evaluate the new model on 16,384 light curves with t E < 40 days. 5 The results are shown in Figure 6. The prediction accuracy on the timescale increases by a factor of ∼ 2, reaching a fractional uncertainty on t E of 6.2%. The performance on the ratio t 0 /t E remains largely unchanged.
The above investigation has two implications. First, it outlines a general strategy to obtain more accurate t E . For example, one may divide the overall range of t E into multiple smaller intervals and train separate locator models. The preliminary estimate of t E from the default run of locator can then be "polished" by the locator model that is optimized for the narrower t E interval including the preliminary estimate. Additionally, the comparison in locator performance between two k values also sheds light on how locator works. In the single-lens case, the most constraining power on t 0 comes from the central region with width of ∼ u 0 t E , whereas the event timescale t E is best constrained by     Figure 8. The performance of locator, evaluated on RM-SEs of lg tE, t0/tE and t0, as a function of k values.
10 1 10 2 10 3 S/N 10 1 RMSE of t0 Trained separately Trained on S/N=10 Trained on S/N=1000 Figure 9. The performance of locator, evaluated on RMSE of t0, as functions of the photometric uncertainty (S/N) of individual data point. Three different training strategies are used. In blue points, locator models are trained and evaluated on datasets with given S/N values separately. In orange points, locator is trained on the dataset with the lowest S/N and then applied to all datasets. In green points, locator is trained on the dataset with the highest S/N and then applied to all datasets. data points extending into the light curve wings. This seems to apply to locator as well. As shown in Figure 8, the timescale parameter can be better constrained when a wider mask window (i.e., a larger k) is used, whereas the event peak time t 0 is largely unaffected by the choice of k. More works are needed to better understand the mechanism of locator (and U-Net in general).
We have fixed the magnitude of the baseline object at 18 and used S/N=33 for all data points in all simulated events. In reality, the baseline magnitude varies from event to event, and the photometric uncertainty varies from data point to data point. However, as discussed previously, the choice of the baseline magnitude does not really matter to the performance of locator, as long as the photometric observations cover a long enough time baseline. We therefore only focus on the impact of varying S/N on the performance of locator. Rather than simulating the realistic noise conditions, 6 we study the impact of different training strategies. Several datasets with different S/N values are generated following the method in Section 3.2. We then train locator models on light curves with a fixed S/N value and then evaluate the performance on light curves from the same dataset. The result of this training strategy is shown as the blue points in Figure 9. Here we only show the RMSE values of t 0 as the trend is very similar for t E . Next, we train one locator model on events with the lowest quality (S/N=10) and apply it to all datasets, the result of which is shown in orange in Figure 9. In the third test, we train locator on events with the highest quality (S/N=1000) and apply it to all datasets. This result is shown in green in the same figure. As expected, locator performs the best if it is trained and evaluated on light curves with the same (or very similar) data quality. However, the performance of locator is only slightly worse if it is trained on the noisest dataset. This is probably because locator, similar to other convolutional neural networks (Goodfellow et al. 2016), is effectively composed of a series of filters learned from the data. When faced with noisy data, it learns to denoise through training, so the relatively large noises in the training data make locator more robust. For the same reason, the accuracy of the predicted t 0 is substantially worse on noisy light curves if locator only saw during the training light curves of the highest quality. For applications to real events, it seems that training one locator on relatively noisy data may be a good strategy that balances efficiency and accuracy.
Finally, we discuss the issue related to wide (s > 1) binaries. As pointed out earlier, the match between the ground truth labels and the predictions of locator for wide binaries is not as good as for close binaries. This larger mismatch, however, does not mean that the inference power of locator is reduced for wide binaries. After visually inspecting a large number of wide binary light curves, we find that the mismatch is related to the use of the reference system. When generating the light curves, we have followed the common convention and adopted the center of mass of the binary as the origin of the coordinate system. However, for wide binaries, the external shear due to a component at wide separation is better modeled in the reference frame defined on the center of magnification. For s > 1, the offset between the center of mass and the closer center of magnification is (e.g., Di Stefano & Mao 1996;Chung et al. 2005;Penny 2014 For very wide binaries, the event timescale given by the mass of the closest lens also takes over the timescale given by the total mass of the binary. Both can lead to large offsets between the ground truth window and the predicted window. The most extreme example that has the largest deviation in t 0 and t E with q = 0.99 and s = 2.8 is shown in Figure 7. Once the reference frame defined on the center of magnification is adopted, the match between the ground truth and the prediction of locator is substantially improved (see also Zhang et al. 2021).

Method
We develop a separate method to infer the timeindependent microlensing parameters, which we reparameterize as ω ≡ (u 0 , lg q, lg s, α/180 • , lg f S ). (13) As shown in Figure 1, this estimator consists of two main parts. A neural controlled differential equation (neural CDE) is used to extract features from the microlensing light curve, which is irregularly sampled and may contain large data gaps. The second part, mixture density network (MDN), infers the probability distribution of the microlensing parameters.
To improve the training speed and reduce the memory requirement, we first apply the log-signature transform (Morrill et al. 2021). This transform outputs a shorter, steadier, and higher dimensional sequence that summarizes the sub-step information. For the simulated light curves with 500 time stamps (see Section 4.2), we set the log-signature depth n D = 3 and the shortened length l = 100. The output thus has v = 5 dimensions.
Neural CDE approximates the underlying process of the time series as a differential equation (Kidger et al. 2020;Kidger 2022) Here the latent state z t is defined on the time interval of the input signal, (t 0 , t n ), with w = 32 dimensions. The continuous signal X t , defined on the same time interval, is obtained from applying natural cubic splines to the input, discrete signal. 7 The initial value and signalderivative of the latent state, ξ θ and f θ , are two different neural networks, where the subscript θ denotes the dependence on the learnable machine learning parameters. By Equation 14, neural CDE gradually extracts features from the signal X t according to a learned policy f θ , records its knowledge in the latent state z t , and output a terminal value z tn that summarizes the useful information of the input signal. We refer interested readers to Kidger et al. (2020) and Kidger (2022) for further details of neural CDE. The summary features from neural CDE is then fed into the MDN (Bishop & Nasrabadi 2006) to estimate the probability distribution of the microlensing parameters Here φω ,Σ denotes the density of a multivariate Gaussian with meanω and covariance matrix Σ. The normalized weight π i , meanω i , and covariance matrix Σ i of each of the n G multivariate Gaussians are given by the neural network of MDN. Compared to the masked autoregressive flow method of Zhang et al. (2021), MDN provides a straightforward and efficient alternative with explicit density estimation, which is useful for the identification of degenerate solutions and the further optimization (see below). We adopt n G = 12 Gaussians with diagonal covariance matrices. This preserves universality as long as n G is large enough (Bishop & Nasrabadi 2006). In the joint framework, estimator, the initial value network ξ θ , the evolution network f θ , and MDN are constructed as residual networks (ResNet, He et al. 2016) with three residual blocks and each block consisting of two fully connected layers of width 1024. The evolution network has an additional activation of hyperbolic tangent to rectify the dynamic. To train the model, we minimize the negative log-likelihood (NLL) loss function To further optimize the microlensing parameters, we start from mean positions of the multivariate Gaussians and minimize the light curve χ 2 values via the Nelder & 7 Note that the use of cubic splines does not mean that neural CDE is just another interpolation scheme. This is solely for creating a continuous embedding of the discrete input. In contrast, standard interpolation schemes only make use of the values at fixed time stamps. A detailed discussion on this issue can be found in Kidger (2022). Mead (1965) simplex algorithm. In this step, the source flux parameter f S is derived analytically for a given set of (q, s, u 0 , α) by maximizing the model likelihood (e.g., Poleski & Yee 2019).

Data generation
For the purpose of learning the time-independent parameters, we generate light curves as functions of the parameter τ . This is equivalent as setting t 0 = 0 and t E = 1 d. The other parameters are sampled in the same way as in Secion 3.2.
Each light curve contains 500 data points randomly drawn within the interval −2 ≤ τ ≤ 2. For a typical timescale t E = 20 d, this corresponds to an averaged cadence of ∼4 hr with a range from 24 min (lower 10%) to 8.8 hr (upper 10%). To mimic the missing data points due to bad weather and/or instrumental failures, we throw away data points within one randomly chosen light curve segment out of the 25 equal-duration segments that constitute the whole light curve. This leaves a data gap of 3.2 d for the aforementioned timescale. Light curves and noises are generated in the same way as in Section 3.2. Finally, we only keep light curves with significant binary features, defined as χ 2 pspl /d.o.f. > 2. Six data sets are generated, each having 10 5 light curves. We use five of them for training estimator Cumulative Distribution π 0 π 0 + π 1 π 0 + π 1 + π 2 Figure 11. Cumulative distributions of the normalized weights in MDN. Here π0, π1, and π2 are the top three highest weights in our MDN of 12 multivariate Gaussians. The blue curve indicates that there is one dominating Gaussian in about half of the cases, and the green curve shows that the combination of the top three Gaussians dominates the density distribution in almost all cases. and the last one for validation and test. To further understand the performance of estimator on light curves with different cadences, we have also generated data sets with various numbers of data points within the interval −2 ≤ τ ≤ 2. For these additional data sets, we have assumed irregular sampling but no gap.

Performance of estimator
We train estimator with the ADAM optimizer (Kingma & Ba 2014). The learning rate, initially set to 10 −4 , drops by ten percent after each epoch and reaches a minimum of 10 −6 . Training samples are labeled by their input parameters ω and divided into mini-batches of size 128. The training took ∼ 1 day for 50 epochs on one NVIDIA Tesla V100 GPU, and the average loss drops to −8.4 on the validation set of size 1024.
We apply the optimized estimator to 16,384 light curves from the test set and obtain the probability distribution of the microlensing parameters given by the summation of 12 multivariate Gaussians. We then apply the Nelder & Mead (1965)  input position in the parameter space. The median χ 2 values are ∼ 500 for both solutions, comparable to the number of data points in the simulated light curve and thus approximately the degrees of freedom. This suggests that the found solutions can indeed describe the data very well. Even though in ∼20%-30% of cases the model χ 2 values are substantially large compared to the number of data points, the returned parameters are still good representatives of the input parameters, as will be discussed later. Appendix A shows an example event whose model prediction gives large χ 2 and yet reasonably accurate binary parameters. We have included a mixture of n G = 12 Gaussians in the density estimate. Is it flexible enough? To answer this question, we first look at the distribution of the Gaussian weights given by estimator. As shown in Figure 11, the probability distribution is dominated by one Gaussian in ∼50% of the tested events and by the top three Gaussians with the highest weights in nearly all tested events. As for the optimized solutions, 74% and 85% of the minimum χ 2 and minimum RMSE solutions are given by the top four Gaussians, respectively. Only in 1% of the tested events does the Gaussian with the lowest weight lead to either solution. Therefore, we conclude that 12 Gaussians can satisfactorily describe the probability distribution of microlensing model parameters.
We compare the estimator predictions before and after the simplex optimization with the input parameters in Figures 12 and 13. For each tested event, we choose as the model prediction the set of parameters that is closest to the input position in the parameter space, in order to minimize the impact of degenerate solutions. The results change only marginally if we pick the Gaussian with the highest weight as the model prediction, as seen in Figure 14. We use the median absolute deviation (MAD) to evaluate the agreement between input and prediction. Compared to RMSE, MAD is more robust against outliers. As shown in Figure 12, the direct output of estimator already contains a set of param- eters that is fairly close to the input truth. Among all parameters, binary mass ratio q and separation s are of special interest to the physical interpretation of the event. The direct output of estimator yields a typical deviation of 0.057 (0.0155) for lg q (lg s), corresponding to a fractional uncertainty of 13% (3.6%) on q (s). The accuracy of the prediction is further improved after the optimization, as shown in Figure 13. Specifically, the optimized prediction can achieve median deviations of ∼ 0.03 and ∼ 0.006 for lg q and lg s, corresponding to fractional errors of ∼ 7% and ∼ 1.4% in q and s, respectively. Overall, the agreement between the input and predicted values is very well, demonstrating the high accuracy of our method on light curves with realistic sampling conditions. We have also applied the optimized estimator to light curves in the additional data sets, which have different averaged cadence but no data gap. 8 The results are shown in the lower right panel of Figure 13 for the optimized prediction. There are a few interesting trends to report. First, the performance of estimator is better, although fairly marginally, on events with 500 points and no data gap. Second, the prediction accuracy decreases with decreasing number of data points as expected, but the degradation is typically small. Specifically, with 200 data points the prediction accuracy is only worse by ∼ 10 0.13 − 1 ≈ 35% compared to the default case with 500 data points. The only exception is seen in the case with 100 data points, for which the performance is reduced by a factor of 2. This is probably because the averaged cadence becomes too low to con- strain binaries with mass ratios down to 10 −3 . 9 Nevertheless, the overall uncertainty in the case of 100 points remains fairly small and acceptable. These results may provide useful guidance to the future application of our method to real light curves with variant cadences and/or the detection of low mass-ratio planets. We show in Figure 15 an example event whose anomalous feature is only partially covered by data points. The probability distribution from estimator, represented by the multivariate Gaussians, is already at the neighborhood of the ground truth position. The optimization with Nelder & Mead (1965) algorithm leads to even better match in microlensing parameters as well as model the light curve. Events like this one, which have irregular samplings and large data gaps, pose challenges to machine learning methods that base on the standard RNN approach. 9 The averaged cadence in terms of τ is 4/100, and thus the smallest binary mass ratio that it can detect is (4/100) 2 = 1.6 × 10 −3 .
Our method also identify multiple solutions in events that suffer model degeneracy. We defer to Section 5.2 for further discussions on this issue.

Joint pipeline
Combining the methods described in Sections 3 and 4, the joint pipeline of MAGIC works as follows. First, we feed the light curve into locator, which outputs the predicted t 0 and t E . If t E is small, we can further refine it using locator with a larger k. Then we use the predicted t 0 and t E to shift and rescale the light curve. The resulting light curve is processed by estimator, which in return gives the positions and shapes of n G = 12 Gaussians in the parameter space. If needed, these Gaussians can be further optimized to obtain more accurate prediction or predictions (see Section 4.1).
The performance of the combined pipeline is largely determined by the propagation of inaccurate parameters from locator to estimator. We therefore test the joint pipeline on 16,384 light curves from the test set of locator. After applying locator, we transform Figure 15. An example event in which the binary caustic-crossing feature is partially missing due to the injected data gap, as shown by the light curve in the top right corner. Red color is for the input model, whereas blue color is for the prediction given by estimator. The corner plot shows the probability distributions of the four MDN Gaussians with the highest weights, with the colors indicating the relative weights. Contours denote the 1-σ boundaries of the multivariate Gaussians with diagonal covariance matrices. As we have chosen axis limits to optimize the overall appearance, not all four contours are visible in every panel. Values and positions of the ground truth and the prediction of estimator on the corner plot are also indicated. Even without the final optimization, the Gaussian distribution with the highest weight gives a fairly good approximation of the input parameters.
the light curves using the predicted values of t 0 and t E . The transformed light curves are then sent to estimator to infer the time-independent microlensing parameters. The set of parameters closest to the position of the input values in the parameter space is adopted as the final prediction.
The comparisons between the input and the predicted values of the time-independent parameters are shown in Figure 16. The typical, fractional deviation in q and s are 25% and 9%, respectively. These deviations remain small, even though they are worse than the estimatoronly case by a factor of ∼ 2. Overall, the joint pipeline achieves fairly high prediction accuracy in the timeindependent parameters even with the use of the inaccurate time-dependent parameters. We expect the accuracy to be further improved once the predictions are optimized via the Nelder & Mead (1965) algorithm, as seen in the estimator-only case.
A substantial fraction of binary events can be well described by a Paczynski (1986) model plus relatively short binary perturbation. For such events, one may mask out the anomalous feature and fit a single-lens model to obtain precise and accurate peak time t 0 and event timescale t E . Combining these values with estimator, we can achieve nearly as accurate time-independent parameters as we do in the estimator-alone case.
It is worth noting that the light curves generated for locator are similar to but not the same as those used for training estimator. First, locator uses binary light curves with χ 2 pspl /d.o.f. > 1.25, whereas estimator uses binaries with χ 2 pspl /d.o.f. > 2. Both correspond to a ∆χ 2 > 500 between single and binary models, under the reasonable assumption that the binary signals are usually local features. 10 Furthermore, the time window is fixed and the timescale parameters varies in the data set for locator. This means the number of data points within −2 ≤ τ ≤ 2 varies too. Nevertheless, we have shown in Section 4.3 that the performance of estimator is only mildly affected with the various length of input time series.
In addition to its high accuracy, our method demonstrates high efficiency in the parameter estimation. Specifically, locator and estimator can process ∼ 1000 and ∼ 400 light curves per second on a single GPU, respectively. As shown in Figure 13, the output of estimator is already accurate enough for purposes like statistical studies. If needed, the optimization step can be performed. Compared to the two machine learning steps, this traditional optimization step is slower, but there are only up to 12 initial positions to be optimized, which altogether may require up to a few thousand light curve evaluations. For comparison, the traditional approach usually requires optimizations at individual grid points on a fine grid in the (lg q, lg s, α) space (e.g., Albrow et al. 2000;Wang et al. 2022). Our approach therefore represents a speedup by factors of ∼ 10 2 -10 5 .

A Real Event Application
We apply MAGIC to KMT-2019-BLG-0371 (hereafter KB190371 for short), which is a short-timescale (t E = 6.5 d) binary event with q ≈ 0.08 or 0.12 (Kim et al. 2021). KB190371 is chosen for this demonstration based on several considerations. First, its causticcrossing feature is prominent and U-shaped, which is typical for binary events. Second, this event suffers from the close/wide degeneracy (Griest & Safizadeh 1998;Dominik 1999) 11 , thus enabling additional tests of MAGIC. Furthermore, except for the relatively short timescale, the other microlensing parameters are fairly typical in comparison to those used in our simulation.
The available data sets from Kim et al. (2021) are not aligned to the same magnitude system. In order to align the data sets and then apply MAGIC, we mask out the caustic-crossing feature, fit a Paczynski (1986) model to the data, and use the best-fit flux parameters to align all other data sets to the OGLE one. We note that, although here we have made use of the special property of KB190371 to align the data sets, such a data alignment task is the same as deriving the source color and can be performed without reference to specific microlensing models (Yoo et al. 2004;Gould et al. 2010a), as has been practiced in the analysis of many real microlensing events (e.g., Gould et al. 2010b;Zhu et al. 2017).
To apply our joint pipeline, we first preprocess the data and make the light curve similar to those used in the training procedure. To be more specific, we only include data points falling in the time interval of about 100 days centered at the anomaly, and filter out data points with large (> 0.1 mag) error bars. KB190371 was observed with high cadence, so we bin every five consecutive points into one. Finally, the light curve is shifted to the same base magnitude as those used in locator.
Then we feed the preprocessed light curve into locator. Because the timescale is relatively small, we set k = 2. The output mask is rectified into binary values {0, 1} with a threshold of 0.5. This step gives the prediction t 0 (HJD ) = 8592.391 and t E = 6.68 d, both of which are close to the values given in Kim et al. (2021). We transform the light curve with the predicted t 0 and t E and only keep the segment within −2 ≤ τ ≤ 2. The rest of application follows closely the procedure given in Section 5.1. With the positions of the 12 multivariate Gaussians given by estimator, we perform Nelder & Mead (1965) optimization to further refine the microlensing parameters. Because the light curve shows prominent finite-source effect, we also include in the optimization the scaled source size ρ.
The results are shown in Figure 17. The top two Gaussians with the highest weights match closely to the close (s < 1) and wide (s > 1) solutions, respectively. The corresponding model light curves from the optimized solutions both match the data reasonably well, with almost indistinguishable χ 2 values. According to these solutions, the event can be explained by a binary with (lg q, lg s) = (−1.16, −0.057) or (−0.95, 0.180). Given the imperfect alignment as well as the slightly differ-ent preprocessing procedures, we consider these values in good match with those given by the discovery paper. Kim et al. (2021) found the source scaled parameter ρ to be (6.46±0.17)×10 −3 and (6.99±0.21)×10 −3 for the close and wide solutions, respectively. These are quite different from the value (10 −3 ) that we have used to train MAGIC. Nevertheless, the optimized parameters yield ρ = 6.6 × 10 −3 and 7.3 × 10 −3 for the close and wide solutions, respectively, in good agreement with the values from the discovery paper. In other words, the presence of the finite-source effect in KB190371 does not seem to confuse MAGIC, even though MAGIC has not learned to infer the source scaled parameter ρ. This is probably because the finite-source effect in KB190371 only affects a limited number of data points around the caustic entrance and exit. While the majority of binary microlensing events bear the same feature, some may have their binary signals substantially modified by the finite-source effect (e.g., Hwang et al. 2018). Therefore, it may be necessary to include in MAGIC the ability to infer ρ. See Section 6 for further discussion.
To further demonstrate the capability of MAGIC to handle large data gaps, we exclude observations from the South Africa site of KMTNet. This creates a substantial gap during the exit of the caustic crossing, as shown in the top right panel of Figure 18. The same procedure is applied as before, and the results are shown in Figure  18. The two degenerate solutions are again identified at high accuracy, indicating the robustness of MAGIC in handling large data gaps as well as model degeneracy.

DISCUSSION
We present Microlensing Analysis Guided by Intelligent Computation (MAGIC), a machine learning method that is readily applicable to the analysis of binary microlensing light curves from ongoing and future photometric surveys. We divide the microlensing model parameters into two groups and develop separate deep learning neural networks. For the inference of timedependent parameters, namely the peak time t 0 and the event timescale t E , we construct locator, which is a U-Net architecture (Ronneberger et al. 2015). For the inference of time-independent parameters, namely (u 0 , q, s, α, f S ), we construct estimator, which uses neural CDE (Kidger et al. 2020;Kidger 2022) to extract features and mixture density network (MDN; Bishop & Nasrabadi 2006) to infer the probability distribution. The output parameters can be further refined through the standard optimization method. Depending on properties of the event under investigation, the two components can be used separately or jointly to achieve the best performance.
A little MAGIC can take you a long way. 12 We train and validate both locator and estimator as well as the joint pipeline on simulated light curves. To mimic the sampling condition in the realistic world, we introduce irregular sampling and large data gap in the light curve simulation. As demonstrated in previous sections, MAGIC can recover the input parameters at high accuracy for the majority of tested events. In particular, estimator alone can achieve typical fractional errors of ∼ 7% and ∼ 1.4% for the binary mass ratio q and the separation s, respectively. As demonstrated in the analysis of a real microlensing event (KB190371; Kim et al. 2021), MAGIC is also capable of identifying degenerate models that explain the same data. The proposed method is fast to perform, with an estimated speedup of ∼ 10 2 -10 5 depending on whether or not the optimization is needed.
While MAGIC is readily applicable to the analysis of real world data sets, a few improvements can be made 12 From Roald Dahl.
to make MAGIC more powerful and more robust. First, this work focuses on binaries with relatively large mass ratios (q > 10 −3 ), and extension into lower values may be necessary in order to apply MAGIC to the analysis of low-mass planetary events. Second, the source scaled size ρ is not inferred in the current version of MAGIC for reasons given in Section 3.2, although it can still be constrained in the optimization stage (see Section 5). The inclusion of ρ into MAGIC involves effectively both classification and regression tasks, as one will need to determine whether or not the finite source effect is detected and if yes, infer the value of ρ. Furthermore, the inclusion of additional effects-especially the microlensing parallax effect (Gould 1992) and the binary source effect (Griest & Hu 1992)-may be necessary in order for MAGIC to provide robust estimates on binaries with relatively weak signals. Last but not least, MAGIC can also be applied to the analysis of microlensing events with more than two lenses, for which the sampling-based method becomes much more timeconsuming (e.g., Gaudi et al. 2008;Kuang et al. 2021).
MAGIC uses neural CDE to handle irregular sampling and data gaps. Compared to other existing methods that have been used to handle the same issues in astronomical time series (e.g., Charnock & Moss 2017;Naul et al. 2018), neural CDE takes directly the irregularly sampled time series and models the underlying process with a differential equation. This approach has been demonstrated to achieve good performances in other disciplines that also involve time series data (e.g., Kidger et al. 2020). As irregular sampling and data gaps are common in astronomical observations, we expect that neural CDE will be useful to many other fields in time domain astronomy. Figure 17. The application of MAGIC to KB190371 and the output. The light curve and the MAGIC predictions are shown in the top right panel. Two solutions, labeled "close" and "wide," have comparable χ 2 values, and the corresponding microlensing geometries (i.e., caustic curve and source trajectory) are shown below the light curve plot. The corner plot shows the positions of the optimized solutions as well as the top two Gaussians that correspond to the "close" and "wide" solutions.

APPENDIX
A. EVENTS WITH LARGE χ 2 As described in Section 4.3, estimator achieves reasonably high accuracy in microlensing parameters but fails to obtain good light curve fits (i.e. small χ 2 values) in ∼ 20% − 30% of the test cases. This phenomenon arises because a small deviation in the parameter estimation may still lead to a large deviation in the light curve behavior, thus resulting in a large χ 2 value. Here we examine a specific example, which is the event with the largest χ 2 value in the test set. As a result, the exit portion of the caustic crossing feature is not covered with observations. Even so, MAGIC is able to identify the two degenerate solutions at high accuracy.
The predicted parameters only differ from the ground truth values by (0.02, 0.04, 0.02, 0.04) in (u 0 , lg q, lg s, α/180 • ). However, the corresponding light curve deviates from the ground truth drastically, as shown by the red curve in the top panel of Figure 19. To illustrate the impact on the light curve behaviour caused by minor changes in the parameter values, we change one-by-one the values of α, u 0 , and lg s from the model prediction to the ground truth values, and plot the corresponding light curves in the same figure. We observe that the light curve fit remains bad until all the three parameters are fine tuned. This confirms that even when the parameter estimation is accurate enough, the light curve fit may still looks bad. Nevertheless, an accurate estimation of the parameters is enough for most applications.