Reconstructing Missing Variables for Multivariate Time Series Forecasting via Conditional Generative Flows

1 Introduction

Multivariate time series forecasting is critical in many real world applications, such as traffic prediction [Bai et al., 2020], tourism demand analysis [Du Preez and Witt, 2003] and air quality estimation [Du et al., 2019]. However, real world multivariate time series data are usually imperfect and incomplete. For instance, in wind speed forecasting, when there are malfuncitoning sensors, long term data unavailability arises; when there are dynamic fluctuations in resource availability, resource domain shifts arise [Chauhan et al.,

* Corresponding authors

2022]. This can be generalized as the task of Variable Subset time-series Forecast (VSF). Figure 1 shows the VSF task has two steps: 1) in the training step, a multivariate time series forecasting model can access to all variables for training; 2) in the inference step, only a subset of variables are available and the other variables are missing. The objective is to forecast the future values of the available variable subset.

The VSF task introduce two unique issues for existing multivariate time series forecasting literature: 1) Issue 1: incomplete inter-variable correlation awareness: Due to the unavailability of variables, existing multivariate time series forecasting models is incapable of modeling the interdependency between missing variables and available variables. The incomplete correlation awareness hinder forecasting performance [Guo et al., 2019; Bai et al., 2018]; 2) Issue 2: variable subset agnostic generalization: In a dynamic open environment, it is usually unlikely to know which variables will be absent during inference/testing in advance. Once available variable subset is changed, we have to retrain forecasting models, which is inefficient and untimely. Moreover, training many subset-specific forecasting models over all possible variable combinations is computationally costly and impractical in real world deployment.

To tackle the two issues of VSF, an intuitive idea is to impute missing data. Prior imputation techniques capture both inter-variable and intra-variable correlation to impute missing values [Luo et al., 2018; Tang et al., 2020]. However, the complete absence (not partial) of variables impedes their ability to capture inter-variable information (Issue 1), making such techniques unsuitable for the VSF setting. VSF remained understudied until the work in [Chauhan et al., 2022] introduced a wrapper approach, i.e., Forecast Distance Weighting (FDW), to recover missing variables from the nearest neighbors selected by ensemble weighting. However, top-K neighbors based recovery only consider partial information and could suffer if any top-K neighbor is biased, and the hyper parameter K is not generalizable. Moreover, FDW needs to re-search neighbors whenever the missing variable subset changes (Issue 1 & 2).

To address the issues, we develop deep time series generative reconstruction flow for VSF. Our approach builds upon two major contributions:

First, Time-series Reconstruction Flows (TRF) for recovering missing subset. To address the first issue, we propose to reformulate the missing variable reconstruction problem as deep collective generative task that model the dynamic temporal pattern of time series data. This generative reformulation inspires a density estimation perspective to generate the missing variables. Specifically, based on the available variable subset, we intend to estimate the unknown conditional density of unavailable variables via generative networks and rebuild the missing variables by leveraging the estimated distribution. Inspired by human reasoning [Byrne, 2019; Hu et al., 2024], which seeks to convey complex ideas through the use of familiar and simple objects, we can express the intricate density of inaccessible variables by transforming a simple distribution via a set of mapping function. To implement the mapping process, we introduce invertible Time-series Reconstruction Flows (TRF), a conditional flows-based generative framework, to reconstruct the missing variables for multivariate forecasting. The fully invertible flow structure of TRF can ensure accurate reconstruction by mapping any distribution of missing variable subset to the Gaussian distribution through forward transformation and rebuilding it with the inverse transformation. In addition, we incorporate the available subset as the conditional input, enabling us to effectively constrain the generation of missing variables by leveraging the inter-variable correlation between the missing and available variable subset. To update TRF, we utilize a log-likelihood training schema. This proves advantageous for our framework as it facilitates learning the generation process without introducing extensive constraints or losses, making our framework easy to train.

Second, Variable-Agnostic learning paradigm for TRF. Retraining forecasting models for a dynamic available variable subset and iterating through all possible variable subset combinations are both impractical and costly. In response, we propose a meta-learning based training framework, VariableAgnostic Meta Learning, to help TRF to quickly adapt to new available variables subsets. Specifically, we regard diverse variable subsets as different meta-learning tasks. For each epoch, we sample a batch of variable subsets and update the parameters of TRF for each variable subset. We aggregate the subset-specific parameters and update the TRF by this aggregated parameters. Hence, the parameters of TRF can be more generalized for all possible variable subset. For a new missing variable subset, a small number of gradient steps will yield good generalization performance, allowing TRF to quickly adapt to this new scenario and tackle Issue 2.

In summary, we reformulate the VSF problem from a novel density estimation perspective. We propose a novel flowbased framework TRF to recover the information of missing variable subset. To improve the generalization ability of TRF, we design a training framework, namely VariableAgnostic Meta Learning. Ultimately, we conduct extensive experiments to indicate the effectiveness of our framework.

2 RelatedWorks

Variable Subset Time-series Forecasting. Multivariate time series forecasting has been an abiding research topic since last century [Yi et al., 2024b; Yi et al., 2024a; Hu et al., 2023; Fan et al., 2023; Ren et al., 2022; Fan et al., 2022]. To model temporal patterns and pairwise dependencies among variables, researchers has proposed several GNN-based models, such as DCRNN [Li et al., 2017], MTGNN [Wu et al., 2020], and ASTGCN [Guo et al., 2019]. However, Variable Subset Time-series Forecast is underexplored. Missing several variables in the inference phase will considerably hinder the ability of GNN-based models to capture the inter-dependencies among variables. To tackle this problem, [Chauhan et al., 2022] propose their wrapper technique by borrowing data from training set. However, they fail to consider the distribution shift of time series and make an imprecise estimation of the missing variables. In this work, we utilize a deep learning method to estimate the density of missing variables and the empirical results demonstrate the effectiveness and superiority of our proposed solution.

Normalizing Flows. Normalizing flows are a set of bijective mapping functions which are designed to model an unknown distribution based on the given data [Dinh et al., 2014; Dinh et al., 2016]. For example, based on the coupling layer proposed by [Dinh et al., 2016], Glow [Kingma and Dhariwal, 2018] is constructed with more complex structure and achieves superior performance. MAF, proposed by [Papamakarios et al., 2017], utilizes the autoregressive mechanism for better density estimation.

Meta-learning. Dissimilar to traditional machine learning methods, meta-learning algorithms aim to acquire generalizable information and transfer the knowledge to new tasks [Vanschoren, 2018]. For instance, MAML algorithm proposed by [Finn et al., 2017] can prompt faster training in new tasks by a small number of gradient steps with a small amount of training data. Reptile [Nichol and Schulman, 2018] initials the parameters of a neural network via stochastic gradient descent on each task.

3 ProblemFormulations

Let X(N) t 2 RN denote the recording of N distinct variables at time-step t, we define historical L values as the lookback window X(N) t-L:t = ¶X(N) t-L+1 ↪ · · · ↪X(N) t ♦, and define the future H values as the horizon window X(N) t:t+H = ¶X(N) t+1 ↪ · · · ↪X(N) t+H ♦. For the classical setting of multivariate time series forecasting, we intend to find a mapping function F✓ to project lookback window into horizon window, which can be formulated as: X(N) t:t+H = F✓(X(N) t-L:t ).

However, in real life, the data collection for time series can meet a common issue: due to the long-term unavailability of data [Yick et al., 2008] and resource domain shift in time series [Hamid and Wibowo, 2018], chances are that only partial variables are available in the future [Chauhan et al., 2022]. Accordingly, this problem of time series forecasting with partial variables during inference is formulated as Variable Subset Forecasting (VSF) problem. Under such settings, we denote the available variable subset as S and the universe of variables as N, i.e., S ⇢ N and ♣S♣ ⌧ ♣N♣, where N is denoted as the set of N integers, [[1▷▷▷N]]. In the training phase, each input X(N) t is a vector in R♣N♣⇥D while in the inference phase, each data point can be represented as X(S) t 2 R♣S♣⇥D . In other words, though the models can take all variables in training, only a small subset of variables is considered in inference; formally, the VSF problem can be written as:

where Ttrain and Ttest are sets of the time-step of training data and test data.

4 Methodology

We introduce our proposed framework for VSF problem. First, we formulate VSF problem from a novel generative density estimation perspective in Section 4.1. We reconstruct the missing variables via density estimation and generation. Then, we concretely propose a novel generative network, Time-series Reconstruction Flows (TRF) to recover the information of absent variables and elaborate how to integrate our framework and forecasting models in Section 4.2. To further improve the generalization ability of the model, we propose meta training framework, Variable-Agnostic Meta Learning, to TRF in Section 4.3.

4.1 Density Estimation and Generation for Variable Subset of Time Series

In VSF problem, the unavailability of the majority variables (i.e., ♣N♣ - ♣S♣) would lead to the missing values and absent inter-dependencies between variables in S and N∩S during the inference, which largely hinder their forecasting performances. For this problem, previous VSF work [Chauhan et al., 2022] used fixed assumption of similarity and retrieved neighbors of available variables for imputation. However, it encounters constraints in completely recovering information that is missing from the unavailable subset, along with significant computational costs of measuring pair-wise similarity. In this paper, we reconsider the VSF problem from a novel density estimation perspective. We can fully recover the missing information by leveraging generative models to simultaneously reconstruct all the unavailable variables, without introducing additional computational cost.

Formally, given the fully observable training data, we try to estimate the unknown conditional density pX (XN∩S ♣XS ) of the missing variables XN∩S and available subset data XS , where X is the input space. Inspired by the human's logic of explaining intricate things via understandable objects ,we express the density of missing variables via Gaussian distribution pZ . As a result, we utilize a set of bijective functions f✓ to learn the dependencies of time series data and transform the original density into pZ, i.e., f✓ : X ! Z. Hence, we can express the conditional density pX(XN∩S ♣XS ) by the change of variables formula [Dinh et al., 2014]: pX(XN∩S ♣XS ) = pZ (f✓ (XN ∩S ; XS )) ∣∣∣det @f✓ @XN∩S (XN ∩S ; XS )∣∣∣. 4.2 Time-series Reconstruction Flows

To estimate the intricate density of high dimensional data, there are a lot of generative networks such as GAN [Goodfellow et al., 2014], VAE [Kingma and Welling, 2013]. However, GAN-based models suffer from model collapse and VAE-based model cannot estimate temporal density accurately by maximizing the lower bound of likelihood, so they are not suitable for time series data with complex temporal dependencies. Inspired by current normalizing flows [Dinh et al., 2016; Papamakarios et al., 2017], we propose Time-series Reconstruction Flows (TRF), a conditional flows-based generative model to estimate the density of missing variables of time series. Specifically, along the conditional generation formulation as proposed in Section 4.1, we consider observed variables as conditional subset embedding to reconstruct the missing variables, so that our generation framework can better access the information of available subset.

Algorithm 1 Variable-Agnostic Meta Learning

Input: Dvariables: Domain of the variables subsets, ↵, β: step size hyperparameters

• 1: Randomly initialize the parameter ✓ of the TRFlows f✓
• 2: while not done do
• 3: Sample batch of variables subset Si ⇠ Dvariables
• 4: for all Si do
• 5: Evaluate r✓LSi (f✓) with respect to N variables.
• 7: end for
• 6: Compute adapted parameters with the gradient descent: ✓0i = ✓ -↵r✓LSi (f✓)
• 8: Update ✓ ✓ -βr✓ PSi⇠Dvariables LSi (f✓0i ) 9: end while

Conditional Subset Embedding

To model the complex density of time series, the generative network needs to capture the temporal dependency among data. While directly modeling the dependencies of the whole long time series is unrealistic, we cooperate generation with the forecasting models, and design window based subset embedding for the variables. Formally, we utilize linear layers to obtain the embedding conditional vector HS as: HS = Projection(X(S) t-L:t ), where HS 2 RT⇥S⇥D and X(S) t-L:t is lookback window of available subset S.

Conditional TRF

We construct our window-based TRF consisting of conditional autoregressive layers and batch normalization [Ioffe and Szegedy, 2015]. Conditional autoregressive layers aim to estimate the distribution pX of missing variables via a couple of invertible transformation functions. Formally, we input the horizon window of missing variables X(N-S) t-L:t and accessible variables X(S) t-L:t , where X(N-S) t-L:t 2 N∩S, and X(S) t-L:t 2 S. We can rewrite the expression the original density with horizon windows as:

We can further express the conditional probability by the chain rule of probability as pX(X(N-S) t-L:t ♣X(S) t-L:t ↪ ✓) = Qi pXi (X(N-S) i ♣X(N-S) 1:i-1 ↪X(S) t-L:t ↪ ✓). Due to the dynamic pattern of time series data, the transformation process is difficult and requires a couple of layers, i.e., f✓ = f 0 ✓ ◦ f 1 ✓ ◦ · · · ◦ f N ✓ . So for the forward calculation, we define the input of the first flows layer as h0 := X(N-S) t-L:t , and the output of the last layer as hN . We can utilize the following recursion to construct the nth layer:

where, f n ✓↪si and f n ✓↪bi are the unconstrained functions to obtain scale and bias. The log determinant can be calculated as log ∣∣∣det @f n ✓ @hn (hn ; HS )∣∣∣ = - Pi f✓↪si (hn1:i-1 ; HS ).

To implement the function f n ✓↪si and f n ✓↪bi , we utilize the masking technique of MADE [Germain et al., 2015]. For the nth , the general form of masking matrix can be formulated as M W n k0 ↪k = 1mn (k0 )≥mn-1 (k), where k and k0 are the hidden unit index, and the values of mn (k) for each hidden layer n 2 ¶1↪ · · · ↪N♦ are sampled uniformly. And for the last layer, the masking rule is M V k↪d = 1d>mN (k). With the help of MADE, we can estimate the density without the sequential loop, and accelerate the training process of TRFlows.

We utilize batch normalization to further improve the propagation of signal. With batch normalization, we can not only stabilize the training process, but also train a deeper stack of coupling layers to model more complex conditional distribution. The procedure of normalization of nth layer can be written as hn 7! hn-e θ pe σ2+✏ . Here, e θ and e σ2 are the estimated batch mean and variance. We can calculate the Jacobian determinant as (Qi (e σ2 i + ✏))- 12 . To update TRF, we can derive the negative log-likelihood as the loss function by applying chain rules to Equation 2:

Forecasting Process

In the inference stage of VSF, we randomly sample noisy vectors from Gaussian distribution and employ the inverse function of equation 3 to recursively reconstruct missing variables, which can be written as: hni = exp(-f n ✓↪si (hn+1 1:i-1 ;HS )) · (hn+1 i - f n ✓↪bi (hn+1 1:i-1 ;HS )). We denote the final reconstructed missing variables as ˆX(N-S) t-L:t . We concatenate the available variables and reconstructed missing variables in the inference stage of the multivariate time series forecasting models for prediction, which can be formulated as: X(N) t:t+H = F✓(CONCAT(X(S) t-L:t ↪ ˆX(N-S) t-L:t )). It should be noted that our proposed method is applicable to any forecasting models that leverage inter-correlation.

4.3 Variable-Agnostic Meta Learning

Reconstructing missing variables set N∩S is challenging since the subset S is randomly chosen, and training one specific generative model for every possible combination of variable subsets is unrealistic. As a result, we need to improve the generalization ability of TRFlows f✓ to make a precise estimation of all possible combination of missing variables. Inspired by the meta-learning technique [Finn et al., 2017], we design Variable-Agnostic Meta Learning as the training framework to learn the common information and transfer the knowledge to new variable subset.

Formally, we define the domain of the variables subsets as Dvariables = ¶S1↪ ▷▷▷↪ Sz♦, where z is the number of variable subsets, and each variable subset Si contains N variables with different missing variables combinations. It should be noted that the selection of variable subset should be consistent with that of inference phase. We sample a batch of variable subsets Si, and when adapting to different subsets, we update the parameters ✓ of flows f✓ to ✓0i : ✓0i = ✓-↵r✓LSi (f✓) where ↵ is the learning rate of the inner loop optimization, and LSi (f✓) is the loss following equation 4. To train the model parameters, we optimize the performance of the adapted generation model f✓0i via diverse tasks sampled fromDvariables. We formulate the meta-objective as:

We target to update the model parameters in such a manner that optimal performance can be achieved via minimum gradient steps on a new task. So we employ the stochastic gradient descent (SGD) and update the parameters ✓ by computing the gradient of LSi (f✓0i ), which can be written as: ✓ ✓ - βr✓ PSi⇠Dvariables LSi (f✓0i ), where β is the meta learning rate of the outer loop optimization. The training process of the Meta-flows is shown in Algorithm 1.

5 Experimental Results

5.1 Experimental Setup

Dataset. We conducted our experiments on four real world datasets: (1) METR-LA [Li et al., 2017] includes statistics on traffic speed gathered from 207 loop detectors on county highways during Mar 1st 2012 and Jun 30th 2012 in Los Angeles. (2) SOLAR contains the solar power output of the 137 plants in Alabama in 2007. (3) TRAFFIC describes road occupancy rates collected by 862 sensors in the San Francisco Bay area from 2015 to 2016. (4) ECG5000 contains 140 electrocardiograms (ECG) with a length of 5000 each.

Multivariate Time Series Forecasting Model Backbone Setting. After missing variable reconstruction, we employed the widely-used models as backbones for our downstream multivariate time series forecasting task: (1) T-GCN [Zhao et al., 2019]: utilizes graph convolutional network and gated recurrent unit tomodel spatial-temporal correlation. (2) ASTGCN [Guo et al., 2019] exploits spatial and temporal attentions to model correlations among dynamic time series.

(3) MTGNN [Wu et al., 2020] utilizes a graph learning to describe uni-directed relations among variables.

To evaluate the average performance of our method, we randomly set the number of available variables ♣S♣ for 100 times, in order to indicate various missing percentages,. During training, all the time series backbone models were trained with all variables available in lookback windows and horizon windows. During inference (testing), we created three settings: 1) Partial: partial variables were observed and available. We set missing variables to 0. This setting aims to show missing variables will lower time series forcasting performance. 2) Reconstruct: Not just partial observed variables were available, but also the other missing variables were reconstructed by our method in both in lookback and horizon windows; in other words, our reconstruction method plays a role. This setting aims to show reconstruction can advance forecasting. By default, we adopted the Reconstruct setting during inference (testing), unless specified otherwise. 3) Oracle: All variables were observed and there are no missing variables. Yet, we evaluate the forecasting performance only on the partially available variable subset (S) in Setting 2). This setting aims show effective reconstruction of missing variables can approach the performances of all variables available. The lengths of horizon windows and lookback windows were set to 12.

Evaluation Metrics. Since the Oracle setting use data with fully observed and available variables, the performance of the Oracle setting can be seen as the upper bound of forecasting. Thereafter, we measure the performance difference between Reconstruct and the Oracle settings with respect to MAE and RMSE, which is formulated as ∆MAE = Ereconstruct MAE -Eoracle MAE Eoracle MAE ⇥ 100%, and ∆RMSE = Ereconstruct RMSE -Eoracle RMSE Eoracle RMSE ⇥100%, where, Ereconstruct MAE and Ereconstruct RMSE represent the average MAE and RMSE of the Reconstruct setting, and Eoracle MAE and Eoracle RMSE are the av-

erage MAE and RMSE of the oracle setting.

5.2 Overall Performance

We evaluated the performance under the reconstruct setting (i.e., partial variables available with missing variable reconstruction) with respect to different available variable subset sizes (♣S♣), i.e., 5%, 15%, 25% of total variables. Table 1 presents a comparative analysis of three time series forecasting backbone models under the partial setting (without reconstruction) and the reconstruct setting (with our TRF reconstruction). The results show that our TRF reconstruction of missing variables advance the three backbone forecasting models. In the Appendix D.1, Table 5 details the results under an oracle setting. Notably, when operating under a reconstruct setting with a subset of variables, there is a marked decrease in the performance of the backbone models. However, the integration of our reconstruction method leads to significant enhancements across all backbone models. For instance, in the TRAFFIC dataset, with only 5% variables available, integrating our TRF method can boost the performance of MTGNN by 47.86% (i.e., improving the Mean Absolute Error (MAE) from 22.21 to 11.58). Similarly, for the SOLAR dataset with 5% variables available, ASTGCN's performance is elevated by 29.97% (reducing the MAE from 6.257 to 4.382). These results support the effectiveness of our generative flow as missing variable reconstruction strategy.

5.3 Comparison with Variable Subset Time-Series Forecasting under the Reconstruct Setting

Under the reconstruct setting (i.e., reconstruct missing variables then forecast), we compared our method with the stateof-the-art VSF model, namely the Forecast Distance Weighting (FDW) ensemble method [Chauhan et al., 2022] Table 2 shows that our method outperform FDW, achieving an average performance improvement of over 5%. One possible explanation is that our deep generative flow learning method is capable to model temporal patterns to estimate time series density distribution of missing variables for recovery. This feature is particularly significant as non-parametric methods often fail to effectively recover missing values.

5.4 Comparison with ImputationMethods

Imputation can complete missing data as well. We compared our method with six popular imputation baseline methods, i.e., CSDI [Tashiro et al., 2021], kNNE [Domeniconi and Yan, 2004], MICE [Van Buuren and Groothuis-Oudshoorn, 2011], TRMF [Yu et al., 2016], SAITS [Du et al., 2023] and US-GAN [Miao et al., 2021] with respect to various forecasting backbone models under the reconstruct setting. Due to page limits, we only present the results for MTGNN while other results are detailed in the Appendix D.2. Table 3 shows that our framework can outperform all of the imputation baselines, except the RMSE of traffic dataset. This can be explained as imputation methods are insufficient in modeling inter-variable correlation under time series dynamics.

5.5 Ablation Study

We evaluated the effectiveness of each component of our method on four datasets. We denoted TRF w/o VAML as the variant of our method without Variables-Agnostic Meta Learning. TRF w/o BZ is another variant without batch normalization. Figure 3 shows that removing different components reduces forecasting accuracy for diverse backbones. This illustrates that all components of our method works as one for estimating the unknown density of missing variables and reconstruct the missing relationship. Moreover, we observed that TRF w/o VAML generally perform better than TRF w/o BZ. A potential explanation is that stabilizing the training process is more crucial for recovering missing information compared to emphasizing generalization ability.

5.6 Model Analysis

Horizon Analysis. We evaluated the influence of different horizon lengths on model performance. Figure 5 shows that when dealing with a smaller horizon window, the performance difference between the reconstruct setting and the oracle setting is smaller, indicating a better forecasting accuracy. On the contrary, we observed that under the oracle setting, forecasting accuracy is not influenced by the horizon window size. This is because, all variables are observed and available under the oracle setting.

Analysis of the Number of Coupling Layers. We evaluated the impact of the number of coupling layers. Figure 4 shows that the performance of our proposed framework exhibit limited correlation with the number of coupling layers under T-GCN. In contrast, the performance is sensitive to the number of coupling layers under ASTGCN and MTGNN. For instance, ∆MAE and ∆RMSE decrease when the number of coupling layers increases, while ∆MAE and ∆RMSE become larger when the number of coupling layers increases. This shows there exist an optimal coupling layer number. The potential reason for this trend is that increasing the number of coupling layers in a certain level can improve the performance of framework since flows without deep stack cannot model intricate density, but when the number of coupling layers gets large enough, it adds additional training challenges to impede model convergence.

Analysis of the Number of Hidden Layers. Figure 6 and Figure 7 show how the performance of our method changes over the number of hidden layers in each MADE. We observed that the lager the number of hidden layers is, the poorer the forecasting accuracy is. It is highly possible that more hidden layers will bring more training obstacles and more complex masking transformation for the TRF so that the performances of TRF are greatly hindered.

6 Conclusion

In this work, we reformulate the VSF problem from the perspective of density estimation. To assist the forecasting model in reconstructing the missing values and capturing inter-dependencies among variables, we first estimate the unknown density of inaccessible variable subset, and then reconstruct the missing variables by generation. Specifically, we propose Time-series Reconstruction Flows (TRF), a conditional flows-based generative framework consisting of conditional autoregressive layers based on subset embedding, to estimate the density and rebuild missing variables. To further improve the generalization ability, we design the VariablesAgnostic Meta Learning as the training framework such that TRF can adapt to estimate the density new variables subset more accurately. Finally, we conduct experiments on realworld datasets to indicate our proposed method can recover the most of the performance of forecasting models.

Acknowledgments

This research is funded by the Science and Technology Development Fund (FDCT), Macau SAR (file no. 0123/2023/RIA2, 001/2024/SKL), the Start-up Research Grant of University of Macau (File no. SRG2021-00017IOTSC), and National Science Foundation (NSF) via the grant numbers: 2040950, 2006889, 2045567.

Contribution Statement

In this work, Xuanming Hu and Wei Fan contributed equally by leading the project, offering theoretical support, and handling the overall model design, code implementation, experimental design, and paper writing. Haifeng Chen, Pengyang Wang, and Yanjie Fu gave valuable feedback on the paper drafts. All authors reviewed and approved the manuscript.