#include <iostream>
#include <CosmoTool/miniargs.hpp>
#include <cmath>
#include "libqhull/qhull_a.h"
#include "voz.h"
#include "voz_io.hpp"

using namespace std;

#define DL for (d=0;d<3;d++)

bool checkParameters(int *numdiv, int *b)
{
  for (int i = 0; i < 3; i++)
    {
      if (numdiv[i] < 0 || b[i] < 0 || b[i] >= numdiv[i])
        return false;
    }
  return true;
}


int main(int argc, char *argv[]) {
  const float BF = std::numeric_limits<float>::max();

  int exitcode;
  pid_t i, j, np;
  PositionData pdata;
  coordT rtemp[3], *parts;
  coordT deladjs[3*MAXVERVER], points[3*MAXVERVER];
  pointT intpoints[3*MAXVERVER];
  FILE *pos, *out;
  char *posfile, outfile[200], *suffix, *outDir;
  PARTADJ *adjs;
  float *vols;
  realT predict, xmin,xmax,ymin,ymax,zmin,zmax;
  pid_t *orig;
 
  int isitinbuf;
  char isitinmain, d;
  pid_t nvp, nvpall, nvpbuf;
  realT width, width2, totwidth, totwidth2, bf, s, g;
  double border, boxsize[3];
  realT c[3];
  int b[3];
  int numdiv[3];
  double totalvol;
  
  CosmoTool::MiniArgDesc args[] = {
    { "POSITION FILE", MINIARG_STRING, &posfile },
    { "BORDER SIZE", MINIARG_DOUBLE, &border },
    { "BOX_X", MINIARG_DOUBLE, &boxsize[0] },
    { "BOX_Y", MINIARG_DOUBLE, &boxsize[1] },
    { "BOX_Z", MINIARG_DOUBLE, &boxsize[2] },
    { "SUFFIX", MINIARG_STRING, &suffix },
    { "NUM_DIVISION_X", MINIARG_INT, &numdiv[0] },
    { "NUM_DIVISION_Y", MINIARG_INT, &numdiv[1] },
    { "NUM_DIVISION_Z", MINIARG_INT, &numdiv[2] },
    { "B0", MINIARG_INT, &b[0] },
    { "B1", MINIARG_INT, &b[1] },
    { "B2", MINIARG_INT, &b[2] },
    { "OUTPUT DIRECTORY", MINIARG_STRING, &outDir },
    { 0, MINIARG_NULL, 0 }
  };

  if (!CosmoTool::parseMiniArgs(argc, argv, args))
    return 1;

  if (!checkParameters(numdiv, b))
    return 2;
  
  /* Boxsize should be the range in r, yielding a range 0-1 */
  if (!pdata.readFrom(posfile))
    return 3;

  (cout << pdata.np << " particles" << endl).flush();

  pdata.findExtrema();

  (cout << boost::format("np: %d, x: %f,%f; y: %f,%f; z: %f,%f") 
    % pdata.np
    % xyz_min[0] % xyz_max[0]
    % xyz_min[1] % xyz_max[1]
    % xyz_min[2] % xyz_max[2]).flush();


  width = 1./(float)numdiv;
  width2 = 0.5*width;
  if (border > 0.) bf = border;
  else bf = 0.1;
      /* In units of 0-1, the thickness of each subregion's buffer*/
  totwidth = width+2.*bf;
  totwidth2 = width2 + bf;
  
  s = width/(float)NGUARD;
  if ((bf*bf - 2.*s*s) < 0.) {
    printf("bf = %f, s = %f.\n",bf,s);
    printf("Not enough guard points for given border.\nIncrease guards to >= %f\n.",
	   sqrt(2.)*width/bf);
    exit(0);
  }
  g = (bf / 2.)*(1. + sqrt(1 - 2.*s*s/(bf*bf)));
  printf("s = %f, bf = %f, g = %f.\n",s,bf,g);

  fflush(stdout);

  adjs = (PARTADJ *)malloc(np*sizeof(PARTADJ));
  if (adjs == NULL) {
    printf("Unable to allocate adjs\n");
    exit(0);
  }
  
  DL c[d] = ((float)b[d])*width;
  printf("c: %f,%f,%f\n",c[0],c[1],c[2]);
  /* Assign temporary array*/
  nvpbuf = 0; /* Number of particles to tesselate, including
		 buffer */
  nvp = 0; /* Without the buffer */
  for (i=0; i<np; i++) {
    isitinbuf = 1;
    isitinmain = 1;
    DL {
      rtemp[d] = (double)r[i][d] - (double)c[d];
      if (rtemp[d] > 0.5) rtemp[d] --;
      if (rtemp[d] < -0.5) rtemp[d] ++;
      isitinbuf = isitinbuf && (fabs(rtemp[d]) < totwidth2);
      isitinmain = isitinmain && (fabs(rtemp[d]) <= width2);
    }
  
    if (isitinbuf) nvpbuf++;
    if (isitinmain) nvp++;
  }
  
  nvpbuf += 6*(NGUARD+1)*(NGUARD+1); /* number of guard
					points */

  parts = (coordT *)malloc(3*nvpbuf*sizeof(coordT));
  orig = (pid_t *)malloc(nvpbuf*sizeof(pid_t));

  if (parts == NULL) {
    printf("Unable to allocate parts\n");
    fflush(stdout);
  }
  if (orig == NULL) {
    printf("Unable to allocate orig\n");
    fflush(stdout);
  }

  nvp = 0; nvpall = 0; /* nvp = number of particles without buffer */
  xmin = BF; xmax = -BF; ymin = BF; ymax = -BF; zmin = BF; zmax = -BF;
  for (i=0; i<np; i++) {
    isitinmain = 1;
    DL {
      rtemp[d] = (realT)r[i][d] - (realT)c[d];
      if (rtemp[d] > 0.5) rtemp[d] --;
      if (rtemp[d] < -0.5) rtemp[d] ++;
      isitinmain = isitinmain && (fabs(rtemp[d]) <= width2);
    }
    if (isitinmain) {
      parts[3*nvp] = rtemp[0];
      parts[3*nvp+1] = rtemp[1];
      parts[3*nvp+2] = rtemp[2];
      orig[nvp] = i;
      nvp++;
      if (rtemp[0] < xmin) xmin = rtemp[0];
      if (rtemp[0] > xmax) xmax = rtemp[0];
      if (rtemp[1] < ymin) ymin = rtemp[1];
      if (rtemp[1] > ymax) ymax = rtemp[1];
      if (rtemp[2] < zmin) zmin = rtemp[2];
      if (rtemp[2] > zmax) zmax = rtemp[2];
    }
  }
  printf("nvp = %d\n",nvp);
  printf("x: %f,%f; y: %f,%f; z:%f,%f\n",xmin,xmax,ymin,ymax,zmin,zmax);
  nvpbuf = nvp;
  for (i=0; i<np; i++) {
    isitinbuf = 1;
    isitinmain = 1;

    DL {
      rtemp[d] = (realT)r[i][d] - (realT)c[d];
      if (rtemp[d] > 0.5) rtemp[d] --;
      if (rtemp[d] < -0.5) rtemp[d] ++;
      isitinbuf = isitinbuf && (fabs(rtemp[d])<totwidth2);
      isitinmain = isitinmain && (fabs(rtemp[d]) <= width2);
    }
    if (isitinbuf && !isitinmain) {
      /*printf("%3.3f ",sqrt(rtemp[0]*rtemp[0] + rtemp[1]*rtemp[1] +
	rtemp[2]*rtemp[2]));
	printf("|%2.2f,%2.2f,%2.2f,%f,%f",r[i][0],r[i][1],r[i][2],width2,totwidth2);*/
      parts[3*nvpbuf] = rtemp[0];
      parts[3*nvpbuf+1] = rtemp[1];
      parts[3*nvpbuf+2] = rtemp[2];
      orig[nvpbuf] = i;

      nvpbuf++;
      if (rtemp[0] < xmin) xmin = rtemp[0];
      if (rtemp[0] > xmax) xmax = rtemp[0];
      if (rtemp[1] < ymin) ymin = rtemp[1];
      if (rtemp[1] > ymax) ymax = rtemp[1];
      if (rtemp[2] < zmin) zmin = rtemp[2];
      if (rtemp[2] > zmax) zmax = rtemp[2];
    }
  }
  printf("nvpbuf = %d\n",nvpbuf);
  printf("x: %f,%f; y: %f,%f; z:%f,%f\n",xmin,xmax,ymin,ymax,zmin,zmax);
  nvpall = nvpbuf;
  predict = pow(width+2.*bf,3)*(float)np;
  printf("There should be ~ %f points; there are %d\n",predict,nvpbuf);

  pdata.destroy();
  
  /* Add guard points */
  for (i=0; i<NGUARD+1; i++) {
    for (j=0; j<NGUARD+1; j++) {
      /* Bottom */
      parts[3*nvpall]   = -width2 + (realT)i * s;
      parts[3*nvpall+1] = -width2 + (realT)j * s;
      parts[3*nvpall+2] = -width2 - g;
      nvpall++;
      /* Top */
      parts[3*nvpall]   = -width2 + (realT)i * s;
      parts[3*nvpall+1] = -width2 + (realT)j * s;
      parts[3*nvpall+2] = width2 + g;
      nvpall++;
    }
  }
  for (i=0; i<NGUARD+1; i++) { /* Don't want to overdo the corners*/
    for (j=0; j<NGUARD+1; j++) {
      parts[3*nvpall]   = -width2 + (realT)i * s;
      parts[3*nvpall+1] = -width2 - g;
      parts[3*nvpall+2] = -width2 + (realT)j * s;
      nvpall++;
      
      parts[3*nvpall]   = -width2 + (realT)i * s;
      parts[3*nvpall+1] = width2 + g;
      parts[3*nvpall+2] = -width2 + (realT)j * s;
      nvpall++;
    }
  }
  for (i=0; i<NGUARD+1; i++) {
    for (j=0; j<NGUARD+1; j++) {
      parts[3*nvpall]   = -width2 - g;
      parts[3*nvpall+1] = -width2 + (realT)i * s;
      parts[3*nvpall+2] = -width2 + (realT)j * s;
      nvpall++;
      
      parts[3*nvpall]   = width2 + g;
      parts[3*nvpall+1] = -width2 + (realT)i * s;
      parts[3*nvpall+2] = -width2 + (realT)j * s;
      nvpall++;
    }
  }
  xmin = BF; xmax = -BF; ymin = BF; ymax = -BF; zmin = BF; zmax = -BF;
  for (i=nvpbuf;i<nvpall;i++) {
    if (parts[3*i] < xmin) xmin = parts[3*i];
    if (parts[3*i] > xmax) xmax = parts[3*i];
    if (parts[3*i+1] < ymin) ymin = parts[3*i+1];
    if (parts[3*i+1] > ymax) ymax = parts[3*i+1];
    if (parts[3*i+2] < zmin) zmin = parts[3*i+2];
    if (parts[3*i+2] > zmax) zmax = parts[3*i+2];
  }
  
  printf("Added guard points to total %d points (should be %d)\n",nvpall,
	 nvpbuf + 6*(NGUARD+1)*(NGUARD+1));
  printf("x: %f,%f; y: %f,%f; z:%f,%f\n",xmin,xmax,ymin,ymax,zmin,zmax);
  
  /* Do tesselation*/
  printf("File read.  Tessellating ...\n"); fflush(stdout);
  exitcode = delaunadj(parts, nvp, nvpbuf, nvpall, &adjs);
  if (exitcode != 0)
   {
     printf("Error while tesselating. Stopping here."); fflush(stdout);
     exit(1);
   }
  
  /* Now calculate volumes*/
  printf("Now finding volumes ...\n"); fflush(stdout);
  vols = (float *)malloc(nvp*sizeof(float));
  
  for (i=0; i<nvp; i++) { /* Just the original particles
			     Assign adjacency coordinate array*/
    /* Volumes */
    for (j = 0; j < adjs[i].nadj; j++)
      DL {
	deladjs[3*j + d] = parts[3*adjs[i].adj[j]+d] - parts[3*i+d];
	if (deladjs[3*j+d] < -0.5) deladjs[3*j+d]++;
	if (deladjs[3*j+d] > 0.5) deladjs[3*j+d]--;
      }
    
    exitcode = vorvol(deladjs, points, intpoints, adjs[i].nadj, &(vols[i]));
    vols[i] *= (float)np;
    if (i % 1000 == 0)
      printf("%d: %d, %f\n",i,adjs[i].nadj,vols[i]);
  }

  /* Get the adjacencies back to their original values */

  for (i=0; i<nvp; i++)
    for (j = 0; j < adjs[i].nadj; j++)
      adjs[i].adj[j] = orig[adjs[i].adj[j]];
  
  totalvol = 0.;
  for (i=0;i<nvp; i++) {
    totalvol += (double)vols[i];
  }
  printf("Average volume = %g\n",totalvol/(float)nvp);
 
  /* Now the output!
     First number of particles */
  sprintf(outfile,"%s/part.%s.%02d.%02d.%02d",outDir,suffix,b[0],b[1],b[2]);

  printf("Output to %s\n\n",outfile);
  out = fopen(outfile,"w");
  if (out == NULL) {
    printf("Unable to open %s\n",outfile);
    exit(0);
  }
  fwrite(&np,1, sizeof(int),out);
  fwrite(&nvp,1, sizeof(int),out);
  printf("nvp = %d\n",nvp);

  /* Tell us where the original particles were */
  fwrite(orig,sizeof(pid_t),nvp,out);
  /* Volumes*/
  fwrite(vols,sizeof(float),nvp,out);
  /* Adjacencies */
  for (i=0;i<nvp;i++) {
    fwrite(&(adjs[i].nadj),1,sizeof(pid_t),out);
    if (adjs[i].nadj > 0)
      fwrite(adjs[i].adj,adjs[i].nadj,sizeof(pid_t),out);
    else printf("0");
  }
  fclose(out);
  
  return(0);
}