#include <cstdlib>
#include <ctime>
#include <cmath>
#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
#include <assert.h>

using namespace std;
//
// Input for the program DMA part
//
int N = 8;			// The size of the cube in power of 2
double Hurst_exponent = 0.1;	// The input of the exponent for the sponge
double theta = 0.5;             // The values of theta in order to calculate the average
int n_max = 23;                 // Max values in the sum n_max << N_max
int n_min = 3;                  // Starting value for the sum >= 3
int m = 0;                      // Default values for the m for the average case to use

int N_0_p_glob = int(pow(2.,N) + 1);
int N_1_p_glob = int(pow(2.,N) + 1);
int N_2_p_glob = int(pow(2.,N) + 1);					

int N_tot_p_glob = int(N_0_p_glob*N_1_p_glob*N_2_p_glob);
//
// Functions for the program
//
int main(int argc, char *argv[]); // The main program function

int main(int argc, char *argv[]){
//
// Initialize a random seed
//
  srand(time(NULL));
//
// Program variables
//
  // RMD data file - Start
  ifstream rmd_data_file;
  stringstream N_number_rmd, H_number_rmd;
  
  string rmd_data = "rmd_data_";
  string buffer;

  N_number_rmd << N_0_p_glob;
  H_number_rmd << Hurst_exponent;
  rmd_data = rmd_data + N_number_rmd.str() + "_";
  rmd_data = rmd_data + H_number_rmd.str() + ".dat";

  rmd_data_file.open(rmd_data.c_str());

  if (rmd_data_file.fail()){
	cout << "THE DATA FILE " + rmd_data + " DOES NOT EXIST" << endl;
	exit(1);
  }
  // RMD data file - End

  // DMA data file - Start
  ofstream dma_data_file;
  stringstream N_number_dma, H_number_dma;
  
  string dma_data = "dma_data_";

  N_number_dma << N_0_p_glob;
  H_number_dma << Hurst_exponent;
  dma_data = dma_data + N_number_dma.str() + "_";
  dma_data = dma_data + H_number_dma.str() + ".dat";

  dma_data_file.open(dma_data.c_str());
    
  // DMA data file - End

//
// The 3d DMA algorithm
//
  cout << "---- DMA ALGORITHM - ESTIMATION OF THE HURST EXPONENT ----" << endl;
  cout << "                                                          " << endl;
  cout << "                 THREE DIMENSIONAL CASE                   " << endl;
  cout << "----                                                  ----" << endl;
  cout << "                                                          " << endl;
  
  cout << "Size of cube: (" << N_0_p_glob << "," << N_1_p_glob << "," << N_2_p_glob << ")" << endl;  
  cout << "Theta: (" << theta << "," << theta << "," << theta << ")" << endl;

  dma_data_file << "# Hurst Exponent: " << Hurst_exponent << endl;
  dma_data_file << "# Size of cube: (" << N_0_p_glob << "," << N_1_p_glob << "," << N_2_p_glob << ")" << endl;
  dma_data_file << "# Theta: (" << theta << "," << theta << "," << theta << ")" << endl;

  // Reading the data - Start
  double *point = new double[N_tot_p_glob]; // The global structure of the grid
  int count = 0;

  cout << "READING THE DATA FILE: " + rmd_data << endl;

  while (getline(rmd_data_file, buffer)){
    if (buffer == "" || buffer[0] == '#') continue;
	
    stringstream bufferstr;
	
    double databuff;

	bufferstr << buffer;
	
    bufferstr >> databuff;
	
	*(point+count) = databuff;
	count++;
  }
  
  rmd_data_file.close();  
  // Reading the data - End
  
  // --- START: Algorithm
  int n_1_min = n_min, n_2_min = n_min, n_3_min = n_min; 
  int n_1_max = n_max, n_2_max = n_max, n_3_max = n_max;
  int N_1_glob = N_0_p_glob, N_2_glob = N_1_p_glob, N_3_glob = N_2_p_glob;

  int m_1 = m, m_2 = m, m_3 = m;
  double theta_1 = theta, theta_2 = theta, theta_3 = theta;

  for (int n_1 = n_1_min; n_1 < (n_1_max + 1); n_1++){
    for (int n_2 = n_2_min; n_2 < (n_2_max + 1); n_2++){
      for (int n_3 = n_3_min; n_3 < (n_3_max + 1); n_3++){
		if ((n_1 == n_2) && (n_2 == n_3)){  // The line of choice for the calculations
		  if ( (n_1 % 2) == 1 ){
		  
			bool average_sum_start = true;
			
            double sigma2_DMA_tot = 0.; 
			
            double sigma2_DMA = 0.; 
            double sigma2_sum_var = 0.;

            double average_sum = 0., average_sum_var = 0.;

            double sigma2_DMA_1 = 0., sigma2_DMA_2 = 0., sigma2_DMA_3 = 0.;

            double sigma2_sum_var_1 = 0., sigma2_sum_var_2 = 0., sigma2_sum_var_3 = 0.;

            double average_sum_1 = 0., average_sum_2 = 0., average_sum_3 = 0. ;
            double average_sum_var_1 = 0., average_sum_var_2 = 0., average_sum_var_3 = 0.;

            sigma2_DMA = 1./(N_1_glob - n_1_max);  // First factor
            sigma2_DMA *= 1./(N_2_glob - n_2_max); // Second factor
            sigma2_DMA *= 1./(N_3_glob - n_3_max); // Third factor
			
			sigma2_DMA_1 = sigma2_DMA_2 = sigma2_DMA_3 = sigma2_DMA;

            m_1 = int(n_1*theta_1); // Setting m_1
            m_2 = int(n_2*theta_2);	// Setting m_2
            m_3 = int(n_3*theta_3);	// Setting m_3

            for (int i_1 = (n_1 - m_1); i_1 <= (N_1_glob - m_1); i_1++){
			  
              int i_2, i_3;
			  
			  i_3 = i_2 = i_1;

			  int f_nn = int((i_3 - 1)*N_2_p_glob*N_1_p_glob + (i_2 - 1)*N_1_p_glob + (i_1 - 1));

			  average_sum = 1./(n_1*n_2*n_3);
			  average_sum_var = 0.;
			  
			  for (int k_1 = (-m_1); k_1 <= (n_1 - 1 - m_1); k_1++){
				for (int k_2 = (-m_2); k_2 <= (n_2 - 1 - m_2); k_2++){
				  for (int k_3 = (-m_3); k_3 <= (n_3 - 1 - m_3); k_3++){
					
					int f_average_nn = int((i_3 - k_3 - 1)*N_2_p_glob*N_1_p_glob + (i_2 - k_2 - 1)*N_1_p_glob + (i_1 - k_1 - 1));
					average_sum_var += *(point+f_average_nn);
					
				  }
				}
			  }
			  
			  average_sum_var_1 = average_sum_var_2 = average_sum_var_3 = average_sum_var;
			  
			  average_sum *= average_sum_var;
			  sigma2_sum_var += pow(((*(point+f_nn)) - average_sum),2.);

			  // Direction positive - Start
			  for (int i_i = (i_1 + 1); i_i <= (N_1_glob - m_1); i_i++){
				
				double average_sum_var_decrease_1 = 0., average_sum_var_increase_1 = 0.;
				int f_average_nn_1;
				
				double average_sum_var_decrease_2 = 0., average_sum_var_increase_2 = 0.;
				int f_average_nn_2;
				
				double average_sum_var_decrease_3 = 0., average_sum_var_increase_3 = 0.;
				int f_average_nn_3;
				
				int f_nn_1 = int((i_3 - 1)*N_2_p_glob*N_1_p_glob + (i_2 - 1)*N_1_p_glob + (i_i - 1));
				int f_nn_2 = int((i_3 - 1)*N_2_p_glob*N_1_p_glob + (i_i - 1)*N_1_p_glob + (i_1 - 1));
				int f_nn_3 = int((i_i - 1)*N_2_p_glob*N_1_p_glob + (i_2 - 1)*N_1_p_glob + (i_1 - 1));
				
				for (int k_face_1 = (-m_1); k_face_1 <= (n_1 - 1 - m_1); k_face_1++){
				  for (int k_face_2 = (-m_2); k_face_2 <= (n_2 - 1 - m_2); k_face_2++){
					
					int k_1, k_2, k_3;
					
					k_2 = k_face_1; k_3 = k_face_2;
					
					k_1 = -m_1; 
					f_average_nn_1 = int((i_3 - k_3 - 1)*N_2_p_glob*N_1_p_glob + (i_2 - k_2 - 1)*N_1_p_glob + (i_i - k_1 - 1));
					average_sum_var_increase_1 += *(point+f_average_nn_1);
					
					k_1 = n_1 - 1 - m_1;
					f_average_nn_1 = int((i_3 - k_3 - 1)*N_2_p_glob*N_1_p_glob + (i_2 - k_2 - 1)*N_1_p_glob + ((i_i - 1) - k_1 - 1));
					average_sum_var_decrease_1 += *(point+f_average_nn_1);
					
					k_1 = k_face_1; k_3 = k_face_2;
					
					k_2 = -m_2; 
					f_average_nn_2 = int((i_3 - k_3 - 1)*N_2_p_glob*N_1_p_glob + (i_i - k_2 - 1)*N_1_p_glob + (i_1 - k_1 - 1));
					average_sum_var_increase_2 += *(point+f_average_nn_2);
					
					k_2 = n_2 - 1 - m_2;
					f_average_nn_2 = int((i_3 - k_3 - 1)*N_2_p_glob*N_1_p_glob + ((i_i - 1) - k_2 - 1)*N_1_p_glob + (i_1 - k_1 - 1));
					average_sum_var_decrease_2 += *(point+f_average_nn_2);
					
					k_1 = k_face_1; k_2 = k_face_2;
					
					k_3 = -m_3; 
					f_average_nn_3 = int((i_i - k_3 - 1)*N_2_p_glob*N_1_p_glob + (i_2 - k_2 - 1)*N_1_p_glob + (i_1 - k_1 - 1));
					average_sum_var_increase_3 += *(point+f_average_nn_3);
					
					k_3 = n_3 - 1 - m_3;
					f_average_nn_3 = int(((i_i - 1) - k_3 - 1)*N_2_p_glob*N_1_p_glob + (i_2 - k_2 - 1)*N_1_p_glob + (i_1 - k_1 - 1));
					average_sum_var_decrease_3 += *(point+f_average_nn_3);
					
				  }
				}
				
				average_sum_var_1 += average_sum_var_increase_1;
				average_sum_var_1 -= average_sum_var_decrease_1;
				
				average_sum_var_2 += average_sum_var_increase_2;
				average_sum_var_2 -= average_sum_var_decrease_2;
				
				average_sum_var_3 += average_sum_var_increase_3;
				average_sum_var_3 -= average_sum_var_decrease_3;
				
				average_sum_1 = 1./(n_1*n_2*n_3);
				average_sum_1 *= average_sum_var_1;
				sigma2_sum_var_1 += pow(((*(point+f_nn_1)) - average_sum_1),2.);
				
				average_sum_2 = 1./(n_1*n_2*n_3);
				average_sum_2 *= average_sum_var_2;
				sigma2_sum_var_2 += pow(((*(point+f_nn_2)) - average_sum_2),2.);
				
				average_sum_3 = 1./(n_1*n_2*n_3);
				average_sum_3 *= average_sum_var_3;
				sigma2_sum_var_3 += pow(((*(point+f_nn_3)) - average_sum_3),2.);
				
			  }
			  // Direction positive - End
			  
			  average_sum_var_1 = average_sum_var_2 = average_sum_var_3 = average_sum_var;

			  // Direction negative - Start
			  for (int i_i_1 = (i_1 - 1); i_i_1 >= (n_1 - m_1); i_i_1--){
				
				double average_sum_var_decrease_1 = 0., average_sum_var_increase_1 = 0.;
				int f_average_nn_1;
				
				double average_sum_var_decrease_2 = 0., average_sum_var_increase_2 = 0.;
				int f_average_nn_2;
				
				double average_sum_var_decrease_3 = 0., average_sum_var_increase_3 = 0.;
				int f_average_nn_3;
				
				int f_nn_1 = int((i_3 - 1)*N_2_p_glob*N_1_p_glob + (i_2 - 1)*N_1_p_glob + (i_i_1 - 1));
				int f_nn_2 = int((i_3 - 1)*N_2_p_glob*N_1_p_glob + (i_i_1 - 1)*N_1_p_glob + (i_1 - 1));
				int f_nn_3 = int((i_i_1 - 1)*N_2_p_glob*N_1_p_glob + (i_2 - 1)*N_1_p_glob + (i_1 - 1));
				
				for (int k_face_1 = (-m_1); k_face_1 <= (n_1 - 1 - m_1); k_face_1++){
				  for (int k_face_2 = (-m_2); k_face_2 <= (n_2 - 1 - m_2); k_face_2++){
					
					int k_1, k_2, k_3;
					
					k_2 = k_face_1; k_3 = k_face_2;
					
					k_1 = n_1 - 1 - m_1;
					f_average_nn_1 = int((i_3 - k_3 - 1)*N_2_p_glob*N_1_p_glob + (i_2 - k_2 - 1)*N_1_p_glob + (i_i_1 - k_1 - 1));
					average_sum_var_increase_1 += *(point+f_average_nn_1);
					
					k_1 = -m_1; 
					f_average_nn_1 = int((i_3 - k_3 - 1)*N_2_p_glob*N_1_p_glob + (i_2 - k_2 - 1)*N_1_p_glob + ((i_i_1 + 1) - k_1 - 1));
					average_sum_var_decrease_1 += *(point+f_average_nn_1);
					
					k_1 = k_face_1; k_3 = k_face_2;
					
					k_2 = n_2 - 1 - m_2;
					f_average_nn_2 = int((i_3 - k_3 - 1)*N_2_p_glob*N_1_p_glob + (i_i_1 - k_2 - 1)*N_1_p_glob + (i_1 - k_1 - 1));
					average_sum_var_increase_2 += *(point+f_average_nn_2);
					
					k_2 = -m_2; 
					f_average_nn_2 = int((i_3 - k_3 - 1)*N_2_p_glob*N_1_p_glob + ((i_i_1 + 1) - k_2 - 1)*N_1_p_glob + (i_1 - k_1 - 1));
					average_sum_var_decrease_2 += *(point+f_average_nn_2);
					
					k_1 = k_face_1; k_2 = k_face_2;
					
					k_3 = n_3 - 1 - m_3;
					f_average_nn_3 = int((i_i_1 - k_3 - 1)*N_2_p_glob*N_1_p_glob + (i_2 - k_2 - 1)*N_1_p_glob + (i_1 - k_1 - 1));
					average_sum_var_increase_3 += *(point+f_average_nn_3);
					
					k_3 = -m_3; 
					f_average_nn_3 = int(((i_i_1 + 1) - k_3 - 1)*N_2_p_glob*N_1_p_glob + (i_2 - k_2 - 1)*N_1_p_glob + (i_1 - k_1 - 1));
					average_sum_var_decrease_3 += *(point+f_average_nn_3);
					
				  }
				}
				
				average_sum_var_1 += average_sum_var_increase_1;
				average_sum_var_1 -= average_sum_var_decrease_1;
				
				average_sum_var_2 += average_sum_var_increase_2;
				average_sum_var_2 -= average_sum_var_decrease_2;
				
				average_sum_var_3 += average_sum_var_increase_3;
				average_sum_var_3 -= average_sum_var_decrease_3;
				
				average_sum_1 = 1./(n_1*n_2*n_3);
				average_sum_1 *= average_sum_var_1;
				sigma2_sum_var_1 += pow(((*(point+f_nn_1)) - average_sum_1),2.);
				
				average_sum_2 = 1./(n_1*n_2*n_3);
				average_sum_2 *= average_sum_var_2;
				sigma2_sum_var_2 += pow(((*(point+f_nn_2)) - average_sum_2),2.);
				
				average_sum_3 = 1./(n_1*n_2*n_3);
				average_sum_3 *= average_sum_var_3;
				sigma2_sum_var_3 += pow(((*(point+f_nn_3)) - average_sum_3),2.);
				
			  }
			  // Direction negative - End
			  
			}

			// Sum up the results - Start
			sigma2_DMA *= sigma2_sum_var;
			
			sigma2_DMA_1 *= sigma2_sum_var_1;
			
			sigma2_DMA_2 *= sigma2_sum_var_2;

			sigma2_DMA_3 *= sigma2_sum_var_3;
			
			sigma2_DMA_tot = sigma2_DMA + sigma2_DMA_1 + sigma2_DMA_2 + sigma2_DMA_3;
			// Sum up the results - End

			dma_data_file << n_1 << "\t" << sigma2_DMA_tot << endl; // Output of the results to a data file
						
          }
        }
      }
    }
  }
  // --- END: Algorithm

  dma_data_file.close();
  
  return 0;
}