Masksembles for Uncertainty Estimation

Deep neural networks have amply demonstrated their prowess but estimating the reliability of their predictions remains challenging. Deep Ensembles are widely considered as being one of the best methods for generating uncertainty estimates but are very expensive to train and evaluate. MC-Dropout is another popular alternative, which is less expensive, but also less reliable. Our central intuition is that there is a continuous spectrum of ensemble-like models of which MC-Dropout and Deep Ensembles are extreme examples. The first uses an effectively infinite number of highly correlated models while the second relies on a finite number of independent models. To combine the benefits of both, we introduce Masksembles. Instead of randomly dropping parts of the network as in MC-dropout, Masksemble relies on a fixed number of binary masks, which are parameterized in a way that allows to change correlations between individual models. Namely, by controlling the overlap between the masks and their density one can choose the optimal configuration for the task at hand. This leads to a simple and easy to implement method with performance on par with Ensembles at a fraction of the cost. We experimentally validate Masksembles on two widely used datasets, CIFAR10 and ImageNet.

Introduction
The ability of deep neural networks to produce useful predictions is now abundantly clear [50,4,44,47,8,9,25,31] but assessing the reliability of these predictions remains a challenge. Among all the methods that can be used to this end, MC-Dropout [12] and Deep Ensembles [32] have emerged as two of the most popular ones. Both of those methods exploit the concept of ensembles to produce uncertainty estimates. MC-Dropout does so implicitly -by train-* This work was supported in part by the Swiss National Science Foundation ing a single stochastic network [14,51,29], where randomness is achieved by dropping different subsets of weights for each observed sample. At test time, one can obtain multiple predictions by running this network multiple times, each time with a different weight configuration. Deep Ensembles, on the other hand, build an explicit ensemble of models, where each model is randomly initialized and trained independently using stochastic gradient descent [46,28]. As importantly, both methods rely on the fact that the outputs of individual models are diverse, which is achieved by introducing stochasticity into the training or testing process.
However, simply adding randomness does not always lead to diverse predictions. In practice, MC-Dropout often performs significantly worse than Deep Ensembles on uncertainty estimation tasks [32,40,17]. We argue that this can be attributed to the fact that each weight is dropped randomly and independently from the others. For many configurations of hyperparameters, in particular the dropout rate, this results in similar weight configurations and, consequently, less diverse predictions [10]. Deep Ensembles seem to not share this weakness for reasons that are not yet fully understood and usually produce significantly more reliable uncertainty estimates. Unfortunately, this comes at a price, both at training and inference time: Building an ensemble requires training multiple models, and using it means loading all of them simultaneously in memory, typically a very scarce resource.
In this work, we introduce Masksembles, an approach to uncertainty estimation that tackles these challenges and produces reliable uncertainty estimates on par with Deep Ensembles at a significantly lower computational cost. The main idea behind the method is simple -to introduce a more structured way to drop model parameters than that of MC-Dropout.
Masksembles produces a fixed number of binary masks which specify the network parameters to be dropped. The properties of the masks effectively define the properties of the final ensemble in terms of capacity and correlations.
During training, for each sample, we randomly choose one of the masks and drop the corresponding parts of the model just like a standard dropout. During inference, we run the model multiple times, once per mask, to obtain a set of predictions and, ultimately, an uncertainty estimate. Our method has three key hyperparameters: the total number of masks, the average overlap between masks and the number of ones and zeros in each mask. Intuitively, using a large number of masks approximates MC-Dropout, while using a set of non-overlapping, completely disjoint masks yields an Ensemble-like behavior. In other words, as illustrated in Fig. 5 Masksembles defines a spectrum of model configurations of which MC-Dropout and Deep Ensembles can be seen as extreme cases. We evaluate our method on several synthetic and real datasets and demonstrate that it outperforms MC-Dropout in terms of accuracy, calibration, and out-of-distribution (OOD) detection performance. When compared to Deep Ensembles, our method has similar performance but is more favorable in terms of training time and memory requirements. Furthermore, Masksembles is simple to implement and can serve as a drop-in replacement for MC-Dropout in virtually any model. To summarize, the main contributions of this work are as follows: • We propose an easy to implement, drop-in replacement method that performs better than MC-Dropout, at a similar computational cost, and that matches Ensembles at a fraction of the cost.
• We provide theoretical insight into two popular uncertainty estimation methods -MC-Dropout and Ensembles, by creating a continuum between the two.
• We provide a comprehensive evaluation of our method on several public datasets. This validates our claims and guides the choice of Masksemble parameters in practical applications.

Related Work
MC-Dropout [12] and Deep Ensembles [32] have emerged as two of the most prominent and practical uncertainty estimation methods for deep neural networks [1,40,17]. In what follows, we provide a brief background on uncertainty estimation, review the best-known methods, and discuss potential alternatives [53,20,2,37,36].

Background
The goal of Uncertainty Estimation (UE) is to produce a measure of confidence for model predictions. Two major types of uncertainty can be modeled within the Bayesian framework [7]. Aleatoric uncertainty, also known as Data uncertainty [11], captures noise that arises from data and therefore is irreducible, because the natural complexity of the data directly causes it. Sensor noise, labels noise, and classes overlapping are among the most common sources of aleatoric uncertainty. Epistemic uncertainty or model uncertainty [11] accounts for uncertainty in the parameters of the model that we learned, that is, it represents our ignorance about parameters of the data generation process. This type of uncertainty is reducible and therefore, given enough data, could be eliminated entirely. In addition to the aforementioned types of uncertainty, there is also Distributional uncertainty. Also known as dataset shift [43], this type of uncertainty is useful when there is a mismatch between the training and testing data distributions, and a model is confronted with unfamiliar samples.
For all the aforementioned types of uncertainty, ideally we want our model to output an additional signal indicating how reliable its prediction is for a particular sample. More formally, given an input, we want our model to output the actual prediction target as well as the corresponding uncertainty measure. For regression tasks, the uncertainty measure can be the variance or confidence intervals around the prediction made by the network [41]. For classification tasks, popular uncertainty measures are the entropy and max-probability [37,33]. Ideally, the measure of uncertainty should be high when the model is incapable of producing accurate prediction for the given input and low otherwise. Generally speaking, uncertainty can be either learned from the data directly [33,27], in particular for big data regimes and aleatoric uncertainty, or could be acquired from a diverse set of predictions [12,33], which are obtained from stochastic regularization techniques or ensembles. In this work, we focus on the latter family of methods.

Ensembles
Deep Ensembles [32] involve training an ensemble of deep neural networks on the same data by initializing each network randomly. This approach is originally inspired by a classical ensembling method, bagging [19]. An uncertainty estimate is obtained by running all the elements in the ensemble and aggregating their predictions. The major drawback of Deep Ensembles is the computational overhead: it requires training multiple independent networks, and during inference it is desirable to keep all these networks in memory. More generally, ensembling has been a popular way to boost the performance of individual neural network models [18,30,33,42,54,56], with the quality improvements often being associated with diverse predictions. In [10] authors provide a comprehensive comparison of several popular ensemble learning methods as well as several approximate Bayesian inference techniques. Extensive experiments on vision tasks suggest that simple ensembles have the lowest correlation between individual models. This is an important advantage of these methods, as ultimately this leads to highly diversified predictions on the samples which are very distinct from those observed in training.
Snapshot-based methods take a slightly different approach towards creating an ensemble, and try to collect a diverse set of weight configurations over the course of a single training. Snapshot ensembles [24] build upon the idea of using cyclical learning rates to ensure that the optimization process can efficiently explore multiple local optima. Similarly, Fast Geometric Ensembling [13,26], uses a method to identify low-error paths between local optima to obtain a collection of diverse representations. Although these approaches partially tackle training time issues, they typically perform worse than the standard ensembles [1] and consume the same amount of memory. Furthermore, because of the specifics of their training procedure, in particular, shared initialization with a pre-trained network, individual snapshots are not guaranteed to be decorrelated, which might hurt the performance [1].

Dropout
Dropout [48] is an extremely popular and simple approach to stochastically regularizing deep neural networks. During training, neurons are randomly dropped with fixed probability, which effectively produces an approximation of several slightly different model architectures and allows a single model to mimic ensemble behavior. Originally Dropout was used during training as a stochastic regularization technique, but [12,11] proposed using it during inference, providing sound mathematical justification for ensembling and Bayesian model averaging analogies. MC-Dropout is extremely easy to use and has zero memory overhead compared to a single model but shows consistently worse performance compared to Deep Ensembles [40,17]. In particular, different predictions made by several forward passes with randomly generated masks seem to be overly correlated and strongly underestimate the variance. In this paper, we postulate that the training procedure, which makes any combination of the weights possible, will tend to create uniformity between predictions. As opposed to Ensemble techniques, MC-Dropout techniques do not let multiple versions of the network actually be created, and the weights are forced to converge jointly.

Other Methods
Although Bayesian inference provides a natural way to assess uncertainty of model predictions, it is prohibitively expensive to apply these techniques to deep neural networks directly. A number of approximate inference techniques were developed for Bayesian Neural Networks (BNN) [34,38,35] that try to alleviate these problems, to name a few: Bayes by Backprop [3] uses Variational Inference to learn a parametric distribution over model weights that approximates the true posterior, Laplace Approximation [45] tries to find a Gaussian approximation to the posterior modes, MC-Dropout [12] could be seen as an approximate Bayesian procedure as well. Although these techniques are theoretically grounded, Deep Ensembles often show significantly better performance in practice [40,1], in terms of both accuracy and quality of the resulting uncertainty estimates.

Method
In this section, we introduce Masksembles, an approach to uncertainty estimation that creates a continuum between single networks, deep ensembles, and MC-Dropout. It builds on the intuition that dropout-based methods can be seen as stochastic variants of ensembles [33], albeit ensembles with infinite cardinality. We will argue that the main reason for MC-Dropout's poor performance is the high correlation between the ensemble elements that makes the overall predictions insufficiently diverse. Masksembles is designed to give control over the correlation between the elements of an ensemble and, hence, achieve a satisfactory compromise between reliable uncertainty estimation and acceptable computational cost.

Masksembles as Structured Dropout
Both Ensembles and MC-Dropout are ensemble methods, that is, they produce multiple predictions for a given input and then use an aggregated measure such as variance or entropy as an uncertainty estimator. Formally, consider a dataset {x i , y i } i , where x i ∈ R F are input features, and y i are scalars or labels. For classification problems, we model conditional distribution p(y|x; θ) as a mixture where θ k are the weights of element k of the ensemble, and N is the size of the ensemble. For Deep Ensembles, individual models are independent and do not share any weights, and hence each θ k is a separate set of weights, trained entirely independently. By contrast, for MC-Dropout there is a single shared θ, and individual models are obtained as θ k =b t k · θ, wherẽ b t k ∈ {0, 1} |θ| are random binary masks which are resampled for each iteration t. The elements ofb t k follow Bernoulli distribution parameterized by dropout probability p, and the parameters corresponding to the same network activations are dropped simultaneously. At inference time, a small number of predictions is produced by randomly sampling the masks. In practice, for any reasonable choice of p, the masksb t k will overlap significantly, which will tend to make predictions p(y|x; θ k ) highly correlated, and thus could lead to underestimated uncertainty. Furthermore, because masks are re-sampled at each training iteration, each network unit needs to be able to form a coherent response with any other unit, which creates a mixing effect and leads to uniformity between the predictions from different masks.
To tackle these issues, we propose using a finite set of pre-defined binary masks, which are generated so that their overlap can be controlled. They are then used to drop corresponding network activations so that the resulting models are sufficiently decorrelated and do not suffer from the mixing effect described above. In what follows, we describe the algorithm for mask generation, explain the way we apply generated masks to models, and discuss the trade-offs it incurs.

Generating Masks
Generating masks with a controllable amount of overlap is central to Masksembles. Let us first introduce the key parameters of our mask generation process: • N -number of masks.
• M -number of ones in each mask.
• S -scale that controls the amount of overlap given N and M.
The masks generation algorithm is summarized below. It starts by generating N vectors of all zeros of size M × S. Then, in each of those vectors, it randomly sets M of these elements to be 1 in each vector. As depicted by Fig. 1, it then looks for positions that are all zero in the resulting vectors. These correspond to features that will be dropped. For S = 1.0, the algorithm then generates N masks of size M that will contain only ones and therefore fully overlapping. Conversely, setting S to large positive values will generate N masks with mostly zero values and therefore very few overlapping one values.

Inserting Masksembles Layer into Networks
We start with the very simple case depicted by Fig. 2. It shows an MLP with one hidden layer of size 3 and 2dimensional inputs and outputs. We insert Masksembles layer right after the hidden layer and increase its size to M × S − D where D is the number of dropped features so that the sizes of the masks and the size of the layer match. Then at run-time, we apply the mask values and ignore all components of the hidden layer for which they are zero. Fig. 2 features 3 different configurations obtained by fixing both M and N -here, N is sufficiently large so that D = 0 -and changing the value of S to 1.0, 1.7, and 2.3, therefore increasing masks sizes and hidden size in range [3,5,7].
Applying such modification to the model creates a larger one that allows for multiple relatively uncorrelated predictions. This process extends naturally to MLPs of any dimension. It also extends to convolutional architectures without any significant modification. As in MC-dropout, we  can drop entire channels [49] instead of individual activations. In other words, treating the whole channel as a single feature is equivalent to the fully connected case discussed above.
In general, this means that we may have to adapt the size of network layers each time we change M, N, or S, which can be cumbersome. Fortunately, keeping S × M equal to the original channels number, 3 in the MLP example, makes it possible to apply Masksembles without any such change of the network configuration.

Parameters and Model Transitions
The chosen parameters N, M and S define several key model properties: • Total model size: The total number of model parameters. In the above MLP example, this is equal to where D is the number of features trimmed during the mask generation procedure. This functional relationship is illustrated in Fig. 3.
• Model capacity: The number of activations that are used during one forward pass of the model, which is equal to M.
• Masks IoU average Intersection over Union (IoU) between generated masks. We derive approximation for this quantity that appears to depend only on S : IoU(S) = 1/(2S − 1) and therefore N, M do not influence it, as show in Fig. 3. In fact, Deep Ensembles and MC-Dropout can be both seen as extreme cases and Masksembles can transition smoothly from one to the other (Fig. 5).
Ensembles: When S goes to infinity, Masks IoU goes to zero and Masksembles starts behaving like Ensembles that uses non-overlapping masks.
Dropout: As N increases, each individual mask configuration becomes less and less likely to be picked during training. In the limit where N goes to infinity, we reproduce the situation where each mask is seen only once, which is equivalent to MC-Dropout. In this case, the dropout rate would be 1 − 1 S . Transition: In Fig. 14, we use a simple classification problem to illustrate the span of behaviors Masksembles can cover. Specifically, we show that reducing mask overlap, that is, decreasing the correlation between binary masks, yields progressively more diverse predictions on out-of-training examples and generates more ensemble-like behavior.

Implementation details
The only extra computation performed in Masksembles over a single network is cheap tensor-to-mask multiplication. Therefore, compared to convolutional and fully connected layers, this product induces negligible computational overhead. Furthermore, we rely on similar optimizations as in [52] to make our implementation more efficient.

Experiments
In this section, we first introduce the evaluation metrics which we use to quantitatively compare Masksembles against MC-Dropout and Ensembles. We then compare our method with the baselines on two broadly used benchmark datasets CIFAR10 [30] and ImageNet [6].
Note that we did not use any complex data augmentation schemes, such as AugMix [23], Mixup [57] or adver-sarial training [15], in order to avoid entangling their effects with those of Masksembles. Ultimately, these techniques are complementary to ours and could easily be used together.

Metrics
Given a classifier with parameters θ, let p(y|x; θ) be the predicted probability for a sample represented by features x to have a label y, and let p gt (y|x) denote the true conditional distribution over the labels given the features. We use several different metrics to analyze how close p(y|x; θ) is to p gt (y|x) and to compare our method to others, in terms of both accuracy and quality of uncertainty estimation.
Accuracy. For classification tasks, we use the standard classification accuracy, i.e. the percentage of correctly classified samples on the test set.
Entropy (ENT). It is one of the most popular measures used to quantify uncertainty [33,37]. For classification tasks it is defined as where N c is the number of classes.
Expected Calibration Error (ECE). ECE quantifies the quality of uncertainty estimates specifically for in-domain uncertainty. A model is considered well-calibrated if its predicted distribution over labels p ( y|x; θ) is close to the true distribution p gt (y|x). As in [16], we define ECE to be the average discrepancy between the predicted and real data distribution: which is a common metric for calibration evaluation.
Out-of-Distribution Detection (OOD ROC / PR). Domain shift occurs when the features in the training data do not follow the same probability distribution as those in the test distribution. Detecting out-of-distribution samples is an important task, and uncertainty estimation can be used for this purpose under the assumption that our model returns higher uncertainty estimates for such samples. To this end, we use the uncertainty measure produced by our models as a classification score that determines whether a sample is in-domain or not. We then apply standard detection metrics, ROC and PR AUCs, to quantify the performance of our model.
Model Size This corresponds to the Total model size defined in Section 3.4. A major motivation for Masksembles is to match the performance of Deep Ensembles while using a smaller model that requires significantly less memory. We use a total number of weights that parameterize our models as a proxy for that. In addition to the model size, we also report the total GPU time used to train any particular model.

CIFAR10
CIFAR10 is among the most popular benchmarks for uncertainty estimation. The dataset is relatively small-scale, and very large models tend to overfit the data easily. We, therefore, picked the Wide-ResNet-16-4 [55] for our experiments and used training procedures similar to the original paper. Following [12], we place dropout layers before all the convolutional and fully connected layers in the MC-Dropout model. Masksembles layers are added exactly the same way.
Accuracy and ECE. One of the standard ways to assess the robustness of a model is to evaluate the dependency of its performance with respect to inputs perturbations [22,40]. In our experiments, we investigate the behaviour of all the considered models under the influence of domain shifts induced by typical sources of noise.
Following [52,40], we evaluate model accuracy and ECE on a corrupted version of CIFAR10 [22]. Namely, we consider 19 different ways to artificially corrupted the images and 5 different levels of severity for each of those corruptions.
We train all the models on the original uncorrupted CI-FAR10 training set and test them both on the original images and their corrupted versions. In Table 1, we report our results on the original uncorrupted images and in Fig. 6 the results on the noisy ones. Each box represents the first and third quantiles of accuracy and ECE while the error bars denote the minimum and maximum values. We tested four different approaches, that is, a single network, MC-Dropout, Ensembles, and our Masksembles. Unsurprisingly, as the severity of the perturbations increases, using a single network is the least robust approach and degrades fastest. Masksembles performs on par with Ensembles and consistently outperforms MC-dropout, even though we extensively tuned MC-Dropout for the best possible performance by doing a grid search on the dropout rate (which turned out to be p ≈ 0.1).
In these experiments, Masksembles was trained with a fixed number of N = 4 masks. We varied both S in the range [1.0, 5.0] and M in the range [0.3C, 2.0C], where C is the number of channels in the original model. Fig. 7 depicts the resulting range of behaviors. The 2D coordinates of the markers depict their accuracy and ECE, their colorsthe corresponding model size, and the size of the star markers denotes the value of the mask-overlap. For comparison purposes, we also display MC-dropout and Ensembles results in a similar manner, simply replacing the star with a square and a circle, respectively. As can be seen, the optimal Masksemble configuration depicted by the green star in the lower right corner delivers performances similar to those of Ensembles for a model size that is almost half the size. This is the one we also used in the experiments depicted by Fig. 6.
In Fig. 8, we chose the Masksembles parameters so that the model size always remains the same: We fix the number of masks N = 4 and vary S and M so that the total model size remains constant, which means that a lower mask overlap results in a lower capacity for each individual model. This procedure is important in practice because it enables us to take a large network, such as ResNet, and add Masksembles layers as we would add Dropout layer without having to change the rest of the architecture.
For this set of experiments we obtain a consistent improvement in calibration quality by reducing mask overlap and achieve ensembles-like performance with S = 5.
OOD ROC and PR. For OOD task, we followed the evaluation protocol of [5]. Namely, we trained our models on CIFAR10, and then used CIFAR10 test images as  our in-distribution samples and images from the SVHN dataset [39] (which belong to classes that are not in CI-FAR10), as our out-of-distribution samples. We report our OOD results in the two bottom rows of Table 1. Again, performance of our method is similar to that of Ensembles at a fraction of the computational cost, and significantly better than that of a single network and MC-dropout.

ImageNet
The ImageNet [6] is another dataset that is widely used to assess the ability of neural networks to produce uncertainty estimates. For this dataset, we rely on the well-known ResNet-50 [21] architecture as a base model, using the setup similar to that in [1]. Similarly to our CIFAR10 experiments, for the MC-Dropout and Masksembles models, we introduce dropout (masking) layers before every convolutional layer starting from stage3 layers.
We followed exactly the same evaluation protocol as in Section 4.2 and report our results on the original images in Table 2 and on the corrupted ones in Fig. 11. In Table 2 and Fig. 10 we use the same procedure as for CIFAR10 to select Masksembles parameters so that the total model size remains fixed.
For this dataset, Masksembles demonstrates its best performance for a choice of parameters corresponding to a Masks IoU of 0.5. Similarly to our results on CIFAR10, the performance of Masksembles both in terms of accuracy and ECE is very close to that of Ensembles and significantly better than that of MC-Dropout. Note that our method achieves these results with a training time and memory consumption that is 4 (four) times smaller than that of Ensembles, and nearly the same as that of a single network. In terms of calibration quality on corrupted images, we achieve Ensemblelike performance starting from S = 2 and up.

Conclusion
In this work, we introduced Masksembles, a new approach to uncertainty estimation in deep neural networks. Instead of using a fixed number of independently trained models as in Ensembles or drawing a new set of random binary masks at each training step as in MC-Dropout, Masksembles predefines a pool of binary masks with a con-   trollable overlap and stochastically iterates through them. By changing the parameters that control the mask generation process, we can span a range of behaviors between those of MC-Dropout and Ensembles. This allows us to identify model configurations that provide a useful tradeoff between the high-quality uncertainty estimates of Deep Ensembles at a high computational cost and the lower performance of MC-Dropout at a lower computational cost. In fact, our experiments demonstrate that we can achieve the performance on par with that of Deep Ensembles at a fraction of the cost. Because Masksembles is easy to implement, we believe it can serve as a drop-in replacement for MC-Dropout. In future work, we will explore ways to accelerate Masksembles by exploiting the structure of the masks, and apply our method to tasks that require a scalable uncertainty estimation method, in particular reinforcement learning and Bayesian optimization.

A. Supplementary Material
In this appendix, we study in more detail the relationship between the choice for N, M, and S parameters of Section 3.2 that control the mask generation process and the resulting model size, mask correlation, and models' diversity.

A.1. Expected models size
We first provide a formal derivation for the expected size of a generated model given the N, M, S parameters. Let us consider N vectors of size M × S filled with zeros, and then in each of those vectors, we randomly set M of these elements to be 1. Let ζ ij be a random variable that represents whether or not j-th position in i-th vector is equal to one. such configurations where j-th position is fixed to be one, therefore ζ ij > 0] be the random variable representing whether at least one 1 appears in the j-th position among all generated vectors. Given that P(ζ ij = 0) = 1− 1 S and P(X) = 1−P(X), then probability of j can be written as To compute the expected model size, we compute the expectation of j sum, since effectively it represents expected number of features that one will acquire after a trimming procedure:

A.2. Average IoU
We now justify our 1 2S−1 approximation for the average mask IoU. Let us consider two vectors produced by our mask generation algorithm. As above, let ζ 1j and ζ 2j be random variables that represent the value in the j-th position in the first and second masks, respectively. Starting from the standard definition of the IoU, we estimate the average intersection I of two such masks as the number of common ones: Given two masks with M ones each and intersection I, their IoU is I 2M−I . Therefore, a simple approximation for expected IoU is In Fig. 12, we plot this value as a function of S along with empirical values and the agreement is excellent.

A.3. Diversity analysis
It is known that less correlated ensembles of models deliver better performance, produce more accurate predictions [33,42,18], and demonstrate lower calibration error [40]. Hence, a strength of Masksembles is that it provides the ability to control how correlated models within Masksembles are by adjusting the N, M, and S parameters.
We now perform a diversity analysis of Masksembles models using the metric of [10]. It involves comparing two different models trained on the same data in terms of how different their predictions are. Measuring fraction of the test data points on which predictions of models disagree, the diversity, and normalizing it by the models error rate, we write Diversity = fraction of disagreed labels 1.0 − accuracy Plotting this diversity against accuracy allows us to look at our models from a bias-variance trade-off perspective. Since we want our models to be less correlated-larger diversityand at the same time to be accurate, the upperright corner of Fig. 13 is where the best models should be. For this experiment we used a single trained model; several Masksembles models with N = 4, fixed M to have the same capacity as the single one and varying S in [2,3,4,5]; ensembles model with 4 members. Fig. 13 depicts the results on CIFAR10 dataset. Ensembles have the largest diversity but Masksembles gives us the ability to achieve very similar results by controlling its parameters. The single model has 0 diversity by the definition.

A.4. Calibration Plots
Since Expected Calibration Error (ECE) provides only limited and aggregated information about model's calibration therefore in this section we present full calibration plots for experiments performed on CIFAR and ImageNet datasets in experiments sections.  For both datasets, the calibration plots support our claim that lower model correlation yields better calibration. Furthermore, reducing masks overlap for Masksembles enables us to match Ensembles behavior very closely.


Comparative Analysis of Recurrent Neural Network Architectures for Reservoir Inflow Forecasting

Due to the stochastic nature and complexity of flow, as well as the existence of hydrological uncertainties, predicting streamflow in dam reservoirs, especially in semi-arid and arid areas, is essential for the optimal and timely use of surface water resources. In this research, daily streamflow to the Ermenek hydroelectric dam reservoir located in Turkey is simulated using deep recurrent neural network (RNN) architectures, including bidirectional long short-term memory (Bi-LSTM), gated recurrent unit (GRU), long short-term memory (LSTM), and simple recurrent neural networks (simple RNN). For this purpose, daily observational flow data are used during the period 2012–2018, and all models are coded in Python software programming language. Only delays of streamflow time series are used as the input of models. Then, based on the correlation coefficient (CC), mean absolute error (MAE), root mean square error (RMSE), and Nash–Sutcliffe efficiency coefficient (NS), results of deep-learning architectures are compared with one another and with an artificial neural network (ANN) with two hidden layers. Results indicate that the accuracy of deep-learning RNN methods are better and more accurate than ANN. Among methods used in deep learning, the LSTM method has the best accuracy, namely, the simulated streamflow to the dam reservoir with 90% accuracy in the training stage and 87% accuracy in the testing stage. However, the accuracies of ANN in training and testing stages are 86% and 85%, respectively. Considering that the Ermenek Dam is used for hydroelectric purposes and energy production, modeling inflow in the most realistic way may lead to an increase in energy production and income by optimizing water management. Hence, multi-percentage improvements can be extremely useful. According to results, deep-learning methods of RNNs can be used for estimating streamflow to the Ermenek Dam reservoir due to their accuracy.

Introduction
Large dam structures are built to reserve water for different supply objectives, such as drinking water irrigation, hydroelectric generation, and flood control. Generating clean energy, especially in developing countries, is a major focus of managers of hydroelectric dams. Accurate predictions of daily several clusters, and used the deep network of SAE-BP (integrating stacked autoencoders with BP neural networks) for modeling. Then, they compared the results with other methods, such as SVM. According to the results, the deep-learning algorithm presented a relatively better performance in flood estimation. Esmaeilzadeh et al. [14] estimated daily flow to the Sattarkhan Dam in Iran by using data-mining methods, including ANNs, the M5 tree method, support vector regression, and hybrid Wavelet-ANN methods. Results indicated that the wavelet artificial neural network (WANN) method estimated flow better than other methods. Amnatsan et al. [15] used the variation analog method (VAM), WANN, and weighted-mean analog methods to estimate input flow to the Sirikit Dam in Thailand. Results revealed that the VAM approach was better than other methods for high-intensity flows. Chiang et al. [16] integrated ensemble techniques into artificial neural networks to reduce model uncertainty in hourly streamflow predictions in Longquan Creek and Jinhua River watersheds in China. Results demonstrated that the ensemble neural networks improved about 19-37% of the accuracy of streamflow predictions as compared to a single neural network. Zhang et al. [17] used ANN, SVM, and long short-term memory (LSTM) methods to operate the Gezhouba Dam reservoir in China in hourly, daily, and monthly time scales. Results showed that the LSTM method, with less computational time, was more accurate in peak periods, and the number of maximum iterations was effective in model performance. Zhou et al. [18] used three ANN architectures (a radial basis function network, an extreme learning machine, and the Elman network) for monthly streamflow forecasting in the Jinsha River Basin in China. The best estimate was made by Elman architectures with an R of 0.906. Kratzert et al. [19] used the LSTM method with a freely available data set for rainfall-runoff modeling. Results showed that this method works well and could be used in hydrological modeling. Zhang et al. [20] used three deep learning algorithms (recurrent neural network (RNN), long short-term memory (LSTM), and gated recurrent unit (GRU)) to predict outflows in the Xiluodu reservoir. They determined that all three models can predict reservoir outflows accurately and efficiently and could be used to control floods and generate power; the number of iterations and hidden nodes mainly influenced the model precision. Kao et al. [21] used the LSTM-based Encoder-Decoder (LSTM-ED) model for multi-step-ahead flood forecasting in the Shihmen Reservoir catchment in Taiwan. Results showed that the proposed model that translated and linked the rainfall sequence with the runoff sequence could improve the reliability of flood forecasting and increase the interpretability of model internals. Zhou et al. [22] adopted an unscented Kalman filter (UKF) post-processing technique to forecast point flood by RNN in the three Gorges Reservoir in China. Results indicated that the proposed method extracted the complex non-linear dependence structure between the model's outputs and observed inflows and overcame the systematic error.
Fossil fuels throughout the world but also in Turkey are inadequate and electricity production from fossil fuels is expensive. Therefore, there is a potential importance of hydroelectric power plants for the production of electricity and the country relies on them for stabilizing the economic conditions. Inflow to the reservoir at the hydroelectric power plant is converted into much cheaper electricity. Increasing prediction accuracy can clearly lead to more useful planning for timely production and increase its generation to further continue improving profits and the efficiency of the dam. However, streamflows have a totally random and stochastic structure. Deep learning can be utilized for the evaluation of the complex and nonlinear association between streamflows, basin properties, and meteorological variables. In a stream, the streamflows are a time series with a seasonal recurrence, and the RNN is a technique that is implemented efficiently under recurrent conditions. RNN architectures have not been sufficiently examined and are not compared in streamflows. Given the significance of flow prediction in hydroelectric dam power production, this study attempts to predict daily inflow to the Ermenek Dam reservoir located in Turkey by utilizing deep learning techniques. Streamflow is precisely a time series, and RNN is the most suitable deep learning approach for the time series. Therefore, this research uses four RNN architectures, including Bi-LSTM, gated recurrent unit (GRU), LSTM, and a simple recurrent neural network (simple RNN). The efficiency of deep learning techniques in the flow prediction is assessed, and the results acquired from these techniques are compared with the ANN technique. The deep-learning focused method suggested in this study can help specialist engineers make more accurate and effective predictions. This proposed research may be considered as a part of creative prediction approaches and stochastic views in hydrology.

Study Area and Data
In the study, data measured in the Ermenek Dam, which is within the boundaries of the Ermenek District of Karaman Province, were used. The district is located between 36 • 58 N and 32 • 53 E. The average height of the district from sea level is 1250 m. The town of Ermenek contains many river sources, plateaus, promenades, as well as sources of historical and natural beauty. The Ermenek River, one of the crucial streams in the region, collects almost all waters of the region. Ermenek District, with an area of 122,297 ha, covers 13.86% of Karaman and 0.16% of Turkey. Of the total area, 21% of the district is farmland, 10% is meadow-pasture, and 38% is forest area [23][24][25][26].
Data used in the research belong to the Ermenek Dam, which was opened in 2009 with an installed capacity of 308.88 MW and an annual electricity generation capacity of 1187 GWh. The body of the dam is a concrete arch, and the crest length is 123 m. The filling volume is 305,000 m 3 , and the storage volume is 4582 billion m 3 [27].
To simulate streamflow to the Ermenek Dam ( Figure 1), daily mean inflow data are used during the 2012-2018 period. The basic characteristics of the data are provided in Table 1, and the time-series plots are displayed in Figure 2. This proposed research may be considered as a part of creative prediction approaches and stochastic views in hydrology.

Study Area and Data
In the study, data measured in the Ermenek Dam, which is within the boundaries of the Ermenek District of Karaman Province, were used. The district is located between 36°58′N and 32°53′E. The average height of the district from sea level is 1250 m. The town of Ermenek contains many river sources, plateaus, promenades, as well as sources of historical and natural beauty. The Ermenek River, one of the crucial streams in the region, collects almost all waters of the region. Ermenek District, with an area of 122,297 ha, covers 13.86% of Karaman and 0.16% of Turkey. Of the total area, 21% of the district is farmland, 10% is meadow-pasture, and 38% is forest area [23][24][25][26].
Data used in the research belong to the Ermenek Dam, which was opened in 2009 with an installed capacity of 308.88 MW and an annual electricity generation capacity of 1187 GWh. The body of the dam is a concrete arch, and the crest length is 123 m. The filling volume is 305,000 m 3 , and the storage volume is 4582 billion m 3 [27].
To simulate streamflow to the Ermenek Dam ( Figure 1), daily mean inflow data are used during the 2012-2018 period. The basic characteristics of the data are provided in Table 1, and the time-series plots are displayed in Figure 2.     For modeling, ANN methods and four different RNN deep-network methods were employed. In total, 70% and 30% of the data were used in training and testing stages, respectively, as indicated in Figure 2. Besides this, several delayed daily mean streamflow are used as input data. For this purpose, the correlation between the time series of streamflows with its delays is obtained, and seven delays (SFt-1, SFt-2, SFt-3, SFt-4, SFt-5, SFt-6, and SFt-7) are selected as inputs due to high correlation ( Figure  3). The correlation coefficients are found to vary between 0.63 and 0.85 in seven-day lag conditions. Similar studies [7,8] have also used seven delays in flow time series as inputs. For modeling, ANN methods and four different RNN deep-network methods were employed. In total, 70% and 30% of the data were used in training and testing stages, respectively, as indicated in Figure 2. Besides this, several delayed daily mean streamflow are used as input data. For this purpose, the correlation between the time series of streamflows with its delays is obtained, and seven delays (SF t-1 , SF t-2 , SF t-3 , SF t-4 , SF t-5 , SF t-6 , and SF t-7 ) are selected as inputs due to high correlation ( Figure 3). The correlation coefficients are found to vary between 0.63 and 0.85 in seven-day lag conditions. Similar studies [7,8] have also used seven delays in flow time series as inputs. For modeling, ANN methods and four different RNN deep-network methods were employed. In total, 70% and 30% of the data were used in training and testing stages, respectively, as indicated in Figure 2. Besides this, several delayed daily mean streamflow are used as input data. For this purpose, the correlation between the time series of streamflows with its delays is obtained, and seven delays (SFt-1, SFt-2, SFt-3, SFt-4, SFt-5, SFt-6, and SFt-7) are selected as inputs due to high correlation ( Figure  3). The correlation coefficients are found to vary between 0.63 and 0.85 in seven-day lag conditions. Similar studies [7,8] have also used seven delays in flow time series as inputs.

Artificial Neural Network (ANN)
An ANN is a distributed knowledge treatment system in which performance essentials are alike to the human brain, and is based on a simulated biological neural network [28]. Each neural network has three layers, namely, input, hidden, and output. The input layer is a layer for providing data provided as inputs to the network. The output layer contains values predicted by the network. The hidden layer is the data analysis location. Usually, the number of selected neurons of the layers is obtained by trial and error. The general architecture of the ANN is displayed in Figure 4, where X (x 1, x 2, ..., x n) = inputs vector, W = connecting weights to the next layer, b k = bias, and y k is the ANN final output. The activation function converts input signals into output signals.
Water 2020, 12, x FOR PEER REVIEW 6 of 18 Figure 3. Correlation plot between streamflow and its delays.

Artificial Neural Network (ANN)
An ANN is a distributed knowledge treatment system in which performance essentials are alike to the human brain, and is based on a simulated biological neural network [28]. Each neural network has three layers, namely, input, hidden, and output. The input layer is a layer for providing data provided as inputs to the network. The output layer contains values predicted by the network. The hidden layer is the data analysis location. Usually, the number of selected neurons of the layers is obtained by trial and error. The general architecture of the ANN is displayed in Figure 4, where X (x1, x2, ..., xn) = inputs vector, W = connecting weights to the next layer, bk = bias, and yk is the ANN final output. The activation function converts input signals into output signals. In Figure 4, N inputs are given from x1 to xn to the counterpart weights Wk1 to Wkn. Initially, the weights are multiplied by their inputs, and then they are summed with the amount of bias to obtain u (Equation (1)): Then, the activation function is adapted on u, i.e., f (u); ultimately, the final output value is obtained as = f(u) of the neuron. The most popular activation functions are Sigmoid, ReLU, and Softmax. In this study, feed-forward neural networks were used.

Recurrent Neural Networks (RNN)
Deep-learning algorithms are a sample of machine-learning algorithms where the purpose is to discover multiple levels of representation of input data. Developed in the 1980s, multilayer RNNs are among the most commonly used models for deep learning [30]. These types of networks have a memory that records the information they have seen so far, and many types exist. Moreover, RNNs are powerful models for sequential data (time series) [31], and they use the previous output to predict the next output. In this case, the networks themselves have repetitive loops. These loops, which are in the hidden neurons, allow the storing of previous input information for a while so that the system can predict future outputs. The hidden layer output is retransmitted t times to the hidden layer. The output of a recursive neuron is only sent to the next layer when the number of iterations is completed. In this case, the output is more comprehensive, and the previous information is kept for longer. Finally, the errors are returned backward to update the weights. In this study, four available RNN architectures are used. In order to improve the readability of this study, a research flow chart is given in Figure 5. In Figure 4, N inputs are given from x 1 to x n to the counterpart weights W k1 to W kn . Initially, the weights are multiplied by their inputs, and then they are summed with the amount of bias to obtain u (Equation (1) Then, the activation function is adapted on u, i.e., f (u); ultimately, the final output value is obtained as y k = f(u) of the neuron. The most popular activation functions are Sigmoid, ReLU, and Softmax. In this study, feed-forward neural networks were used.

Recurrent Neural Networks (RNN)
Deep-learning algorithms are a sample of machine-learning algorithms where the purpose is to discover multiple levels of representation of input data. Developed in the 1980s, multilayer RNNs are among the most commonly used models for deep learning [30]. These types of networks have a memory that records the information they have seen so far, and many types exist. Moreover, RNNs are powerful models for sequential data (time series) [31], and they use the previous output to predict the next output. In this case, the networks themselves have repetitive loops. These loops, which are in the hidden neurons, allow the storing of previous input information for a while so that the system can predict future outputs. The hidden layer output is retransmitted t times to the hidden layer. The output of a recursive neuron is only sent to the next layer when the number of iterations is completed. In this case, the output is more comprehensive, and the previous information is kept for longer. Finally, the errors are returned backward to update the weights. In this study, four available RNN architectures are used. In order to improve the readability of this study, a research flow chart is given in Figure 5.

Simple Recurrent Neural Network (Simple RNN)
A simple RNN is essentially a collection of common neural networks arranged together, each of them transmitting a message to another. In other words, these networks have a memory that stores knowledge about the data seen, but their memory is short term and cannot maintain long-term time series [32]. Figure 6a displays a simple RNN. A simple recurrent network has only one internal memory-ℎ -which is computed from Equation (2): where g() denotes an activation function, U and W are flexible weight matrices of the h layer, b is a bias, and X is an input vector [19].

Simple Recurrent Neural Network (Simple RNN)
A simple RNN is essentially a collection of common neural networks arranged together, each of them transmitting a message to another. In other words, these networks have a memory that stores knowledge about the data seen, but their memory is short term and cannot maintain long-term time series [32]. Figure 6a displays a simple RNN. A simple recurrent network has only one internal memory-h t -which is computed from Equation (2): where g() denotes an activation function, U and W are flexible weight matrices of the h layer, b is a bias, and X is an input vector [19].

Long Short-Term Memory (LSTM)
LSTM is a kind of model or structure for sequential data developed by [33] for the advancement of RNN. It uses a special combination of hidden units, elementwise products, and sums between units to implement gates that control "memory cells." These cells are designed to retain information without modification for long periods [34]. The greatest feature of LSTM is in its capability to learn long-term dependency, which is not possible with simple RNNs. To predict the next step, the weight values on the network have to be updated, which requires the maintenance of information from the initial steps. A simple RNN can only learn a limited number of short-term relationships and it cannot learn long-term series. However, LSTM can learn these long-term dependencies properly, and LSTM has three gates: input, forget, and output ( Figure 6b). The forget gate is embedded to indicate how much the previous memory remembers and how much it has forgotten. For LSTM, the hidden state ℎ is computed as follows: where , , and are the input, forget, and output gates at time t, respectively; , , , and are weights that map the hidden layer input to the three gates of input, forget, and output while , , , and weights matrices map the hidden layer output to gates; , , , and are vectors. Moreover, and ℎ are the outcome of the cell and the outcome of the layer, respectively [35]

Long Short-Term Memory (LSTM)
LSTM is a kind of model or structure for sequential data developed by [33] for the advancement of RNN. It uses a special combination of hidden units, elementwise products, and sums between units to implement gates that control "memory cells." These cells are designed to retain information without modification for long periods [34]. The greatest feature of LSTM is in its capability to learn long-term dependency, which is not possible with simple RNNs. To predict the next step, the weight values on the network have to be updated, which requires the maintenance of information from the initial steps. A simple RNN can only learn a limited number of short-term relationships and it cannot learn long-term series. However, LSTM can learn these long-term dependencies properly, and LSTM has three gates: input, forget, and output ( Figure 6b). The forget gate is embedded to indicate how much the previous memory remembers and how much it has forgotten. For LSTM, the hidden state h t is computed as follows: where i t , f t , and O t are the input, forget, and output gates at time t, respectively; W i , W f , W o , and W c are weights that map the hidden layer input to the three gates of input, forget, and output while U i , U f , U o , and U c weights matrices map the hidden layer output to gates; b i , b f , b o , and b c are vectors. Moreover, C t and h t are the outcome of the cell and the outcome of the layer, respectively [35].

Gated Recurrent Unit (GRU)
GRU is a simple type of LSTM suggested by Cho et al. [36]. The variance with LSTM is that GRU merges the input and forget gates and converts them with an update gate. Therefore, GRU has fewer parameters than LSTM which makes training easier. For GRU, the output value of h t is computed as follows: where r is a reset gate, and z denotes an update gate. The reset gate indicates how the new input can be combined with earlier memory. The update gate also indicates how much previous memory is kept. If the update gate is 1, the previous memory is fully preserved, and if it is 0, the previous memory is completely forgotten. There is a forget gate in the LSTM that automatically determines how much of the previous memory is maintained whereas, in GRU, all previous memories are maintained or completely forgotten (Figure 6c). For some problems, GRUs can provide a performance comparable with LSTMs, but with a lower memory requirement [34,36].

Bidirectional LSTM (Bi-LSTM)
In fact, instead of one direction, the network is processed on two sides: backward and forward with two separate hidden layers. Graves and Schmidhuber [37] demonstrated that bidirectional networks responded better than unidirectional networks in cases such as phonemic clustering. The structure of a bidirectional network is displayed in Figure 6d. Given this figure, these networks have a structure with a forward and backward LSTM layer. The forward layer output order, → h t , is repeatedly computed via inputs in a positive order from time T − n to time T − 1, while the backward layer outcome order, ← h t , is computed using the inverted inputs from time T − n to T − 1 [35]. The output of both forward and backward layers is computed in a similar manner to the unidirectional LSTM. In the Bi-LSTM layer, Y t is computed from Equation (13): where σ is a function used to merge two → h t and ← h t outputs [35]. In all networks, the dropout function is used among layers, which is a technique to prevent network overfitting. This means that learning on different architectures occurs with a variety of neurons. Figure 7a,b illustrate networks before and after dropout, respectively.
where σ is a function used to merge two ℎ ⃗⃗⃗ and ℎ ⃖⃗⃗⃗ outputs [35]. In all networks, the dropout function is used among layers, which is a technique to prevent network overfitting. This means that learning on different architectures occurs with a variety of neurons. Figure 7a and b illustrate networks before and after dropout, respectively. The learning rate is one of the hyperparameters necessary to find the optimal value and usually takes values between 1 and 1 × 10 −7 . In fact, the learning rate expresses the size of the move steps by the network. Figure 8 compares the changes in the loss function versus the epoch based on the learning rate. Working with a low learning rate, it takes a long time to find the best solution. Besides, if the learning rate is too high, it rejects the optimal mode. Since a high learning rate is advantageous in early iterations, and a low one is advantageous in later iterations, a learning rate that slows down as the algorithm progresses is preferred [39]. The learning rate is one of the hyperparameters necessary to find the optimal value and usually takes values between 1 and 1 × 10 −7 . In fact, the learning rate expresses the size of the move steps by the network. Figure 8 compares the changes in the loss function versus the epoch based on the learning rate. Working with a low learning rate, it takes a long time to find the best solution. Besides, if the learning rate is too high, it rejects the optimal mode. Since a high learning rate is advantageous in early iterations, and a low one is advantageous in later iterations, a learning rate that slows down as the algorithm progresses is preferred [39].

Open-Source Software and Codes
This research relies heavily on open-source software. Python 3.6 is used as a programming language. The libraries are used for preprocessing the data and for data management, including Pandas, Numpy, and Scikit-Learn. Deep-learning frameworks employed are TensorFlow and Keras.

Open-Source Software and Codes
This research relies heavily on open-source software. Python 3.6 is used as a programming language. The libraries are used for preprocessing the data and for data management, including Pandas, Numpy, and Scikit-Learn. Deep-learning frameworks employed are TensorFlow and Keras. Different codes are used for each RNN and ANN. All figures are created using Matplotlib.

Evaluation Criteria
In order to evaluate the models, correlation coefficients (CC), Nash-Sutcliff coefficient (NS), root mean square error (RMSE), and mean absolute error (MAE) are used, as presented in Equations (14) to (17), respectively.
where X i is the observation parameter with a mean denoted by X; Y i is the prediction parameter with a mean denoted by Y; N is a number of instances. The more the two first criteria are closer to 1 and the next three values are closer to 0 show the better performance of the model. According to Chiew et al. [41], if NS > 0.90, the simulation is very acceptable; if 0.60 < NS < 0.90, the simulation is acceptable; and if NS < 0.60 as in this case, simulation is unacceptable.

Results and Discussion
In this research, daily streamflow to the Ermenek Dam reservoir in Turkey was simulated using various deep-learning models. In this section, the historical observation streamflow data are compared with the computed streamflow from artificial neural networks and RNN, such as Bi-LSTM, GRU, LSTM, and simple RNN using seven lag days. A network attempts to predict outcomes as accurately as possible. The value of this precision in the network is obtained by the cost function, which tries to penalize the network when it fails. The optimal output is the one with the lowest cost. In this study, for the applied networks of MSE, the cost function is used. A repetition step in training generally works with a division of training data named a batch size. The number of samples for each batch is a hyperparameter, which is normally obtained by trial and error. In this study, the value of this parameter in all models is 512 in the best mode. In each repetition step, the cost function is computed as the mean MSE of these 512 samples of observed and predicted streamflow. The number of iteration steps for neural networks is named an epoch; in each epoch, the streamflow time series is simulated by the network once. Like other networks, neurons or network layers can be selected arbitrarily in recurrent networks. For the purpose of the comparison of models with each other, the structures of all recurrent network models are created identically. In each network, a double hidden layer is used so that there are 200 units in the first layer and 150 units of the neuron in the second layer. The last layer output of the network at the final time step is linked to a dense layer with a single output neuron. Between the layers, a dropout equal to 10% is used. The structure of the neural network is also used in two hidden layers. The first and second layers have 200 and 150 neurons, respectively. In all networks, the ReLU activation function is applied for the hidden layer. The main advantage of using ReLU is that, for all inputs greater than 0, there is a fixed derivative. This constant derivative speeds up network learning. Each method is run with epoch numbers of 100, 200, and 500. Results of different methods based on the evaluation criteria and their running times are presented in Table 2 (Computations are performed using a Pentium B960, 2.20 GHz laptop computer with Windows 8 and 4 GB RAM.). Among the methods, ANN is the fastest; among RNN methods, the speed sequence is as follows: simple RNN, GRU, LSTM, and Bi-LSTM. Different models are run with different iterations. As seen in Table 2, the number of epochs as one of the influential parameters plays a basic role; so in fewer iterations, the accuracy of the model is low, and the error is greater. By increasing the number of iterations, the model gradually converges so that there is not a significant difference between the results of 300 and 500 epochs. For this reason, it is decided that a maximum of 500 epochs is sufficient. Another parameter that influences the accuracy of the models is the number of neurons in the hidden layers. If it is low, the model will not be able to simulate correctly, and if it is too high, there is a risk of overfitting. This problem is solved by using the dropout method, which in fact deactivates several neurons. The values of the learning rate (LR) and decay are also presented in Table 2. Decay occurs due to how much the learning rate is reduced in each step. As seen in Table 2, five different evaluation criteria are used to compare five different prediction methods. Each method has the best result with 300 or 500 epochs. The performance of the testing stage is generally 1-5% lower than that of the training stage. According to Table 2, among the methods used, the LSTM method performs best in 500 epochs with an accuracy of CC = 87% in the testing stage. This network has long-term memory, and its forget gate specifies how much previous memory is kept. The first step is the multiplication of input data in weights and then its summation with the bias, followed by output. In this step, the output is likely to be very different from the actual output. Therefore, errors are returned backward to update the weights and also biases.
Comparison of some metrics by all methods used in the study are given in Table 3. In this table, measured streamflow values (those set aside for the test and the whole streamflow data set) and values predicted by the methods can be compared. The GRU and LSTM models had the closest results to measured values in max, mean, and standard deviation (StDev) statistics in terms of similarity to the measured values. Simple RNN, ANN, and BLSTM methods gave the worst results, respectively.
The loss function, observational streamflow values and predicted values for training and testing stages computed by various methods are displayed in Figure 9. When the loss function and time-series graphs are examined, the most successful methods are LSTM and GRU. The loss charts of testing and training stages of both methods overlap at the highest epoch value. According to the loss function graphs (Figure 9), with the reduction of the modeling error for the training stage, the error of the testing stage decreases, and the distance between the two graph lines reduces. Therefore, it can be said that the dropout function helps to prevent network overfitting. Besides, considering the changes in the loss function value versus epochs (and according to Figure 8) can ensure that the learning rate is selected correctly.  Figure 10 shows the scatter of observed values versus the predicted values for both training and testing stages. Regarding the streamflow graphs (Figures 9 and 10), the networks yield poor results in peak periods and cannot simulate these periods well, except LSTM. One major reason of this is that streamflow abundance has very few high values, and the network cannot properly learn them. This happens even though the streamflow abundance has several low values, so the model can correctly and accurately be learned in the training stage. According to the data, the difference between the minimum and the maximum flow values is high, and the validity of these values is investigated during the study period based on precipitation in the region. In periods where peak streamflow is observed, the maximum rainfall is observed with a significant amount that indicates the correctness of peak flow values. However, among the methods applied, the LSTM network is better than other methods, and is able to simulate peak flow periods fairly. As noted above, the greatest feature of this network is its skill to learn long-term dependencies, and the forget gate makes the network keep or forget the desired amount of previous memory and thus helps to improve the modeling. The GRU method is similar to the LSTM method, and its results are close to those of LSTM. In contrast, the other methods failed to simulate peak flow values. The scatter plots in Figure 10 indicate a relatively small dispersion of observed and predicted data in training and testing stages compared to other methods. Kratzert et al. [19] used deep LSTM networks in similar research and demonstrated that these networks perform better than other networks and can be used in hydrological simulations as an application method.   Figure 10 shows the scatter of observed values versus the predicted values for both training and testing stages. Regarding the streamflow graphs (Figures 9 and 10), the networks yield poor results in peak periods and cannot simulate these periods well, except LSTM. One major reason of this is that streamflow abundance has very few high values, and the network cannot properly learn them. This happens even though the streamflow abundance has several low values, so the model can correctly and accurately be learned in the training stage. According to the data, the difference between the minimum and the maximum flow values is high, and the validity of these values is investigated during the study period based on precipitation in the region. In periods where peak streamflow is observed, the maximum rainfall is observed with a significant amount that indicates the correctness of peak flow values. However, among the methods applied, the LSTM network is better than other methods, and is able to simulate peak flow periods fairly. As noted above, the greatest feature of this network is its skill to learn long-term dependencies, and the forget gate makes the network keep or forget the desired amount of previous memory and thus helps to improve the modeling. The GRU method is similar to the LSTM method, and its results are close to those of LSTM. In contrast, the other methods failed to simulate peak flow values. The scatter plots in Figure 10 indicate a relatively small dispersion of observed and predicted data in training and testing stages compared to other methods. Kratzert et al. [19] used deep LSTM networks in similar research and demonstrated that these networks perform better than other networks and can be used in hydrological simulations as an application method.

Conclusions
In this paper, a novel approach to streamflow simulation based on recurrent neural networks is presented. For this purpose, 70% of the inflow data are trained to ANN and RNNs, including Bi-LSTM, LSTM, GRU, and simple RNN. Then, the 30% remaining data were predicted. For each of the networks, the double hidden layer was used. The number of neurons in the hidden layers, the learning rate, and the number of iterations have a basic role in modeling accuracy, the optimal values of which were obtained using trial and error. Using the dropout function in the network structure can well prevent network overfitting. To attain this purpose, all models were coded in Python programming software. The efficiency of the proposed approach is evaluated for the simulation of daily inflow to the Ermenek Dam reservoir in Turkey. Results of the study point to the better performance of recurrent networks (NS: 0.74, RMSE:17.92) compared to ANN (NS: 0.72, RMSE:18.71). Among different RNN architectures, the LSTM network performed the best and can estimate streamflow with fairly good accuracy (CC: 0.87), and it can be said that these networks can be used as an effective method of streamflow modeling and simulation.
Thermal power plants are not environmentally friendly and cause air pollution and global warming. Due to limited fossil resources in Turkey the importance of the hydroelectric power plant is increased. Considering the purpose of the reservoir construction, which is the production of the hydroelectric power or agricultural irrigation, the prediction of streamflow may help to increase efficiency and achieve the desired purpose. In this study, for the first time, deep learning methods are used to model inflows to the Ermenek Dam reservoir located in Turkey. Results of this study are consistent with those of studies by Nagesh Kumar et al. [5] and Zhang et al. [17] Given that water resource management is essential, it can be acknowledged that the use of these methods can help managers and officials of the Ermenek Dam in future planning. As a result, increasing the accuracy of the model can help researchers and experts to achieve this goal.