#include <iostream>
#include <fstream>
#include <cstdlib>
#include <cstdio>
#include <iomanip>
#include <cmath>

#include "Daubechies.h"
#include "Coefficients.h"

char          Daubechies::TheP;
char          Daubechies::diese;
int           Daubechies::type;
int           Daubechies::dimX;
int           Daubechies::maxIntensity;
std::string   Daubechies::line_of_comments;

double       *Daubechies::v;
double       *Daubechies::workAAA;
double       *Daubechies::workBBB; 
double       *Daubechies::workCCC;
double       *Daubechies::h;
double       *Daubechies::g;
int           Daubechies::max_iter;

void Daubechies::Do_Discrete_Matrix_Transformation(char File_in[], 
           char File_out[], int Total_Iterations, int Dorder)
{

   ReadSignal(File_in);
   AssignCoefficients(Dorder);
   VerifyMaxIterations(Total_Iterations);

   DWT1Diter(Dorder);

   int type;

   std::cout << "Data Compression: press 1\n";
   std::cout << "Edge Detection: press 0\n";
   std::cin  >> type;

   AlterTransform(type);

   WriteSignal(File_out);
   delete[] Daubechies::v;

}

void Daubechies::Do_Inverse_Matrix_Transformation(char File_in[], 
           char File_out[], int Total_Iterations, int Dorder)
{

   ReadSignal(File_in);
   AssignCoefficients(Dorder);
   VerifyMaxIterations(Total_Iterations);

   IDWT1Diter(Dorder);

   WriteSignal(File_out);
   delete[] Daubechies::v;

}

// *****************************************************************************
// ******************************** DWT ****************************************
// *****************************************************************************
  

void Daubechies::DWT1Diter(int Dorder)
{

   int dim_matrix;

   for (int Jtrans=0; Jtrans < Daubechies::max_iter; Jtrans++) {

      dim_matrix = int(Daubechies::dimX/pow(2,Jtrans));

      std::cout << "dim_matrix=" << dim_matrix << "\n";

      AllocateWork(dim_matrix);

      GetCorner(Daubechies::workAAA,                      dim_matrix        );
      DWT1D1   (Daubechies::workAAA, Daubechies::workBBB, dim_matrix, Dorder);
      PutCorner(Daubechies::workBBB,                      dim_matrix        );

      DeallocateWork(dim_matrix);
   }

}

void Daubechies::DWT1D1(double *A, double *B, int N, int Dorder)
{
// Reference: Van Fleet, Discrete Wavelet Transformations, page 270-271

   int p, L;

   L = 2*Dorder - 1;
   p = (L - 1)/2; // number of wrapping rows.

// Compute lowpass and highpass values for the non-wrapping rows.
// To index highpass values, add N/2 to the subscript

   for (int k=1; k <= N/2 -p; k++) {
      B[      k] = 0.0;
      B[N/2 + k] = 0.0;
      for (int j=0; j <= L; j++) {
         B[      k] += Daubechies::h[L - j]*A[2*k - 1 + j];
         B[N/2 + k] += Daubechies::g[L - j]*A[2*k - 1 + j];
      }
   }

// compute lowpass and highpass values for the wrapping rows.
   for (int k=1; k <= p; k++) {
      B[N/2 - p + k] = 0.0;
      B[N   - p + k] = 0.0;
      // This loop is the end sum (7.94)
      for (int j=0; j <= L-2*k; j++) {
         B[N/2 - p + k] += Daubechies::h[L - j]*A[N - L + 2*k + j];
         B[N   - p + k] += Daubechies::g[L - j]*A[N - L + 2*k + j];
      }
      // This loop adds the beginning sum (7.93)
      for (int j=1; j <= 2*k; j++) {
         B[N/2 - p + k] += Daubechies::h[2*k - j]*A[j];
         B[N   - p + k] += Daubechies::g[2*k - j]*A[j];
      }
   }

}

// *****************************************************************************
// ******************************* IDWT ****************************************
// *****************************************************************************
  
void Daubechies::IDWT1Diter(int Dorder)
{
// Reference: Van Fleet, "Discrete Wavelet Transformations", 2007,
// page 169 *********************************
// Careful, there are several errors in the algorithm given in the book by
// Van Fleet
// Follow the argument in the book and some index become clear

   int dim_matrix;

   for (int Jtrans=1; Jtrans <= Daubechies::max_iter; Jtrans++) {

      dim_matrix = int(Daubechies::dimX/pow(2,Daubechies::max_iter - Jtrans));

      AllocateWork(dim_matrix);

      GetCorner(Daubechies::workAAA,                      dim_matrix        );
      IDWT1D1  (Daubechies::workAAA, Daubechies::workBBB, dim_matrix, Dorder);
      PutCorner(Daubechies::workBBB,                      dim_matrix        );

      DeallocateWork(dim_matrix);
   }
}

void Daubechies::IDWT1D1(double *A, double *B, int N, int Dorder)
{
// algorithm page 277-278 Van Fleet

   double  lowPass[N/2 + 1];
   double highPass[N/2 + 1];
   double  InverseTransform[N + 1];

// Separate the matrix vector into low and high pass
   for (int i = 1; i <= N/2; i++) {
       lowPass[i] = A[i      ];
      highPass[i] = A[i + N/2];
   }

// nullify the vector, just in case.
   for (int i = 1; i <= N ; i++) {
      InverseTransform[i] = 0.0;
   }

// Perform the inverse transform on the low pass section
   IDWTht( lowPass, InverseTransform, Daubechies::h, N, Dorder);

// copy result in work space B
   for (int i = 1; i <= N ; i++) {
      B[i] = InverseTransform[i];
   }

// Perform the inverse transform on the high pass section
   IDWTht(highPass, InverseTransform, Daubechies::g, N, Dorder);

// add value to work space B
   for (int i = 1; i <= N ; i++) {
      B[i] += InverseTransform[i];
   }

// B is return

}

void Daubechies::IDWTht(double *signal, double *InverseTransform,
                        double *filter, int N, int Dorder)
{

   int L, p, r, n;

   L = 2*Dorder - 1;
   p = (L - 1)/2; // number of wrapping rows.
   r = p + 1;

   double even[r], odd[r];

   n = N/2;

// Initialize the vectors e and o.
   for (int k=0; k <= p; k++) {
      even[k+1] = filter[2*k  ];
       odd[k+1] = filter[2*k+1];
   }

// Compute the nonwrapping row values
   for (int k=1; k <= n-p; k++) {
      InverseTransform[L + 2*k - 2] = 0.0;
      InverseTransform[L + 2*k - 1] = 0.0;
//------------------------------------------------------------------
//    incorrect index in algorithm 7.2, it should be p+k           |
//    and equation (7.99) is wrong, it is (L+1)/2+k-1 or (L-1)/2+k |
//------------------------------------------------------------------
      for (int j=0; j<=p; j++) {
//       odd-even index start at 1, not 0
         InverseTransform[L + 2*k - 2] +=  odd[j + 1]*signal[j + k];
         InverseTransform[L + 2*k - 1] += even[j + 1]*signal[j + k];
      }
   }

// Compute the wrapping row values
   for (int k=1; k <= p; k++) {
      InverseTransform[2*k - 1] = 0.0;
      InverseTransform[2*k  ] = 0.0;

      // This loop is the beginning sum (7.101), page 275
      for (int j=1; j<= k ; j++) {
         InverseTransform[2*k - 1] +=  odd[r - k + j]*signal[j];
         InverseTransform[2*k    ] += even[r - k + j]*signal[j];
      }

      // This loop adds the end sm (7.102), page 275
//------------------------------------------------------------------
// End sum is incorrect in the textbook, it should be O[j]         |
//------------------------------------------------------------------
      for (int j=1; j<= r-k ; j++) {
         InverseTransform[2*k - 1] +=  odd[j]*signal[n - r + k + j];
         InverseTransform[2*k    ] += even[j]*signal[n - r + k + j];
      }
   }

}

// *****************************************************************************
// *************************** UTIL ********************************************
// *****************************************************************************
  
void Daubechies::ReadSignal(char File_in[])
{

// Notice that this version accepts only square matrices
   int dimY, dummy;
   std::ifstream file_in_A;

   file_in_A.open(File_in);
   if (file_in_A.fail() ) {
      std::cerr << "Unable to open input files\n";
      exit(8);
   }

   file_in_A >> Daubechies::TheP >> Daubechies::type;
   if (Daubechies::type == 888) {
      std::cout << "Real signal\n";
      file_in_A >> Daubechies::dimX;
   }
   Daubechies::v = new double[dimX+1];

   for (int i=1; i<= Daubechies::dimX; i++) {
      file_in_A  >> Daubechies::v[i];
   }

   file_in_A.close();
}

void Daubechies::WriteSignal(char File_out[])
{
   std::ofstream file_out_A;

   file_out_A.open(File_out);
   if (file_out_A.fail() ) {
      std::cerr << "Unable to open output files\n";
      exit(8);
   }

   file_out_A  << Daubechies::TheP << 888 << "\n";
   file_out_A << Daubechies::dimX << "\n";

   for (int i=1; i<= Daubechies::dimX; i++) {
      file_out_A  << std::setw(25) << std::fixed << std::setprecision(15) 
                  << Daubechies::v[i] << "\n";
  }

  file_out_A.close();
} 

void Daubechies::AssignCoefficients(int Dorder)
{

   Daubechies::h = new double[2*Dorder];
   Daubechies::g = new double[2*Dorder];

   if (2*Dorder > Daubechies::dimX) {
      std::cout << "\n\nReduce order to a value less than or equal to "
                << Daubechies::dimX/2 << "; Program is STOPPING\n";
      exit(8);
   }

   for (int i=0; i< 2*Dorder; i++) {
      switch (Dorder) {
         case  1: Daubechies::h[i] = Daub01[i]; break;
         case  2: Daubechies::h[i] = Daub02[i]; break;
         case  3: Daubechies::h[i] = Daub03[i]; break;
         case  4: Daubechies::h[i] = Daub04[i]; break;
         case  5: Daubechies::h[i] = Daub05[i]; break;
         case  6: Daubechies::h[i] = Daub06[i]; break;
         case  7: Daubechies::h[i] = Daub07[i]; break;
         case  8: Daubechies::h[i] = Daub08[i]; break;
         case  9: Daubechies::h[i] = Daub09[i]; break;
         case 10: Daubechies::h[i] = Daub10[i]; break;
         case 11: Daubechies::h[i] = Daub11[i]; break;
         case 12: Daubechies::h[i] = Daub12[i]; break;
         case 13: Daubechies::h[i] = Daub13[i]; break;
         case 14: Daubechies::h[i] = Daub14[i]; break;
         case 15: Daubechies::h[i] = Daub15[i]; break;
         case 16: Daubechies::h[i] = Daub16[i]; break;
         case 17: Daubechies::h[i] = Daub17[i]; break;
         case 18: Daubechies::h[i] = Daub18[i]; break;
         case 19: Daubechies::h[i] = Daub19[i]; break;
         case 20: Daubechies::h[i] = Daub20[i]; break;
         case 21: Daubechies::h[i] = Daub21[i]; break;
         case 22: Daubechies::h[i] = Daub22[i]; break;
         case 23: Daubechies::h[i] = Daub23[i]; break;
         case 24: Daubechies::h[i] = Daub24[i]; break;
         case 25: Daubechies::h[i] = Daub25[i]; break;
         case 26: Daubechies::h[i] = Daub26[i]; break;
         case 27: Daubechies::h[i] = Daub27[i]; break;
         case 28: Daubechies::h[i] = Daub28[i]; break;
         case 29: Daubechies::h[i] = Daub29[i]; break;
         case 30: Daubechies::h[i] = Daub30[i]; break;
         case 31: Daubechies::h[i] = Daub31[i]; break;
         case 32: Daubechies::h[i] = Daub32[i]; break;
         case 33: Daubechies::h[i] = Daub33[i]; break;
         case 34: Daubechies::h[i] = Daub34[i]; break;
         case 35: Daubechies::h[i] = Daub35[i]; break;
         case 36: Daubechies::h[i] = Daub36[i]; break;
         case 37: Daubechies::h[i] = Daub37[i]; break;
         case 38: Daubechies::h[i] = Daub38[i]; break;

         default: std::cout << "order must be between 1 and 38, inclusive\n";
                  exit(8); break;

         }
   }
   for (int k = 0; k <  2*Dorder; k++) {
      g[k] = pow(-1, k)*Daubechies::h[(2*Dorder-1) - k];
//    std::cout << "h[" << k << "]=" << Daubechies::h[k] <<"\n";
//    std::cout << "g[" << k << "]=" << Daubechies::g[k] <<"\n";
   }

}

void Daubechies::VerifyMaxIterations(int Total_Iterations)
{
   int ratio_of_matrix;
   int test_modulo;
   int i;

   Daubechies::max_iter = Total_Iterations;

   i = 1;
   while (i <= Total_Iterations) {
      ratio_of_matrix = int(pow(2,i));
      test_modulo = Daubechies::dimX % ratio_of_matrix;
      if (test_modulo != 0) {
         Daubechies::max_iter = i-1;
         break;
      }
      i++;
   }
   std::cout << "# iterations =" << Daubechies::max_iter << "\n";
}

void Daubechies::AllocateWork(int dim_matrix)
{
   Daubechies::workAAA = new double[dim_matrix+1];
   Daubechies::workBBB = new double[dim_matrix+1];
}

void Daubechies::DeallocateWork(int dim_matrix)
{
   delete[] Daubechies::workAAA;
   delete[] Daubechies::workBBB;
}

void Daubechies::GetCorner(double *workAAA, int dim_matrix)
{

   for (int i=1; i<= dim_matrix; i++) {
      workAAA[i] = v[i];
   }

}

void Daubechies::PutCorner(double *workBBB, int dim_matrix)
{

   for (int i=1; i<= dim_matrix; i++) {
      v[i] = workBBB[i];
   }
}

void Daubechies::AlterTransform(int type)
{

   int dim_blur;
   int zeroed, total_edge;

   double higher_bound, lower_bound;

   dim_blur = int(Daubechies::dimX/pow(2,Daubechies::max_iter));

   if (type == 1) {
      std::cout << "WARNING: Compression: zeroing part of complement{blur}"
                << " of dim="
                << dimX - dim_blur << "\n";
      zeroed = 0;
      total_edge = 0;

      double min =  1.0e+10;
      double max = -1.0e+10;
      for (int i=1; i<= Daubechies::dimX; i++) {
         if (!(i <= dim_blur)) {
            if (Daubechies::v[i] < min) min = v[i];
            if (Daubechies::v[i] > max) max = v[i];
         }
      }
      std::cout << "Edge section: min=" << min << " max=" << max << "\n";
      std::cout << "Data Compression:  type in the higher_bound:";
      std::cin >> higher_bound;
      std::cout << "Data Compression:  type in the lower_bound:";
      std::cin >> lower_bound;

      for (int i=1; i<= dimX; i++) {
         if (!(i <= dim_blur)) {
            total_edge += 1;
            if (lower_bound <= Daubechies::v[i] &&
                               Daubechies::v[i] <= higher_bound) {
               Daubechies::v[i] = 0.0;
               zeroed += 1;
            }
         }
      }
      std::cout << "Number of elements set to zero:" << zeroed
                << " (" << double(zeroed) /double(total_edge)*100. << "\%)\n";

   }
   if (type == 0) {
      std::cout << "WARNING: Edge-detection: Zeroing the blur sub-matrix of dim="
                << dim_blur << "\n";
      for (int i=1; i<= dimX; i++) {
         if ( (i <= dim_blur)) {
            Daubechies::v[i] = 0.0;
         }
      }
   }

}