text:
Research Article

Journal of Korea Water Resources Association. June 2020. 395-408
https://doi.org/10.3741/JKWRA.2020.53.6.395

# MAIN

• 1. INTRODUCTION

• 2. MATERIAL

•   2.1 Data description

•   2.2 Artificial neural networks (ANN)

•   2.3 Adaptive neuro-fuzzy inference system (ANFIS)

•   2.4 Wavelet-based neural networks (WNN)

•   2.5 Genetic algorithm (GA)

•   2.6 K-fold cross validation (CV)

• 3. METHODOLOGY

•   3.1 Models development

•   3.2 Coding for the ANN, ANFIS, and WNN models

•   3.3 Measures of accuracy

• 4. Results and Discussion

• 5. Conclusion

1. INTRODUCTION

Streamflow forecasting based on accurate measurements can be used to design flood mitigation structures for agricultural and urban basins, and build water allocation systems for agricultural, industrial, and commercial purposes (Rezaie-Balf et al., 2019; Samsudin et al., 2011). The complex and significant variability of streamflow can be explained using spatial and temporal characteristics. It has guided the evolution and application of different approaches for estimation, modeling, forecasting, and prediction (Martins et al., 2011). Forecasting of hydrological time series (e.g., streamflow, water stage, evaporation, and groundwater stage etc.) is important for understanding the internal relationship of natural processes (Krishna et al., 2011). Since the complexity of hydrological time series requires specific tools of non-linear and non-stationary dynamic systems, the diverse forecasting methods have been proposed for hydrological forecasting (Sivakumar et al., 2001).

Data-driven approaches, although credible for hydrological forecasting, don’t have the capability to depict physical processes, because they only consider an adequate selection of hydrological variables with temporal and input-output modification. Therefore, these approaches can be categorized as two types (i.e., classicial and computational intelligence approaches) typically. The classical approaches can be exemplified as linear regression (LR), auto regressive integrated moving average (ARIMA), and ARIMA with exogenous input (ARIMAX) etc. The machine learning approaches can be classified as artificial neural networks (ANN), adaptive neuro-fuzzy inference systems (ANFIS), genetic programming (GP), gene expression programming (GEP), model tree (MT), extreme learning machines (ELM), support vector machines (SVM), and multivariate adaptive regression spline (MARS) etc. The field of machine learning approaches has undergone innovative changes for novel techniques of data simulation and processing (Chandwani et al., 2015).

For three decades, the ANN model based on neuron systems has been used for different types of hydrological forecasting. Various approaches have been applied to validate the accuracy and effectiveness of the ANN model for streamflow forecasting (Biswas and Jayawardena, 2014; Badrzadeh et al., 2013; Moradkhani et al., 2004; Cigizoglu, 2003; Abrahart and See, 2000). The ANFIS model (Jang, 1993), which has the merits of the ANN model and the fuzzy system, has been employed for streamflow forecasting (Talei et al., 2010; Nasr and Bruen, 2008; Keskin et al., 2006; Lohani et al., 2006).

Methodologies combining wavelet and machine learning approaches have been utilized for streamflow forecasting. Combined approaches have outstripped conventional models (Nourani et al., 2014). The wavelet-based machine learning approaches, including wavelet-based neural networks (WNN), wavelet-based support vector regression (WSVR), and wavelet-based adaptive neuro fuzzy inference system (WANFIS), have been effectively implemented for hydrological forecasting, including streamflow, water stage, runoff, and groundwater etc. (Seo et al., 2015; Kamruzzaman et al., 2013; Partal, 2009; Partal and Kisi, 2007). Many researchers have attempted to forecast streamflow using wavelet-based machine learning approaches (Zakhrouf et al., 2016, 2020; Shoaib et al., 2014; Badrzadeh et al., 2013; Nourani et al., 2013; Guo et al., 2011; Tiwari and Chatterjee, 2010).

The optimal structure of wavelet-based machine learning approaches can be constructed as a search method. A method for designing machine learning approaches using evolutionary optimization algorithm is proposed to format the best models. Evolutionary machine learning approaches for modeling hydrologic systems have been suggested by Zakhrouf et al. (2018), Kalteh (2015), Sahay and Srivastava (2014), Asadi et al. (2013).

K-fold cross validation (CV) method is one of the methods to assess the algorithmic generalization. The out-of-sample cross validation (OOS-CV) method is a frequently used method for hydrological modeling. This paper attempts to develop a new approach combining wavelet transformation, machine learning approach, evolutionary optimization algorithm, and k-fold cross validation method for multi-step (days) (i.e., t+1, t+2, and t+3 days) streamflow forecasting in the Seybous River, North Algeria. The paper is divided into five chapters. The first part provides a brief introduction. The second part discusses data and application tools, including ANN, ANFIS, WNN, GA, and k-fold CV, respectively. The third part applies the methodology, and results and discussion are presented in the fourth part. Conclusions are stated in the concluding part.

2. MATERIAL

2.1 Data description

The data for training and testing phases of the developed models were obtained from the Seybous River basin, Algeria. The basin is positioned between 36.0 N and 38.0 N (latitudes), and between 7.0 E and 8.0 E (longitudes), Algeria, North Africa, and covers a total catchment area of 6,471.0 km2 (Fig. 1). Daily streamflow data for 14 years (September, 1982-August, 1996) were obtained from National Agency of Water Resources of Mirebeck (14 06 01) gauging station and were utilized for multi-step (days) streamflow forecasting. These selected multi-step (days) can be adequate, considering the rapid surface streamflow in the Mediterranean basin and the watershed size. For this purpose, the first ten years (70% of data) were utilized for the training phase, and the second four years (30% of data) were utilized for the testing phase as shown in Fig. 2.

##### Fig. 1.

The Seybous River basin and Mirebeck station (Aichouri et al., 2015)

##### Fig. 2.

The observed streamflow hydrograph (14 years)

2.2 Artificial neural networks (ANN)

ANN is a parallel computational system based on the architectural and functional standard of biological networks (Imrie et al., 2000). The feedforward neural networks (FFNN) is one of the commonly applied neural network architectures. It consists of an input layer (i.e., first layer) which gets any outside news directly; an output layer (i.e., final layer); and one or more transitional layers (i.e., hidden layers) which divide the input and output layers. The neurons in the layers are usually fully connected using signs from the first layer to the final layer (Zhang et al., 1998), where the signs are handled by a consolidation function combining first all incoming signs, and second an activation function converting output neuron (Günther and Fritsch, 2010). Fig. 3 shows a conventional FFNN architecture.

##### Fig. 3.

The conventional FFNN model architecture

2.3 Adaptive neuro-fuzzy inference system (ANFIS)

ANFIS is a complex function combining an adaptive neural network and a fuzzy inference system (FIS). The basic formation of a FIS consists of a structure that draws input to membership functions (Moosavi et al., 2013). Since ANFIS is based on FIS, it can be described using fuzzy IF-THEN rules (Jang et al., 1997). There are two methods for FIS: Mamdani (Mamdani and Assilian, 1975) and Sugeno (Takagi and Sugeno, 1985). The differences between the two methods can be explained using the consequent part. Mamdani’s method utilizes fuzzy membership functions (MFs), whereas Sugeno’s method applies linear or constant functions. In this study, the procedures presented by Jang et al. (1997) using Sugeno’s approach were utilized for multi-step (days) streamflow forecasting. Fig. 4 shows the conventional ANFIS architecture. The ANFIS was trained using the least-squares method (LSM) and backpropagation gradient descent method (BP-GDM). Detailed explanation for the ANFIS model can be found in Jang (1993).

##### Fig. 4.

The conventional ANFIS model architecture (Seo et al., 2015)

2.4 Wavelet-based neural networks (WNN)

The fundamental objective of the wavelet transform (WT) is to finalize a complete time scale formation of localized and passing events happening at different time scales (Labat et al., 2000). There are two categories of WT: continuous wavelet transform (CWT) and discrete wavelet transform (DWT). The CWT method can be applied for continuous functions or time series (Nourani et al., 2009; Mallat, 1989). Since the CWT method is time-consuming, it requires large resources, whereas the DWT method can be more easily applied (Nourani et al., 2009; Mallat, 1989). The DWT requires less calculation time, is not complicated (Ghaemi et al., 2019; Quilty and Adamowski, 2018; Seo et al., 2015, 2018; Kisi, 2011; Nourani et al., 2009; Mallat, 1989), and involves choosing scales and positions (Adamowski and Sun, 2010).

A WNN model is a combination of WT and ANN, where the original time series for each input neuron of the ANN model is decomposed to some multi-frequency time series using the WT algorithm. The decomposed details (D) and approximation (A) are formatted as new inputs to the ANN model architecture, as shown in Fig. 5. The WNN model consists of a two-stage algorithm. The first stage can be specialized with the multi-level WT where the original time series is decomposed utilizing the WT after choosing the decomposition level. The second stage can be classified as the training and testing phases using the ANN model. The details and approximation, which are obtained by the WT in the first stage, are utilized as input of the ANN model. In the WNN, the number of decomposition levels (details and approximation) of each input is calculated according to full data length (Nourani et al., 2009).

##### Fig. 5.

The schematic diagram of WNN model with one input model

 $$L=int\lbrack log(N)\rbrack$$ (1)

where L defines the decomposition level, N denotes the number of time series data, and int[·] depicts the integer-part function. In this study, the data comprised N = 5110 samples; therefore, the level decomposition was classified as L = 3.

2.5 Genetic algorithm (GA)

GA, developed by Holland (1975), is a heuristic and evolutionary optimization technique based on evolutionary and genetic theory (Kim and Kim, 2008). The procedure of GA starts with random strings representing design or decision variables. Each string is then evaluated to allocate the fitness value, based on the objective and constraint conditions. Then, the termination condition is validated in the algorithm (Gowda and Mayya, 2014). This study focuses on optimizing the structure of three machine learning approaches (i.e., ANN, ANFIS, and WNN) using the evolutionary optimization algorithm (i.e., GA).

The first step is to choose the GA chromosome. Each individual in the population represents a possible configuration for architecture of the machine learning approach. Based on the concept of evolutionary optimization theory, GA starts with a population of chromosomes, which evolve towards optimum solutions by GA operators, including selection, crossover, and mutation. These steps are reiterated from one generation to the next with the aim of arriving at the general optimal solution (Kim and Kim, 2008).

2.6 K-fold cross validation (CV)

The CV is a statistical methodology for comparing and evaluating training algorithms by separating data into two parts. One is utilized to train a model and the other is utilized to test it (Stone, 1974). K-fold CV assesses the generalization of algorithms in the evolutionary machine learning approach. Based on the category of k-fold CV, the data is first separated into k equally measured parts. Subsequently, k iterations of training and testing phases are performed such that a different part of the data is held-out for testing phase within each iteration, while the remaining (k-1) parts are used for training phase (Zhao et al., 2018)

Suppose we have a model with one or more unknown parameters f ($\alpha$), and a data set $\stackrel{^}{y}$ to which the model can be fitted. The primary method to estimate the tuning parameter $\alpha$ using k-fold CV divides the data into rougher parts (Fig. 6). Since the model computes the mean squared error (MSE) for each i= 1, 2, …, k, the cross validation error (CVE) can be calculated as Eq. (2).

##### Fig. 6.

Representation of k-fold cross validation (CV)

 $$CVE(\alpha)=\frac1K\sum_{i=1}^kMSE_i(\alpha)$$ (2)

where CVE defines the cross validation error.

3. METHODOLOGY

3.1 Models development

The machine learning approaches (i.e., ANN, ANFIS, and WNN) were utilized for multi-step (days) streamflow forecasting using previous time series values. The training phase of machine learning approaches provides a non-linear matching between inputs and outputs, and is useful in identifying patterns utilizing complicated data (Liu and Chung, 2014). Since the backpropagation (BP) algorithm cannot be guaranteed to find the minimum error margin, the convergence cannot occur with fast track. The solution for this problem can be applied with a fitness function that tests how well an architecture learns from the data. The fitness function can be expressed as Eq. (3).

 $$F=Min\lbrack CVE\rbrack$$ (3)

where F denotes the fitness function. In general, GA can significantly reduce the weakness of BP algorithm. The data utilized for the training phase (the 75% of the data) were sub-divided into 10 subsets (10-fold cross validation). Real coding was utilized to find the favorable topology for the ANN, ANFIS, and WNN models. The coding to find the best architecture and parameters of machine learning approaches (i.e., ANN, ANFIS, and WNN) is described as follows.

3.2 Coding for the ANN, ANFIS, and WNN models

Feedforward neural networks (FFNN) models with two hidden layers and one neuron in final layer were utilized. A chromosome was built from a series of genes (Fig. 7) to find: the input delay (D), the number of neuron in the hidden layers (NHL1 and NHL2), the activation functions in the hidden and final layers (AFHL1, AFHL2, and AFOL) including linear transfer function (Purelin), symmetric saturating linear transfer function (Satlins), log-sigmoid transfer function (Logsig), and hyperbolic tangent sigmoid transfer function (Tansig), respectively, and the initial connection coefficients of weights and bias (IWB).

##### Fig. 7.

Chromosome encoding for the ANN model

A chromosome of ANFIS model using a series of genes was created to find different parameters (Fig. 8): the input delay (D); the number of membership functions (NMF); the type of membership functions (TMF) including Π-shaped (Pimf), Trapezoidal-shaped (Trapmf), Triangular-shaped (Trimf), Gaussian curve (Gaussmf), and Built-in Gaussian function (Dsigmf); the definition of if-then (AND operation/OR operation) rules type (DRT); and the firing strength of a rule (FSR).

##### Fig. 8.

Chromosome encoding for the ANFIS model

In this study, two methods were used as firing strength of AND rule, including Minimum (Min) and Product (Prod). Also, two methods were used as firing strength OR rule, including Maximum (Max) and the probabilisticOr method (Probor). For the membership functions type, there are many kinds of membership functions (MFs).

A chromosome of different parameters for wavelet-based neural networks (WNN) model was created from a series of genes (Fig. 9) to find: the input delay (D); the type of mother wavelet (TMW) based on the five most frequently used wavelet families, including Haar (Har), Daubechies (Db), Coiflets (Coif), Symlets (Sym), and Biorthogonal (Bior); the number of neuron in hidden layers (NHL1 and NHL2); the activation functions in hidden and output layers (AFHL1, AFHL2, and AFOL), including linear transfer function (Purelin), symmetric saturating linear transfer function (Satlins), log-sigmoid transfer function (Logsig) and hyperbolic tangent sigmoid transfer function (Tansig); and the initial connection weights and bias coefficients (IWB).

##### Fig. 9.

Chromosome encoding for the WNN model

3.3 Measures of accuracy

To assess the performance of three different machine learning approaches (i.e., ANN, ANFIS, and WNN) to forecast multi-step (days) streamflow during the testing phase, four statistical indices (measures of accuracy) were applied, including root mean square error (RMSE), Nash-Sutcliffe efficiency (NSE), correlation coefficient (R), and peak flow criteria (PFC).

RMSE would vary from zero to a large value which presents prefect forecasting by the difference between observed and forecasted streamflow. NSE is considered for evaluating the ability of hydrological models (Rezaie-Balf et al., 2019; Nash and Sutcliffe, 1970). R is an assessment of the precision of hydrologic modeling and is used for comparisons of alternative models. A perfect matching produces R = 1.0 (Kim and Kim, 2008).

These statistical indices (i.e., RMSE, NSE, and R) may not provide the model performance for extreme streamflow events (e.g., floods and drought). Therefore, it is fundamental to evaluate the extreme distributions using PFC for forecasting extreme events (Rezaie-Balf et al., 2019). PFC plays an important role in following the extreme events to achieve the efficient model. The RMSE, NSE, R, and PFC indices can be represented as Eqs. (4)~(7).

 $$RMSE=\sqrt{\frac1N\sum_{i=1}^N\left[Qt_i-\widehat Qt_i\right]^2}$$ (4)
 $$NSE=\left(1-\frac{{\displaystyle\sum_{i=1}^n}\left[Qt_i-\widehat Qt_i\right]^2}{{\displaystyle\sum_{i=1}^N}\left[Qt_i-\overline Qt_i\right]^2}\right)\times100$$ (5)
 $$R=\frac{{\displaystyle\sum_{i=1}^N}\left[Qt_i-\overline Qt\right]\;\left[\widehat Qt_i-\widetilde Qt\right]}{\sqrt{{\displaystyle\sum_{i=1}^N}\left[Qt_i-\overline Qt\right]^2\;{\displaystyle\sum_{i=1}^N}\left[\widehat Qt_i-\widetilde Qt\right]^2}}$$ (6)
 $$PFC=\frac{\sqrt[4]{{\displaystyle\sum_{i=1}^{T_p}}\left[Qt_i-\widehat Qt_i\right]^2\cdot Qt_i^2\;}}{\sqrt{{\displaystyle\sum_{i=1}^{T_p}}Qt_i^2}}$$ (7)

where $Q{t}_{i}$ = the observed streamflow value, $\stackrel{^}{Q}{t}_{i}$ = the forecasted streamflow value, $\overline{)Q}t$ = the mean observed streamflow, $\stackrel{~}{Q}t$ = the mean forecasted streamflow, N = the number of data, Tp = the number of peak streamflow greater than one third of the observed mean peak streamflow.

4. Results and Discussion

Table 1 presents the optimal structure of ANN model using GA which can be found from table 1 for (t+1) day forecasting as follows; input delay = 3, number of neurons (1st hidden layer) = 9, number of neurons (2nd hidden layer) = 8, activation function (1st hidden layer) = log-sigmoid transfer function, activation function (2nd hidden layer) = linear transfer function, and activation function (final layer) = linear transfer function. For (t+2) day forecasting, the following can be suggested; input delay = 3, number of neurons (1st hidden layer) = 9, number of neurons (2nd hidden layer) = 8, activation function (1st hidden layer) = log-sigmoid transfer function, activation function (2nd hidden layer) = log-sigmoid transfer function, and activation function (final layer) = linear transfer function. For (t+3) day forecasting, the following can be provided; input delay = 5, number of neurons (1st hidden layer) = 5, number of neurons (2nd hidden layer) = 8, activation function (1st hidden layer) = symmetric saturating linear transfer function, activation function (2nd hidden layer) = log-sigmoid transfer function, and activation function (final layer) = linear transfer function.

#### The optimal structure for ANN model using GA

 ANN parameters t+1 t+2 t+3 D 3 3 5 NHL1 9 9 5 NHL2 8 8 8 AFHL1 Logsig Logsig Satlins AFHL2 Purelin Logsig Logsig AFOL Purelin Purelin Purelin

Table 2 shows the optimal structure of ANFIS model using GA which can be expressed from table 2 for (t+1) day forecasting as follows; input delay = 1, number of membership functions = 3, type of membership functions = Gaussian curve, firing strength of a rule = product, and definition of if-then rules type = and. For (t+2) day forecasting, the following can be suggested; input delay = 2, number of membership functions = 3, type of membership functions = trapezoidal-shaped, firing strength of a rule = product, and definition of if-then rules type = and. For (t+3) day forecasting, the following can be provided; input delay = 2, number of membership functions = 2, type of membership functions = Gaussian curve, firing strength of a rule = probabilisticOr, and definition of if-then rules type = or.

#### The optimal structure for ANFIS model using GA

 ANFIS parameters t+1 t+2 t+3 D 1 2 2 NMF 3 3 2 TMF Gaussmf Trapmf Gaussmf FSR Prod Prod Probor DRT And And Or

Table 3 proposes the optimal structure of WNN model using GA which can be supposed from table 3 for (t+1) day forecasting as follows; input delay = 4, type of mother wavelet = Symlets 7, number of neurons (1st hidden layer) = 5, number of neurons (2nd hidden layer) = 10, activation function (1st hidden layer) = log-sigmoid transfer function, activation function (2nd hidden layer) = log-sigmoid transfer function, and activation function (final layer) = linear transfer function. For (t+2) day forecasting, the following can be suggested; input delay = 4, type of mother wavelet = Daubechies 9, number of neurons (1st hidden layer) = 8, number of neurons (2nd hidden layer) = 10, activation function (1st hidden layer) = symmetric saturating linear transfer function, activation function (2nd hidden layer) = log-sigmoid transfer function, and activation function (final layer) = linear transfer function. For (t+3) day forecasting, the following can be provided; input delay = 5, type of mother wavelet = Symlets 6, number of neurons (1st hidden layer) = 8, number of neurons (2nd hidden layer) = 7, activation function (1st hidden layer) =linear transfer function, activation function (2nd hidden layer) = symmetric saturating linear transfer function, and activation function (final layer) = linear transfer function.

#### The optimal structure for WNN model using GA

 WNN parameters t+1 t+2 t+3 D 4 4 5 MWT sym7 db9 Sym6 NHL1 5 8 8 NHL2 10 10 7 AFHL1 Logsig Satlins Purelin AFHL2 Logsig Logsig Satlins AFOL Purelin Purelin Purelin

Table 4 suggests the performances of ANN, ANFIS, and WNN models for different multi-step (days) (i.e., t+1, t+2, and t+3 days) forecasting. The statistical results of standalone (i.e., ANN and ANFIS) models yielded similar performances based on RMSE, NSE, R, and PFC indices for training and testing phases. The performances of standalone models were not better than those of hybrid model (i.e., WNN). For example, the values of RMSE and PFC for the WNN model (RMSE = 8.590 m3/sec and PFC = 0.252, (t+1) day forecasting) were lower than those of ANN (RMSE = 19.120 m3/sec and PFC = 0.446, (t+1) day forecasting) and ANFIS (RMSE = 18.520 m3/sec and PFC = 0.444, (t+1) day forecasting) models in the testing phase. In addition, the values of NSE and R for the WNN model (NSE = 92.000% and R = 0.969, (t+1) day forecasting) were higher than those of ANN (NSE = 60.400% and R = 0.783, (t+1) day forecasting) and ANFIS (NSE = 62.860% and R = 0.793, (t+1) day forecasting) models in the testing phase. The performances of hybrid model were superior to those of standalone models.

As the multi-step (days) for the three models increased from (t+1) to (t+3) days, the model accuracy decreased. The hybrid model applied sub-time series by using WT as model input, while the standalone models utilized the original time series as model input without WT. It can be suggested from table 4 that the hybrid model using sub-time series as input data of the standalone models can improve the performance of conventional standalone models. Fig. 10 shows the scatter diagrams for the ANN, ANFIS, and WNN models in the testing phase. Fig. 11 presents the relative errors of peak flow for the ANN, ANFIS, and WNN models in the testing phase. Figs. 10 and 11 show that the performance of hybrid model was superior to that of standalone models. From the performance comparison, the performance of three models with (t+1) day was superior to that of (t+2) and (t+3) days. In addition, Figs. 12(a)~12(c) display the daily streamflow hydrographs for the ANN, ANFIS, and WNN models based on (t+1) day forecasting in the testing phase.

#### Comparison between the performance results obtained by the models: ANN, ANFIS, and WNN in the training and testing phases

 Training Testing RMSE (m3/sec) NSE (%) R PFC RMSE (m3/sec) NSE (%) R PFC ANN t+1 23.850 70.100 0.837 0.463 19.120 60.400 0.783 0.446 t+2 30.990 49.520 0.704 0.513 25.390 30.200 0.560 0.501 t+3 37.270 27.000 0.520 0.645 30.800 -2.710 0.352 0.552 ANFIS t+1 25.610 65.530 0.810 0.479 18.520 62.860 0.793 0.444 t+2 35.910 32.240 0.568 0.626 26.090 26.310 0.535 0.510 t+3 38.470 22.210 0.472 0.653 27.140 20.260 0.457 0.546 WNN t+1 3.530 99.350 0.997 0.103 8.590 92.000 0.969 0.252 t+2 5.230 98.560 0.993 0.128 13.360 80.680 0.917 0.318 t+3 7.950 96.680 0.983 0.103 12.950 81.850 0.907 0.391

The machine learning approaches were combined a data pre-processing technique and an evolutionary optimization algorithm based on k-fold cross validation. In addition, an attempt to do multi-step (days) streamflow forecasting is different compared to the previous studies (Rezaie-Balf et al., 2019; Seo et al., 2018; Zakhrouf et al., 2018; Abdollahi et al., 2017; Ravansalar et al., 2017; Zakhrouf et al., 2016). The forecasts confirmed the accuracy and effectiveness of developed models, based on the statistical indices and scatter diagrams. Different machine learning approaches (e.g., gene expression programming (GEP), support vector machines (SVM), extreme learning machines (ELM), model tree (MT), and multivariate adaptive regression spline (MARS) etc.) combined data pre-processing techniques (e.g., principal component analysis (PCA), maximum entropy spectral analysis (MESA), empirical mode decomposition (EMD), and ensemble empirical mode decomposition (EEMD) etc.) and evolutionary optimization algorithms (e.g., particle swarm optimization (PSO), ant colony optimization (ACO), harmony search (HS), tabu search (TS), and crow search algorithm (CSA) etc.) can also be pursued in further studies.

##### Fig. 10.

Scatter diagrams for the ANN, ANFIS, and WNN models in the testing phase

##### Fig. 11.

Relative errors of the peak flows for the ANN, ANFIS, and WNN models in the testing phase

##### Fig. 12.

Daily streamflow hydrographs for the ANN, ANFIS, and WNN models based on the (t+1) day forecasting in the testing phase

5. Conclusion

This study suggests three machine learning approaches, including artificial neural networks (ANN), adaptive neuro-fuzzy inference system (ANFIS), and wavelet-based neural networks (WNN) models for multi-step (days) streamflow forecasting in Seybous River basin, Algeria. The accuracy and effectiveness of developed models (i.e., ANN, ANFIS, and WNN) are evaluated based on four statistical indices, including root mean square error (RMSE), Nash-Sutcliffe efficiency (NSE), correlation coefficient (R), and peak flow criteria (PFC).

Based on the combination of evolutionary optimization algorithm and k-fold cross validation, the statistical results of hybrid (i.e., WNN) model are superior to those of standalone (i.e., ANN and ANFIS) models for different multi-step (days). Also, the results of (t+1) day streamflow forecasting are superior to those of (t+2) and (t+3) days, based on statistical indices and scatter diagrams. The combination of machine learning approach and data pre-processing technique based on the evolutionary optimization algorithm and k-fold cross validation can be a potential implement for accurate multi-step (days) streamflow forecasting. Hybrid methodologies combining diverse machine learning approaches and data pre-processing technique based on other evolutionary optimization algorithms and cross validations can be recommended for further studies.

## References

1
Abdollahi, S., Raeisi, J., Khalilianpour, M., Ahmadi, F., and Kisi, O. (2017). "Daily mean streamflow prediction in perennial and non-perennial rivers using four data driven techniques." Water Resources Management, Vol. 31, No. 15, pp. 4855-4874.
10.1007/s11269-017-1782-7
2
Abrahart, R.J., and See, L. (2000). "Comparing neural network and autoregressive moving average techniques for the provision of continuous river flow forecasts in two contrasting catchments." Hydrological Processes, Vol. 14, No. 11‐12, pp. 2157-2172.
10.1002/1099-1085(20000815/30)14:11/12<2157::AID-HYP57>3.0.CO;2-S
3
Adamowski, J., and Sun, K. (2010). "Development of a coupled wavelet transform and neural network method for flow forecasting of non-perennial rivers in semi-arid watersheds." Journal of Hydrology, Vol. 390, No. 1, pp. 85-91.
10.1016/j.jhydrol.2010.06.033
4
Aichouri, I., Hani A., Bougherira, N., Djabri, L., Chaffai, H., and Lallahem, S. (2015). "River flow model using artificial neural networks." Energy Procedia, Vol. 74, pp. 1007-1014.
10.1016/j.egypro.2015.07.832
5
Asadi, S., Shahrabi, J., Abbaszadeh, P., and Tabanmehr, S. (2013). "A new hybrid artificial neural networks for rainfall-runoff process modeling." Neurocomputing, Vol. 121, pp. 470-480.
10.1016/j.neucom.2013.05.023
6
Badrzadeh, H., Sarukkalige, R., and Jayawardena, A.W. (2013). "Impact of multi-resolution analysis of artificial intelligence models inputs on multi-step ahead river flow forecasting." Journal of Hydrology, Vol. 507, pp. 75-85.
10.1016/j.jhydrol.2013.10.017
7
Biswas, R.K., and Jayawardena, A.W. (2014). "Water level prediction by artificial neural network in a flashy transboundary river of Bangladesh." Global Nest Journal, Vol. 16, No. 2, pp. 432-444.
10.30955/gnj.001226
8
Chandwani, V., Vyas, S.K., Agrawal, V., and Sharma, G. (2015). "Soft computing approach for rainfall-runoff modelling: A review." Aquatic Procedia, Vol. 4, No. 2015, pp. 1054-1061.
10.1016/j.aqpro.2015.02.133
9
Cigizoglu, H.K. (2003). "Estimation, forecasting and extrapolation of river flows by artificial neural networks." Hydrological Sciences Journal, Vol. 48, No. 3, pp. 349-361.
10.1623/hysj.48.3.349.45288
10
Ghaemi, A., Rezaie-Balf, M., Adamowski, J., Kisi, O., and Quilty, J. (2019). "On the applicability of maximum overlap discrete wavelet transform integrated with MARS and M5 model tree for monthly pan evaporation prediction." Agricultural and Forest Meteorology, Vol. 278, p. 107647.
10.1016/j.agrformet.2019.107647
11
Gowda, C.C., and Mayya, S.G. (2014). "Comparison of back propagation neural network and genetic algorithm neural network for stream flow prediction." Journal of Computational Environmental Sciences, Vol. 2014, p. 290127.
10.1155/2014/290127
12
Günther, F., and Fritsch, S. (2010). "Neuralnet: Training of neural networks." The R journal, Vol. 2, No. 1, pp. 30-38.
10.32614/RJ-2010-006
13
Guo, J., Zhou, J., Qin, H., Zou, Q., and Li, Q. (2011). "Monthly streamflow forecasting based on improved support vector machine model." Expert Systems with Applications, Vol. 38, No. 10, pp. 13073-13081.
10.1016/j.eswa.2011.04.114
14
Holland, J.H. (1975). Adaptation in natural and artificial systems. University of Michigan Press, Ann Arbor, M.I., U.S.
15
Imrie, C.E., Durucan, S., and Korre, A. (2000). "River flow prediction using artificial neural networks: Generalisation beyond the calibration range." Journal of Hydrology, Vol. 233, No. 1, pp. 138-153.
10.1016/S0022-1694(00)00228-6
16
Jang, J.S.R. (1993). "ANFIS: Adaptive-network-based fuzzy inference system." IEEE Transactions on Systems, Man, and Cybernetics, Vol. 23, No. 3, pp. 665-685.
10.1109/21.256541
17
Jang, J.S.R., Sun, C.T., and Mizutani, E. (1997). Neuro-fuzzy and soft computing: A computational approach to learning and machine intelligence. Prentice-Hall, New Jersey, U.S.
10.1109/TAC.1997.633847
18
Kalteh, A.M. (2015). "Wavelet genetic algorithm-support vector regression (wavelet GA-SVR) for monthly flow forecasting." Water Rsources Mnagement, Vol. 29, pp. 1283-1293.
10.1007/s11269-014-0873-y
19
Kamruzzaman, M., Metcalfe, A.V., and Beecham, S. (2013). "Wavelet-based rainfall-stream flow models for the southeast Murray Darling basin." Journal of Hydrologic Engineering, Vol. 19, No. 7, pp. 1283-1293.
10.1061/(ASCE)HE.1943-5584.0000894
20
Keskin, M.E., Taylan, D., and Terzi, Ö. (2006). "Adaptive neural-based fuzzy inference system (ANFIS) approach for modelling hydrological time series." Hydrological Sciences Journal, Vol. 51, No. 4, pp. 588-598.
10.1623/hysj.51.4.588
21
Kim, S., and Kim, H.S. (2008). "Neural networks and genetic algorithm approach for nonlinear evaporation and evapotranspiration modeling." Journal of Hydrology, Vol. 351, No. 3, pp. 299-317.
10.1016/j.jhydrol.2007.12.014
22
Kisi, O. (2011). "A combined generalized regression neural network wavelet model for monthly streamflow prediction." KSCE Journal of Civil Engineering, Vol. 15, No. 8, pp. 1469-1479.
10.1007/s12205-011-1004-4
23
Krishna, B., Satyaji Rao, Y.R., and Nayak, P.C. (2011). "Time series modeling of river flow using wavelet neural networks." Journal of Water Resource and Protection, Vol. 3, pp. 50-59.
10.4236/jwarp.2011.31006
24
Labat, D., Ababou, R., and Mangin, A. (2000). "Rainfall-runoff relations for karstic springs. Part II: Continuous wavelet and discrete orthogonal multiresolution analyses." Journal of Hydrology, Vol. 238, pp. 149-178.
10.1016/S0022-1694(00)00322-X
25
Liu, W.C., and Chung, C.E. (2014). "Enhancing the predicting accuracy of the water stage using a physical-based model and an artificial neural network-genetic algorithm in a river system." Water, Vol. 6, No. 6, pp. 1642-1661.
10.3390/w6061642
26
Lohani, A.K., Goel, N.K., and Bhatia, K.K.S. (2006). "Takagi-Sugeno fuzzy inference system for modeling stage-discharge relationship." Journal of Hydrology, Vol. 331, pp. 146-160.
10.1016/j.jhydrol.2006.05.007
27
Mallat, S.G. (1989). "A theory for multiresolution signal decomposition: The wavelet representation." IEEE Transactions on Pattern Analysis and Machine Intelligence, Vol. 11, No. 7, pp. 674-693.
10.1109/34.192463
28
Mamdani, E.H., and Assilian, S. (1975). "An experiment in linguistic synthesis with a fuzzy logic controller." International Journal of Man-machine Studies, Vol. 7, No. 1, pp. 1-13.
10.1016/S0020-7373(75)80002-2
29
Martins, O.Y., Sadeeq, M.A., and Ahaneku, I.E. (2011). "ARMA modelling of Benue river flow dynamics: Comparative study of PAR model." Open Journal of Modern Hydrology, Vol. 1, pp. 1-9.
10.4236/ojmh.2011.11001
30
Moosavi, V., Vafakhah, M., Shirmohammadi, B., and Behnia, N. (2013). "A wavelet-ANFIS hybrid model for groundwater level forecasting for different prediction periods." Water Resources Management, Vol. 27, pp. 1301-1321.
10.1007/s11269-012-0239-2
31
Moradkhani, H., Hsu, K.L., Gupta, H.V., and Sorooshian, S. (2004). "Improved streamflow forecasting using self-organizing radial basis function artificial neural networks." Journal of Hydrology, Vol. 295, No. 1, pp. 246-262.
10.1016/j.jhydrol.2004.03.027
32
Nash, J.E., and Sutcliffe, J.V. (1970). "River flow forecasting through conceptual models I: A discussion of principles." Journal of Hydrology, Vol. 10, pp. 282-290.
10.1016/0022-1694(70)90255-6
33
Nasr, A., and Bruen, M. (2008). "Development of neuro-fuzzy models to account for temporal and spatial variations in a lumped rainfall-runoff model." Journal of Hydrology, Vol. 349, pp. 277-290.
10.1016/j.jhydrol.2007.10.060
34
Nourani, V., Alami, M.T., and Aminfar, M.H. (2009). "A combined neural-wavelet model for prediction of Ligvanchai watershed precipitation." Engineering Applications of Artificial Intelligence, Vol. 22, No. 3, pp. 466-472.
10.1016/j.engappai.2008.09.003
35
Nourani, V., Baghanam, A.H., Adamowski, J., and Kisi, O. (2014). "Applications of hybrid wavelet-artificial intelligence models in hydrology: A review." Journal of Hydrology, Vol. 514, pp. 358-377.
10.1016/j.jhydrol.2014.03.057
36
Nourani, V., Hosseini, B.A., Adamowski, J., and Gebremicheal, M. (2013). "Using self-organizing maps and wavelet transforms for space-time pre-processing of satellite precipitation and runoff data in neural network based rainfall-runoff modeling." Journal of Hydrology, Vol. 476, pp. 228-243.
10.1016/j.jhydrol.2012.10.054
37
Partal, T. (2009). "Modelling evapotranspiration using discrete wavelet transform and neural networks." Hydrological Processes, Vol. 23, No. 25, pp. 3545-3555.
10.1002/hyp.7448
38
Partal, T., and Kisi, O. (2007). "Wavelet and neuro-fuzzy conjunction model for precipitation forecasting." Journal of Hydrology, Vol. 342, pp. 199-212.
10.1016/j.jhydrol.2007.05.026
39
Quilty, J., and Adamowski, J. (2018). "Addressing the incorrect usage of wavelet-based hydrological and water resources forecasting models for real-world applications with best practices and a new forecasting framework." Journal of Hydrology, Vol. 563, pp. 336-353.
10.1016/j.jhydrol.2018.05.003
40
Ravansalar, M., Rajaee, T., and Kisi, O. (2017). "Wavelet-linear genetic programming: A new approach for modeling monthly streamflow." Journal of Hydrology, Vol. 549, pp. 461-475.
10.1016/j.jhydrol.2017.04.018
41
Rezaie-Balf, M., Kim, S., Fallah, H., and Alaghmand, S. (2019). "Daily river flow forecasting using ensemble empirical mode decomposition based heuristic regression models: Application on the perennial rivers in Iran and South Korea." Journal of Hydrology, Vol. 572, pp. 470-485.
10.1016/j.jhydrol.2019.03.046
42
Sahay, R.R., and Srivastava, A. (2014). "Predicting monsoon floods in rivers embedding wavelet transform, genetic algorithm and neural network." Water Resources Management, Vol. 28, pp. 301-317.
10.1007/s11269-013-0446-5
43
Samsudin, R., Saad, P., and Shabri, A. (2011). "River flow time series using least squares support vector machines." Hydrology and Earth System Sciences, Vol. 15, pp. 1835-1852.
10.5194/hess-15-1835-2011
44
Seo, Y., Kim, S., Kisi, O., and Singh, V.P. (2015). "Daily water level forecasting using wavelet decomposition and artificial intelligence techniques." Journal of Hydrology, Vol. 520, pp. 224-243.
10.1016/j.jhydrol.2014.11.050
45
Seo, Y., Kim, S., and Singh, V. (2018). "Machine learning models coupled with variational mode decomposition: A new approach for modeling daily rainfall-runoff." Atmosphere, Vol. 9, No. 7, p. 251.
10.3390/atmos9070251
46
Shoaib, M., Shamseldin, A,Y., and Melville, B,W. (2014). "Comparative study of different wavelet based neural network models for rainfall-runoff modeling." Journal of Hydrology, Vol. 515, pp. 47-58.
10.1016/j.jhydrol.2014.04.055
47
Sivakumar, B., Berndtsson, R., Olsson, J., and Jinno, K. (2001). "Evidence of chaos in the rainfall-runoff process." Hydrological Sciences Journal, Vol. 46, No. 1, pp. 131-145.
10.1080/02626660109492805
48
Stone, M. (1974). "Cross-validatory choice and assessment of statistical predictions." Journal of the royal statistical society. Series B (Methodological), Vol. 36, No. 2, pp. 111-147.
10.1111/j.2517-6161.1974.tb00994.x
49
Takagi, T., and Sugeno, M. (1985). "Fuzzy identification of systems and its applications to modeling and control." IEEE Transactions on Systems, Man, and Cybernetics, Vol. 15, No. 1, pp. 116-132.
10.1109/TSMC.1985.6313399
50
Talei, A., Chye, C.L.H., and Wong, T.S.W. (2010). "Evaluation of rainfall and discharge inputs used by adaptive network-based Fuzzy Inference Systems (ANFIS) in rainfall-runoff modeling." Journal of Hydrology, Vol. 391, pp. 248-262.
10.1016/j.jhydrol.2010.07.023
51
Tiwari, M.K., and Chatterjee, C. (2010). "Development of an accurate and reliable hourly flood forecasting model using wavelet-bootstrap-ANN (WBANN) hybrid approach." Journal of Hydrology, Vol. 394, pp. 458-470.
10.1016/j.jhydrol.2010.10.001
52
Zakhrouf, M., Bouchelkia, H., and Stamboul, M. (2016). "Neuro-wavelet (WNN) and neuro-fuzzy (ANFIS) systems for modeling hydrological time series in arid areas. A case study: The catchment of Aïn Hadjadj (Algeria)." Desalination and Water Treatment, Vol. 57, No. 37, pp. 17182-17194.
10.1080/19443994.2015.1085908
53
Zakhrouf, M., Bouchelkia, H., Stamboul, M., Kim, S., and Heddam, S. (2018). "Time series forecasting of river flow using an integrated approach of wavelet multi-resolution analysis and evolutionary data-driven models. A case study: Sebaou River (Algeria)." Physical Geography, Vol. 39, No. 6, pp. 506-522.
10.1080/02723646.2018.1429245
54
Zakhrouf, M., Bouchelkia, H., Stamboul, M., and Kim, S. (2020). "Novel hybrid approaches based on evolutionary strategy for streamflow forecasting in the Chellif River, Algeria." Acta Geophysica, Vol. 68, No. 1, pp. 167-180.
10.1007/s11600-019-00380-5
55
Zhang, G., Patuwo, B.E., and Hu, M.Y. (1998). "Forecasting with artificial neural networks: The state of the art." International Journal of Forecasting, Vol. 14, No. 1, pp. 35-62.
10.1016/S0169-2070(97)00044-7
56
Zhao, X., Guo, X., Luo, J., and Tan, X. (2018). "Efficient detection method for foreign fibers in cotton." Information Processing in Agriculture, Vol. 5, No. 3, pp. 320-328.
10.1016/j.inpa.2018.04.002


text:
## Self-sustained emission in semi-infinite non-Hermitian systems at the exceptional point

X. Z. Zhang, L. Jin, Z. Song

Complex potential and non-Hermitian hopping amplitude are building blocks of a non-Hermitian quantum network. Appropriate configuration, such as PT-symmetric distribution, can lead to the full real spectrum. To investigate the underlying mechanism of this phenomenon, we study the phase diagram of a semi-infinite non-Hermitian system. It consists of a finite non-Hermitian cluster and a semi-infinite lead. Based on the analysis of the solution of the concrete systems, it is shown that it can have the full real spectrum without any requirements on the symmetry and the wave function within the lead becomes a unidirectional plane wave at the exceptional point. This universal dynamical behavior is demonstrated as the persistent emission and reflectionless absorption of wave packets in the typical non-Hermitian systems containing the complex on-site potential and non-Hermitian hopping amplitude.

http://arxiv.org/abs/1212.0086

Quantum Physics (quant-ph)

## Hermitian scattering behavior for the non-Hermitian scattering center

L. Jin, Z. Song

We study the scattering problem for the non-Hermitian scattering center, which consists of two Hermitian clusters with anti-Hermitian couplings between them. Counterintuitively, it is shown that it acts as a Hermitian scattering center, satisfying $$|r| ^{2}+|t| ^{2}=1$$, i.e., the Dirac probability current is conserved, when one of two clusters is embedded in the waveguides. This conclusion can be applied to an arbitrary parity-symmetric real Hermitian graph with additional PT-symmetric potentials, which is more feasible in experiment. Exactly solvable model is presented to illustrate the theory. Bethe ansatz solution indicates that the transmission spectrum of such a cluster displays peculiar feature arising from the non-Hermiticity of the scattering center.

http://arxiv.org/abs/1109.2187
Quantum Physics (quant-ph)

## Hermitian dynamics in a class of pseudo-Hermitian networks

L. Jin, Z. Song

We investigate the connection between pseudo-Hermitian and Hermitian descriptions for a lattice, which consists of a set of isomorphic pseudo-Hermitian clusters. We show that such non-Hermitian systems can act as Hermitian systems. This is made possible by considering the dynamics of a state involving the superposition of a single eigenmode of each isomorphic clusters. It still holds when multiple eigenmodes are involved when additional restriction on the state is imposed. This Hermitian dynamics is demonstrated for the case of an exactly solvable PT-symmetric ladder system.

http://arxiv.org/abs/1107.4311
Quantum Physics (quant-ph)

## Perfect State Transfer in PT-symmetric Non-Hermitian Networks

X. Z. Zhang, L. Jin, Z. Song

We systematically study the parity- and time-reversal- (PT) symmetric non-Hermitian version of a quantum network proposed in the work of Christandl et al. [Phys. Rev. Lett. 92, 187902 (2004)]. We show that such an extended model still allows conditional perfect state transfer within the unbroken PT-symmetric region, but not arbitrary. This is due to the fact that the evolution operator at certain period is equivalent to the PT operator for the real-valued wavefunction in the elaborate PT-symmetric Hilbert space. The critical behavior of the eigenstates is also discussed at the exceptional point.

http://arxiv.org/abs/1106.0087
Quantum Physics (quant-ph)

## Phase diagram and scaling behavior of a PT symmetric non-Hermitian Bose-Hubbard system

L. Jin, Z. Song

We study the phase diagram of a PT symmetric non-Hermitian Bose-Hubbard model. Using both analytical and numerical approaches, we find that even small on-site interaction can break the PT symmetry drastically. It has been demonstrated that the scaling law can be established for the exceptional point in small U limit. Based on numerical approach, we also find that the phase diagram shows rich structure for medium U: there exist multiple regions in which the PT symmetry is unbroken.

http://arxiv.org/abs/1105.6186
Statistical Mechanics (cond-mat.stat-mech); Quantum Physics (quant-ph)

## A physical interpretation for the non-Hermitian Hamiltonian

L. Jin, Z. Song

We explore a way of finding the link between a non-Hermitian Hamiltonian and a Hermitian one. Based on the analysis of Bethe Ansatz solutions for a class of non-Hermitian Hamiltonians and the scattering problems for the corresponding Hermitian Hamiltonians. It is shown that a scattering state of an arbitrary Hermitian lattice embedded in a chain as the scattering center shares the same wave function with the corresponding non-Hermitian tight binding lattice, which consists of the Hermitian lattice with two additional on-site complex potentials, no matter the non-Hermitian is broken PT symmetry or even non-PT. An exactly solvable model is presented to demonstrate the main points of this article.

http://arxiv.org/abs/1101.0351
Quantum Physics (quant-ph)

format: plain


text:
# Frequency modes of the rectangular shell

12 views

The shapes of three natural modes having the frequencies $\omega_1, \omega_2, \omega_3$ of the rectangular shell are presented in the figure. The exciting pressure $p(t)$ applied uniformly all over the one side of the shell has the form $p(t) = Pe^{jωt}$.
Make a sketch of the normal displacement of the gravity point of the shell against frequency, if the excitation frequency varies within bounds $0.5\omega_1< \omega <2\omega_1$ and static displacement of that point equals to $u_0$.

Links to photo of frequency nodes (sorry for low quality)->1 .

Can somebody help me and tell me how i should get started?

* The image is also not clear.  I guess for $\omega_1$ the only nodes are at the edges of the shell, for $\omega_2$ there is a node at the centre of the shell but I cannot make out the mode shape for $\omega_3$.


text:
# Can conservation of angular momentum be derived from energy consideration? [closed]

About the context of this question:
A dramatic illustration of back and forth conversion of kinetic energy and gravitational potential energy is when a celestial body has an highly eccentric orbit. Vivid example: Halley's comet.

As we know, at perihelion Halley's comet has far more velocity than at aphelion. As we know, the reason for that is that at perihelion Halley's comet has been falling towards the Sun for 38 years. For 38 years the gravity of the Sun has been doing work, increasing the kinetic energy of Halley's comet.

The conversion of gravitational potential energy to kinetic energy completely determines the shape and velocity profile of the comet's trajectory.

I have always wondered: is it possible to use the work-energy theorem to derive conservation of angular momentum? As I mentioned: in the case of a central force just the energy conversion completely determines the trajectory.

Of course, implicit in such an attempt is that the setup has rotational symmetry. The force is a central force, the force law is the same for every orientation.

Here is an attempt at such a derivation, but it falls short. The expression for kinetic energy that it uses is in terms of the angular velocity $$\omega$$, but in order for a celestial object to move from aphelion to perihelion there must be some radial velocity, that is inconsistent. The attempt does arrive at conservation of angular momentum, but well... maybe that's just a fluke.

So my question is: can the outline that is below be improved to the level of an actual derivation?

Rotational kinetic energy:

$$E_k = \frac{1}{2}mr^2\omega^2 \quad (1)$$

Required centripetal force at angular velocity $$\omega$$ and radial distance $$r$$:

$$F = -m\omega^2r \quad (2)$$

The Work-Energy theorem in differential form:

$$\frac{dE}{dr} = F \quad (3)$$

Substituting the expressions for $$E$$ and $$F$$ in the differential equation:

$$\frac{d(\frac{1}{2}r^2\omega^2)}{dr} = -\omega^2r \quad (4)$$

To show how this implies a relation between $$\omega$$ and $$r$$ the factor $$\omega^2$$ is substituted with $$\frac{C}{r^4}$$, ($$C$$ is a constant), and where applicable factors of $$r$$ are allowed to drop away against each other.

$$\frac{d(\frac{1}{2}\frac{C}{r^2})}{dr} = -\frac{C}{r^3} \quad (5)$$

So we have the following relation:

$$\omega^2 = \frac{C}{r^4} \quad (6)$$

Hence:

$$\omega^2r^4 = C \quad (7)$$

Which is of course a squared version of the expression that says that angular momentum $$\omega r^2$$ is a constant.

Pushing the attempt further:
I assume using polar coordinates is best in this case, because that aligns with the rotational symmetry of the setup.

Assumption: the total kinetic energy is the following sum:

$$E_{kin.total} = \tfrac{1}{2}m\left(\tfrac{dr}{dt}\right)^2 + \tfrac{1}{2}mr^2\omega^2 \qquad (8)$$

The component velocities of these contributions are perpendicular, so I expect these component velocities to sum as the sides of a right triangle.

Assumption:
The amount of centripetal force that is available for radial acceleration is the total centripetal force minus the required centripetal force.

The thing that is tripping me up here is the assignment of minus signs. What is needed is an expression with the total centripetal force on one side, and the radial acceleration and required centripetal force on the other side. Of course everything in the radial direction must be treated like a vector; the minus signs must be correct.

Assumption on how to apply the work-energy theorem:
The motion of a celestial body from aphelion to perihelion can be thought of as spiraling in. Compare that gradual motion to moving down a ramp. To apply the work-energy theorem it isn't necessary to know the shape of the ramp. The total change in energy is always equal to the total difference in energy potential from top to bottom of the ramp, independent of its shape.

So I think the derivation only needs to consider the change in potential energy as a function of the radial distance, no need for separate radial and tangential component there.

Pushing the attempt further:

The image illustrates the kind of case that this question is about.

Image credit:
University of Toronto, Department of Physics, Angular Momentum

Of course, the case will be treated in an idealized form, so that exclusively the principle is examined.

The idealization:
All motion is treated as frictionless.

Initial condition: the radial distance of the circumnavigating mass is kept constant. As we know, in the case of circular motion the magnitude of the required centripetal force is $$m\omega^2r$$ From here on I will use the expression 'required centripetal force' for 'magnitude of centripetal force required to sustain circular motion'.

In the next phase the centripetal force is increased. So in this fase there is a surplus of centripetal force, and the circular motion contracts.

As long as the contraction continues the centripetal force is doing work, increasing the rotational kinetic energy.

In the final state the increase of centripetal force is leveled off.
As we know, in the final state the angular velocity is much higher than in the initial state. The angular velocity is higher because the centripetal force has been doing work.

What is interesting here is that even though the centripetal force acts purely in radial direction, the end result is an increase of angular velocity. That is the nature of circumnavigating motion: any change of radial velocity transforms to angular velocity.

So: initial state is unambiguious; circular motion at a particular angular velocity, and the end state is unambiguous too: circular motion at a decreased radial distance, with an increased angular velocity.

We can calculate the increase of angular velocity in two ways:

• Apply the principle of conservation of angular momentum

• Solve analytically for the amount of work that the centripetal force has done when the system has contracted. The change in energy from initial state to end state is equal to the amount of work done, hence the end state is uniquely determined.

Those two independent calculations must end up with the same answer.

So: how to arrive at an analytical solution?

As noted: when the total force is larger than the required centripetal force there will be an acceleration in the radial direction. For the magnitudes of the total force, the radial acceleration and the required centripetal force respectively:

$$F_{total} = m\tfrac{d^2r}{dt^2} + m\omega^2r \qquad (9)$$

Repeating equation (3):

$$\frac{dE}{dr} = F \qquad (10)$$

Equation (10) expresses that in this case the force is at all times acting in radial direction, which means that all the work that is done is in the radial direction. Therefore the derivative of the total energy with respect to $$r$$ should be equal to the force.

The glaring problem with this approach: In equation (8) a decomposition in two perpendicular components is proposed, but in the case of the force the two contributions on the right side of equation (9) are not perpendicular components, those two contributions are aligned. This way of approaching the problem is neither purely scalar, nor purely vectorial, but something in between, which may well invalidate the whole approach.

Subsitituting (8) and (9) into (10) $$\frac{d(\tfrac{1}{2}m\left(\tfrac{dr}{dt}\right)^2 + \tfrac{1}{2}mr^2\omega^2)}{dr} = m\tfrac{d^2r}{dt^2} + m\omega^2r \qquad (11)$$

the following two equations are a separation: equation (11) is separated into two statements:

$$\frac{d(\tfrac{1}{2}m\left(\tfrac{dr}{dt}\right)^2)}{dr} = m\tfrac{d^2r}{dt^2} \qquad (12)$$

$$\frac{d(\tfrac{1}{2}mr^2\omega^2)}{dr} = m\omega^2r \qquad (13)$$

Equation (12) has on the left the derivative with respect to $$r$$ of the kinetic energy in radial direction, and on the right the acceleration in radial direction.
Assertion: the term on the left side and the term on the right side will be the same at every point in time; it's the work-energy theorem.
Assertion: equation (12) is satisfied at every point in time, that implies that equation (13) is constrained to be satisfied at every point in time too.

Equation (13) is very similar in form to equation (4), except in equation (4) the $$m\omega^2r$$ term has a minus sign.

Therefore in the case of equation (13) to accommodate the non-presence of a minus sign you would have to relax the constraint to a statement about magnitude; asserting that the magnitude of the term on the left side must be equal to the magnitude of the term on the right side. Presumably this minus-sign-problem arises because the whole approach is neither purely scalar, nor purely vectorial

Now to what would count as a useful stackexchange answer:
In this form this is a problem in mathematical physics. When a rotating system contracts, find a way to calculate the increase of rotational kinetic energy. The initial state and the end state are both unambiguous, therefore an unambiguous calculation must be possible. I have outlined my own attempt, and I have outlined the glaring problem with it: neither purely scalar nor purely vectorial. Can the outline be improved?

More generally:
Discussion of necessary and sufficient condition.
First the case of conservation of linear momentum, in the case of two masses that exert a force upon each other, which of course causes mutual acceleration. In that case a necessary condition for conservation of momentum is that the force law supports it. $$F=ma$$ supports conservation of linear momentum. By contrast: if instead of $$F=ma$$ the force law would be this:

$$a = \frac{F}{m^2}$$

Then you don't have conservation of linear momentum.

In the case of angular momentum it is of course a necessary condition is that the force that the masses exert upon each other is a central force, but that is not a sufficient condition. In addition the force law must support it, just as in the case of linear momentum.

## closed as unclear what you're asking by Aaron Stevens, GiorgioP, Jon Custer, Buzz, ZeroTheHeroApr 15 at 3:57

Please clarify your specific problem or add additional details to highlight exactly what you need. As it's currently written, it’s hard to tell exactly what you're asking. See the How to Ask page for help clarifying this question. If this question can be reworded to fit the rules in the help center, please edit the question.

• The work-kinetic energy theorem is basically a slightly different form of $F=ma$, so you're asking if you can derive conservation laws given the equation of motion. And the answer is, of course you can! – knzhou Apr 14 at 0:04
• @knzhou Yeah, it's about $F = ma$. So, in the case of conservation of linear momentum when two objects collide in a perfectly elastic collision: then conservation of linear momentum in that forceful interaction can be regarded as a consequence of $F = ma$. If the two objects have unequal mass then the least massive object will undergo the largest acceleration. Like collision, orbital mechanics is two objects, both accelerating each other. – Cleonis Apr 14 at 0:30
• Perhaps another way to say this is: conservation laws are just intuitive and useful ways of summarizing some of the information that comes out of $F=ma$. If you have all the $F=ma$ (and similar) equations, then it's not surprising you get everything. – knzhou Apr 14 at 1:18
• Your talk about total, required, and available centripetal force doesn't make sense to me. Honestly I'm trying hard to put something together here, but anything I try just ends up being circular or conservation of angular momentum pops up without energy being directly tied in. – Aaron Stevens Apr 14 at 4:36
• The question has been put on hold. Independent of that I had already come to the conclusion that this question is not suitable for stackexchange. I still think the question is an interesting one to raise, but stackexchange isn't a good place for it. So I intend to delete it, in a day or so. My thanks to knzhou, for pointing out that it is about $F=ma$, the work-energy theorem coming in because the work-energy theorem is equivalent with $F=ma$ – Cleonis Apr 15 at 16:26

I did not read the whole of the derivation, but

The force is a central force

implies conservation of angular momentum with respect to the center the force is directed to. So a derivation that starts from that assumption is not a derivation of conservation of angular momentum from the work-energy theorem.

• I misread it too, but the OP isn't asking about energy conservation. – Aaron Stevens Apr 13 at 22:26
• It seems to me that the condition of rotational symmetry in itself only constrains that something angular must be conserved. Additional steps are required to arrive at conservation of $\omega r^2$. It seems to me there are multiple ways to arrive at conservation of $\omega r^2$, that is why I think exploring the a route via energy conversion can be interesting. If possible then that would be an interesting illustration of the interconnectedness of energy considerations and momentum considerations – Cleonis Apr 13 at 22:31
• As far as I understand, he is asking whether it is possible to use energy conservation to derive the conservation of angular momentum: starts with "is it possible to use the work-energy theorem to derive conservation of angular momentum" and ends with "Which is of course a squared version of the expression that says that angular momentum $ωr^2$ is a constant." and asks "can the outline that is below be improved to the level of an actual derivation?". The answer is no, with a central force angular momentum is conserved, independently of anything else. – JTS Apr 13 at 22:32
• @Cleonis: here you are a case in which energy is not conserved and angular momentum is conserved: the force has a component proportional to the component of velocity directed towards the center. – JTS Apr 13 at 22:33
• @JTS The work-energy theorem isn't energy conservation – Aaron Stevens Apr 13 at 22:35

In general this is not possible: Conservation of energy and conservation of angular momentum pertain to two distinctly different symmetries in the Lagrangians describing physical systems. And thanks to Noether's theorem, we know that there is a 1-to-1 correspondence between symmetries and conserved quantities.

• There seems to be a confusion here. My question hinges on energy conversion (in this case between potential energy and kinetic energy). This question does not necessarily involve energy conservation – Cleonis Apr 14 at 1:29
• Your question title speaks about conservation though. – AtmosphericPrisonEscape Apr 14 at 1:30
• Also you get an important thing wrong in your question: The energy landscape does not determine the trajectory in the Kepler problem fully. The other free parameter is the angular momentum. – AtmosphericPrisonEscape Apr 14 at 1:32
• The title says 'energy consideration', which unfortunately is easily misread as energy conservation. Too easy to misread; in retrospect not a good choice of word, in the title. – Cleonis Apr 14 at 1:43
• About how much is determined by the energy landscape: given a particular radial distance at aphelion, and a particular amount of kinetic energy at aphelion, it seems to me the subsequent motion is completely determined. Of course, if the starting specification is only a total amount of energy then additional information is needed to narrow down to one particular trajectory. Anyway, these details do not materially affect the purpose of the question – Cleonis Apr 14 at 2:24

Your "proof" is circular. You substitute in $$\omega^2=C/r^4$$ and show that $$\omega^2r^4=C$$

You assume angular momentum is conserved in your derivation then.

• I have edited the question to add enumeration to the expressions. About the subsitution from exp. (4) to exp. (5): it is not the case that any substitution will result in a satisfied equation. The point is that eq. (4) imposes a constraint. There is only one subsitution that results in a satisfied equation; in that sense that particular substitution is enforced. So the suggestion of circular reasoning is incorrect. (Of course, the outline does fall short, but for the reason given in the original question.) – Cleonis Apr 13 at 22:58
• @Cleonis No it is circular. You will always get a true expression no matter what you replace $\omega$ with. This is because the left and right sides of equation 4 are exactly the same thing. Just take the derivative and you will see this. Therefore you are in fact assuming angular momentum conservation and then getting angular momentum conservation. The angular momentum conservation doesn't arise from something else. I could argue $\omega^2r^5$ is also conserved using your arguments – Aaron Stevens Apr 13 at 23:06
• Well, expression (4) is not one where a partial derivative is taken. When you don't do a subsitution then writing out the derivative means also writing out the derivative of the $\omega$ term with respect to $r$ – Cleonis Apr 13 at 23:16
• @Cleonis Oh. By saying $F_{rad}=-m\omega^2r$ you are assuming a constant $r$ and therefore constant $\omega$. In general the radial component of the force is given by $F_{rad}=m\ddot r-m\omega^2r$ where $\ddot r$ is the second rate of change of $r$. I thought you were assuming uniform circular motion by using this expression. – Aaron Stevens Apr 13 at 23:39
• @Cleonis And of course you have to add in the energy the radial part of the motion too. I think it will not be possible to deduce how to split the energy into the parts that are due to radial and tangential velocity using the work-energy theorem alone, as in the work-energy theorem (scalar) you lose the information about the vectorial character of motion. The correct splitting is what will give you the conservation of angular momentum. This is for the moment in part guessing but it might be correct and also the thing that I could add to my answer. Need to think about it a bit more. – JTS Apr 13 at 23:49


text:
# How did the Kepler Telescope reach its orbit?

The Kepler Telescope was (is) in a heliocentric Earth trailing orbit, but I want to figure out how it was able to reach this orbit exactly.

It was launched on a Delta II with 3 stages and I studied the different burns and the parking orbit (picture below). After coasting in the parking orbit, the second and third stage put the telescope in its appropriate orbit. Still I find it difficult to reconstruct how exactly it went from LEO Parking orbit to a heliocentric orbit approximately the same as Earth and trailing it. How exactly did those 2 burns put it in this orbit?

A professor helped me by stating it made several widening ellipses from LEO to its orbit (she gave me this picture). But how does this work with just 2 burns?

Can somebody show me how the trajectory from LEO went and where the burns took place?

Thanks!

## 2 Answers

Going from the diagram you presented, it seem that after the first of those two burns (the second burn of the second stage) if was in an elliptical Earth orbit ("94 x 1180 n.mi."). A minute or so later, still pretty much at the perigee of that orbit, the third stage lit up and burned for 90 seconds. At the end of that it was in a hyperbolic Earth orbit with perigee 99 n.mi. On that orbit it escapes from Earth, slowing as it gets further away. Had Earth been all alone in the Universe, it would have approached a velocity of $$\sqrt{0.61} = 0.78 km/s$$ as the distance to Earth approached infinity. In fact though, the influence of the Sun becomes the dominant consideration once you get more than a million km or so from Earth, and from that perspective, it was now in a similar orbit to Earth's, modified by that 0.78 km/s, and the direction was chosen so that this was an ellipse a little larger than Earth's orbit, which takes a little longer, so that it slowly dropped behind Earth.

• Thank you very much, it really helped! Quick question, the 94 nmi is the distance from Earths surface to the perigee and 1180 nmi the distance from Earths surface to the apogee, right? Dec 5 '20 at 15:44
• And is the second burn of the second stage, which creates an elliptical orbit, just used to bring the launch vehicle up to a higher speed? And afterwards the third stage gives it its escape velocity? Isn't this possible in one burn from LEO to escape velocity? Dec 5 '20 at 19:13
• @veeke After the second burn they had to pause to separate the second stage. Those two burns are with different engines. Otherwise one longer burn would be optimal Dec 5 '20 at 22:56
• If I'm not mistaken. The last burn (third stage) put the spacecraft in it's hyperbolic escape trajectory from Earth. But doesn't their need to be another burn to establish the heliocentric orbit with it's specific parameters? Dec 8 '20 at 17:20
• @Veeke There may have been some tiny tweaks, but the speed and direction of the Earth escape orbit are chosen to lead to the desired heliocentric orbit. The deeper in Earth's gravity well you can do your burns, the more efficient it is. Dec 8 '20 at 18:46

If they were indeed quoted correctly here, then in that case your professor is not exactly right.

A plot of the first few days (not shown) shows a featureless departure from Earth to heliocentric orbit, no spiraling.

### Nonetheless, there be spirals!

Depending on what you choose for the center of your plot and whether your frame rotates once per year or is "inertial1" it can absolutely look like a spiral around the Earth!

And I can see why one might assume, having seen the spirally pattern and not realize it always looks like that in some frame, that one might think that it looks like a spacecraft spiraling away from Earth, perhaps under ion propulsion.

But it's just the way to bodies in similar orbits look if their periods are very slightly different and one or both are slightly elliptical.

Plots of the two STEREO spacecraft will look something like this too, except there's a pair of counter-spiraling orbits since one leads and one lags.

In the plots distances are in kilometers, and the data is from JPL's Horizons every 6 hours since launch 2009-03-07 04:57.

1not a rotating frame but still accelerating because origin moves with the Earth

import numpy as np
import matplotlib.pyplot as plt

class Body(object):
def __init__(self, name):
self.name = name

def rotz(vecs, th):
x, y, z = vecs
cth, sth = np.cos(th), np.sin(th)
xrot = x * cth - y * sth
yrot = y * cth + x * sth
return np.vstack((xrot, yrot, z))

fnames = ('Kepler Sun full horizons_results.txt',
'Kepler Earth full horizons_results.txt',
'Kepler Kepler full horizons_results.txt')

bodies = []
for fname in fnames:

with open(fname, 'r') as infile:
lines = infile.read().splitlines()

iSOE = [i for i, line in enumerate(lines) if "$$SOE" in line][0] iEOE = [i for i, line in enumerate(lines) if "$$EOE" in line][0]

print(iSOE, iEOE, lines[iSOE], lines[iEOE])
lines = [line.split(',') for line in lines[iSOE+1:iEOE]]
JD  = np.array([float(line[0]) for line in lines])
pos = np.array([[float(item) for item in line[2:5]] for line in lines])
vel = np.array([[float(item) for item in line[5:8]] for line in lines])

b = Body(fname.split()[0])
b.JD = JD
b.pos = pos.T.copy()
b.vel = vel.T.copy()

bodies.append(b)

sun, earth, kepler = bodies

days = kepler.JD - kepler.JD[0]

x, y, z = earth.pos - sun.pos
theta = np.arctan2(y, x)

for body in bodies:
body.posr = rotz(body.pos, -theta)

if True:
plt.figure()
plt.subplot(2, 2, 1)
x, y, z = kepler.pos - earth.pos
plt.plot(x, y, '-r', linewidth=1.0)
plt.plot([0], [0], 'ob')
plt.gca().set_aspect('equal')
plt.title('Kepler minus Earth "inertial"')

plt.subplot(2, 2, 2)
x, y, z = kepler.posr - earth.posr
plt.plot(x, y, '-r', linewidth=1.0)
plt.plot([0], [0], 'ob')
plt.gca().set_aspect('equal')
plt.title('Kepler minus Earth rotating')

plt.subplot(2, 2, 3)
x, y, z = earth.pos - sun.pos
plt.plot(x, y, '-b', linewidth=1.5)
x, y, z = kepler.pos - sun.pos
plt.plot(x, y, '-r', linewidth=1.0)
plt.plot([0], [0], 'oy')
plt.gca().set_aspect('equal')
plt.title('Kepler & Earth  minus Sun "inertial"')

plt.subplot(2, 2, 4)
x, y, z = earth.posr - sun.posr
plt.plot(x, y, '-b', linewidth=2.5)
x, y, z = kepler.posr - sun.posr
plt.plot(x, y, '-r', linewidth=1.0)
plt.plot([0], [0], 'oy')
plt.gca().set_aspect('equal')
plt.title('Kepler & Earth  minus Sun rotating')
plt.show()

• Thank you very much! I'm trying to reconstruct your python code in matlab. How did you retrieve the data from the sun, in what reference frame? Dec 6 '20 at 11:22
• @Veeke I've included a screen shot of the setup in Horizons, but I'll type that out in plaintext when I get to a real keyboard. I think it's the Horizon defaults; barycenter origin, ecliptic J2000.0 etc... You could also consider using this as an opportunity to try out Python a little bit; if you do a regular Anaconda installation (not mini-conda) you will have an incredible array of scientific and technical libraries ready to go, and of course everything is free forever! :-)
– uhoh
Dec 6 '20 at 14:00
• Sorry if I'm asking stupid things. But your screenshot setup only gives data from Kepler relative to the Sun, right? In your code you use sun.pos to calculate theta if I'm not mistaken. Dec 6 '20 at 14:00
• @Veeke no problem! The setup for coordinate origin shows Kepler's state vectors relative to to Solar System Barycenter, which spends about half its time inside the Sun and half just outside it. I could have specified coordinates relative to the Sun itself but I always use SSB to avoid getting them mixed up. For these plots it wouldn't make much difference.
– uhoh
Dec 6 '20 at 14:08
• Got it! big thanks for the help!! Dec 6 '20 at 15:52


text:
# New on CTAN: frege

• From: CTAN Announcements <ctan-ann@xxxxxxxx>
• Date: Sun, 15 Jul 2012 17:16:37 +0200 (CEST)

On Sat, 14 Jul 2012 Quirin Pamp submitted the

frege

package.

Summary description: Latex package for typesetting Frege's Begriffsschrift

Announcement text:
----------------------------------------------------------------------

The package frege defines a number of new commands for typesetting fregean
Begriffsschrift in LaTeX.

It is loosely based on the package begriff with the same purpose. However it
offers a number of improvements including better relative lengths of the content
stroke with respect to other strokes, content strokes that point at the middle
of lines rather than the bottom, a greater width for the assertion stroke as
compared to the content stroke, a more intuitive structure for the conditional,
greater care taken to allow for the linewidth in the spacing of formulas.

----------------------------------------------------------------------

This package is located at
http://mirror.ctan.org/macros/latex/contrib/frege
http://www.ctan.org/pkg/frege
(if the package is new it may take a day for that information to
appear). We are supported by the TeX Users Group http://www.tug.org .
Please join a users group; see http://www.tug.org/usergroups.html .

------------------------------------------------------------------------

For the CTAN Team
Rainer Schöpf

## Relevant Pages

... You need a fresh cvs update from today since I just fixed one bug. ... Updated the MacOSX build but not the other. ... package require tkpng ... You can also style the text rendering with stroke and fill colors, ...
(comp.lang.tcl)
• Re: Using a Canon EOS-30D for deepsky?
... If you think that Pete's posting simply to stroke his ego, ... It cames in the same package. ...
(uk.sci.astronomy)
• Re: A LaTeX challenge
... HO> * package `collect' ... HO> * package `environ' ... Good luck ... Any man can work when every stroke of his hand brings down the fruit ...
(comp.text.tex)


text:
# Wydawnictwa / Czasopisma IMPAN / Colloquium Mathematicum / Wszystkie zeszyty

## The spectrum of characters of ultrafilters on $\omega$

### Tom 111 / 2008

Colloquium Mathematicum 111 (2008), 213-220 MSC: Primary 03E05, 03E17. DOI: 10.4064/cm111-2-5

#### Streszczenie

We show the consistency of the statement: “the set of regular cardinals which are the characters of ultrafilters on $\omega$ is not convex”. We also deal with the set of $\pi$-characters of ultrafilters on $\omega$.

#### Autorzy

• Saharon ShelahThe Hebrew University of Jerusalem
Einstein Institute of Mathematics
Edmond J. Safra Campus, Givat Ram
Jerusalem 91904, Israel
and
Department of Mathematics
Hill Center, Busch Campus
Rutgers, The State University of New Jersey


text:
## Vertex operators in conformal field theory on $${\mathbb{P}}^ 1$$ and monodromy representations of braid group.(English)Zbl 0661.17021

Conformal field theory and solvable lattice models, Symp., Kyoto/Jap. 1986, Adv. Stud. Pure Math. 16, 297-372 (1988).
[For the entire collection see Zbl 0646.00016.]
The paper is a very clear exposition of the theory of vertex operators in two-dimensional conformal field theory, and the relation with monodromy representations of the braid group. The ideas presented here originate from a paper by V. G. Knizhnik and A. B. Zamolodchikov [Nucl. Phys. B 247, 83-103 (1984; Zbl 0661.17020)], but the present paper goes much deeper in the mathematical formulation and the foundations of the theory. For people who wish not to read such a long paper, there is a 10- page introduction listing the notation and all the basic theorems proved here.
The paper gives a summary of representation theory of the affine sl(2,$${\mathbb{C}})$$ Lie algebra, denoted by $$A_ 1^{(1)}$$, including the natural action of the operators of the Virasoro algebra. Then, the existence and uniqueness of vertex operators of spin j is shown: vertex operators are defined as operators satisfying the gauge condition and the equations of motion. It is shown that a triplet $$\left( \begin{matrix} j\\ j_ 2\quad j_ 1\end{matrix} \right)$$, satisfying the $$\ell$$-constrained Clebsch-Gordan condition, determines uniquely a vertex operator. Next, N- point functions are defined as compositions of vertex operators. Differential equations (called the fundamental equations) satisfied by these N-point functions, which have only regular singularities, give rise to properties of the vertex operators. First, convergence of compositions of vertex operators is proven. The commutation relation of vertex operators is equivalently rephrased in terms of the connection matrix of the fundamental equations, and is explicitly calculated in the case where one of the j-values is equal to $$1/2$$. Then, the monodromies of the fundamental equations give rise to representations of the braid group $$B_ N$$. The monodromy representation is explicitly determined in a special case (where implicitly $$j=$$ vertex operators are involved). This representation is in fact an irreducible representation of the Hecke algebra $$H_ N(q)$$ of type $$A_{N-1}$$, with $$q=\exp (2\pi \sqrt{- 1}/(\ell +2))$$. It is remarkable that the representation obtained here is with q a root of unity, since the Hecke algebra $$H_ N(q)$$ is not semi- simple in that case.
Reviewer: J.Van der Jeugt

### MSC:

 17B65 Infinite-dimensional Lie (super)algebras 81T40 Two-dimensional field theories, conformal field theories, etc. in quantum mechanics 20F36 Braid groups; Artin groups

### Citations:

Zbl 0646.00016; Zbl 0661.17020


text:
# 6.5. Statistical significance¶

Note

This page is very concise and avoids derivations. The focus here is on a working exposure to the topic. The corresponding lecture will add intuition.

Previously, we estimated

$\log(\text{price})= 8.2+ 1.53 *\log(\text{carat}) + 0.33* \text{Ideal} + 0.18* \log(\text{carat})\cdot \text{Ideal}$

Those coefficients are estimates, not gospel. They come from the sample of data we have. There is some uncertainty about what the “true” value of the coefficients is in the unseen “population.”

So, when we run the regression, it would be nice to get some extra info about those estimates. sm_ols does just that. For that regression (click the + sign to see the code) we get the following info:

# load some data to practice regressions
import seaborn as sns
import numpy as np
from statsmodels.formula.api import ols as sm_ols # need this

# this alteration is not strictly necessary to practice a regression
# but we use this in livecoding
diamonds2 = (diamonds.query('carat < 2.5')               # censor/remove outliers
.assign(lprice = np.log(diamonds['price']))  # log transform price
.assign(lcarat = np.log(diamonds['carat']))  # log transform carats
.assign(ideal = diamonds['cut'] == 'Ideal')

# some regression packages want you to explicitly provide
# a variable for the constant
.assign(const = 1)
)

(
sm_ols('lprice ~ lcarat + ideal + lcarat*ideal',
data=diamonds2.query('cut in ["Fair","Ideal"]'))
).fit().summary().tables[1] # the summary() is multiple tables stitched together. I only care to print the params here.

coef std err t P>|t| [0.025 0.975] 8.1954 0.007 1232.871 0.000 8.182 8.208 0.3302 0.007 46.677 0.000 0.316 0.344 1.5282 0.015 103.832 0.000 1.499 1.557 0.1822 0.015 12.101 0.000 0.153 0.212

The information next to each coefficient in this table is about the level of precision associated with the estimated coefficients.

Column

Meaning

“std err”

The standard error of the coefficient estimate. This gives an indication of how much the estimated coefficient likely varies from the population coefficient. There are several ways to compute a standard error, and the choice is important! However, it’s beyond the scope of this class, and we will use the default option throughout.

“t”

The “t-stat” = $$\beta$$ divided by its standard error.

“P>t”

The “p-value” is the probability that the coefficient is different than zero by random chance.

“[0.025 0.975]”

The 95% confidence interval for the coefficient.

Important

We use these columns, particularly the “t-stat” and “p-value”, to assess the probability that the coefficient is different from zero by random chance.

• A t-stat of 1.645 corresponds to a p-value of 0.10; meaning only 10% of the time would you get that coefficient randomly

• A t-stat of 1.96 corresponds to a p-value of 0.05; this is a common “threshold” to say a “relationship is statistically significant” and that “the relationship between X and y is not zero”

• A t-stat of 2.58 corresponds to a p-value of 0.01

## 6.5.1. Some practical guidance¶

1. You can focus on the p-values rather than the t-stats. However, knowing the threshold values of the t-stat (those above) is useful, as many people distribute research by discussing t-stats instead.

2. If a p-value for the coefficient in a regression is below 0.05,

• We say that the relationship between that variable and Y is “statistically significant” at the 5% level.

• Now, consider the direction: Is the coefficient positive or negative? Does this align with your intuition and economic theory?

• Now, consider the “size” of the coefficient: Is a “reasonable” change in X leading to a “small” or “big” change in Y? If the relationship is small, it may not be important enough to care about even if true.

• Now, take a step back: Statistically significant does NOT mean X causes Y. You need additional information to make these claims. That’s what the next page is about.

3. Practically, if the p-value is above 0.05, most researchers consider completely disregard the coefficient (and ignore the sign and the value). Because you can’t say the coefficient is statistically distinguishable from zero, they basically the interpret the coefficient for that variable as being zero. Meaning: Ignore the sign, ignore the value, assume the coefficient is zero.


text:
0
101 Sep 08, 2005 at 22:23

Hi to all,

I need to introduce HDR textures in my DX9 engine. Actually, IMHO the best texture format for HDR images is A16R16G16B16. The problem is that every texel occupies 8 bytes! An extremely huge quantity of memory.
Compression techniques are necessary. I looked at DX9 demos and I see a nice compression method called RGBE8. Basically, it rewrite every channel (r, g, b) in the form: m * 2\^e. m represent the mantissa of the channel and e is the exponent (stored in the alpha channel). Exponent is common to all the three channels.
This method requires 4 bytes per texel, that is an improvement over the naive one.
Actually I can extend this compression up to 1,5 bytes per texel in the following way: I compress RGB channels with DXT1 compression (0,5 bytes per texel). Then, I use another texture to store the exponent uncompressed (this increase the precision because exponent precision is crucial). Total: 1,5 bytes per texel.
A quite good compression.
However, there is a problem: this method works only if you use point filtering. Other filtering methods, like the bilinearing filter, make the things go wrong. :sad:
I have tried to apply point filtering to exponent texture (to fetch it correctly) and bilinearing in the RGB channels… but this still expose too artifacts.
I have another idea: splitting the 16 bit per channel texture in two 8 bit per channel textures (low byte & high byte texture) and apply DXT1 compression to both. Total: 1 byte per texel! Problem: an error in the “low byte” texture is not very noticeable. But an error in the “high byte” texture (because DXT1 is a lossy compression as you know) is quite devastating… HELL! :angry:
A solution to the problem is to store the low byte texture compressed with DXT1 and the high byte texture uncompressed. Total: 3,5 bytes per texel. :huh:
Here bilinearing filters don’t make problems.
Good quality, but It’s not a great compression… I will be happy if I can compress it a bit more (up to 2 bytes per texel or less :).
I’m very happy if you cool guys can suggest me others methods.

• AGPX

#### 22 Replies

0
101 Sep 09, 2005 at 07:40

It should work to use lower resolution for the exponent. I am not shure about the artifacs, but is should only be noticeable in high contrast edges.

0
101 Sep 09, 2005 at 20:56

@geon

It should work to use lower resolution for the exponent. I am not shure about the artifacs, but is should only be noticeable in high contrast edges. [snapback]21034[/snapback]

Actually I do the RGBE encoding with the following code:

int encode(float r, float g, float b,
float &encodedR,
float &encodedG,
float &encodedB)
{
const float maxComponent = max(max(r, g), b);
const int exp = (int)ceilf(log2f(maxComponent));
const float divisor = powf(2.0f, exp);
encodedR = r / divisor;
encodedG = g / divisor;
encodedB = b / divisor;
return exp;
}

So: r = encodedR * 2\^exp, g = encodedG * 2\^exp and b = encodedB * 2\^exp;

I found that it works better if I return exp * 16 (instead of exp), but in high contrast edge some Halo effects appear.

Anyway, I dunno if is convenient to do bilinear filtering manually in the pixel shader…
This sound too nasty and slow… :dry:

The following image is a shot from my editor. I use RGBE normally (without multiply by 16):

The following is with the exponent multiplication:

The last show halo and artifacts (again with multiplication by 16):

Exponent is clipped to stay in the range [-8, 7] when multiplied by 16. So it lie in the range [-128, 112]. I also sum 128 to remove the sign ([0, 240]). So it can stay in a single byte. Yes, dynamic range is less than previous, but images look really better.

I don’t understand the nature of the artifacts, at all. Halo should be due to bilinear filtering. The problem is that filtering should happen AFTER decoding and not before.

Is this what you mean for “use lower resolution for exponent”?

Thanks,

• AGPX
0
101 Sep 09, 2005 at 21:52

Hi again,

I have fixed the artifact problem. It’s due to DXT1 compression of the RGB channels (not the exponent). However, halo still occurs… (thanks to mikez for his hints).

• AGPX
0
101 Sep 09, 2005 at 21:53

@AGPX

Is this what you mean for “use lower resolution for exponent”?

Not at all…

As you showed in your code, you can encode several values with the same exponent (hence RGBE). There is nothing stoping you from encoding even more values to a single exponent.

A neat usage of this would be to store the color / exponent in a separate texture, as mentioned, BUT! why store so many of those exponents? If you take a block of say 4x4 pixels, the exponents should generally be virtually the same. (Except in very high contrast.) So; just store the exponent once and reuse it.

If you use this solution, you must of course use no interpolation for the E channel, since the RGB values “expect” exactly the original exponent. The RGB channels should be fine to interpolate, though, expect for seams between 4x4 blocks. (This could be really ugly.)

One could, however, take this into account. Let’s say you use linear interpolation. (Can we assume the interpolation will return exactly the same values for every card/driver? I think this is very well defined by D3D… ?) Then, sample down the original exponent texture first, and encode the RGB values based on the interpolated values. Theoretically, this should work perfectly in most cases. (You might need to run an extra pass on the downsampled exponent texture to give each pixel the highest value of its neighbours, so no interpolation gives less than needed max-RGB-value. Then the final color would allways stay too dark.)

Also, if you match the block-size to the DXT compression size (4x4), the RGB values should still compress nicely.

0
101 Sep 09, 2005 at 22:03

thanks to all repliers.

0
101 Sep 09, 2005 at 22:32

Actually, I have two textures: one for RGB and the other for exponent.
I have tried to apply bilinear filtering only on the RGB, disabling it in the exponent. Here the result:

The halo is more evident.

The problem is that mantissas related to two quite different exponent, cannot be interpolated linearly. :sad:

My research still continue…

Helps are always appreciated, thanks.

• AGPX
0
101 Sep 10, 2005 at 16:24

the only way working out correctly would be to do the linear sampling manually by doing 4 pointsamples…

on ps3.0 you could optimise it by determine if you actually are at an edge of the exponent, and do the full bilinear only as needed.

this way, you could gain much speed to get it near to pointsample performance.

0
101 Sep 10, 2005 at 20:29

Hi again,

I will become crazy… without doubts. :wacko:

I have changed my method. I split my 16 bit integer texture (A16R16G16B16) in two 8-bit per channel texture (one for the low byte, and the other for the high byte).
The low texture is compressed with DXT3 compression. The other one is stored uncompressed (A8R8G8B8). According to my calculations, this method should be immune from bilinear filtering issue. And it is… but only with Reference Rasterizer! With HAL device, the thing don’t works. :blink:
WHY??????

Take a small look to the pixel shader I used to reassemble the two 8 bit textures in a 16 bit one:

sampler2D lightMapLow;
sampler2D lightMapHigh;

float4 Expose(in float4 rgba, float exposure)
{
return 1 - exp(-rgba * exposure);
}

float4 texHDR2D(in sampler2D texLow, in sampler2D texHigh, in float2 tex)
{
return tex2D(texLow, tex) + tex2D(texHigh, tex) * 256.0;
}

float4 ps_main(float2 inTex: TEXCOORD0) : COLOR0
{
return Expose(texHDR2D(lightMapLow, lightMapHigh, inTex), 1.2);
}

And now take a look at the results using Reference Rasterizer:

And the following is with HAL device on Sapphire Radeon 9800 Pro (128 Mb, 256 bit):

HOW THE HELL IS THIS POSSIBLE?

P.S.: I have tried to use A8R8G8B8 for low texture instead of DXT3, but NONE change!

P.S.2: Works well on HAL only if I switch from bilinear filtering to point filtering!

P.S.3: Here you can see the low & high textures generated by my lightmapper:

LOW:

HIGH:

The high is not totally 0, they have some pixel set to 1,2 and 3. The textures appears to be ok.

0
101 Sep 10, 2005 at 23:27

Hi,

all the problems are now solved.

I post here my solution to support other peoples with a similar problem.
Basically, the problem is due to the high interpolators imprecision.
Reference rasterizer is preciser than HAL device.
So basically, there are nothing to do: you have to perform POINT filtering and make bilinear filtering via pixel shader. Here the code to perform the filtering (gently posted to me by a guy on #flipcode. Thanks to you!)

float4 bilinear(in sampler2D s, in float2 t)
{
const float2 wh = {256.0, 256.0};
const float2 dwh = 1 / wh;
float2 dxy = t * wh - floor(t * wh);
float4 a = lerp(tex2D(s, t + float2(0,   0)), tex2D(s, t + float2(dwh.x,   0)), dxy.x);
float4 b = lerp(tex2D(s, t + float2(0, dwh.y)), tex2D(s, t + float2(dwh.y, dwh.y)), dxy.x);
return lerp(a, b, dxy.y);
}

This way, also the RGBE encoding works well, so I switch to it, because it requires less storage memory. Oh, well, at least… only if I can compress it with DXT1-5 without introduction of artifacts… I will investigate on that tomorrow and, finally, I think that I’ll write a tutorial on HDR compression… could be helpful, until hardware manufacturer will introduce compression for 16 bit textures! (to say all the truth some 16-bit FourCC formats exists, but are largely not implemented).

Thanks to all for help. :D

0
101 Sep 11, 2005 at 21:25

thanks a lot for the posts, really nice to see that someone else is struggeling with bilinear interpolation of hdr-textures ;)

0
101 Sep 12, 2005 at 00:55

Reference rasterizer is preciser than HAL device.

Should read: Reference rasterizer is preciser than ATi HAL device. :rolleyes:

0
165 Sep 12, 2005 at 01:59

Yup. ATI skimps on precision, using only 24-bit floats for internal pipeline. But nVidia uses the full 32 bits.

0
101 Sep 13, 2005 at 15:27

Hello. Really nice topic :). And what with RGBE cubemaps?

0
101 Sep 13, 2005 at 20:47

Yup. ATI skimps on precision, using only 24-bit floats for internal pipeline. But nVidia uses the full 32 bits.

AFAIK the interpolators are FX8 only.

0
165 Sep 13, 2005 at 22:37

Yeah, the bilinear interpolation in the texture fetch may well be only 8-bit. But even in a pixel shader, the temp registers are only 24-bit rather than 32-bit.

0
101 Sep 13, 2005 at 22:49

No, I mean I’m pretty sure that the interpolators on Radeons are not floating point. The bilinear interpolation is only 6 bit btw…

0
101 Sep 13, 2005 at 22:57

yes but what with bilinear filtering on cubemaps? How to fetch the neightbours??

0
101 Sep 14, 2005 at 00:02

ok I’ve got it - It’s on ATI’s page.

0
101 Sep 14, 2005 at 09:58

@bonzaj

ok I’ve got it - It’s on ATI’s page.

0
165 Sep 14, 2005 at 16:09

Radeon 9500 and up supports floating point pbuffers, so I’m pretty sure they must have floating-point temp registers. Unless they are doing some kind of fancy bit-packing stuff…

0
101 Sep 15, 2005 at 15:55

@Reedbeta

Radeon 9500 and up supports floating point pbuffers, so I’m pretty sure they must have floating-point temp registers. Unless they are doing some kind of fancy bit-packing stuff…

They have FP24 temp registers, but the interpolators are not FP.

0
101 Sep 15, 2005 at 18:45

kusma: here You go :)


text:
A recent line of research has investigated a new profiling technique based on deep learning as an alternative to the well-known template attack.
The advantage of this new profiling approach is twofold: $(1)$ the approximation of the information leakage by a multivariate Gaussian distribution is relaxed (leading to a more generic approach) and $(2)$ the pre-processing phases such as the traces realignment or the selection of the Points of Interest (PoI) are no longer mandatory, in some cases, to succeed the key recovery (leading to a less complex security evaluation roadmap).
The related published works have demonstrated that Deep Learning based Side-Channel Attacks (DL-SCA) are very efficient when targeting cryptographic implementations protected with the common side-channel countermeasures such as masking, jitter and random delays insertion.
In this paper, we assess the efficiency of this new profiling attack under different realistic and practical scenarios. First, we study the impact of the intrinsic characteristics of the manipulated data-set (\emph{i.e.} distance in time samples between the PoI, the dimensionality of the area of interest and the pre-processing of the data) on the robustness of the attack. We demonstrate that the deep learning techniques are sensitive to these parameters and we suggest some practical recommendations that can be followed to enhance the profiling and the key recovery phases. Second, we discuss the tolerance of DL-SCA with respect to a deviation from the idealized leakage models and provide a comparison with the well-known stochastic attack. Our results show that DL-SCA are still efficient in such a context.
Then, we target a more complex masking scheme based on Shamir's secret sharing and prove that this new profiling approach is still performing well. Finally, we conduct a security evaluation of a batch of several combinations of side-channel protections using simulations and real traces captured on the ChipWhisperer board. The experimental results obtained confirm that DL-SCA are very efficient even when a cryptographic implementation combines several side-channel countermeasures.

All News


text:
# Q.33. The variation of potential energy with intermolecular separation is given by For r > ${r}_{1}$, force is (1) Attractive (2) Repulsive (3) Zero (4) May be attractive, repulsive or zero

Dear Student,

so option (4) is correct.
If you have any more doubts just ask here on the forum and our experts will try to help you out as soon as possible.

Regards

• 0
I don't know
• 0
What are you looking for?