#include <vector> 
#include <fstream> 
#include <cmath>
#include <iostream> 
#include <cstdlib>
#include <sstream>
#include <cfloat>

using namespace std; 


int main (int argc, const char* argv[]) 
{
  if (system("clear")) {}

  cout << "======= The Lagrangian Zel'dovich Void Finder =======" << endl;
  
  string dir_input = "/home/andrii/voronoi/voids/code/federico_upd/input/"; //test
  string dir_output = "/home/andrii/voronoi/voids/code/federico_upd/output/"; //test
  

  // ----- get the parameters of the method from the parameter file ----- 

  string par_file = (argc>1) ? argv[1] : dir_input+"input_param.ini";
  
  string file_in = dir_input+"coord_LCDM_0.5.txt";
  string type = "a";
  string output = "coord_LCDM_0.5";
  int n_max = 1000;
  string mode = "p";
  int n_rnd = 3;

  ifstream finput (par_file.c_str());
  if (!finput) { cerr <<"Error in opening the parameter file "<<par_file<<"!\n\a"; exit(1); }

  string line, par_name, par_value, eq;
  
  while (getline(finput, line)) {
   
    stringstream line_string (line);
    line_string >> par_name;
      
    if (par_name!="#" && !par_name.empty()) {
    
      line_string >> eq;
      line_string >> par_value;
      
      if (eq != "=" || par_value.empty()) {
    cerr << "Error in the parameters file at " << par_name << endl;
    cerr << "Using the default value for " << par_name << endl;
      }

      if (par_name == "input") file_in = dir_input+par_value;
      if (par_name == "type") type = par_value;
      if (par_name == "output") output = par_value;
      if (par_name == "n_max") n_max = atoi(par_value.c_str());
      if (par_name == "mode") mode = par_value;
      if (par_name == "n_rnd") n_rnd = atoi(par_value.c_str());

    }   
  }
  
  finput.clear(); finput.close();

  cout  << endl << "file_in = " << file_in << endl;
  cout  << endl << "type = " << type << endl;
  cout  << "output = " << output << endl;
  cout  << "n_max = " << n_max << endl;
  cout  << "mode = " << mode << endl;
  cout  << "n_rnd = " << n_rnd << endl << endl;
 

  // --------------------------------------------------------------------

  
  srand(time(NULL)); // seed for randomization


  vector<double> xc, yc, zc;
  double shift_x,shift_z,shift_y,min_k,dist_k,dist_min,dist,n1,n2,n3,sum1,sum2,dist_new4,dist_new2,min_x,min_y,max_x,max_y,min_z,max_z;
  vector<int> num(20);

  if (type=="a") {

      ifstream fin (file_in.c_str());   
      if (!fin) { cerr <<"Error in opening ASCII file "<<file_in<<"!\n\a"; exit(1); } 

      while (fin >>n1>>n2>>n3) {
      xc.push_back(n1);
      yc.push_back(n2);
      zc.push_back(n3);
    } 
  fin.clear(); fin.close();
  }

  if (type=="b") { 
      short num;
      float val;
      ifstream fin (file_in.c_str(), ios::in|ios::binary|ios::ate);
       if (!fin) { cerr <<"Error in opening binary file "<<file_in<<"!\n\a"; exit(2); } 
       if (fin.is_open())
           fin.seekg (0, ios::beg);
           fin.read ((char*)(&num), 2);
           int  n_blocks=num;
           for (int i=1; i<=n_blocks; ++i) {
               fin.read ((char*)(&num), 4);
               int n_objs=num;
               for (int j=1; j<=n_objs; ++j){ 
                  fin.read ((char*)(&val), 4);
                  xc.push_back(val);
                  fin.read ((char*)(&val), 4);
                  yc.push_back(val);
                  fin.read ((char*)(&val), 4);
                  zc.push_back(val);
               }
               fin.read ((char*)(&num), 4);
               int n_objs2=num;
               if (n_objs2!=n_objs) {cout<<"Wrong reading of input binary file"<<endl; exit(3);}
            }
          fin.close();
        }
  //else {cout << "Unable to open input binary file"<<endl; exit(3);}

  
  int imax=int(xc.size());

  vector<double> dist_link(imax),xr(imax),yr(imax),zr(imax);
  vector<unsigned int> link(imax),used(imax), used_real(imax);
 
  min_x=DBL_MAX;max_x=-DBL_MAX;
  min_y=DBL_MAX;max_y=-DBL_MAX;
  min_z=DBL_MAX;max_z=-DBL_MAX;


  for (int i=0; i<imax; ++i) {
    used[i]=0;
    if (min_x>xc[i]) {min_x=xc[i];}
    if (min_y>yc[i]) {min_y=yc[i];}
    if (min_z>zc[i]) {min_z=zc[i];}

    if (max_x<xc[i]) {max_x=xc[i];}
    if (max_y<yc[i]) {max_y=yc[i];}
    if (max_z<zc[i]) {max_z=zc[i];}
  }

  double size_x=max_x-min_x;
  double size_y=max_y-min_y;
  double size_z=max_z-min_z;


  for (int o=1; o<=10; ++o) 
    num[o]=0;
  
  for (int i=0; i<imax; ++i) {
    for (int o=1; o<=10; ++o) {
      if (xc[i]>min_x+(max_x-min_x)*0.1*(o-1) && xc[i]<min_x+(max_x-min_x)*0.1*o){
    num[o]=num[o]+1; 
      }
    }
  }
  
  double c = 0., d = 0., e = 0., f = 0.; 
  for (int o=1; o<=10; ++o) {
    c += min_x+(max_x-min_x)*0.1*(o-0.5);
    d += num[o];
    e += pow(min_x+(max_x-min_x)*0.1*(o-0.5),2);
    f += (min_x+(max_x-min_x)*0.1*(o-0.5))*num[o];
  }
  double koef_k=(f-0.1*c*d)/(e-0.1*c*c);
  double koef_b=0.1*(d-koef_k*c);
 
  double integ=1/(koef_k*0.5*(max_x*max_x-min_x*min_x)+koef_b*(max_x-min_x));

  double mx,sum,rnd; 

  int rndd;

  for (int rndd=1; rndd<=n_rnd; ++rndd) { /* generation of random field */

    cout<<rndd<<" of "<<n_rnd<<" random realizations"<<endl;

    for (int i=0; i<imax; ++i) { /* generation of random field */

      if (mode=="t") {
    rnd=double(rand())/double(RAND_MAX);
    mx=min_x;
    sum=0;
    while (sum<rnd) {
      mx=mx+0.05;
      sum=sum+(koef_k*mx+koef_b)*integ*0.05;
    }
    xr[i]=mx;
      }

      if (mode=="n"||mode=="p")  {xr[i]=min_x+(max_x-min_x)*double(rand())/double(RAND_MAX);}
      yr[i]=min_y+(max_y-min_y)*double(rand())/double(RAND_MAX);
      zr[i]=min_z+(max_z-min_z)*double(rand())/double(RAND_MAX);
      used[i]=0;
      used_real[i]=0;

    }



    stringstream ss;
    ss << rndd;
    string str = ss.str();

    string name = dir_output+"displ";
    name.append("_");
    name.append(str);
    name.append("_");
    name.append(output);
    ofstream fout1(name.c_str());
    fout1.precision(7);

    //test

    name = dir_output+"test";
    name.append("_");
    name.append(str);
    name.append("_");
    name.append(output);
    ofstream ftest(name.c_str());
    ftest.precision(7);


    //test


    int jmax = imax;

    // looking for the closest pairs

    sum=0;
    while (sum<0.95*jmax) {
      
      int i=int(jmax*double(rand())/double(RAND_MAX)); // random number of real 
      int j=int(jmax*double(rand())/double(RAND_MAX)); // random number of random
      if (used[j]!=1 && used_real[i]!=1) {

    if (mode!="p") {
      min_k=DBL_MAX;
      dist_k=(xc[i]-xr[j])*(xc[i]-xr[j])+(yc[i]-yr[j])*(yc[i]-yr[j])+(zc[i]-zr[j])*(zc[i]-zr[j]);
      if (dist_k<min_k) min_k=dist_k; 
    }

    else {
      min_k=DBL_MAX;
      for (int xk=-1; xk<=1; ++xk) {
        for (int yk=-1; yk<=1; ++yk) {
          for (int zk=-1; zk<=1; ++zk) {
        dist_k=(xc[i]-xr[j]-xk*size_x)*(xc[i]-xr[j]-xk*size_x)+(yc[i]-yr[j]-yk*size_y)*(yc[i]-yr[j]-yk*size_y)+(zc[i]-zr[j]-zk*size_z)*(zc[i]-zr[j]-zk*size_z);
        if (dist_k<min_k) min_k=dist_k;
          }
        }
      }
    }

        dist=min_k;
        link[i]=j; dist_link[i]=dist; 
        used[j]=1; used_real[i]=1; sum=sum+1;
      }
    }
    while (sum<jmax) { 
      for (int i=0; i<jmax; ++i) {
        for (int j=0; j<jmax; ++j) {
      if (used[j]!=1 && used_real[i]!=1) {
        
        if (mode!="p") {
          min_k=DBL_MAX;
          dist_k=(xc[i]-xr[j])*(xc[i]-xr[j])+(yc[i]-yr[j])*(yc[i]-yr[j])+(zc[i]-zr[j])*(zc[i]-zr[j]);
          if (dist_k<min_k) {min_k=dist_k;}

        }

        else {
          min_k=DBL_MAX;
          for (int xk=-1; xk<=1; ++xk) {
        for (int yk=-1; yk<=1; ++yk) {
          for (int zk=-1; zk<=1; ++zk) {
            dist_k=(xc[i]-xr[j]-xk*size_x)*(xc[i]-xr[j]-xk*size_x)+(yc[i]-yr[j]-yk*size_y)*(yc[i]-yr[j]-yk*size_y)+(zc[i]-zr[j]-zk*size_z)*(zc[i]-zr[j]-zk*size_z);
            if (dist_k<min_k) {min_k=dist_k;}
          }
        }
          }
        }

        dist=min_k;
        link[i]=j; dist_link[i]=dist;
        used[j]=1; used_real[i]=1; sum=sum+1;
      }
    }
      }
    }



    //swaping

    for (int k=1; k<=n_max; ++k) { //iteration steps
      
      if (k==int(n_max*0.1)) { cout <<"\r.....10% completed\r"; cout.flush(); }
      if (k==int(n_max*0.15)) { cout <<"\r.....15% completed\r"; cout.flush(); }
      if (k==int(n_max*0.2)) { cout <<"\r.....20% completed\r"; cout.flush(); }
      if (k==int(n_max*0.25)) { cout <<"\r.....25% completed\r"; cout.flush(); }
      if (k==int(n_max*0.3)) { cout <<"\r.....30% completed\r"; cout.flush(); }
      if (k==int(n_max*0.4)) { cout <<"\r.....40% completed\r"; cout.flush(); }
      if (k==int(n_max*0.5)) { cout <<"\r.....50% completed\r"; cout.flush(); }
      if (k==int(n_max*0.6)) { cout <<"\r.....60% completed\r"; cout.flush(); }
      if (k==int(n_max*0.7)) { cout <<"\r.....70% completed\r"; cout.flush(); }
      if (k==int(n_max*0.8)) { cout <<"\r.....80% completed\r"; cout.flush(); }
      if (k==int(n_max*0.9)) { cout <<"\r.....90% completed\r"; cout.flush(); }
      
      for (int o=0; o<jmax; ++o) { //over all data points
    int ii=int(jmax*double(rand())/double(RAND_MAX));  //random choosing of two pairs real-random
    int i=int(jmax*double(rand())/double(RAND_MAX)); 
    sum1=dist_link[i]+dist_link[ii]; // sum of two lenghts pair size

    if (mode!="p") {
      min_k=DBL_MAX;
      dist_k=(xc[i]-xr[link[ii]])*(xc[i]-xr[link[ii]])+(yc[i]-yr[link[ii]])*(yc[i]-yr[link[ii]])+(zc[i]-zr[link[ii]])*(zc[i]-zr[link[ii]]); // new distance 1 after possible swaping
      if (dist_k<min_k) {min_k=dist_k;} 
    }

    else {
      min_k=DBL_MAX; 
      for (int xk=-1; xk<=1; ++xk) {
        for (int yk=-1; yk<=1; ++yk) {
          for (int zk=-1; zk<=1; ++zk) {
                dist_k=(xc[i]-xr[link[ii]]-xk*size_x)*(xc[i]-xr[link[ii]]-xk*size_x)+(yc[i]-yr[link[ii]]-yk*size_y)*(yc[i]-yr[link[ii]]-yk*size_y)+(zc[i]-zr[link[ii]]-zk*size_z)*(zc[i]-zr[link[ii]]-zk*size_z);
                if (dist_k<min_k) {min_k=dist_k;}
          }
        }
      }
    }

    dist_new4=min_k;

    if (mode!="p") {
      min_k=DBL_MAX;
      dist_k=(xc[ii]-xr[link[i]])*(xc[ii]-xr[link[i]])+(yc[ii]-yr[link[i]])*(yc[ii]-yr[link[i]])+(zc[ii]-zr[link[i]])*(zc[ii]-zr[link[i]]); // new distance 2 after possible swaping
      if (dist_k<min_k) {min_k=dist_k;} 
    }

    else {
      min_k=DBL_MAX;
      for (int xk=-1; xk<=1; ++xk) {
        for (int yk=-1; yk<=1; ++yk) {
          for (int zk=-1; zk<=1; ++zk) {
                dist_k=(xc[ii]-xr[link[i]]-xk*size_x)*(xc[ii]-xr[link[i]]-xk*size_x)+(yc[ii]-yr[link[i]]-yk*size_y)*(yc[ii]-yr[link[i]]-yk*size_y)+(zc[ii]-zr[link[i]]-zk*size_z)*(zc[ii]-zr[link[i]]-zk*size_z);
        if (dist_k<min_k) {min_k=dist_k;}
          }
        }
      }
    }  

    dist_new2=min_k;
    sum2=dist_new4+dist_new2; // new sum of distances after swapping

    if (sum2<sum1) { // case of decreasing of the distance sum after swaping
      int mih=link[i];
      link[i]=link[ii];
      link[ii]=mih;
      dist_link[i]=dist_new4;
      dist_link[ii]=dist_new2;
    }
      }

//test begin
double total=0;
  for (int o=0; o<jmax; ++o) { //over all data points
 total=total+dist_link[o];}
ftest<<total<<endl;
//test end
    }

    for (int i=0; i<jmax; ++i) { 
      shift_x=xr[link[i]]-xc[i];
      shift_y=yr[link[i]]-yc[i];
      shift_z=zr[link[i]]-zc[i];


      if (shift_x>size_x*0.5) {shift_x=xr[link[i]]-size_x-xc[i];}
      if (shift_x<-size_x*0.5) {shift_x=xr[link[i]]+size_x-xc[i];}

      if (shift_y>size_y*0.5) {shift_y=yr[link[i]]-size_y-yc[i];}
      if (shift_y<-size_y*0.5) {shift_y=yr[link[i]]+size_y-yc[i];}

      if (shift_z>size_z*0.5) {shift_z=zr[link[i]]-size_z-zc[i];}
      if (shift_z<-size_z*0.5) {shift_z=zr[link[i]]+size_z-zc[i];}

      fout1<<xr[link[i]]<<" "<<yr[link[i]]<<" "<<zr[link[i]]<<" "<<shift_x<<" "<<shift_y<<" "<<shift_z<<endl; 
      // output random position x,y,z, shift at randomization dx,dy,dz
    }

    fout1.clear(); fout1.close(); 
    ftest.clear(); ftest.close(); //for test 
  }

  cout <<"\r                  \r"<<endl;
      
  return 0;
}